Prealignment.cpp 8.79 KB
Newer Older
Simon Spannagel's avatar
Simon Spannagel committed
1
2
3
4
5
6
7
8
9
10
/**
 * @file
 * @brief Implementation of module Prealignment
 *
 * @copyright Copyright (c) 2018-2020 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.
 */

Morag Williams's avatar
Morag Williams committed
11
12
13
14
15
#include "Prealignment.h"

using namespace corryvreckan;
using namespace std;

16
17
18
Prealignment::Prealignment(Configuration config, std::shared_ptr<Detector> detector)
    : Module(std::move(config), detector), m_detector(detector) {

19
    max_correlation_rms = m_config.get<double>("max_correlation_rms", Units::get<double>(6, "mm"));
Morag Williams's avatar
Morag Williams committed
20
    damping_factor = m_config.get<double>("damping_factor", 1.0);
21

22
    method = m_config.get<std::string>("method", "mean");
23
    fit_range_rel = m_config.get<int>("fit_range_rel", 500);
24

25
26
27
    // Backwards compatibilty: also allow timing_cut to be used for time_cut_abs
    m_config.setAlias("time_cut_abs", "timing_cut", true);

Jens Kroeger's avatar
Jens Kroeger committed
28
29
30
31
32
33
34
35
    if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
        throw InvalidCombinationError(
            m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
    } else if(m_config.has("time_cut_abs")) {
        timeCut = m_config.get<double>("time_cut_abs");
    } else {
        timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution();
    }
36
37
    LOG(DEBUG) << "Setting max_correlation_rms to : " << max_correlation_rms;
    LOG(DEBUG) << "Setting damping_factor to : " << damping_factor;
Morag Williams's avatar
Morag Williams committed
38
39
40
41
}

void Prealignment::initialise() {

42
43
    // get the reference detector:
    std::shared_ptr<Detector> reference = get_reference();
Morag Williams's avatar
Morag Williams committed
44

45
46
47
48
49
50
51
52
53
    // Correlation plots
    std::string title = m_detector->name() + ": correlation X;x_{ref}-x [mm];events";
    correlationX = new TH1F("correlationX", title.c_str(), 1000, -10., 10.);
    title = m_detector->name() + ": correlation Y;y_{ref}-y [mm];events";
    correlationY = new TH1F("correlationY", title.c_str(), 1000, -10., 10.);
    // 2D correlation plots (pixel-by-pixel, local coordinates):
    title = m_detector->name() + ": 2D correlation X (local);x [px];x_{ref} [px];events";
    correlationX2Dlocal = new TH2F("correlationX_2Dlocal",
                                   title.c_str(),
54
                                   m_detector->nPixels().X(),
55
56
                                   -0.5,
                                   m_detector->nPixels().X() - 0.5,
57
                                   reference->nPixels().X(),
58
59
                                   -0.5,
                                   reference->nPixels().X() - 0.5);
60
61
62
    title = m_detector->name() + ": 2D correlation Y (local);y [px];y_{ref} [px];events";
    correlationY2Dlocal = new TH2F("correlationY_2Dlocal",
                                   title.c_str(),
63
                                   m_detector->nPixels().Y(),
64
65
                                   -0.5,
                                   m_detector->nPixels().Y() - 0.5,
66
                                   reference->nPixels().Y(),
67
68
                                   -0.5,
                                   reference->nPixels().Y() - 0.5);
69
70
71
72
    title = m_detector->name() + ": 2D correlation X (global);x [mm];x_{ref} [mm];events";
    correlationX2D = new TH2F("correlationX_2D", title.c_str(), 100, -10., 10., 100, -10., 10.);
    title = m_detector->name() + ": 2D correlation Y (global);y [mm];y_{ref} [mm];events";
    correlationY2D = new TH2F("correlationY_2D", title.c_str(), 100, -10., 10., 100, -10., 10.);
Morag Williams's avatar
Morag Williams committed
73
74
}

75
StatusCode Prealignment::run(std::shared_ptr<Clipboard> clipboard) {
Morag Williams's avatar
Morag Williams committed
76

77
    // Get the clusters
78
    auto clusters = clipboard->getData<Cluster>(m_detector->name());
79
80
    if(clusters == nullptr) {
        LOG(DEBUG) << "Detector " << m_detector->name() << " does not have any clusters on the clipboard";
81
        return StatusCode::Success;
82
    }
Morag Williams's avatar
Morag Williams committed
83

84
85
    // Get clusters from reference detector
    auto reference = get_reference();
86
    auto referenceClusters = clipboard->getData<Cluster>(reference->name());
87
88
    if(referenceClusters == nullptr) {
        LOG(DEBUG) << "Reference detector " << reference->name() << " does not have any clusters on the clipboard";
89
        return StatusCode::Success;
90
    }
Morag Williams's avatar
Morag Williams committed
91

92
93
94
95
96
    // Loop over all clusters and fill histograms
    for(auto& cluster : (*clusters)) {
        // Loop over reference plane pixels to make correlation plots
        for(auto& refCluster : (*referenceClusters)) {
            double timeDifference = refCluster->timestamp() - cluster->timestamp();
Morag Williams's avatar
Morag Williams committed
97

98
            // Correlation plots
Jens Kroeger's avatar
Jens Kroeger committed
99
            if(abs(timeDifference) < timeCut) {
100
101
                correlationX->Fill(refCluster->global().x() - cluster->global().x());
                correlationX2D->Fill(cluster->global().x(), refCluster->global().x());
102
                correlationX2Dlocal->Fill(cluster->column(), refCluster->column());
103
104
                correlationY->Fill(refCluster->global().y() - cluster->global().y());
                correlationY2D->Fill(cluster->global().y(), refCluster->global().y());
105
                correlationY2Dlocal->Fill(cluster->row(), refCluster->row());
106
            }
Morag Williams's avatar
Morag Williams committed
107
108
109
        }
    }

110
    return StatusCode::Success;
Morag Williams's avatar
Morag Williams committed
111
112
113
}

