ClicpixAnalysis.C 24.4 KB
Newer Older
1
2
3
4
5
6
7
8
#include "ClicpixAnalysis.h"
#include "Pixel.h"
#include "Cluster.h"
#include "Track.h"

ClicpixAnalysis::ClicpixAnalysis(bool debugging)
: Algorithm("ClicpixAnalysis"){
  debug = debugging;
9
10
  m_associationCut = 0.05; // 100 um
  m_proximityCut = 0.0005; // 125 um
11
  timepix3Telescope = false;
12
13
}

14
15
16
17
18
19
template <typename T>
std::string convertToString(T number) {
  std::ostringstream ss;
  ss << number;
  return ss.str();
}
20
21
22
23
24

void ClicpixAnalysis::initialise(Parameters* par){
 
  parameters = par;
  
25
26
27
  // Initialise member variables
  m_eventNumber = 0; m_triggerNumber = 0; dutID = parameters->DUT; m_lostHits = 0.;

28
29
30
31
32
33
34
  // Cluster/pixel histograms
  hHitPixels = new TH2F("hHitPixels","hHitPixels",64,0,64,64,0,64);
  hColumnHits = new TH1F("hColumnHits","hColumnHits",64,0,64);
  hRowHits = new TH1F("hRowHits","hRowHits",64,0,64);
  
  hClusterSizeAll = new TH1F("hClusterSizeAll","hClusterSizeAll",20,0,20);
  hClusterTOTAll = new TH1F("hClusterTOTAll","hClusterTOTAll",50,0,50);
35
  hClustersPerEvent = new TH1F("hClustersPerEvent","hClustersPerEvent",500,0,500);
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
  hClustersVersusEventNo = new TH1F("hClustersVersusEventNo","hClustersVersusEventNo",60000,0,60000);
  hGlobalClusterPositions = new TH2F("hGlobalClusterPositions","hGlobalClusterPositions",200,-2.0,2.0,300,-1.,2);
  
  hClusterWidthRow = new TH1F("hClusterWidthRow","hClusterWidthRow",20,0,20);
  hClusterWidthCol = new TH1F("hClusterWidthCol","hClusterWidthCol",20,0,20);
  
  // Track histograms
  hGlobalTrackDifferenceX = new TH1F("hGlobalTrackDifferenceX","hGlobalTrackDifferenceX",1000,-10.,10.);
  hGlobalTrackDifferenceY = new TH1F("hGlobalTrackDifferenceY","hGlobalTrackDifferenceY",1000,-10.,10.);
  
  hGlobalResidualsX = new TH1F("hGlobalResidualsX","hGlobalResidualsX",600,-0.3,0.3);
  hGlobalResidualsY = new TH1F("hGlobalResidualsY","hGlobalResidualsY",600,-0.3,0.3);
  hAbsoluteResiduals = new TH1F("hAbsoluteResiduals","hAbsoluteResiduals",600,0.,0.600);
  hGlobalResidualsXversusX = new TH2F("hGlobalResidualsXversusX","hGlobalResidualsXversusX",600,-3.,3.,600,-0.3,0.3);
  hGlobalResidualsXversusY = new TH2F("hGlobalResidualsXversusY","hGlobalResidualsXversusY",600,-3.,3.,600,-0.3,0.3);
  hGlobalResidualsYversusY = new TH2F("hGlobalResidualsYversusY","hGlobalResidualsYversusY",600,-3.,3.,600,-0.3,0.3);
  hGlobalResidualsYversusX = new TH2F("hGlobalResidualsYversusX","hGlobalResidualsYversusX",600,-3.,3.,600,-0.3,0.3);
  
  hGlobalResidualsXversusColWidth = new TH2F("hGlobalResidualsXversusColWidth","hGlobalResidualsXversusColWidth",30,0,30,600,-0.3,0.3);
  hGlobalResidualsXversusRowWidth = new TH2F("hGlobalResidualsXversusRowWidth","hGlobalResidualsXversusRowWidth",30,0,30,600,-0.3,0.3);
  hGlobalResidualsYversusColWidth = new TH2F("hGlobalResidualsYversusColWidth","hGlobalResidualsYversusColWidth",30,0,30,600,-0.3,0.3);
  hGlobalResidualsYversusRowWidth = new TH2F("hGlobalResidualsYversusRowWidth","hGlobalResidualsYversusRowWidth",30,0,30,600,-0.3,0.3);
  
59
60
61
  hTrackInterceptRow = new TH1F("hTrackInterceptRow","hTrackInterceptRow",660,-1.,65.);
  hTrackInterceptCol = new TH1F("hTrackInterceptCol","hTrackInterceptCol",660,-1.,65.);

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
  hAbsoluteResidualMap = new TH2F("hAbsoluteResidualMap","hAbsoluteResidualMap",50,0,50,25,0,25);
  hXresidualVersusYresidual = new TH2F("hXresidualVersusYresidual","hXresidualVersusYresidual",600,-0.3,0.3,600,-0.3,0.3);
  
  hAssociatedClustersPerEvent = new TH1F("hAssociatedClustersPerEvent","hAssociatedClustersPerEvent",50,0,50);;
  hAssociatedClustersVersusEventNo = new TH1F("hAssociatedClustersVersusEventNo","hAssociatedClustersVersusEventNo",60000,0,60000);
  hAssociatedClustersVersusTriggerNo = new TH1F("hAssociatedClustersVersusTriggerNo","hAssociatedClustersVersusTriggerNo",50,0,50);
  hAssociatedClusterRow = new TH1F("hAssociatedClusterRow","hAssociatedClusterRow",64,0,64);
  hAssociatedClusterColumn = new TH1F("hAssociatedClusterColumn","hAssociatedClusterColumn",64,0,64);
  hFrameEfficiency = new TH1F("hFrameEfficiency","hFrameEfficiency",6000,0,6000);
  hFrameTracks = new TH1F("hFrameTracks","hFrameTracks",6000,0,6000);
  hFrameTracksAssociated = new TH1F("hFrameTracksAssociated","hFrameTracksAssociated",6000,0,6000);
  
  hClusterSizeAssociated = new TH1F("hClusterSizeAssociated","hClusterSizeAssociated",20,0,20);
  hClusterWidthRowAssociated = new TH1F("hClusterWidthRowAssociated","hClusterWidthRowAssociated",20,0,20);
  hClusterWidthColAssociated = new TH1F("hClusterWidthColAssociated","hClusterWidthColAssociated",20,0,20);
  
  hClusterTOTAssociated = new TH1F("hClusterTOTAssociated","hClusterTOTAssociated",50,0,50);
  hClusterTOTAssociated1pix = new TH1F("hClusterTOTAssociated1pix","hClusterTOTAssociated1pix",50,0,50);
  hClusterTOTAssociated2pix = new TH1F("hClusterTOTAssociated2pix","hClusterTOTAssociated2pix",50,0,50);
  hClusterTOTAssociated3pix = new TH1F("hClusterTOTAssociated3pix","hClusterTOTAssociated3pix",50,0,50);
  hClusterTOTAssociated4pix = new TH1F("hClusterTOTAssociated4pix","hClusterTOTAssociated4pix",50,0,50);
  
  hPixelResponseX = new TH1F("hPixelResponseX","hPixelResponseX",600,-0.3,0.3);
  hPixelResponseY = new TH1F("hPixelResponseY","hPixelResponseY",600,-0.3,0.3);
  
  hPixelResponseGlobalX = new TH1F("hPixelResponseGlobalX","hPixelResponseGlobalX",600,-0.3,0.3);
  hPixelResponseGlobalY = new TH1F("hPixelResponseGlobalY","hPixelResponseGlobalY",600,-0.3,0.3);
  
  hPixelResponseXOddCol = new TH1F("hPixelResponseXOddCol","hPixelResponseXOddCol",600,-0.3,0.3);
  hPixelResponseXEvenCol = new TH1F("hPixelResponseXEvenCol","hPixelResponseXEvenCol",600,-0.3,0.3);
  hPixelResponseYOddCol = new TH1F("hPixelResponseYOddCol","hPixelResponseYOddCol",600,-0.3,0.3);
  hPixelResponseYEvenCol = new TH1F("hPixelResponseYEvenCol","hPixelResponseYEvenCol",600,-0.3,0.3);
  
  hEtaDistributionX = new TH2F("hEtaDistributionX","hEtaDistributionX",25,0,25,25,0,25);
  hEtaDistributionY = new TH2F("hEtaDistributionY","hEtaDistributionY",25,0,25,25,0,25);
  
  hResidualsLocalRow2pix = new TH1F("hResidualsLocalRow2pix","hResidualsLocalRow2pix",600,-0.3,0.3);
  hResidualsLocalCol2pix = new TH1F("hResidualsLocalCol2pix","hResidualsLocalCol2pix",600,-0.3,0.3);
  hClusterTOTRow2pix = new TH1F("hClusterTOTRow2pix","hClusterTOTRow2pix",50,0,50);
  hClusterTOTCol2pix = new TH1F("hClusterTOTCol2pix","hClusterTOTCol2pix",50,0,50);
  hClusterTOTRatioRow2pix = new TH1F("hClusterTOTRatioRow2pix","hClusterTOTRatioRow2pix",100,0,1);
  hClusterTOTRatioCol2pix = new TH1F("hClusterTOTRatioCol2pix","hClusterTOTRatioCol2pix",100,0,1);
  hPixelTOTRow2pix = new TH1F("hPixelTOTRow2pix","hPixelTOTRow2pix",50,0,50);
  hPixelTOTCol2pix = new TH1F("hPixelTOTCol2pix","hPixelTOTCol2pix",50,0,50);
  
  hResidualsLocalRow2pixClusterTOT = new TH2F("hResidualsLocalRow2pixClusterTOT","hResidualsLocalRow2pixClusterTOT",50,0,50,600,-0.3,0.3);
  hResidualsLocalRow2pixPixelIntercept = new TH2F("hResidualsLocalRow2pixPixelIntercept","hResidualsLocalRow2pixPixelIntercept",50,0,25,600,-0.3,0.3);
  
  // Maps
  hTrackIntercepts = new TH2F("hTrackIntercepts","hTrackIntercepts",200,-2.0,2.0,300,-1.,2);
  hTrackInterceptsAssociated = new TH2F("hTrackInterceptsAssociated","hTrackInterceptsAssociated",200,-2.0,2.0,300,-1.,2);
  hGlobalAssociatedClusterPositions = new TH2F("hGlobalAssociatedClusterPositions","hGlobalAssociatedClusterPositions",200,-2.0,2.0,300,-1.,2);
  hTrackInterceptsPixel = new TH2F("hTrackInterceptsPixel","hTrackInterceptsPixel",50,0,50,25,0,25);
  hTrackInterceptsPixelAssociated = new TH2F("hTrackInterceptsPixelAssociated","hTrackInterceptsPixelAssociated",50,0,50,25,0,25);
  hTrackInterceptsChip = new TH2F("hTrackInterceptsChip","hTrackInterceptsChip",65,-0.5,64.5,65,-0.5,64.5);
  hTrackInterceptsChipAssociated = new TH2F("hTrackInterceptsChipAssociated","hTrackInterceptsChipAssociated",65,-0.5,64.5,65,-0.5,64.5);
118
  hTrackInterceptsChipUnassociated = new TH2F("hTrackInterceptsChipUnassociated","hTrackInterceptsChipUnassociated",65,-0.5,64.5,65,-0.5,64.5);
119
120
121
122
123
124
  hTrackInterceptsChipLost = new TH2F("hTrackInterceptsChipLost","hTrackInterceptsChipLost",65,-0.5,64.5,65,-0.5,64.5);
  
  hPixelEfficiencyMap = new TH2F("hPixelEfficiencyMap","hPixelEfficiencyMap",50,0,50,25,0,25);
  hChipEfficiencyMap = new TH2F("hChipEfficiencyMap","hChipEfficiencyMap",65,-0.5,64.5,65,-0.5,64.5);
  hGlobalEfficiencyMap = new TH2F("hGlobalEfficiencyMap","hGlobalEfficiencyMap",200,-2.0,2.0,300,-1.,2);
  
125
126
127
128
  hInterceptClusterSize1 = new TH2F("hInterceptClusterSize1","hInterceptClusterSize1",25,0,25,25,0,25);
  hInterceptClusterSize2 = new TH2F("hInterceptClusterSize2","hInterceptClusterSize2",25,0,25,25,0,25);
  hInterceptClusterSize3 = new TH2F("hInterceptClusterSize3","hInterceptClusterSize3",25,0,25,25,0,25);
  hInterceptClusterSize4 = new TH2F("hInterceptClusterSize4","hInterceptClusterSize4",25,0,25,25,0,25);
129
  
130
131
132
133
134
135
136
137
138
139
140
141
142
143
  m_nBinsX=32; m_nBinsY=32;
  hMapClusterSizeAssociated = new TH2F("hMapClusterSizeAssociated","hMapClusterSizeAssociated",m_nBinsX,0,parameters->detector[dutID]->nPixelsX(),m_nBinsY,0,parameters->detector[dutID]->nPixelsY());
  
  for(int x=0;x<m_nBinsX;x++){
    for(int y=0;y<m_nBinsY;y++){
      int id = x+y*m_nBinsX;
      std::string name = "hMapClusterTOTAssociated1pix"+convertToString(id);
      TH1F* hMapEntryClusterTOTAssociated1pix = new TH1F(name.c_str(),name.c_str(),50,0,50);
      hMapClusterTOTAssociated1pix[id] = hMapEntryClusterTOTAssociated1pix;
    }
  }
  
  
  
144
145
146
147
148
149
150
151
152
153
154
}

