DUTAnalysis.C 10.6 KB
Newer Older
1
#include "DUTAnalysis.h"
2
3
4
#include "Pixel.h"
#include "Cluster.h"
#include "Track.h"
Daniel Hynds's avatar
Daniel Hynds committed
5
#include "SpidrSignal.h"
6
7
8
9

DUTAnalysis::DUTAnalysis(bool debugging)
: Algorithm("DUTAnalysis"){
  debug = debugging;
Daniel Hynds's avatar
Daniel Hynds committed
10
  m_digitalPowerPulsing = false;
11
12
13
14
15
}


void DUTAnalysis::initialise(Parameters* par){
 
16
  // Pick up a copy of the parameters
17
18
19
  parameters = par;

  // Initialise single histograms
Daniel Hynds's avatar
Daniel Hynds committed
20
21
  tracksVersusTime = new TH1F("tracksVersusTime","tracksVersusTime",300000,0,300);
  associatedTracksVersusTime = new TH1F("associatedTracksVersusTime","associatedTracksVersusTime",300000,0,300);
22
23
  residualsX = new TH1F("residualsX","residualsX",400,-0.2,0.2);
  residualsY = new TH1F("residualsY","residualsY",400,-0.2,0.2);
Daniel Hynds's avatar
Daniel Hynds committed
24
  residualsTime = new TH1F("residualsTime","residualsTime",2000,-0.000001,0.000001);
25
  
Daniel Hynds's avatar
Daniel Hynds committed
26
27
28
29
30
  hTrackCorrelationX = new TH1F("hTrackCorrelationX","hTrackCorrelationX",4000,-10.,10.);
  hTrackCorrelationY = new TH1F("hTrackCorrelationY","hTrackCorrelationY",4000,-10.,10.);
  hTrackCorrelationTime = new TH1F("hTrackCorrelationTime","hTrackCorrelationTime",2000000,-0.005,0.005);
  clusterToTVersusTime = new TH2F("clusterToTVersusTime","clusterToTVersusTime",300000,0.,300.,200,0,1000);
  
31
  residualsTimeVsTime = new TH2F("residualsTimeVsTime","residualsTimeVsTime",20000,0,200,400,-0.0005,0.0005);
Daniel Hynds's avatar
Daniel Hynds committed
32

Daniel Hynds's avatar
Daniel Hynds committed
33
34
35
  tracksVersusPowerOnTime = new TH1F("tracksVersusPowerOnTime","tracksVersusPowerOnTime",1200000,-0.01,0.11);
  associatedTracksVersusPowerOnTime = new TH1F("associatedTracksVersusPowerOnTime","associatedTracksVersusPowerOnTime",1200000,-0.01,0.11);

Daniel Hynds's avatar
Daniel Hynds committed
36
37
38
  hAssociatedTracksGlobalPosition = new TH2F("hAssociatedTracksGlobalPosition","hAssociatedTracksGlobalPosition",200,-10,10,200,-10,10);
  hUnassociatedTracksGlobalPosition = new TH2F("hUnassociatedTracksGlobalPosition","hUnassociatedTracksGlobalPosition",200,-10,10,200,-10,10);

39
40
41
  // Initialise member variables
  m_eventNumber = 0;
  m_nAlignmentClusters = 0;
Daniel Hynds's avatar
Daniel Hynds committed
42
43
  m_powerOnTime = 0;
  m_powerOffTime = 0;
Daniel Hynds's avatar
Daniel Hynds committed
44
45
  m_shutterOpenTime = 0;
  m_shutterCloseTime = 0;
Daniel Hynds's avatar
Daniel Hynds committed
46

47
48
}