void Prealignment::finalise() {
114
115
116
117
118
119
120
121
122

    double rmsX = correlationX->GetRMS();
    double rmsY = correlationY->GetRMS();
    if(rmsX > max_correlation_rms or rmsY > max_correlation_rms) {
        LOG(ERROR) << "Detector " << m_detector->name() << ": RMS is too wide for prealignment shifts";
        LOG(ERROR) << "Detector " << m_detector->name() << ": RMS X = " << Units::display(rmsX, {"mm", "um"})
                   << " , RMS Y = " << Units::display(rmsY, {"mm", "um"});
    }

123
124
    // Move all but the reference:
    if(!m_detector->isReference()) {
125
126
127
128

        double shift_X = 0.;
        double shift_Y = 0.;

129
        LOG(INFO) << "Using prealignment method: " << method;
130
        if(method == "gauss_fit") {
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
            int binMaxX = correlationX->GetMaximumBin();
            double fit_low_x =
                correlationX->GetXaxis()->GetBinCenter(binMaxX) - m_detector->getSpatialResolution().x() * fit_range_rel;
            double fit_high_x =
                correlationX->GetXaxis()->GetBinCenter(binMaxX) + m_detector->getSpatialResolution().x() * fit_range_rel;

            int binMaxY = correlationY->GetMaximumBin();
            double fit_low_y =
                correlationY->GetXaxis()->GetBinCenter(binMaxY) - m_detector->getSpatialResolution().y() * fit_range_rel;
            double fit_high_y =
                correlationY->GetXaxis()->GetBinCenter(binMaxY) + m_detector->getSpatialResolution().y() * fit_range_rel;

            LOG(DEBUG) << "Fit range in x direction from: " << Units::display(fit_low_x, {"mm", "um"}) << " to "
                       << Units::display(fit_high_x, {"mm", "um"});
            LOG(DEBUG) << "Fit range in y direction from: " << Units::display(fit_low_y, {"mm", "um"}) << " to "
                       << Units::display(fit_high_y, {"mm", "um"});

            correlationX->Fit("gaus", "Q", "", fit_low_x, fit_high_x);
            correlationY->Fit("gaus", "Q", "", fit_low_y, fit_high_y);
150
151
            shift_X = correlationX->GetFunction("gaus")->GetParameter(1);
            shift_Y = correlationY->GetFunction("gaus")->GetParameter(1);
152
        } else if(method == "mean") {
153
154
            shift_X = correlationX->GetMean();
            shift_Y = correlationY->GetMean();
155
        } else if(method == "maximum") {
Katharina Dort's avatar
Katharina Dort committed
156
157
158
159
            int binMaxX = correlationX->GetMaximumBin();
            shift_X = correlationX->GetXaxis()->GetBinCenter(binMaxX);
            int binMaxY = correlationY->GetMaximumBin();
            shift_Y = correlationY->GetXaxis()->GetBinCenter(binMaxY);
160
161
        } else {
            throw InvalidValueError(m_config, "method", "Invalid prealignment method");
162
163
        }

164
165
166
        LOG(DEBUG) << "Shift (without damping factor)" << m_detector->name()
                   << ": x = " << Units::display(shift_X, {"mm", "um"})
                   << " , y = " << Units::display(shift_Y, {"mm", "um"});
167
168
        LOG(INFO) << "Move in x by = " << Units::display(shift_X * damping_factor, {"mm", "um"})
                  << " , and in y by = " << Units::display(shift_Y * damping_factor, {"mm", "um"});
169
170
171
172
        LOG(INFO) << "Detector position after shift in x = "
                  << Units::display(m_detector->displacement().X() + damping_factor * shift_X, {"mm", "um"})
                  << " , and in y = "
                  << Units::display(m_detector->displacement().Y() + damping_factor * shift_Y, {"mm", "um"});
173
174
        m_detector->displacement(XYZPoint(m_detector->displacement().X() + damping_factor * shift_X,
                                          m_detector->displacement().Y() + damping_factor * shift_Y,
175
                                          m_detector->displacement().Z()));
Morag Williams's avatar
Morag Williams committed
176
177
    }
}