DUTAnalysis.cpp 15 KB
Newer Older
1
#include "DUTAnalysis.h"
2
#include "objects/Cluster.h"
3
#include "objects/MCParticle.h"
Daniel Hynds's avatar
Daniel Hynds committed
4
#include "objects/Pixel.h"
5
#include "objects/SpidrSignal.h"
Simon Spannagel's avatar
Simon Spannagel committed
6
#include "objects/Track.h"
7

8
9
using namespace corryvreckan;

10
DUTAnalysis::DUTAnalysis(Configuration config, std::vector<Detector*> detectors)
11
    : Module(std::move(config), std::move(detectors)) {
12
13
    m_digitalPowerPulsing = m_config.get<bool>("digitalPowerPulsing", false);
    m_DUT = m_config.get<std::string>("DUT");
14
    m_useMCtruth = m_config.get<bool>("useMCtruth", false);
15
    timingCut = m_config.get<double>("timingCut", Units::convert(200, "ns"));
16
    spatialCut = m_config.get<double>("spatialCut", Units::convert(200, "um"));
17
    chi2ndofCut = m_config.get<double>("chi2ndofCut", 3.);
18
19
}

20
void DUTAnalysis::initialise() {
Simon Spannagel's avatar
Simon Spannagel committed
21
22
23
24
25

    // Initialise single histograms
    tracksVersusTime = new TH1F("tracksVersusTime", "tracksVersusTime", 300000, 0, 300);
    associatedTracksVersusTime = new TH1F("associatedTracksVersusTime", "associatedTracksVersusTime", 300000, 0, 300);
    residualsX = new TH1F("residualsX", "residualsX", 400, -0.2, 0.2);
26
    residualsXfine = new TH1F("residualsXfine", "residualsXfine", 400, -0.05, 0.05);
27
    residualsX1pix = new TH1F("residualsX1pix", "residualsX1pix", 400, -0.2, 0.2);
28
    residualsX2pix = new TH1F("residualsX2pix", "residualsX2pix", 400, -0.2, 0.2);
Simon Spannagel's avatar
Simon Spannagel committed
29
    residualsY = new TH1F("residualsY", "residualsY", 400, -0.2, 0.2);
Simon Spannagel's avatar
Simon Spannagel committed
30

31
    clusterTotAssociated = new TH1F("clusterTotAssociated", "clusterTotAssociated", 10000, 0, 10000);
Simon Spannagel's avatar
Simon Spannagel committed
32
    clusterSizeAssociated = new TH1F("clusterSizeAssociated", "clusterSizeAssociated", 30, 0, 30);
33
34
    clusterSizeAssociated_X = new TH1F("clusterSizeAssociated_X", "clusterSizeAssociated_X", 30, 0, 30);
    clusterSizeAssociated_Y = new TH1F("clusterSizeAssociated_Y", "clusterSizeAssociated_Y", 30, 0, 30);
fpipper's avatar
fpipper committed
35
    residualsTime = new TH1F("residualsTime", "residualsTime", 20000, -1000, +1000);
Simon Spannagel's avatar
Simon Spannagel committed
36
37
38

    hTrackCorrelationX = new TH1F("hTrackCorrelationX", "hTrackCorrelationX", 4000, -10., 10.);
    hTrackCorrelationY = new TH1F("hTrackCorrelationY", "hTrackCorrelationY", 4000, -10., 10.);
fpipper's avatar
fpipper committed
39
    hTrackCorrelationTime = new TH1F("hTrackCorrelationTime", "hTrackCorrelationTime", 2000000, -5000, 5000);
Simon Spannagel's avatar
Simon Spannagel committed
40
41
    clusterToTVersusTime = new TH2F("clusterToTVersusTime", "clusterToTVersusTime", 300000, 0., 300., 200, 0, 1000);

fpipper's avatar
fpipper committed
42
43
    residualsTimeVsTime = new TH2F("residualsTimeVsTime", "residualsTimeVsTime", 20000, 0, 200, 1000, -1000, +1000);
    residualsTimeVsSignal = new TH2F("residualsTimeVsSignal", "residualsTimeVsSignal", 20000, 0, 100000, 1000, -1000, +1000);
Simon Spannagel's avatar
Simon Spannagel committed
44
45
46
47
48
49
50
51
52
53

    tracksVersusPowerOnTime = new TH1F("tracksVersusPowerOnTime", "tracksVersusPowerOnTime", 1200000, -0.01, 0.11);
    associatedTracksVersusPowerOnTime =
        new TH1F("associatedTracksVersusPowerOnTime", "associatedTracksVersusPowerOnTime", 1200000, -0.01, 0.11);

    hAssociatedTracksGlobalPosition =
        new TH2F("hAssociatedTracksGlobalPosition", "hAssociatedTracksGlobalPosition", 200, -10, 10, 200, -10, 10);
    hUnassociatedTracksGlobalPosition =
        new TH2F("hUnassociatedTracksGlobalPosition", "hUnassociatedTracksGlobalPosition", 200, -10, 10, 200, -10, 10);

Daniel Hynds's avatar
Daniel Hynds committed
54
55
    if(m_useMCtruth) {
        residualsXMCtruth = new TH1F("residualsXMCtruth", "residualsXMCtruth", 400, -0.2, 0.2);
56
        telescopeResolution = new TH1F("telescopeResolution", "telescopeResolution", 400, -0.2, 0.2);
Daniel Hynds's avatar
Daniel Hynds committed
57
58
    }

Simon Spannagel's avatar
Simon Spannagel committed
59
60
61
62
63
64
65
    // Initialise member variables
    m_eventNumber = 0;
    m_nAlignmentClusters = 0;
    m_powerOnTime = 0;
    m_powerOffTime = 0;
    m_shutterOpenTime = 0;
    m_shutterCloseTime = 0;
66
67
}

