Detector.cpp 15.2 KB
Newer Older
1
2
3
4
5
6
7
8
/** @file
 *  @brief Implementation of the detector model
 *  @copyright Copyright (c) 2017 CERN and the Corryvreckan 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.
 */

9
10
11
12
#include <fstream>
#include <map>
#include <string>

Simon Spannagel's avatar
Simon Spannagel committed
13
14
15
16
17
#include "Math/RotationX.h"
#include "Math/RotationY.h"
#include "Math/RotationZ.h"
#include "Math/RotationZYX.h"

18
19
#include "Detector.hpp"
#include "core/utils/log.h"
20
21
22
23

using namespace ROOT::Math;
using namespace corryvreckan;

24
25
Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {

Simon Spannagel's avatar
Simon Spannagel committed
26
    // Role of this detector:
27
28
29
30
31
32
33
34
35
36
37
38
    auto roles = config.getArray<std::string>("role", std::vector<std::string>{"none"});
    for(auto& role : roles) {
        std::transform(role.begin(), role.end(), role.begin(), ::tolower);
        if(role == "none") {
            m_role |= DetectorRole::NONE;
        } else if(role == "reference") {
            m_role |= DetectorRole::REFERENCE;
        } else if(role == "dut") {
            m_role |= DetectorRole::DUT;
        } else {
            throw InvalidValueError(config, "role", "Detector role does not exist.");
        }
Simon Spannagel's avatar
Simon Spannagel committed
39
40
41
    }

    // Detector position and orientation
42
43
    m_displacement = config.get<ROOT::Math::XYZPoint>("position", ROOT::Math::XYZPoint());
    m_orientation = config.get<ROOT::Math::XYZVector>("orientation", ROOT::Math::XYZVector());
44
    m_orientation_mode = config.get<std::string>("orientation_mode", "xyz");
45
46
47
    if(m_orientation_mode != "xyz" && m_orientation_mode != "zyx") {
        throw InvalidValueError(config, "orientation_mode", "Invalid detector orientation mode");
    }
48

49
50
51
    // Number of pixels
    auto npixels = config.get<ROOT::Math::DisplacementVector2D<Cartesian2D<int>>>("number_of_pixels");
    // Size of the pixels
52
    m_pitch = config.get<ROOT::Math::XYVector>("pixel_pitch");
53
54
55
56

    // Intrinsic position resolution, defaults to 4um:
    m_resolution = config.get<ROOT::Math::XYVector>("resolution", ROOT::Math::XYVector(0.004, 0.004));

57
    m_detectorName = config.getName();
58
59
60

    if(Units::convert(m_pitch.X(), "mm") >= 1 or Units::convert(m_pitch.Y(), "mm") >= 1 or
       Units::convert(m_pitch.X(), "um") <= 1 or Units::convert(m_pitch.Y(), "um") <= 1) {
61
62
        LOG(WARNING) << "Pixel pitch unphysical for detector " << m_detectorName << ": " << std::endl
                     << Units::display(m_pitch, {"nm", "um", "mm"});
63
64
    }

65
66
67
68
    m_detectorType = config.get<std::string>("type");
    m_nPixelsX = npixels.x();
    m_nPixelsY = npixels.y();
    m_timingOffset = config.get<double>("time_offset", 0.0);
69
    m_roi = config.getMatrix<int>("roi", std::vector<std::vector<int>>());
70
71
72
73

    this->initialise();

    LOG(TRACE) << "Initialized \"" << m_detectorType << "\": " << m_nPixelsX << "x" << m_nPixelsY << " px, pitch of "
74
75
76
               << Units::display(m_pitch, {"mm", "um"});
    LOG(TRACE) << "  Position:    " << Units::display(m_displacement, {"mm", "um"});
    LOG(TRACE) << "  Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")";
77
78
79
80
81
    if(m_timingOffset > 0.) {
        LOG(TRACE) << "Timing offset: " << m_timingOffset;
    }

    if(config.has("mask_file")) {
82
        m_maskfile_name = config.get<std::string>("mask_file");
83
84
85
        std::string mask_file = config.getPath("mask_file");
        LOG(DEBUG) << "Adding mask to detector \"" << config.getName() << "\", reading from " << mask_file;
        setMaskFile(mask_file);
86
        processMaskFile();
87
88
89
    }
}

