Detector.cpp 16.8 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
    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;
32
        } else if(role == "reference" || role == "ref") {
33
34
35
            m_role |= DetectorRole::REFERENCE;
        } else if(role == "dut") {
            m_role |= DetectorRole::DUT;
36
37
        } else if(role == "auxiliary" || role == "aux") {
            m_role |= DetectorRole::AUXILIARY;
38
39
40
        } else {
            throw InvalidValueError(config, "role", "Detector role does not exist.");
        }
Simon Spannagel's avatar
Simon Spannagel committed
41
42
    }

43
44
45
46
47
    // Auxiliary devices cannot hold other roles:
    if(static_cast<bool>(m_role & DetectorRole::AUXILIARY) && m_role != DetectorRole::AUXILIARY) {
        throw InvalidValueError(config, "role", "Auxiliary devices cannot hold any other detector role");
    }

Simon Spannagel's avatar
Simon Spannagel committed
48
    // Detector position and orientation
49
50
    m_displacement = config.get<ROOT::Math::XYZPoint>("position", ROOT::Math::XYZPoint());
    m_orientation = config.get<ROOT::Math::XYZVector>("orientation", ROOT::Math::XYZVector());
51
    m_orientation_mode = config.get<std::string>("orientation_mode", "xyz");
52
53
54
    if(m_orientation_mode != "xyz" && m_orientation_mode != "zyx") {
        throw InvalidValueError(config, "orientation_mode", "Invalid detector orientation mode");
    }
55

56
    // Number of pixels
57
    m_nPixels = config.get<ROOT::Math::DisplacementVector2D<Cartesian2D<int>>>("number_of_pixels");
58
    // Size of the pixels
59
    m_pitch = config.get<ROOT::Math::XYVector>("pixel_pitch");
60
61
62
63

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

64
    m_detectorName = config.getName();
65
66
67

    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) {
68
69
        LOG(WARNING) << "Pixel pitch unphysical for detector " << m_detectorName << ": " << std::endl
                     << Units::display(m_pitch, {"nm", "um", "mm"});
70
71
    }

72
    m_detectorType = config.get<std::string>("type");
73
    std::transform(m_detectorType.begin(), m_detectorType.end(), m_detectorType.begin(), ::tolower);
74
    m_timingOffset = config.get<double>("time_offset", 0.0);
75
    m_roi = config.getMatrix<int>("roi", std::vector<std::vector<int>>());
76
77
78

    this->initialise();

79
    LOG(TRACE) << "Initialized \"" << m_detectorType << "\": " << m_nPixels.X() << "x" << m_nPixels.Y() << " px, pitch of "
80
81
82
               << 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 << ")";
83
84
85
86
87
    if(m_timingOffset > 0.) {
        LOG(TRACE) << "Timing offset: " << m_timingOffset;
    }

    if(config.has("mask_file")) {
88
        m_maskfile_name = config.get<std::string>("mask_file");
89
90
91
        std::string mask_file = config.getPath("mask_file");
        LOG(DEBUG) << "Adding mask to detector \"" << config.getName() << "\", reading from " << mask_file;
        setMaskFile(mask_file);
92
        processMaskFile();
93
94
95
    }
}

96
97
98
99
100
101
102
103
104
105
106
107
std::string Detector::name() const {
    return m_detectorName;
}

std::string Detector::type() const {
    return m_detectorType;
}

XYVector Detector::size() const {
    return XYVector(m_pitch.X() * m_nPixels.X(), m_pitch.Y() * m_nPixels.Y());
}

108
bool Detector::isReference() const {
109
    return static_cast<bool>(m_role & DetectorRole::REFERENCE);
110
111
}

112
bool Detector::isDUT() const {
113
    return static_cast<bool>(m_role & DetectorRole::DUT);
114
115
}

116
117
118
119
bool Detector::isAuxiliary() const {
    return static_cast<bool>(m_role & DetectorRole::AUXILIARY);
}

120
121
122
// Functions to set and check channel masking
void Detector::setMaskFile(std::string file) {
    m_maskfile = file;
123
}
124