Simon Spannagel's avatar
Simon Spannagel committed
68
69
StatusCode DUTAnalysis::run(Clipboard* clipboard) {

70
71
    LOG(TRACE) << "Power on time: " << Units::display(m_powerOnTime, {"ns", "us", "s"});
    LOG(TRACE) << "Power off time: " << Units::display(m_powerOffTime, {"ns", "us", "s"});
Simon Spannagel's avatar
Simon Spannagel committed
72

Simon Spannagel's avatar
Simon Spannagel committed
73
    //    if(clipboard->get_persistent("eventStart") < 13.5)
Simon Spannagel's avatar
Simon Spannagel committed
74
    //        return Success;
Simon Spannagel's avatar
Simon Spannagel committed
75
76
77
78
79
80
81
82
83
84
85
86
87
88

    // Track chi2/ndof cut
    // 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
    if(m_shutterCloseTime != 0 && m_shutterCloseTime > m_shutterOpenTime)
        m_shutterOpenTime = 0;
    if(m_shutterOpenTime != 0 && m_shutterOpenTime > m_shutterCloseTime)
        m_shutterCloseTime = 0;

    // Now update the power pulsing with any new signals
89
    SpidrSignals* spidrData = (SpidrSignals*)clipboard->get(m_DUT, "SpidrSignals");
Simon Spannagel's avatar
Simon Spannagel committed
90
91
92
93
94
95
96
97
98
99
100
101
    // 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];
            // 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
102
                if(std::abs(Units::convert(signal->timestamp() - m_shutterOpenTime, "s")) < 0.001) {
Simon Spannagel's avatar
Simon Spannagel committed
103
                    continue;
Simon Spannagel's avatar
Simon Spannagel committed
104
                }
Simon Spannagel's avatar
Simon Spannagel committed
105
                m_shutterOpenTime = signal->timestamp();
Simon Spannagel's avatar
Simon Spannagel committed
106
                LOG(TRACE) << "Shutter opened at " << Units::display(m_shutterOpenTime, {"ns", "us", "s"});
Simon Spannagel's avatar
Simon Spannagel committed
107
108
109
110
            }
            if(signal->type() == "shutterClosed") {
                // There may be multiple power on/off in 1 time window. At the moment,
                // take earliest if within 1ms
111
                if(std::abs(Units::convert(signal->timestamp() - m_shutterCloseTime, "s")) < 0.001) {
Simon Spannagel's avatar
Simon Spannagel committed
112
                    continue;
Simon Spannagel's avatar
Simon Spannagel committed
113
                }
Simon Spannagel's avatar
Simon Spannagel committed
114
                m_shutterCloseTime = signal->timestamp();
Simon Spannagel's avatar
Simon Spannagel committed
115
                LOG(TRACE) << "Shutter closed at " << Units::display(m_shutterCloseTime, {"ns", "us", "s"});
Simon Spannagel's avatar
Simon Spannagel committed
116
117
            }
        }