90
bool Detector::isReference() {
91
    return static_cast<bool>(m_role & DetectorRole::REFERENCE);
92
93
94
}

bool Detector::isDUT() {
95
    return static_cast<bool>(m_role & DetectorRole::DUT);
96
97
}

98
99
100
// Functions to set and check channel masking
void Detector::setMaskFile(std::string file) {
    m_maskfile = file;
101
}
102

103
void Detector::processMaskFile() {
104
105
106
107
108
109
110
111
112
    // Open the file with masked pixels
    std::ifstream inputMaskFile(m_maskfile, std::ios::in);
    if(!inputMaskFile.is_open()) {
        LOG(ERROR) << "Could not open mask file " << m_maskfile;
    } else {
        int row = 0, col = 0;
        std::string id;
        std::string line;
        // loop over all lines and apply masks
Simon Spannagel's avatar
Simon Spannagel committed
113
        while(inputMaskFile >> id >> col >> row) {
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
            if(id == "c") {
                LOG(TRACE) << "Masking column " << col;
                int nRows = nPixelsY();
                for(int r = 0; r < nRows; r++) {
                    maskChannel(col, r);
                }
            } else if(id == "r") {
                LOG(TRACE) << "Masking row " << row;
                int nColumns = nPixelsX();
                for(int c = 0; c < nColumns; c++) {
                    maskChannel(c, row);
                }
            } else if(id == "p") {
                LOG(TRACE) << "Masking pixel " << col << " " << row;
                maskChannel(col, row); // Flag to mask a pixel
            } else {
                LOG(WARNING) << "Could not parse mask entry (id \"" << id << "\", col " << col << " row " << row << ")";
            }
        }
133
        LOG(INFO) << m_masked.size() << " masked pixels";
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
    }
}

void Detector::maskChannel(int chX, int chY) {
    int channelID = chX + m_nPixelsX * chY;
    m_masked[channelID] = true;
}

bool Detector::masked(int chX, int chY) {
    int channelID = chX + m_nPixelsX * chY;
    if(m_masked.count(channelID) > 0)
        return true;
    return false;
}

// Function to initialise transforms
void Detector::initialise() {

    // Make the local to global transform, built from a displacement and
    // rotation
154
    Translation3D translations = Translation3D(m_displacement.X(), m_displacement.Y(), m_displacement.Z());
155

156
    Rotation3D rotations;
157
    if(m_orientation_mode == "xyz") {
158
        rotations = RotationZ(m_orientation.Z()) * RotationY(m_orientation.Y()) * RotationX(m_orientation.X());
159
    } else if(m_orientation_mode == "zyx") {
160
        rotations = RotationZYX(m_orientation.x(), m_orientation.y(), m_orientation.x());
161
    }
162

163
164
    m_localToGlobal = Transform3D(rotations, translations);
    m_globalToLocal = m_localToGlobal.Inverse();
165

Simon Spannagel's avatar
Simon Spannagel committed
166
167
    // Find the normal to the detector surface. Build two points, the origin and a unit step in z,
    // transform these points to the global co-ordinate frame and then make a vector pointing between them
168
    m_origin = PositionVector3D<Cartesian3D<double>>(0., 0., 0.);
169
    m_origin = m_localToGlobal * m_origin;
170
    PositionVector3D<Cartesian3D<double>> localZ(0., 0., 1.);
171
    localZ = m_localToGlobal * localZ;
172
173
174
175
176
177
178
179
180
    m_normal = PositionVector3D<Cartesian3D<double>>(
        localZ.X() - m_origin.X(), localZ.Y() - m_origin.Y(), localZ.Z() - m_origin.Z());
}

// Function to update transforms (such as during alignment)
void Detector::update() {
    this->initialise();
}