StatusCode ClicpixAnalysis::run(Clipboard* clipboard){
  
  // Get the clicpix clusters in this event
  Clusters* clusters = (Clusters*)clipboard->get(dutID,"clusters");
  if(clusters == NULL){
    if(debug) tcout<<"No clusters for "<<dutID<<" on the clipboard"<<endl;
    return Success;
  }
  
Daniel Hynds's avatar
Daniel Hynds committed
155
  // If this is the first or last trigger don't use the data
156
  if( (m_triggerNumber == 0 || m_triggerNumber == 19) && !timepix3Telescope){
157
158
159
160
161
162
    m_eventNumber++;
    m_triggerNumber++;
    return Success;
  }
  
  // Fill the histograms that only need clusters/pixels
Daniel Hynds's avatar
Daniel Hynds committed
163
  fillClusterHistos(clusters);
164
  
165
166
167
168
169
170
171
  // Get the tracks in this event
  Tracks* tracks = (Tracks*)clipboard->get("tracks");
  if(tracks == NULL){
    if(debug) tcout<<"No tracks on the clipboard"<<endl;
    return Success;
  }
  
172
173
174
175
176
177
178
179
180
181
182
  //Set counters
  double nClustersAssociated=0, nValidTracks=0;
  
  // Loop over tracks, and compare with Clicpix clusters
  Tracks::iterator itTrack;
  for (itTrack = tracks->begin(); itTrack != tracks->end(); itTrack++) {
    
    // Get the track
    Track* track = (*itTrack);
    if (!track) continue;
    
183
184
185
    // Cut on the track chi2/ndof
    if(track->chi2ndof() < 3.0) continue;

186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
    // Get the track intercept with the clicpix plane (global and local co-ordinates)
    PositionVector3D< Cartesian3D<double> > trackIntercept = parameters->detector[dutID]->getIntercept(track);
    PositionVector3D< Cartesian3D<double> > trackInterceptLocal = *(parameters->detector[dutID]->m_globalToLocal) * trackIntercept;
    
    // Plot the difference between track intercepts and all clicpix clusters
    // Also record which cluster is the closest
    Clusters::iterator itCorrelate;
    bool matched=false; Clusters::iterator bestCluster;
    double xresidualBest = 10000.; double yresidualBest = 10000.;
    double absoluteresidualBest = sqrt(xresidualBest*xresidualBest + yresidualBest*yresidualBest);

    for (itCorrelate = clusters->begin(); itCorrelate != clusters->end(); ++itCorrelate) {
     
      // Get the distance between this cluster and the track intercept (global)
      double xcorr = (*itCorrelate)->globalX()-trackIntercept.X();
      double ycorr = (*itCorrelate)->globalY()-trackIntercept.Y();
      
      // Fill histograms on correlations
      hGlobalTrackDifferenceX->Fill(xcorr);
      hGlobalTrackDifferenceY->Fill(ycorr);
      
      // Look if this cluster can be considered associated.
      // Cut on the distance from the track in x and y
Daniel Hynds's avatar
Daniel Hynds committed
209
      if ( fabs(ycorr) > m_associationCut || fabs(xcorr) > m_associationCut) continue;
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
      
      // Found a matching cluster
      matched=true;
      
      // Check if this is a better match than the previously found cluster, and store the information
      double absoluteresidual = sqrt(xcorr*xcorr + ycorr*ycorr);
      if( absoluteresidual <= absoluteresidualBest){
        absoluteresidualBest=absoluteresidual;
        xresidualBest=xcorr;
        yresidualBest=ycorr;
        bestCluster=itCorrelate;
      }

    }
  
    // Get the track intercept position along the chip
    double chipInterceptCol = parameters->detector[dutID]->getColumn(trackInterceptLocal);
    double chipInterceptRow = parameters->detector[dutID]->getRow(trackInterceptLocal);

    // Get the track intercept position along the pixel
    double pixelInterceptX = parameters->detector[dutID]->inPixelX(trackInterceptLocal);
    double pixelInterceptY = parameters->detector[dutID]->inPixelY(trackInterceptLocal);
    
    // Cut on the track intercept - this makes sure that it actually went through the chip
    if( chipInterceptCol < 0.5 ||
       chipInterceptRow < 0.5 ||
       chipInterceptCol > (parameters->detector[dutID]->nPixelsX()-0.5) ||
       chipInterceptRow > (parameters->detector[dutID]->nPixelsY()-0.5) ) continue;
    
    // Check if the hit is near a masked pixel
240
241
    bool hitMasked = checkMasked(chipInterceptRow,chipInterceptCol);
    if(hitMasked) continue;
242
243
    
    // Check there are no other tracks nearby
244
245
    bool proximityCut = checkProximity(track,tracks);
    if(proximityCut) continue;
246
247
248
249
250
251
252
253
254
255
256
    
    // Now have tracks that went through the device
    hTrackIntercepts->Fill(trackIntercept.X(),trackIntercept.Y());
    hFrameTracks->Fill(m_eventNumber);
    nValidTracks++;
    hTrackInterceptsChip->Fill(chipInterceptCol,chipInterceptRow);
    hTrackInterceptsPixel->Fill(pixelInterceptX,pixelInterceptY);
    
    if(matched){
      
      // Some test plots of pixel response function
257
      fillResponseHistos(trackIntercept.X(),trackIntercept.Y(),(*bestCluster));
258
      
259
260
261
      // Add this cluster to the list of associated clusters held by this track. This will allow alignment to take place
      track->addAssociatedCluster(*bestCluster);
      
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
      // Now fill all of the histograms/efficiency counters that we want
      hTrackInterceptsAssociated->Fill(trackIntercept.X(),trackIntercept.Y());
      hGlobalResidualsX->Fill(xresidualBest);
      hGlobalResidualsY->Fill(yresidualBest);
      hGlobalAssociatedClusterPositions->Fill((*bestCluster)->globalX(),(*bestCluster)->globalY());
      nClustersAssociated++;
      hAssociatedClusterRow->Fill((*bestCluster)->row());
      hAssociatedClusterColumn->Fill((*bestCluster)->column());
      hTrackInterceptsChipAssociated->Fill(chipInterceptCol,chipInterceptRow);
      hTrackInterceptsPixelAssociated->Fill(pixelInterceptX,pixelInterceptY);
      hFrameTracksAssociated->Fill(m_eventNumber);
      hGlobalResidualsXversusX->Fill(trackIntercept.X(),xresidualBest);
      hGlobalResidualsXversusY->Fill(trackIntercept.Y(),xresidualBest);
      hGlobalResidualsYversusY->Fill(trackIntercept.Y(),yresidualBest);
      hGlobalResidualsYversusX->Fill(trackIntercept.X(),yresidualBest);
//      hGlobalResidualsXversusColWidth->Fill((*bestCluster)->colWidth(),xresidualBest);
//      hGlobalResidualsXversusRowWidth->Fill((*bestCluster)->rowWidth(),xresidualBest);
//      hGlobalResidualsYversusColWidth->Fill((*bestCluster)->colWidth(),yresidualBest);
//      hGlobalResidualsYversusRowWidth->Fill((*bestCluster)->rowWidth(),yresidualBest);
      hAbsoluteResidualMap->Fill(pixelInterceptX,pixelInterceptY,sqrt(xresidualBest*xresidualBest + yresidualBest*yresidualBest));
      hXresidualVersusYresidual->Fill(xresidualBest,yresidualBest);
      hAbsoluteResiduals->Fill(sqrt(xresidualBest*xresidualBest + yresidualBest*yresidualBest));
      hClusterTOTAssociated->Fill((*bestCluster)->tot());
      hClusterSizeAssociated->Fill((*bestCluster)->size());
//      hClusterWidthColAssociated->Fill((*bestCluster)->colWidth());
//      hClusterWidthRowAssociated->Fill((*bestCluster)->rowWidth());
288
289
            
      hMapClusterSizeAssociated->Fill(chipInterceptCol,chipInterceptRow,(*bestCluster)->tot());
290
291
292
293
      
      if((*bestCluster)->size() == 1){
        hClusterTOTAssociated1pix->Fill((*bestCluster)->tot());
        hInterceptClusterSize1->Fill(pixelInterceptX,pixelInterceptY);
294
295
        int id = floor(chipInterceptCol*m_nBinsX/parameters->detector[dutID]->nPixelsX())+floor(chipInterceptRow*m_nBinsY/parameters->detector[dutID]->nPixelsY())*m_nBinsX;
        hMapClusterTOTAssociated1pix[id]->Fill((*bestCluster)->tot());
296
297
298
299
300
301
302
303
304
305
306
307
308
309
      }
      if((*bestCluster)->size() == 2){
        hClusterTOTAssociated2pix->Fill((*bestCluster)->tot());
        hInterceptClusterSize2->Fill(pixelInterceptX,pixelInterceptY);
      }
      if((*bestCluster)->size() == 3){
        hClusterTOTAssociated3pix->Fill((*bestCluster)->tot());
        hInterceptClusterSize3->Fill(pixelInterceptX,pixelInterceptY);
      }
      if((*bestCluster)->size() == 4){
        hClusterTOTAssociated4pix->Fill((*bestCluster)->tot());
        hInterceptClusterSize4->Fill(pixelInterceptX,pixelInterceptY);
      }
    }else{
310
311
312
      
      // Search for lost hits. Basically loop through all pixels in all clusters and see if any are close.
      // Large clusters (such as from deltas) can pull the cluster centre sufficiently far from the track
313
      bool pixelMatch = false; int size=0;
314
/*      for (itCorrelate = clusters->begin(); itCorrelate != clusters->end(); ++itCorrelate) {
315
        if(pixelMatch) break;
316
        // Loop over pixels
317
318
319
320
        Pixels::const_iterator itPixel;
        for (itPixel = (*itCorrelate)->pixels().begin(); itPixel != (*itCorrelate)->pixels().end(); itPixel++) {
          
          // Get the pixel global position
321
322
          tcout<<"New pixel"<<endl;
          tcout<<"Row = "<<(*itPixel)->m_row<<", column = "<<(*itPixel)->m_column<<endl;
323
324
          PositionVector3D<Cartesian3D<double> > pixelPositionLocal = parameters->detector[dutID]->getLocalPosition((*itPixel)->m_row,(*itPixel)->m_column);
          PositionVector3D<Cartesian3D<double> > pixelPositionGlobal = *(parameters->detector[dutID]->m_localToGlobal) * pixelPositionLocal;
325
          
326
327
328
          // Check if it is close to the track
          if( fabs( pixelPositionGlobal.X() - trackIntercept.X() ) > m_associationCut ||
              fabs( pixelPositionGlobal.Y() - trackIntercept.Y() ) > m_associationCut ) continue;
329
330
          
          pixelMatch=true;
331
          break;
332
333
334
        }
      }
      // Increment counter for number of hits found this way
335
*/      if(pixelMatch){
336
337
        m_lostHits++;
        hTrackInterceptsChipLost->Fill(chipInterceptCol,chipInterceptRow);
338
      }//*/
339
340
      
      if(!pixelMatch) hTrackInterceptsChipUnassociated->Fill(chipInterceptCol,chipInterceptRow);
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
    }
  }
  
  // Histograms relating to this mimosa frame
  hAssociatedClustersPerEvent->Fill(nClustersAssociated);
  hAssociatedClustersVersusEventNo->Fill(m_eventNumber,nClustersAssociated);
  hAssociatedClustersVersusTriggerNo->Fill(m_triggerNumber,nClustersAssociated);
  if(nValidTracks != 0) hFrameEfficiency->Fill(m_eventNumber,nClustersAssociated/nValidTracks);

	// Increment event counter and trigger number
  m_eventNumber++; m_triggerNumber++;
  
  // Return value telling analysis to keep running
  return Success;
}