125
void Detector::processMaskFile() {
126
127
128
129
130
131
132
133
    // 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;
        // loop over all lines and apply masks
Jens Kroeger's avatar
Jens Kroeger committed
134
        while(inputMaskFile >> id) {
135
            if(id == "c") {
Jens Kroeger's avatar
Jens Kroeger committed
136
                inputMaskFile >> col;
137
138
                if(col > nPixels().X() - 1) {
                    LOG(ERROR) << "Column " << col << " outside of pixel matrix, chip has only " << nPixels().X()
Jens Kroeger's avatar
Jens Kroeger committed
139
140
                               << " columns!";
                }
141
                LOG(TRACE) << "Masking column " << col;
142
                for(int r = 0; r < nPixels().Y(); r++) {
143
144
145
                    maskChannel(col, r);
                }
            } else if(id == "r") {
Jens Kroeger's avatar
Jens Kroeger committed
146
                inputMaskFile >> row;
147
148
                if(row > nPixels().Y() - 1) {
                    LOG(ERROR) << "Row " << col << " outside of pixel matrix, chip has only " << nPixels().Y() << " rows!";
Jens Kroeger's avatar
Jens Kroeger committed
149
                }
150
                LOG(TRACE) << "Masking row " << row;
151
                for(int c = 0; c < nPixels().X(); c++) {
152
153
154
                    maskChannel(c, row);
                }
            } else if(id == "p") {
Jens Kroeger's avatar
Jens Kroeger committed
155
                inputMaskFile >> col >> row;
156
157
158
                if(col > nPixels().X() - 1 || row > nPixels().Y() - 1) {
                    LOG(ERROR) << "Pixel " << col << " " << row << " outside of pixel matrix, chip has only "
                               << nPixels().X() << " x " << nPixels().Y() << " pixels!";
Jens Kroeger's avatar
Jens Kroeger committed
159
                }
160
161
162
                LOG(TRACE) << "Masking pixel " << col << " " << row;
                maskChannel(col, row); // Flag to mask a pixel
            } else {
Jens Kroeger's avatar
Jens Kroeger committed
163
                LOG(ERROR) << "Could not parse mask entry (id \"" << id << "\")";
164
165
            }
        }
166
        LOG(INFO) << m_masked.size() << " masked pixels";
167
168
169
170
    }
}

void Detector::maskChannel(int chX, int chY) {
171
    int channelID = chX + m_nPixels.X() * chY;
172
173
174
    m_masked[channelID] = true;
}

175
bool Detector::masked(int chX, int chY) const {
176
    int channelID = chX + m_nPixels.X() * chY;
177
178
179
180
181
182
183
184
185
186
    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
187
    Translation3D translations = Translation3D(m_displacement.X(), m_displacement.Y(), m_displacement.Z());
188

189
    Rotation3D rotations;
190
    if(m_orientation_mode == "xyz") {
191
        rotations = RotationZ(m_orientation.Z()) * RotationY(m_orientation.Y()) * RotationX(m_orientation.X());
192
    } else if(m_orientation_mode == "zyx") {
193
        rotations = RotationZYX(m_orientation.x(), m_orientation.y(), m_orientation.x());
194
    }
195

196
197
    m_localToGlobal = Transform3D(rotations, translations);
    m_globalToLocal = m_localToGlobal.Inverse();
198

Simon Spannagel's avatar
Simon Spannagel committed
199
200
    // 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
201
    m_origin = PositionVector3D<Cartesian3D<double>>(0., 0., 0.);
202
    m_origin = m_localToGlobal * m_origin;
203
    PositionVector3D<Cartesian3D<double>> localZ(0., 0., 1.);
204
    localZ = m_localToGlobal * localZ;
205
206
207
208
209
210
211
212
213
    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();
}

214
Configuration Detector::getConfiguration() const {
215
216
217
218

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

219
    // Store the role of the detector
220
221
222
223
224
225
226
    std::vector<std::string> roles;
    if(this->isDUT()) {
        roles.push_back("dut");
    }
    if(this->isReference()) {
        roles.push_back("reference");
    }
227
228
229
    if(this->isAuxiliary()) {
        roles.push_back("auxiliary");
    }
230
231
232

    if(!roles.empty()) {
        config.setArray("role", roles);
233
234
    }

235
    config.set("position", m_displacement, {"um", "mm"});
236
    config.set("orientation_mode", m_orientation_mode);
237
    config.set("orientation", m_orientation, {"deg"});
238
    config.set("number_of_pixels", m_nPixels);
239
240

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

243
244
245
    // Intrinsic resolution:
    config.set("resolution", m_resolution, {"um"});

246
    if(m_timingOffset != 0.) {
247
        config.set("time_offset", m_timingOffset, {"ns", "us", "ms", "s"});
248
249
250
251
252
253
    }

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

254
255
    config.setMatrix("roi", m_roi);

256
257
258
    return config;
}