181
182
183
184
185
Configuration Detector::getConfiguration() {

    Configuration config(name());
    config.set("type", m_detectorType);

186
    // Store the role of the detector
187
188
189
190
191
192
193
194
195
196
    std::vector<std::string> roles;
    if(this->isDUT()) {
        roles.push_back("dut");
    }
    if(this->isReference()) {
        roles.push_back("reference");
    }

    if(!roles.empty()) {
        config.setArray("role", roles);
197
198
    }

199
    config.set("position", m_displacement, {"um", "mm"});
200
    config.set("orientation_mode", m_orientation_mode);
201
    config.set("orientation", m_orientation, {"deg"});
202
203
204
205
    auto npixels = ROOT::Math::DisplacementVector2D<Cartesian2D<int>>(m_nPixelsX, m_nPixelsY);
    config.set("number_of_pixels", npixels);

    // Size of the pixels
206
    config.set("pixel_pitch", m_pitch, {"um"});
207

208
209
210
    // Intrinsic resolution:
    config.set("resolution", m_resolution, {"um"});

211
    if(m_timingOffset != 0.) {
212
        config.set("time_offset", m_timingOffset, {"ns", "us", "ms", "s"});
213
214
215
216
217
218
    }

    if(!m_maskfile_name.empty()) {
        config.set("mask_file", m_maskfile_name);
    }

219
220
    config.setMatrix("roi", m_roi);

221
222
223
    return config;
}

224
// Function to get global intercept with a track
225
PositionVector3D<Cartesian3D<double>> Detector::getIntercept(const Track* track) {
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

    // Get the distance from the plane to the track initial state
    double distance = (m_origin.X() - track->m_state.X()) * m_normal.X();
    distance += (m_origin.Y() - track->m_state.Y()) * m_normal.Y();
    distance += (m_origin.Z() - track->m_state.Z()) * m_normal.Z();
    distance /= (track->m_direction.X() * m_normal.X() + track->m_direction.Y() * m_normal.Y() +
                 track->m_direction.Z() * m_normal.Z());

    // Propagate the track
    PositionVector3D<Cartesian3D<double>> globalIntercept(track->m_state.X() + distance * track->m_direction.X(),
                                                          track->m_state.Y() + distance * track->m_direction.Y(),
                                                          track->m_state.Z() + distance * track->m_direction.Z());
    return globalIntercept;
}

241
PositionVector3D<Cartesian3D<double>> Detector::getLocalIntercept(const Track* track) {
242
243
244
    return globalToLocal(getIntercept(track));
}

245
// Function to check if a track intercepts with a plane
246
bool Detector::hasIntercept(const Track* track, double pixelTolerance) {
247
248
249
250
251

    // First, get the track intercept in global co-ordinates with the plane
    PositionVector3D<Cartesian3D<double>> globalIntercept = this->getIntercept(track);

    // Convert to local co-ordinates
252
    PositionVector3D<Cartesian3D<double>> localIntercept = this->m_globalToLocal * globalIntercept;
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273

    // Get the row and column numbers
    double row = this->getRow(localIntercept);
    double column = this->getColumn(localIntercept);

    // Check if the row and column are outside of the chip
    bool intercept = true;
    if(row < (pixelTolerance - 0.5) || row > (this->m_nPixelsY - 0.5 - pixelTolerance) || column < (pixelTolerance - 0.5) ||
       column > (this->m_nPixelsX - 0.5 - pixelTolerance))
        intercept = false;

    return intercept;
}

// Function to check if a track goes through/near a masked pixel
bool Detector::hitMasked(Track* track, int tolerance) {

    // First, get the track intercept in global co-ordinates with the plane
    PositionVector3D<Cartesian3D<double>> globalIntercept = this->getIntercept(track);

    // Convert to local co-ordinates
274
    PositionVector3D<Cartesian3D<double>> localIntercept = this->m_globalToLocal * globalIntercept;
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292

    // Get the row and column numbers
    int row = static_cast<int>(floor(this->getRow(localIntercept) + 0.5));
    int column = static_cast<int>(floor(this->getColumn(localIntercept) + 0.5));

    // Check if the pixels around this pixel are masked
    bool hitmasked = false;
    for(int r = (row - tolerance); r <= (row + tolerance); r++) {
        for(int c = (column - tolerance); c <= (column + tolerance); c++) {
            if(this->masked(c, r))
                hitmasked = true;
        }
    }

    return hitmasked;
}

// Functions to get row and column from local position
293
double Detector::getRow(const PositionVector3D<Cartesian3D<double>> localPosition) {
294
    double row = ((localPosition.Y() + m_pitch.Y() / 2.) / m_pitch.Y()) + m_nPixelsY / 2.;
295
296
    return row;
}
297
double Detector::getColumn(const PositionVector3D<Cartesian3D<double>> localPosition) {
298
    double column = ((localPosition.X() + m_pitch.X() / 2.) / m_pitch.X()) + m_nPixelsX / 2.;
299
300
301
302
303
304
    return column;
}