void ClicpixAnalysis::finalise(){
  
  if(debug) tcout<<"Analysed "<<m_eventNumber<<" events"<<endl;
  
Daniel Hynds's avatar
Daniel Hynds committed
361
362
363
364
365
366
367
368
369
370
371
  // Compare number of valid tracks, and number of clusters associated in order to get global efficiency
  double nTracks = hTrackIntercepts->GetEntries();
  double nClusters = hGlobalAssociatedClusterPositions->GetEntries();
  double efficiency = nClusters/nTracks;
  double errorEfficiency = (1./nTracks) * sqrt(nClusters*(1-efficiency));
  
  if(nTracks == 0){ efficiency = 0.; errorEfficiency = 1.;}
  
  tcout<<"***** Clicpix efficiency calculation *****"<<endl;
  tcout<<"***** ntracks: "<<(int)nTracks<<", nclusters "<<(int)nClusters<<endl;
  tcout<<"***** Efficiency: "<<100.*efficiency<<" +/- "<<100.*errorEfficiency<<" %"<<endl;
372
  tcout<<"***** If including the "<<(int)m_lostHits<<" lost pixel hits, this becomes "<<100.*(m_lostHits+nClusters)/nTracks<<" %"<<endl;
Daniel Hynds's avatar
Daniel Hynds committed
373
374
  tcout<<endl;

375
}
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398