Daniel Hynds's avatar
Daniel Hynds committed
49
50
51
52
53
StatusCode DUTAnalysis::run(Clipboard* clipboard){
  
//  tcout<<"Power on time: "<<m_powerOnTime/(4096. * 40000000.)<<endl;
//  tcout<<"Power off time: "<<m_powerOffTime/(4096. * 40000000.)<<endl;
//  tcout<<endl;
54
  
Daniel Hynds's avatar
Daniel Hynds committed
55
56
57
58
  cout<<std::setprecision(10);
  
  if(parameters->currentTime < 13.5) return Success;
  
59
  // Timing cut for association
60
  double timingCut = 200./1000000000.; // 200 ns
61
62
63
  long long int timingCutInt = (timingCut * 4096. * 40000000.);
  
  // Spatial cut
64
65
66
  double spatialCut = 0.2; // 200 um

  // Track chi2/ndof cut
Daniel Hynds's avatar
Daniel Hynds committed
67
  double chi2ndofCut = 3.;
68

Daniel Hynds's avatar
Daniel Hynds committed
69
70
71
72
  // Power pulsing variable initialisation - get signals from SPIDR for this device
  double timeSincePowerOn = 0.;
  
  // If the power was switched off/on in the last event we no longer have a power on/off time
Daniel Hynds's avatar
Daniel Hynds committed
73
74
  if(m_shutterCloseTime != 0 && m_shutterCloseTime > m_shutterOpenTime) m_shutterOpenTime = 0;
  if(m_shutterOpenTime != 0 && m_shutterOpenTime > m_shutterCloseTime) m_shutterCloseTime = 0;
Daniel Hynds's avatar
Daniel Hynds committed
75
76
77
78
79
80
81
82
83
84
  
  // Now update the power pulsing with any new signals
  SpidrSignals* spidrData = (SpidrSignals*)clipboard->get(parameters->DUT,"SpidrSignals");
  // If there are new signals
  if(spidrData != NULL){
    // Loop over all signals registered
    int nSignals = spidrData->size();
    for(int iSig=0;iSig<nSignals;iSig++){
      // Get the signal
      SpidrSignal* signal = (*spidrData)[iSig];
Daniel Hynds's avatar
Daniel Hynds committed
85
86
87
88
89
90
    	// Register the power on or power off time, and whether the shutter is open or not
      if(signal->type() == "shutterOpen"){
        // There may be multiple power on/off in 1 time window. At the moment, take earliest if within 1ms
        if( fabs( double(signal->timestamp()-m_shutterOpenTime)/(4096. * 40000000.)) < 0.001 ) continue;
        m_shutterOpenTime = signal->timestamp();
//        tcout<<"Shutter opened at "<<double(m_shutterOpenTime)/(4096.*40000000.)<<endl;
Daniel Hynds's avatar
Daniel Hynds committed
91
      }
Daniel Hynds's avatar
Daniel Hynds committed
92
93
94
95
96
97
      if(signal->type() == "shutterClosed"){
        // There may be multiple power on/off in 1 time window. At the moment, take earliest if within 1ms
        if( fabs( double(signal->timestamp()-m_shutterCloseTime)/(4096. * 40000000.)) < 0.001 ) continue;
        m_shutterCloseTime = signal->timestamp();
//        cout<<endl;
//        tcout<<"Shutter closed at "<<double(m_shutterCloseTime)/(4096.*40000000.)<<endl;
Daniel Hynds's avatar
Daniel Hynds committed
98
99
100
101
      }
    }
  }
  
102
  // Get the tracks from the clipboard
103
  Tracks* tracks = (Tracks*)clipboard->get("tracks");
104
105
  if(tracks == NULL){
    if(debug) tcout<<"No tracks on the clipboard"<<endl;
Daniel Hynds's avatar
Daniel Hynds committed
106
    return Success;
107
  }
Daniel Hynds's avatar
Daniel Hynds committed
108
  
109
  // Get the DUT clusters from the clipboard
110
  Clusters* clusters = (Clusters*)clipboard->get(parameters->DUT,"clusters");
111
112
113
114
115
116
117
118
  if(clusters == NULL){
    if(debug) tcout<<"No DUT clusters on the clipboard"<<endl;
  }

  // Loop over all tracks
  for(int itTrack=0;itTrack<tracks->size();itTrack++){
    
    // Get the track pointer
119
    Track* track = (*tracks)[itTrack];
120
121
122
    
    // Cut on the chi2/ndof
    if(track->chi2ndof() > chi2ndofCut) continue;
Daniel Hynds's avatar
Daniel Hynds committed
123
124
125
126
127
128
129
130
    
    // Check if it intercepts the DUT
    PositionVector3D<Cartesian3D<double> > globalIntercept = parameters->detector[parameters->DUT]->getIntercept(track);
    if(!parameters->detector[parameters->DUT]->hasIntercept(track,1.)) continue;
    
    // Check that it doesn't go through/near a masked pixel
    if(parameters->detector[parameters->DUT]->hitMasked(track,1.)) continue;

131
132
    tracksVersusTime->Fill( (double)track->timestamp() / (4096.*40000000.) );
    
Daniel Hynds's avatar
Daniel Hynds committed
133
134
135
136
137
138
139
140
    
    timeSincePowerOn = (double)(track->timestamp() - m_shutterOpenTime) / (4096.*40000000.);
    if(timeSincePowerOn > 0. && timeSincePowerOn < 0.0002){
//      cout<<endl;
//      tcout<<"Track at time "<<double(track->timestamp())/(4096.*40000000.)<<" has time shutter open of "<<timeSincePowerOn<<endl;
//      tcout<<"Shutter open time is "<<double(m_shutterOpenTime)/(4096.*40000000.)<<", shutter close time is "<<double(m_shutterCloseTime)/(4096.*40000000.)<<endl;
    }

Daniel Hynds's avatar
Daniel Hynds committed
141
142
    // Check time since power on (if power pulsing).
    // If power off time not known it will be 0. If it is known, then the track should arrive before the power goes off
Daniel Hynds's avatar
Daniel Hynds committed
143
144
145
146
    if( (m_shutterOpenTime != 0 && m_shutterCloseTime == 0) ||
        (m_shutterOpenTime != 0 && ( (m_shutterCloseTime > m_shutterOpenTime && m_shutterCloseTime - track->timestamp() > 0) ||
                                     (m_shutterOpenTime > m_shutterCloseTime && track->timestamp() - m_shutterOpenTime >= 0) ))) {
      timeSincePowerOn = (double)(track->timestamp() - m_shutterOpenTime) / (4096.*40000000.);
Daniel Hynds's avatar
Daniel Hynds committed
147
      tracksVersusPowerOnTime->Fill(timeSincePowerOn);
Daniel Hynds's avatar
Daniel Hynds committed
148
149
150
151
152
//      if(timeSincePowerOn < (0.0002)){
//        cout<<endl;
//        tcout<<"Track at time "<<parameters->currentTime<<" has time shutter open of "<<timeSincePowerOn<<endl;
//        tcout<<"Shutter open time is "<<double(m_shutterOpenTime)/(4096.*40000000.)<<", shutter close time is "<<double(m_shutterCloseTime)/(4096.*40000000.)<<endl;
//      }
Daniel Hynds's avatar
Daniel Hynds committed
153
    }
Daniel Hynds's avatar
Daniel Hynds committed
154

Daniel Hynds's avatar
Daniel Hynds committed
155
156
157
    // If no DUT clusters then continue to the next track
    if(clusters == NULL) continue;
    /*
Daniel Hynds's avatar
Daniel Hynds committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
    // Correlation plot
    for(int itCluster=0;itCluster<clusters->size();itCluster++){
      
      // Get the cluster pointer
      Cluster* cluster = (*clusters)[itCluster];
      
      // Check if the cluster is close in time
//      if( abs(cluster->timestamp() - track->timestamp()) > timingCutInt ) continue;

      // Check distance between track and cluster
      ROOT::Math::XYZPoint intercept = track->intercept(cluster->globalZ());

      // Fill the correlation plot
      hTrackCorrelationX->Fill(intercept.X() - cluster->globalX());
      hTrackCorrelationY->Fill(intercept.Y() - cluster->globalY());
      hTrackCorrelationTime->Fill( (double)(track->timestamp() - cluster->timestamp()) / (4096.*40000000.));
      
      if( fabs(intercept.X() - cluster->globalX()) < 0.1 &&
          fabs(intercept.Y() - cluster->globalY()) < 0.1){
        residualsTime->Fill((double)(track->timestamp() - cluster->timestamp()) / (4096.*40000000.));
        residualsTimeVsTime->Fill( (double)track->timestamp() / (4096.*40000000.), (double)(track->timestamp() - cluster->timestamp()) / (4096.*40000000.));
      }
    }
Daniel Hynds's avatar
Daniel Hynds committed
181
    */
Daniel Hynds's avatar
Daniel Hynds committed
182
    
183
    // Loop over all DUT clusters
Daniel Hynds's avatar
Daniel Hynds committed
184
    bool associated = false;
185
186
187
    for(int itCluster=0;itCluster<clusters->size();itCluster++){

      // Get the cluster pointer
188
189
190
191
      Cluster* cluster = (*clusters)[itCluster];
      
      // Fill the tot histograms on the first run
      if(itTrack == 0) clusterToTVersusTime->Fill((double)cluster->timestamp() / (4096.*40000000.), cluster->tot());
192
193
      
      // Check if the cluster is close in time
Daniel Hynds's avatar
Daniel Hynds committed
194
      if( !m_digitalPowerPulsing && abs(cluster->timestamp() - track->timestamp()) > timingCutInt ) continue;
195
196
197
198
199
200
201
202
203
      
      // Check distance between track and cluster
      ROOT::Math::XYZPoint intercept = track->intercept(cluster->globalZ());
      double xdistance = intercept.X() - cluster->globalX();
      double ydistance = intercept.Y() - cluster->globalY();
      if( abs(xdistance) > spatialCut) continue;
      if( abs(ydistance) > spatialCut) continue;
 
      // We now have an associated cluster! Fill plots
Daniel Hynds's avatar
Daniel Hynds committed
204
205
      associated = true;
//      tcout<<"Found associated cluster"<<endl;
206
207
208
209
210
      associatedTracksVersusTime->Fill( (double)track->timestamp() / (4096.*40000000.) );
      residualsX->Fill(xdistance);
      residualsY->Fill(ydistance);
      track->addAssociatedCluster(cluster);
      m_nAlignmentClusters++;
Daniel Hynds's avatar
Daniel Hynds committed
211
      hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(),globalIntercept.Y());
212
      
Daniel Hynds's avatar
Daniel Hynds committed
213
      // Fill power pulsing response
Daniel Hynds's avatar
Daniel Hynds committed
214
215
216
      if( (m_shutterOpenTime != 0 && m_shutterCloseTime == 0) ||
         (m_shutterOpenTime != 0 && ( (m_shutterCloseTime > m_shutterOpenTime && m_shutterCloseTime - track->timestamp() > 0) ||
                                     (m_shutterOpenTime > m_shutterCloseTime && track->timestamp() - m_shutterOpenTime >= 0) ))) {
Daniel Hynds's avatar
Daniel Hynds committed
217
218
219
        associatedTracksVersusPowerOnTime->Fill(timeSincePowerOn);
      }

220
221
222
223
224
      // Only allow one associated cluster per track
      break;
      
    }
    
Daniel Hynds's avatar
Daniel Hynds committed
225
226
    if(!associated) hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(),globalIntercept.Y());
    
227
228
229
230
231
  }

  // Increment event counter
  m_eventNumber++;
  
Daniel Hynds's avatar
Daniel Hynds committed
232
//  if(m_nAlignmentClusters > 10000) return Failure;
233
  // Return value telling analysis to keep running
Daniel Hynds's avatar
Daniel Hynds committed
234
  return Success;
235
236
237
238
239
240
241
}

void DUTAnalysis::finalise(){
  
  if(debug) tcout<<"Analysed "<<m_eventNumber<<" events"<<endl;
  
}
242
243

// Function to check if a track goes through a given device
244
//bool DUTAnalysis::intercept(Track*, string device){
245
246
247
248
249
250
251
252
253
254
//  
//  // Get the global intercept of the track and the device
//  ROOT::Math::XYZPoint intercept = track->intercept(cluster->globalZ());
//
//  // Transform to the local co-ordinates
//  
//  // Check if the row/column number is outside the acceptable range
//  
//  
//}