259
// Function to get global intercept with a track
260
PositionVector3D<Cartesian3D<double>> Detector::getIntercept(const Track* track) const {
261
262

    // Get the distance from the plane to the track initial state
263
264
265
266
267
    double distance = (m_origin.X() - track->state().X()) * m_normal.X();
    distance += (m_origin.Y() - track->state().Y()) * m_normal.Y();
    distance += (m_origin.Z() - track->state().Z()) * m_normal.Z();
    distance /= (track->direction().X() * m_normal.X() + track->direction().Y() * m_normal.Y() +
                 track->direction().Z() * m_normal.Z());
268
269

    // Propagate the track
270
271
272
    PositionVector3D<Cartesian3D<double>> globalIntercept(track->state().X() + distance * track->direction().X(),
                                                          track->state().Y() + distance * track->direction().Y(),
                                                          track->state().Z() + distance * track->direction().Z());
273
274
275
    return globalIntercept;
}

276
PositionVector3D<Cartesian3D<double>> Detector::getLocalIntercept(const Track* track) const {
277
278
279
    return globalToLocal(getIntercept(track));
}

280
// Function to check if a track intercepts with a plane
281
bool Detector::hasIntercept(const Track* track, double pixelTolerance) const {
282
283
284
285
286

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

    // Convert to local co-ordinates
287
    PositionVector3D<Cartesian3D<double>> localIntercept = this->m_globalToLocal * globalIntercept;
288
289
290
291
292
293
294

    // 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;
Morag Jean Williams's avatar
Morag Jean Williams committed
295
    if(row < pixelTolerance || row > (this->m_nPixels.Y() - pixelTolerance) || column < pixelTolerance ||
296
       column > (this->m_nPixels.X() - pixelTolerance))
297
298
299
300
301
302
        intercept = false;

    return intercept;
}

// Function to check if a track goes through/near a masked pixel
303
bool Detector::hitMasked(Track* track, int tolerance) const {
304
305
306
307
308

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

    // Convert to local co-ordinates
309
    PositionVector3D<Cartesian3D<double>> localIntercept = this->m_globalToLocal * globalIntercept;
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327

    // 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
328
double Detector::getRow(const PositionVector3D<Cartesian3D<double>> localPosition) const {
329
    double row = localPosition.Y() / m_pitch.Y() + static_cast<double>(m_nPixels.Y()) / 2. + (1 - m_nPixels.Y() % 2) / 2.;
330
331
    return row;
}
332
double Detector::getColumn(const PositionVector3D<Cartesian3D<double>> localPosition) const {
333
    double column = localPosition.X() / m_pitch.X() + static_cast<double>(m_nPixels.X()) / 2. + (1 - m_nPixels.X() % 2) / 2.;
334
335
336
337
    return column;
}

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

Simon Spannagel's avatar
Simon Spannagel committed
340
    return PositionVector3D<Cartesian3D<double>>(
341
        m_pitch.X() * (column - (m_nPixels.X()) / 2.), m_pitch.Y() * (row - (m_nPixels.Y()) / 2.), 0.);
342
343
}

344
// Function to get in-pixel position
345
ROOT::Math::XYVector Detector::inPixel(const double column, const double row) const {
346
347
    // a pixel ranges from (col-0.5) to (col+0.5)
    return XYVector(m_pitch.X() * (column - floor(column) - 0.5), m_pitch.Y() * (row - floor(row) - 0.5));
348
}
349

350
ROOT::Math::XYVector Detector::inPixel(const PositionVector3D<Cartesian3D<double>> localPosition) const {
351
352
    double column = getColumn(localPosition);
    double row = getRow(localPosition);
353
    return inPixel(column, row);
354
}
355
356

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

359
360
361
362
363
    // Empty region of interest:
    if(m_roi.empty()) {
        return true;
    }

364
365
366
367
368
369
370
371
372
373
374
375
    // 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:
376
bool Detector::isWithinROI(Cluster* cluster) const {
377

378
379
380
381
382
    // Empty region of interest:
    if(m_roi.empty()) {
        return true;
    }

383
    // Loop over all pixels of the cluster
384
    for(auto& pixel : cluster->pixels()) {
385
386
387
388
389
390
        if(winding_number(pixel->coordinates(), m_roi) == 0) {
            return false;
        }
    }
    return true;
}
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
417
418
419
420
421

/* 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]
422
    for(size_t i = 0; i < polygon.size(); i++) {
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
        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;
}