Daniel Hynds's avatar
Daniel Hynds committed
118
    }
119

Simon Spannagel's avatar
Simon Spannagel committed
120
121
122
123
124
    // Get the tracks from the clipboard
    Tracks* tracks = (Tracks*)clipboard->get("tracks");
    if(tracks == NULL) {
        LOG(DEBUG) << "No tracks on the clipboard";
        return Success;
Daniel Hynds's avatar
Daniel Hynds committed
125
    }
Daniel Hynds's avatar
Daniel Hynds committed
126

Simon Spannagel's avatar
Simon Spannagel committed
127
    // Get the DUT clusters from the clipboard
128
    Clusters* clusters = (Clusters*)clipboard->get(m_DUT, "clusters");
Simon Spannagel's avatar
Simon Spannagel committed
129
130
    if(clusters == NULL) {
        LOG(DEBUG) << "No DUT clusters on the clipboard";
Daniel Hynds's avatar
Daniel Hynds committed
131
    }
132

133
    // Get the MC particles from the clipboard
Daniel Hynds's avatar
Daniel Hynds committed
134
135
136
137
    MCParticles* mcParticles = (MCParticles*)clipboard->get(m_DUT, "mcparticles");
    if(mcParticles == NULL) {
        LOG(DEBUG) << "No DUT MC particles on the clipboard";
    }
138

Simon Spannagel's avatar
Simon Spannagel committed
139
    // Loop over all tracks
140
141
    bool first_track = true;
    for(auto& track : (*tracks)) {
Simon Spannagel's avatar
Simon Spannagel committed
142
143

        // Cut on the chi2/ndof
144
        if(track->chi2ndof() > chi2ndofCut) {
Simon Spannagel's avatar
Simon Spannagel committed
145
            continue;
146
        }
Simon Spannagel's avatar
Simon Spannagel committed
147

148
        // Get the DUT detector:
149
        auto detector = get_detector(m_DUT);
Simon Spannagel's avatar
Simon Spannagel committed
150
        // Check if it intercepts the DUT
151
152
        PositionVector3D<Cartesian3D<double>> globalIntercept = detector->getIntercept(track);
        if(!detector->hasIntercept(track, 1.)) {
Simon Spannagel's avatar
Simon Spannagel committed
153
            continue;
154
        }
Simon Spannagel's avatar
Simon Spannagel committed
155
156

        // Check that it doesn't go through/near a masked pixel
157
        if(detector->hitMasked(track, 1.)) {
Simon Spannagel's avatar
Simon Spannagel committed
158
            continue;
159
        }
Simon Spannagel's avatar
Simon Spannagel committed
160

161
        tracksVersusTime->Fill(Units::convert(track->timestamp(), "s"));
Simon Spannagel's avatar
Simon Spannagel committed
162

163
164
165
166
167
168
        timeSincePowerOn = track->timestamp() - m_shutterOpenTime;
        if(timeSincePowerOn > 0. && timeSincePowerOn < Units::convert(200, "us")) {
            LOG(TRACE) << "Track at time " << Units::display(track->timestamp(), {"ns", "us", "s"})
                       << " has time shutter open of " << Units::display(timeSincePowerOn, {"ns", "us", "s"});
            LOG(TRACE) << "Shutter open time is " << Units::display(m_shutterOpenTime, {"ns", "us", "s"})
                       << ", shutter close time is " << Units::display(m_shutterCloseTime, {"ns", "us", "s"});
Simon Spannagel's avatar
Simon Spannagel committed
169
170
171
172
173
174
175
176
177
        }

        // 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
        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)))) {
178
            timeSincePowerOn = track->timestamp() - m_shutterOpenTime;
Simon Spannagel's avatar
Simon Spannagel committed
179
            tracksVersusPowerOnTime->Fill(timeSincePowerOn);
180
181

            if(timeSincePowerOn < 200000) {
Simon Spannagel's avatar
Simon Spannagel committed
182
                LOG(TRACE) << "Track at time " << Units::display(clipboard->get_persistent("eventStart"), {"ns", "us", "s"})
Simon Spannagel's avatar
Simon Spannagel committed
183
184
185
                           << " has time shutter open of " << Units::display(timeSincePowerOn, {"ns", "us", "s"});
                LOG(TRACE) << "Shutter open time is " << Units::display(m_shutterOpenTime, {"ns", "us", "s"})
                           << ", shutter close time is " << Units::display(m_shutterCloseTime, {"ns", "us", "s"});
186
            }
Simon Spannagel's avatar
Simon Spannagel committed
187
188
189
190
191
        }

        // If no DUT clusters then continue to the next track
        if(clusters == NULL)
            continue;
