AnalysisEfficiency.cpp 7.44 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
/**
 * @file
 * @brief Implementation of [AnalysisEfficiency] module
 * @copyright Copyright (c) 2017 CERN and the Allpix Squared authors.
 * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
 * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
 * Intergovernmental Organization or submit itself to any jurisdiction.
 */

#include "AnalysisEfficiency.h"

12
#include "objects/Cluster.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
13
#include "objects/Pixel.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
14
#include "objects/Track.hpp"
15
16
17
18
19
20
21

using namespace corryvreckan;

AnalysisEfficiency::AnalysisEfficiency(Configuration config, std::shared_ptr<Detector> detector)
    : Module(std::move(config), detector) {
    m_detector = detector;

22
23
    m_timeCutFrameEdge = m_config.get<double>("time_cut_frameedge", static_cast<double>(Units::convert(20, "ns")));
    m_chi2ndofCut = m_config.get<double>("chi2ndof_cut", 3.);
24
25
26
27
28
29
30
}

void AnalysisEfficiency::initialise() {

    auto pitch_x = static_cast<double>(Units::convert(m_detector->pitch().X(), "um"));
    auto pitch_y = static_cast<double>(Units::convert(m_detector->pitch().Y(), "um"));

31
32
33
34
    std::string title = m_detector->name() + " Pixel efficiency map;x_{track} mod " + std::to_string(pitch_x) +
                        "#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;efficiency";
    hPixelEfficiencyMap = new TProfile2D("pixelEfficiencyMap",
                                         title.c_str(),
35
36
37
38
39
40
41
42
                                         static_cast<int>(pitch_x),
                                         0,
                                         pitch_x,
                                         static_cast<int>(pitch_y),
                                         0,
                                         pitch_y,
                                         0,
                                         1);
43
44
45
    title = m_detector->name() + " Chip efficiency map;x [px];y [px];efficiency";
    hChipEfficiencyMap = new TProfile2D("chipEfficiencyMap",
                                        title.c_str(),
46
                                        m_detector->nPixels().X(),
47
                                        0,
48
49
                                        m_detector->nPixels().X(),
                                        m_detector->nPixels().Y(),
50
                                        0,
51
                                        m_detector->nPixels().Y(),
52
53
                                        0,
                                        1);
54
55
56
    title = m_detector->name() + " Global efficiency map;x [mm];y [mm];efficiency";
    hGlobalEfficiencyMap = new TProfile2D("globalEfficiencyMap",
                                          title.c_str(),
57
58
59
60
61
62
63
64
65
66
                                          300,
                                          -1.5 * m_detector->size().X(),
                                          1.5 * m_detector->size().X(),
                                          300,
                                          -1.5 * m_detector->size().Y(),
                                          1.5 * m_detector->size().Y(),
                                          0,
                                          1);
}

67
StatusCode AnalysisEfficiency::run(std::shared_ptr<Clipboard> clipboard) {
68
69
70
71
72

    // Get the telescope tracks from the clipboard
    Tracks* tracks = reinterpret_cast<Tracks*>(clipboard->get("tracks"));
    if(tracks == nullptr) {
        LOG(DEBUG) << "No tracks on the clipboard";
73
        return StatusCode::Success;
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
118
119
120
121
122
123
124
125
126
127
    }

    // Loop over all tracks
    for(auto& track : (*tracks)) {
        bool has_associated_cluster = false;
        bool is_within_roi = true;
        LOG(DEBUG) << "Looking at next track";

        // Cut on the chi2/ndof
        if(track->chi2ndof() > m_chi2ndofCut) {
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
            continue;
        }

        // Check if it intercepts the DUT
        auto globalIntercept = m_detector->getIntercept(track);
        auto localIntercept = m_detector->globalToLocal(globalIntercept);

        if(!m_detector->hasIntercept(track, 1.)) {
            LOG(DEBUG) << " - track outside DUT area";
            continue;
        }

        // Check that track is within region of interest using winding number algorithm
        if(!m_detector->isWithinROI(track)) {
            LOG(DEBUG) << " - track outside ROI";
            is_within_roi = false;
        }

        // Check that it doesn't go through/near a masked pixel
        if(m_detector->hitMasked(track, 1.)) {
            LOG(DEBUG) << " - track close to masked pixel";
            continue;
        }

        // Discard tracks which are very close to the frame edges
        if(fabs(track->timestamp() - clipboard->get_persistent("eventEnd")) < m_timeCutFrameEdge) {
            // Late edge - eventEnd points to the end of the frame`
            LOG(DEBUG) << " - track close to end of readout frame: "
                       << Units::display(fabs(track->timestamp() - clipboard->get_persistent("eventEnd")), {"us", "ns"})
                       << " at " << Units::display(track->timestamp(), {"us"});
            continue;
        } else if(fabs(track->timestamp() - clipboard->get_persistent("eventStart")) < m_timeCutFrameEdge) {
            // Early edge - eventStart points to the beginning of the frame
            LOG(DEBUG) << " - track close to start of readout frame: "
                       << Units::display(fabs(track->timestamp() - clipboard->get_persistent("eventStart")), {"us", "ns"})
                       << " at " << Units::display(track->timestamp(), {"us"});
            continue;
        }

        // Count this as reference track:
        total_tracks++;

        // Calculate in-pixel position of track in microns
128
129
130
        auto inpixel = m_detector->inPixel(localIntercept);
        auto xmod = static_cast<double>(Units::convert(inpixel.X(), "um"));
        auto ymod = static_cast<double>(Units::convert(inpixel.Y(), "um"));
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

        // Get the DUT clusters from the clipboard
        Clusters* clusters = reinterpret_cast<Clusters*>(clipboard->get(m_detector->name(), "clusters"));
        if(clusters == nullptr) {
            LOG(DEBUG) << " - no DUT clusters";
        } else {

            // Loop over all DUT clusters to find matches:
            for(auto* cluster : (*clusters)) {
                LOG(DEBUG) << " - Looking at next DUT cluster";

                auto associated_clusters = track->associatedClusters();
                if(std::find(associated_clusters.begin(), associated_clusters.end(), cluster) != associated_clusters.end()) {
                    LOG(DEBUG) << "Found associated cluster " << (*cluster);
                    has_associated_cluster = true;
                    matched_tracks++;
                    break;
                }
            }
        }

        hGlobalEfficiencyMap->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
        hChipEfficiencyMap->Fill(
            m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept), has_associated_cluster);

        // For pixels, only look at the ROI:
        if(is_within_roi) {
            hPixelEfficiencyMap->Fill(xmod, ymod, has_associated_cluster);
        }
    }

162
    return StatusCode::Success;
163
164
165
166
167
168
169
}

void AnalysisEfficiency::finalise() {
    LOG(STATUS) << "Total efficiency of detector " << m_detector->name() << ": "
                << (100 * matched_tracks / (total_tracks > 0 ? total_tracks : 1)) << "%, measured with " << total_tracks
                << " tracks";
}