AlignmentTrackChi2.cpp 15.4 KB
Newer Older
1
#include "AlignmentTrackChi2.h"
2
3

#include <TVirtualFitter.h>
4
#include <numeric>
5

6
using namespace corryvreckan;
7
using namespace std;
8

9
10
// Global container declarations
Tracks globalTracks;
11
std::shared_ptr<Detector> globalDetector;
12
13
int detNum;

14
AlignmentTrackChi2::AlignmentTrackChi2(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
15
    : Module(std::move(config), std::move(detectors)) {
16
17
    m_numberOfTracksForAlignment = m_config.get<size_t>("number_of_tracks", 20000);
    nIterations = m_config.get<size_t>("iterations", 3);
18

19
    m_pruneTracks = m_config.get<bool>("prune_tracks", false);
20
    m_alignPosition = m_config.get<bool>("alignPosition", true);
21
22
23
    if(m_alignPosition) {
        LOG(INFO) << "Aligning positions";
    }
24
    m_alignOrientation = m_config.get<bool>("alignOrientation", true);
25
26
27
    if(m_alignOrientation) {
        LOG(INFO) << "Aligning orientations";
    }
28
    m_maxAssocClusters = m_config.get<size_t>("max_associated_clusters", 1);
29
    m_maxTrackChi2 = m_config.get<double>("max_track_chi2ndof", 10.);
30
    LOG(INFO) << "Aligning telescope";
31
}
32

33
// During run, just pick up tracks and save them till the end
34
StatusCode AlignmentTrackChi2::run(std::shared_ptr<Clipboard> clipboard) {
35

Simon Spannagel's avatar
Simon Spannagel committed
36
    // Get the tracks
37
38
    Tracks* tracks = reinterpret_cast<Tracks*>(clipboard->get("tracks"));
    if(tracks == nullptr) {
39
        return StatusCode::Success;
Simon Spannagel's avatar
Simon Spannagel committed
40
    }
41

Simon Spannagel's avatar
Simon Spannagel committed
42
    // Make a local copy and store it
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
    for(auto& track : (*tracks)) {

        // Apply selection to tracks for alignment
        if(m_pruneTracks) {
            // Only allow one associated cluster:
            if(track->associatedClusters().size() > m_maxAssocClusters) {
                LOG(DEBUG) << "Discarded track with " << track->associatedClusters().size() << " associated clusters";
                m_discardedtracks++;
                continue;
            }

            // Only allow tracks with certain Chi2/NDoF:
            if(track->chi2ndof() > m_maxTrackChi2) {
                LOG(DEBUG) << "Discarded track with Chi2/NDoF - " << track->chi2ndof();
                m_discardedtracks++;
                continue;
            }
        }

Simon Spannagel's avatar
Simon Spannagel committed
62
63
        Track* alignmentTrack = new Track(track);
        m_alignmenttracks.push_back(alignmentTrack);
64
    }
65

Simon Spannagel's avatar
Simon Spannagel committed
66
    // If we have enough tracks for the alignment, tell the event loop to finish
67
68
    if(m_alignmenttracks.size() >= m_numberOfTracksForAlignment) {
        LOG(STATUS) << "Accumulated " << m_alignmenttracks.size() << " tracks, interrupting processing.";
69
70
71
        if(m_discardedtracks > 0) {
            LOG(INFO) << "Discarded " << m_discardedtracks << " input tracks.";
        }
72
        return StatusCode::EndRun;
73
    }
74

Simon Spannagel's avatar
Simon Spannagel committed
75
    // Otherwise keep going
76
    return StatusCode::Success;
77
78
}

79
80
81
82
// ========================================
//  Minimisation functions for Minuit
// ========================================

83
// METHOD 0
Simon Spannagel's avatar
Simon Spannagel committed
84
85
86
87
// This method will move the detector in question, refit all of the tracks, and
// try to minimise the
// track chi2. If there were no clusters from this detector on any tracks then
// it would do nothing!
88
void AlignmentTrackChi2::MinimiseTrackChi2(Int_t&, Double_t*, Double_t& result, Double_t* par, Int_t) {
Simon Spannagel's avatar
Simon Spannagel committed
89
90

    // Pick up new alignment conditions
91
92
93
94
95
96
    globalDetector->displacementX(par[detNum * 6 + 0]);
    globalDetector->displacementY(par[detNum * 6 + 1]);
    globalDetector->displacementZ(par[detNum * 6 + 2]);
    globalDetector->rotationX(par[detNum * 6 + 3]);
    globalDetector->rotationY(par[detNum * 6 + 4]);
    globalDetector->rotationZ(par[detNum * 6 + 5]);
Simon Spannagel's avatar
Simon Spannagel committed
97
98

    // Apply new alignment conditions
99
    globalDetector->update();
Simon Spannagel's avatar
Simon Spannagel committed
100
101
102
103
104

    // The chi2 value to be returned
    result = 0.;

    // Loop over all tracks
105
    for(size_t iTrack = 0; iTrack < globalTracks.size(); iTrack++) {
Simon Spannagel's avatar
Simon Spannagel committed
106
107
108
109
110
        // Get the track
        Track* track = globalTracks[iTrack];
        // Get all clusters on the track
        Clusters trackClusters = track->clusters();
        // Find the cluster that needs to have its position recalculated
111
        for(size_t iTrackCluster = 0; iTrackCluster < trackClusters.size(); iTrackCluster++) {
Simon Spannagel's avatar
Simon Spannagel committed
112
            Cluster* trackCluster = trackClusters[iTrackCluster];
113
            if(globalDetector->name() != trackCluster->detectorID()) {
Simon Spannagel's avatar
Simon Spannagel committed
114
                continue;
115
            }
116

Simon Spannagel's avatar
Simon Spannagel committed
117
            // Recalculate the global position from the local
118
119
120
            auto positionLocal = trackCluster->local();
            auto positionGlobal = globalDetector->localToGlobal(positionLocal);
            trackCluster->setClusterCentre(positionGlobal);
Simon Spannagel's avatar
Simon Spannagel committed
121
122
123
124
125
126
127
        }

        // Refit the track
        track->fit();

        // Add the new chi2
        result += track->chi2();
128
129
130
    }
}

Daniel Hynds's avatar
Daniel Hynds committed
131
132
133
134
// ==================================================================
//  The finalise function - effectively the brains of the alignment!
// ==================================================================

135
void AlignmentTrackChi2::finalise() {
Simon Spannagel's avatar
Simon Spannagel committed
136
137
138
139
140

    // If not enough tracks were produced, do nothing
    // if(m_alignmenttracks.size() < m_numberOfTracksForAlignment) return;

    // Make the fitting object
141
    TVirtualFitter* residualFitter = TVirtualFitter::Fitter(nullptr, 50);
142
    residualFitter->SetFCN(MinimiseTrackChi2);
Simon Spannagel's avatar
Simon Spannagel committed
143
144
145
146
147
148

    // Set the global parameters
    globalTracks = m_alignmenttracks;

    // Set the printout arguments of the fitter
    Double_t arglist[10];
149
    arglist[0] = -1;
Simon Spannagel's avatar
Simon Spannagel committed
150
151
152
153
154
155
    residualFitter->ExecuteCommand("SET PRINT", arglist, 1);

    // Set some fitter parameters
    arglist[0] = 1000;  // number of function calls
    arglist[1] = 0.001; // tolerance

156
157
158
159
160
161
162
163
    // Store the alignment shifts per detector:
    std::map<std::string, std::vector<double>> shiftsX;
    std::map<std::string, std::vector<double>> shiftsY;
    std::map<std::string, std::vector<double>> shiftsZ;
    std::map<std::string, std::vector<double>> rotX;
    std::map<std::string, std::vector<double>> rotY;
    std::map<std::string, std::vector<double>> rotZ;

164
165
    // Loop over all planes. For each plane, set the plane alignment parameters which will be varied, and then minimise the
    // track chi2 (sum of biased residuals). This means that tracks are refitted with each minimisation step.
166
    for(size_t iteration = 0; iteration < nIterations; iteration++) {
Simon Spannagel's avatar
Simon Spannagel committed
167

Simon Spannagel's avatar
Simon Spannagel committed
168
169
170
        LOG_PROGRESS(STATUS, "alignment_track") << "Alignment iteration " << (iteration + 1) << " of " << nIterations;

        int det = 0;
171
172
        for(auto& detector : get_detectors()) {
            string detectorID = detector->name();
173

Simon Spannagel's avatar
Simon Spannagel committed
174
            // Do not align the reference plane
175
            if(detector->isReference() || detector->isDUT()) {
Simon Spannagel's avatar
Simon Spannagel committed
176
                continue;
177
178
            }

Simon Spannagel's avatar
Simon Spannagel committed
179
            // Say that this is the detector we align
180
181
            globalDetector = detector;

Simon Spannagel's avatar
Simon Spannagel committed
182
183
            detNum = det;
            // Add the parameters to the fitter (z displacement not allowed to move!)
184
185
186
187
188
189
190
191
192
193
194
195
            if(m_alignPosition) {
                residualFitter->SetParameter(
                    det * 6 + 0, (detectorID + "_displacementX").c_str(), detector->displacement().X(), 0.01, -50, 50);
                residualFitter->SetParameter(
                    det * 6 + 1, (detectorID + "_displacementY").c_str(), detector->displacement().Y(), 0.01, -50, 50);

            } else {
                residualFitter->SetParameter(
                    det * 6 + 0, (detectorID + "_displacementX").c_str(), detector->displacement().X(), 0, -50, 50);
                residualFitter->SetParameter(
                    det * 6 + 1, (detectorID + "_displacementY").c_str(), detector->displacement().Y(), 0, -50, 50);
            }
196
            residualFitter->SetParameter(
197
                det * 6 + 2, (detectorID + "_displacementZ").c_str(), detector->displacement().Z(), 0, -10, 500);
Simon Spannagel's avatar
Simon Spannagel committed
198

199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
            if(m_alignOrientation) {
                residualFitter->SetParameter(
                    det * 6 + 3, (detectorID + "_rotationX").c_str(), detector->rotation().X(), 0.001, -6.30, 6.30);
                residualFitter->SetParameter(
                    det * 6 + 4, (detectorID + "_rotationY").c_str(), detector->rotation().Y(), 0.001, -6.30, 6.30);
                residualFitter->SetParameter(
                    det * 6 + 5, (detectorID + "_rotationZ").c_str(), detector->rotation().Z(), 0.001, -6.30, 6.30);
            } else {
                residualFitter->SetParameter(
                    det * 6 + 3, (detectorID + "_rotationX").c_str(), detector->rotation().X(), 0, -6.30, 6.30);
                residualFitter->SetParameter(
                    det * 6 + 4, (detectorID + "_rotationY").c_str(), detector->rotation().Y(), 0, -6.30, 6.30);
                residualFitter->SetParameter(
                    det * 6 + 5, (detectorID + "_rotationZ").c_str(), detector->rotation().Z(), 0, -6.30, 6.30);
            }
214
215
216
            auto old_position = detector->displacement();
            auto old_orientation = detector->rotation();

Simon Spannagel's avatar
Simon Spannagel committed
217
218
219
            // Fit this plane (minimising global track chi2)
            residualFitter->ExecuteCommand("MIGRAD", arglist, 2);

220
221
222
223
224
225
226
227
            // Retrieve fit results:
            auto displacementX = residualFitter->GetParameter(det * 6 + 0);
            auto displacementY = residualFitter->GetParameter(det * 6 + 1);
            auto displacementZ = residualFitter->GetParameter(det * 6 + 2);
            auto rotationX = residualFitter->GetParameter(det * 6 + 3);
            auto rotationY = residualFitter->GetParameter(det * 6 + 4);
            auto rotationZ = residualFitter->GetParameter(det * 6 + 5);

228
            // Store corrections:
229
230
231
232
233
234
235
236
237
238
239
240
            shiftsX[detectorID].push_back(
                static_cast<double>(Units::convert(detector->displacement().X() - old_position.X(), "um")));
            shiftsY[detectorID].push_back(
                static_cast<double>(Units::convert(detector->displacement().Y() - old_position.Y(), "um")));
            shiftsZ[detectorID].push_back(
                static_cast<double>(Units::convert(detector->displacement().Z() - old_position.Z(), "um")));
            rotX[detectorID].push_back(
                static_cast<double>(Units::convert(detector->rotation().X() - old_orientation.X(), "deg")));
            rotY[detectorID].push_back(
                static_cast<double>(Units::convert(detector->rotation().Y() - old_orientation.Y(), "deg")));
            rotZ[detectorID].push_back(
                static_cast<double>(Units::convert(detector->rotation().Z() - old_orientation.Z(), "deg")));
241
242

            LOG(INFO) << detector->name() << "/" << iteration << " dT"
243
244
                      << Units::display(detector->displacement() - old_position, {"mm", "um"}) << " dR"
                      << Units::display(detector->rotation() - old_orientation, {"deg"});
245

Simon Spannagel's avatar
Simon Spannagel committed
246
247
            // Now that this device is fitted, set parameter errors to 0 so that they
            // are not fitted again
248
249
250
251
252
253
            residualFitter->SetParameter(det * 6 + 0, (detectorID + "_displacementX").c_str(), displacementX, 0, -50, 50);
            residualFitter->SetParameter(det * 6 + 1, (detectorID + "_displacementY").c_str(), displacementY, 0, -50, 50);
            residualFitter->SetParameter(det * 6 + 2, (detectorID + "_displacementZ").c_str(), displacementZ, 0, -10, 500);
            residualFitter->SetParameter(det * 6 + 3, (detectorID + "_rotationX").c_str(), rotationX, 0, -6.30, 6.30);
            residualFitter->SetParameter(det * 6 + 4, (detectorID + "_rotationY").c_str(), rotationY, 0, -6.30, 6.30);
            residualFitter->SetParameter(det * 6 + 5, (detectorID + "_rotationZ").c_str(), rotationZ, 0, -6.30, 6.30);
Simon Spannagel's avatar
Simon Spannagel committed
254
255
256

            // Set the alignment parameters of this plane to be the optimised values
            // from the alignment
257
258
259
260
261
262
            detector->displacementX(displacementX);
            detector->displacementY(displacementY);
            detector->displacementZ(displacementZ);
            detector->rotationX(rotationX);
            detector->rotationY(rotationY);
            detector->rotationZ(rotationZ);
263
            detector->update();
Simon Spannagel's avatar
Simon Spannagel committed
264
265
266
            det++;
        }
    }
267

Simon Spannagel's avatar
Simon Spannagel committed
268
    LOG_PROGRESS(STATUS, "alignment_track") << "Alignment finished, " << nIterations << " iteration.";
269

Simon Spannagel's avatar
Simon Spannagel committed
270
    // Now list the new alignment parameters
271
    for(auto& detector : get_detectors()) {
Simon Spannagel's avatar
Simon Spannagel committed
272
        // Do not align the reference plane
273
        if(detector->isReference() || detector->isDUT()) {
Simon Spannagel's avatar
Simon Spannagel committed
274
275
            continue;
        }
Simon Spannagel's avatar
Simon Spannagel committed
276

277
        LOG(STATUS) << detector->name() << " new alignment: " << std::endl
278
279
                    << "T" << Units::display(detector->displacement(), {"mm", "um"}) << " R"
                    << Units::display(detector->rotation(), {"deg"});
280
281
282
283
284
285
286

        // Fill the alignment convergence graphs:
        std::vector<double> iterations(nIterations);
        std::iota(std::begin(iterations), std::end(iterations), 0);

        std::string name = "alignment_correction_displacementX_" + detector->name();
        align_correction_shiftX[detector->name()] =
287
            new TGraph(static_cast<int>(shiftsX[detector->name()].size()), &iterations[0], &shiftsX[detector->name()][0]);
288
289
        align_correction_shiftX[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_shiftX[detector->name()]->GetYaxis()->SetTitle("correction [#mum]");
290
291
292
293
        align_correction_shiftX[detector->name()]->Write(name.c_str());

        name = "alignment_correction_displacementY_" + detector->name();
        align_correction_shiftY[detector->name()] =
294
            new TGraph(static_cast<int>(shiftsY[detector->name()].size()), &iterations[0], &shiftsY[detector->name()][0]);
295
296
        align_correction_shiftY[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_shiftY[detector->name()]->GetYaxis()->SetTitle("correction [#mum]");
297
298
299
300
        align_correction_shiftY[detector->name()]->Write(name.c_str());

        name = "alignment_correction_displacementZ_" + detector->name();
        align_correction_shiftZ[detector->name()] =
301
            new TGraph(static_cast<int>(shiftsZ[detector->name()].size()), &iterations[0], &shiftsZ[detector->name()][0]);
302
303
        align_correction_shiftZ[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_shiftZ[detector->name()]->GetYaxis()->SetTitle("correction [#mum]");
304
305
306
307
        align_correction_shiftZ[detector->name()]->Write(name.c_str());

        name = "alignment_correction_rotationX_" + detector->name();
        align_correction_rotX[detector->name()] =
308
            new TGraph(static_cast<int>(rotX[detector->name()].size()), &iterations[0], &rotX[detector->name()][0]);
309
310
        align_correction_rotX[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_rotX[detector->name()]->GetYaxis()->SetTitle("correction [deg]");
311
312
313
314
        align_correction_rotX[detector->name()]->Write(name.c_str());

        name = "alignment_correction_rotationY_" + detector->name();
        align_correction_rotY[detector->name()] =
315
            new TGraph(static_cast<int>(rotY[detector->name()].size()), &iterations[0], &rotY[detector->name()][0]);
316
317
        align_correction_rotY[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_rotY[detector->name()]->GetYaxis()->SetTitle("correction [deg]");
318
319
320
321
        align_correction_rotY[detector->name()]->Write(name.c_str());

        name = "alignment_correction_rotationZ_" + detector->name();
        align_correction_rotZ[detector->name()] =
322
            new TGraph(static_cast<int>(rotZ[detector->name()].size()), &iterations[0], &rotZ[detector->name()][0]);
323
324
        align_correction_rotZ[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_rotZ[detector->name()]->GetYaxis()->SetTitle("correction [deg]");
325
        align_correction_rotZ[detector->name()]->Write(name.c_str());
Daniel Hynds's avatar
Daniel Hynds committed
326
    }
327
}