fpipper's avatar
fpipper committed
192

Simon Spannagel's avatar
Simon Spannagel committed
193
        // Correlation plot
fpipper's avatar
fpipper committed
194
195
196
197
198
199
        for(int itCluster = 0; itCluster < clusters->size(); itCluster++) {

            // Get the cluster pointer
            Cluster* cluster = (*clusters)[itCluster];

            // Check if the cluster is close in time
200
            // if( std::abs(cluster->timestamp() - track->timestamp()) > timingCut)
201
            //    continue;
fpipper's avatar
fpipper committed
202
203
204
205
206
207
208

            // 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());
fpipper's avatar
fpipper committed
209
            hTrackCorrelationTime->Fill(Units::convert(track->timestamp() - cluster->timestamp(), "ns"));
fpipper's avatar
fpipper committed
210
211

            if(fabs(intercept.X() - cluster->globalX()) < 0.1 && fabs(intercept.Y() - cluster->globalY()) < 0.1) {
fpipper's avatar
fpipper committed
212
213
214
215
                residualsTime->Fill(Units::convert(track->timestamp() - cluster->timestamp(), "ns"));
                residualsTimeVsTime->Fill(Units::convert(track->timestamp(), "ns"),
                                          Units::convert(track->timestamp() - cluster->timestamp(), "ns"));
                residualsTimeVsSignal->Fill(cluster->tot(), Units::convert(track->timestamp() - cluster->timestamp(), "ns"));
fpipper's avatar
fpipper committed
216
            }
Simon Spannagel's avatar
Simon Spannagel committed
217
218
219
220
        }

        // Loop over all DUT clusters
        bool associated = false;
221
        for(auto& cluster : (*clusters)) {
Simon Spannagel's avatar
Simon Spannagel committed
222
223

            // Fill the tot histograms on the first run
224
            if(first_track == 0) {
fpipper's avatar
fpipper committed
225
                clusterToTVersusTime->Fill(Units::convert(cluster->timestamp(), "ns"), cluster->tot());
Simon Spannagel's avatar
Simon Spannagel committed
226
227
            }

228
            // Check distance between track and cluster
Simon Spannagel's avatar
Simon Spannagel committed
229
230
231
            ROOT::Math::XYZPoint intercept = track->intercept(cluster->globalZ());
            double xdistance = intercept.X() - cluster->globalX();
            double ydistance = intercept.Y() - cluster->globalY();
232
233
234
235
236
237
238
            if(abs(xdistance) > spatialCut || abs(ydistance) > spatialCut) {
                LOG(DEBUG) << "Discarding DUT cluster with distance (" << abs(xdistance) << "," << abs(ydistance) << ")";
                continue;
            }

            LOG(DEBUG) << "Found associated cluster with distance (" << abs(xdistance) << "," << abs(ydistance) << ")";
            track->addAssociatedCluster(cluster);
Simon Spannagel's avatar
Simon Spannagel committed
239

Simon Spannagel's avatar
Simon Spannagel committed
240
241
242
            // We now have an associated cluster! Fill plots
            associated = true;
            LOG(TRACE) << "Found associated cluster";
243
            associatedTracksVersusTime->Fill(Units::convert(track->timestamp(), "s"));
Simon Spannagel's avatar
Simon Spannagel committed
244
            residualsX->Fill(xdistance);
245
            residualsXfine->Fill(xdistance);
246
247
            if(cluster->size() == 1)
                residualsX1pix->Fill(xdistance);
248
249
            if(cluster->size() == 2)
                residualsX2pix->Fill(xdistance);
Simon Spannagel's avatar
Simon Spannagel committed
250
            residualsY->Fill(ydistance);
251
252
            clusterTotAssociated->Fill(cluster->tot());
            clusterSizeAssociated->Fill(cluster->size());
253
254
            clusterSizeAssociated_X->Fill(cluster->columnWidth());
            clusterSizeAssociated_Y->Fill(cluster->rowWidth());
Simon Spannagel's avatar
Simon Spannagel committed
255
256
257
258
            track->addAssociatedCluster(cluster);
            m_nAlignmentClusters++;
            hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());