// Check if a track has gone through or near a masked pixel
bool ClicpixAnalysis::checkMasked(double chipInterceptRow, double chipInterceptCol){

  // Get the pixel row and column number
  int rowNumber = ceil(chipInterceptRow);
  int colNumber = ceil(chipInterceptCol);
  
  // Check if that pixel is masked, or the neighbour to the left, right, etc...
  if(parameters->detector[dutID]->masked(colNumber-1,rowNumber-1)) return true;
  if(parameters->detector[dutID]->masked(colNumber,rowNumber-1)) return true;
  if(parameters->detector[dutID]->masked(colNumber+1,rowNumber-1)) return true;
  if(parameters->detector[dutID]->masked(colNumber-1,rowNumber)) return true;
  if(parameters->detector[dutID]->masked(colNumber,rowNumber)) return true;
  if(parameters->detector[dutID]->masked(colNumber+1,rowNumber)) return true;
  if(parameters->detector[dutID]->masked(colNumber-1,rowNumber+1)) return true;
  if(parameters->detector[dutID]->masked(colNumber,rowNumber+1)) return true;
  if(parameters->detector[dutID]->masked(colNumber+1,rowNumber+1)) return true;

  // If not masked
  return false;
}

399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
// Check if there is another track close to the selected track.
// "Close" is defined as the intercept at the clicpix
bool ClicpixAnalysis::checkProximity(Track* track, Tracks* tracks){
  
  // Get the intercept of the interested track at the dut
  bool close = false;
  PositionVector3D< Cartesian3D<double> > trackIntercept = parameters->detector[dutID]->getIntercept(track);
  
  // Loop over all other tracks and check if they intercept close to the track we are considering
  Tracks::iterator itTrack;
  for (itTrack = tracks->begin(); itTrack != tracks->end(); itTrack++) {
    
    // Get the track
    Track* track2 = (*itTrack);
    // Get the track intercept with the clicpix plane (global co-ordinates)
    PositionVector3D< Cartesian3D<double> > trackIntercept2 = parameters->detector[dutID]->getIntercept(track2);
    // If track == track2 do nothing
    if(trackIntercept.X() == trackIntercept2.X() && trackIntercept.Y() == trackIntercept2.Y()) continue;
    if( fabs(trackIntercept.X() - trackIntercept2.X()) <= m_proximityCut ||
       fabs(trackIntercept.Y() - trackIntercept2.Y()) <= m_proximityCut ) close = true;
    
  }
  return close;
}

