AnalysisEfficiency.cpp 10.3 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
    m_timeCutFrameEdge = m_config.get<double>("time_cut_frameedge", Units::get<double>(20, "ns"));
23
    m_pixelTolerance = m_config.get<double>("pixel_tolerance", 1.);
24
    m_chi2ndofCut = m_config.get<double>("chi2ndof_cut", 3.);
25
    m_inpixelBinSize = m_config.get<double>("inpixel_bin_size", Units::get<double>(1.0, "um"));
26
27
28
29
30
31
32
}

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"));

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";
35
36
    hPixelEfficiencyMap_trackPos = new TProfile2D("pixelEfficiencyMap_trackPos",
                                                  title.c_str(),
37
                                                  static_cast<int>(ceil(pitch_x / m_inpixelBinSize)),
38
39
                                                  0,
                                                  pitch_x,
40
                                                  static_cast<int>(ceil(pitch_y / m_inpixelBinSize)),
41
42
43
44
                                                  0,
                                                  pitch_y,
                                                  0,
                                                  1);
45
    title = m_detector->name() + " Chip efficiency map;x [px];y [px];efficiency";
46
47
48
49
50
51
52
53
54
55
    hChipEfficiencyMap_trackPos = new TProfile2D("chipEfficiencyMap_trackPos",
                                                 title.c_str(),
                                                 m_detector->nPixels().X(),
                                                 0,
                                                 m_detector->nPixels().X(),
                                                 m_detector->nPixels().Y(),
                                                 0,
                                                 m_detector->nPixels().Y(),
                                                 0,
                                                 1);
56
    title = m_detector->name() + " Global efficiency map;x [mm];y [mm];efficiency";
57
58
59
60
61
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
    hGlobalEfficiencyMap_trackPos = new TProfile2D("globalEfficiencyMap_trackPos",
                                                   title.c_str(),
                                                   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);
    title = m_detector->name() + " Chip efficiency map;x [px];y [px];efficiency";
    hChipEfficiencyMap_clustPos = new TProfile2D("chipEfficiencyMap_clustPos",
                                                 title.c_str(),
                                                 m_detector->nPixels().X(),
                                                 0,
                                                 m_detector->nPixels().X(),
                                                 m_detector->nPixels().Y(),
                                                 0,
                                                 m_detector->nPixels().Y(),
                                                 0,
                                                 1);
    title = m_detector->name() + " Global efficiency map;x [mm];y [mm];efficiency";
    hGlobalEfficiencyMap_clustPos = new TProfile2D("globalEfficiencyMap_clustPos",
                                                   title.c_str(),
                                                   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);
89
90
}

91
StatusCode AnalysisEfficiency::run(std::shared_ptr<Clipboard> clipboard) {
92
93
94
95
96

    // 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";
97
        return StatusCode::Success;
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
    }

    // 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);

116
        if(!m_detector->hasIntercept(track, m_pixelTolerance)) {
117
            LOG(DEBUG) << " - track outside DUT area: " << localIntercept;
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
            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;
        }

Simon Spannagel's avatar
Simon Spannagel committed
133
134
135
        // Get the event:
        auto event = clipboard->get_event();

136
        // Discard tracks which are very close to the frame edges
Simon Spannagel's avatar
Simon Spannagel committed
137
        if(fabs(track->timestamp() - event->end()) < m_timeCutFrameEdge) {
138
139
            // Late edge - eventEnd points to the end of the frame`
            LOG(DEBUG) << " - track close to end of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
140
141
                       << Units::display(fabs(track->timestamp() - event->end()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
142
            continue;
Simon Spannagel's avatar
Simon Spannagel committed
143
        } else if(fabs(track->timestamp() - event->start()) < m_timeCutFrameEdge) {
144
145
            // Early edge - eventStart points to the beginning of the frame
            LOG(DEBUG) << " - track close to start of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
146
147
                       << Units::display(fabs(track->timestamp() - event->start()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
148
149
150
151
152
153
154
            continue;
        }

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

        // Calculate in-pixel position of track in microns
155
156
157
        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"));
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

        // 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++;
174
175
176
177
178
                    auto clusterLocal = m_detector->globalToLocal(cluster->global());
                    hGlobalEfficiencyMap_clustPos->Fill(
                        cluster->global().x(), cluster->global().y(), has_associated_cluster);
                    hChipEfficiencyMap_clustPos->Fill(
                        m_detector->getColumn(clusterLocal), m_detector->getRow(clusterLocal), has_associated_cluster);
179
180
181
182
                    break;
                }
            }
        }
183
184
        hGlobalEfficiencyMap_trackPos->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
        hChipEfficiencyMap_trackPos->Fill(
185
186
187
            m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept), has_associated_cluster);
        // For pixels, only look at the ROI:
        if(is_within_roi) {
188
189
190
191
192
193
            hPixelEfficiencyMap_trackPos->Fill(xmod, ymod, has_associated_cluster);
        }
        if(has_associated_cluster == false) {
            hGlobalEfficiencyMap_clustPos->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
            hChipEfficiencyMap_clustPos->Fill(
                m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept), has_associated_cluster);
194
195
196
        }
    }

197
    return StatusCode::Success;
198
199
200
}

void AnalysisEfficiency::finalise() {
201
202
    LOG(INFO) << "No. matched tracks=" << matched_tracks;
    LOG(INFO) << "Total no. tracks=" << total_tracks;
203
204
205
206
    LOG(STATUS) << "Total efficiency of detector " << m_detector->name() << ": "
                << (100 * matched_tracks / (total_tracks > 0 ? total_tracks : 1)) << "%, measured with " << total_tracks
                << " tracks";
}