259
            // Get associated MC particle info and plot
Daniel Hynds's avatar
Daniel Hynds committed
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
            if(m_useMCtruth) {
                if(mcParticles == nullptr)
                    continue;
                // Find the closest MC particle
                double smallestDistance(DBL_MAX);
                ROOT::Math::XYZPoint particlePosition;
                for(auto& particle : (*mcParticles)) {
                    ROOT::Math::XYZPoint entry = particle->getLocalStart();
                    ROOT::Math::XYZPoint exit = particle->getLocalEnd();
                    ROOT::Math::XYZPoint centre(
                        (entry.X() + exit.X()) / 2., (entry.Y() + exit.Y()) / 2., (entry.Z() + exit.Z()) / 2.);
                    double distance = sqrt((centre.X() - cluster->localX()) * (centre.X() - cluster->localX()) +
                                           (centre.Y() - cluster->localY()) * (centre.Y() - cluster->localY()));
                    if(distance < smallestDistance) {
                        particlePosition.SetXYZ(centre.X(), centre.Y(), centre.Z());
                    }
276
                }
277
278
                auto detector = get_detector(cluster->detectorID());
                residualsXMCtruth->Fill(cluster->localX() + detector->size().X() / 2 - particlePosition.X());
279
                auto interceptLocal = detector->globalToLocal(intercept);
280
                telescopeResolution->Fill(interceptLocal.X() + detector->size().X() / 2 - particlePosition.X());
Daniel Hynds's avatar
Daniel Hynds committed
281
282
            }

Simon Spannagel's avatar
Simon Spannagel committed
283
284
285
286
287
288
289
290
291
292
293
294
295
296
            // Fill power pulsing response
            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)))) {
                associatedTracksVersusPowerOnTime->Fill(timeSincePowerOn);
            }

            // Only allow one associated cluster per track
            break;
        }

        if(!associated)
            hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
297
298

        first_track = false;
299
    }
300

Simon Spannagel's avatar
Simon Spannagel committed
301
302
    // Increment event counter
    m_eventNumber++;
303

Simon Spannagel's avatar
Simon Spannagel committed
304
305
306
    //  if(m_nAlignmentClusters > 10000) return Failure;
    // Return value telling analysis to keep running
    return Success;
307
308
}

Simon Spannagel's avatar
Simon Spannagel committed
309
void DUTAnalysis::finalise() {
310

Simon Spannagel's avatar
Simon Spannagel committed
311
    LOG(DEBUG) << "Analysed " << m_eventNumber << " events";
312
}
313
314

// Function to check if a track goes through a given device
Simon Spannagel's avatar
Simon Spannagel committed
315
// bool DUTAnalysis::intercept(Track*, string device){
316
//
317
318
319
320
//  // Get the global intercept of the track and the device
//  ROOT::Math::XYZPoint intercept = track->intercept(cluster->globalZ());
//
//  // Transform to the local co-ordinates
321
//
322
//  // Check if the row/column number is outside the acceptable range
323
324
//
//
325
//}