Daniel Hynds's avatar
Daniel Hynds committed
424
425
426
427
428
429
430
431
432
433
434
435
436
// Small sub-routine to fill histograms that only need clusters
void ClicpixAnalysis::fillClusterHistos(Clusters* clusters){
  
	// Pick up column to generate unique pixel id
  int nCols = parameters->detector[dutID]->nPixelsX();
  Clusters::iterator itc;
  
  // Check if this is a new clicpix frame (each frame may be in several events) and
  // fill histograms
  bool newFrame = false;
  for (itc = clusters->begin(); itc != clusters->end(); ++itc) {

    // Loop over pixels and check if there are pixels not known
437
    Pixels* pixels = (*itc)->pixels();
Daniel Hynds's avatar
Daniel Hynds committed
438
    Pixels::iterator itp;
439
    for (itp = pixels->begin(); itp != pixels->end(); itp++) {
Daniel Hynds's avatar
Daniel Hynds committed
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
      // Check if this clicpix frame is still the current
      int pixelID = (*itp)->m_column + nCols*(*itp)->m_row;
      if( m_hitPixels[pixelID] != (*itp)->m_adc){
        // New frame! Reset the stored pixels and trigger number
        if(!newFrame){m_hitPixels.clear(); newFrame=true;}
        m_hitPixels[pixelID] = (*itp)->m_adc;
        m_triggerNumber=0;
      }
      hHitPixels->Fill((*itp)->m_column,(*itp)->m_row);
      hColumnHits->Fill((*itp)->m_column);
      hRowHits->Fill((*itp)->m_row);
    }
    
    // Fill cluster histograms
    hClusterSizeAll->Fill((*itc)->size());
    hClusterTOTAll->Fill((*itc)->tot());
    hGlobalClusterPositions->Fill((*itc)->globalX(),(*itc)->globalY());
    
  }
  
  hClustersPerEvent->Fill(clusters->size());
  hClustersVersusEventNo->Fill(m_eventNumber,clusters->size());
  
  return;
}
465

466
467
468
469
// Sub-routine to look at pixel response, ie. how far from the pixel is the track intercept for the pixel to still see charge
void ClicpixAnalysis::fillResponseHistos(double trackInterceptX, double trackInterceptY, Cluster* cluster){
  
  // Loop over pixels in the cluster and show their distance from the track intercept
470
  Pixels* pixels = cluster->pixels();
471
  Pixels::iterator itp;
472
  for (itp = pixels->begin(); itp != pixels->end(); itp++) {
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
    
    // Get the pixel
    Pixel* pixel = (*itp);
    // Get the pixel local then global position
    PositionVector3D<Cartesian3D<double> > pixelPositionLocal = parameters->detector[dutID]->getLocalPosition(pixel->m_row,pixel->m_column);
    PositionVector3D<Cartesian3D<double> > pixelPositionGlobal = *(parameters->detector[dutID]->m_localToGlobal) * pixelPositionLocal;

    // Fill the response histograms
    hPixelResponseX->Fill( pixelPositionGlobal.X() - trackInterceptX );
    hPixelResponseY->Fill( pixelPositionGlobal.Y() - trackInterceptY );
    
  }
  
  return;
}

489