// Function to get local position from row and column
PositionVector3D<Cartesian3D<double>> Detector::getLocalPosition(double row, double column) {

Simon Spannagel's avatar
Simon Spannagel committed
305
    return PositionVector3D<Cartesian3D<double>>(
306
        m_pitch.X() * (column - m_nPixelsX / 2.), m_pitch.Y() * (row - m_nPixelsY / 2.), 0.);
307
308
}

309
// Function to get in-pixel position
310
double Detector::inPixelX(const PositionVector3D<Cartesian3D<double>> localPosition) {
311
    double column = getColumn(localPosition);
312
    double inPixelX = m_pitch.X() * (column - floor(column));
313
    return inPixelX;
314
}
315
double Detector::inPixelY(const PositionVector3D<Cartesian3D<double>> localPosition) {
316
    double row = getRow(localPosition);
317
    double inPixelY = m_pitch.Y() * (row - floor(row));
318
    return inPixelY;
319
}
320
321
322
323

// Check if track position is within ROI:
bool Detector::isWithinROI(const Track* track) {

324
325
326
327
328
    // Empty region of interest:
    if(m_roi.empty()) {
        return true;
    }

329
330
331
332
333
334
335
336
337
338
339
340
    // Check that track is within region of interest using winding number algorithm
    auto localIntercept = this->getLocalIntercept(track);
    auto coordinates = std::make_pair(this->getColumn(localIntercept), this->getRow(localIntercept));
    if(winding_number(coordinates, m_roi) != 0) {
        return true;
    }

    // Outside ROI:
    return false;
}

// Check if cluster is within ROI and/or touches ROI border:
341
342
bool Detector::isWithinROI(Cluster* cluster) {

343
344
345
346
347
    // Empty region of interest:
    if(m_roi.empty()) {
        return true;
    }

348
349
350
351
352
353
354
355
    // Loop over all pixels of the cluster
    for(auto& pixel : (*cluster->pixels())) {
        if(winding_number(pixel->coordinates(), m_roi) == 0) {
            return false;
        }
    }
    return true;
}
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386

/* isLeft(): tests if a point is Left|On|Right of an infinite line.
 * via: http://geomalgorithms.com/a03-_inclusion.html
 *    Input:  three points P0, P1, and P2
 *    Return: >0 for P2 left of the line through P0 and P1
 *            =0 for P2  on the line
 *            <0 for P2  right of the line
 *    See: Algorithm 1 "Area of Triangles and Polygons"
 */
int Detector::isLeft(std::pair<int, int> pt0, std::pair<int, int> pt1, std::pair<int, int> pt2) {
    return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second));
}

/* Winding number test for a point in a polygon
 * via: http://geomalgorithms.com/a03-_inclusion.html
 *      Input:   x, y = a point,
 *               polygon = vector of vertex points of a polygon V[n+1] with V[n]=V[0]
 *      Return:  wn = the winding number (=0 only when P is outside)
 */
int Detector::winding_number(std::pair<int, int> probe, std::vector<std::vector<int>> polygon) {
    // Two points don't make an area
    if(polygon.size() < 3) {
        LOG(DEBUG) << "No ROI given.";
        return 0;
    }

    int wn = 0; // the  winding number counter

    // loop through all edges of the polygon

    // edge from V[i] to  V[i+1]
387
    for(size_t i = 0; i < polygon.size(); i++) {
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
        auto point_this = std::make_pair(polygon.at(i).at(0), polygon.at(i).at(1));
        auto point_next = (i + 1 < polygon.size() ? std::make_pair(polygon.at(i + 1).at(0), polygon.at(i + 1).at(1))
                                                  : std::make_pair(polygon.at(0).at(0), polygon.at(0).at(1)));

        // start y <= P.y
        if(point_this.second <= probe.second) {
            // an upward crossing
            if(point_next.second > probe.second) {
                // P left of  edge
                if(isLeft(point_this, point_next, probe) > 0) {
                    // have  a valid up intersect
                    ++wn;
                }
            }
        } else {
            // start y > P.y (no test needed)

            // a downward crossing
            if(point_next.second <= probe.second) {
                // P right of  edge
                if(isLeft(point_this, point_next, probe) < 0) {
                    // have  a valid down intersect
                    --wn;
                }
            }
        }
    }
    return wn;
}