AlignmentTrackChi2.cpp 15.2 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>("align_position", true);
21
22
23
    if(m_alignPosition) {
        LOG(INFO) << "Aligning positions";
    }
24
    m_alignOrientation = m_config.get<bool>("align_orientation", 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
    globalDetector->displacement(XYZPoint(par[detNum * 6 + 0], par[detNum * 6 + 1], par[detNum * 6 + 2]));
    globalDetector->rotation(XYZVector(par[detNum * 6 + 3], par[detNum * 6 + 4], par[detNum * 6 + 5]));
Simon Spannagel's avatar
Simon Spannagel committed
93
94

    // Apply new alignment conditions
95
    globalDetector->update();
Simon Spannagel's avatar
Simon Spannagel committed
96
97
98
99
100

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

    // Loop over all tracks
101
    for(size_t iTrack = 0; iTrack < globalTracks.size(); iTrack++) {
Simon Spannagel's avatar
Simon Spannagel committed
102
103
104
105
106
        // 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
107
        for(size_t iTrackCluster = 0; iTrackCluster < trackClusters.size(); iTrackCluster++) {
Simon Spannagel's avatar
Simon Spannagel committed
108
            Cluster* trackCluster = trackClusters[iTrackCluster];
109
            if(globalDetector->name() != trackCluster->detectorID()) {
Simon Spannagel's avatar
Simon Spannagel committed
110
                continue;
111
            }
112

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

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

        // Add the new chi2
        result += track->chi2();
124
125
126
    }
}

Daniel Hynds's avatar
Daniel Hynds committed
127
128
129
130
// ==================================================================
//  The finalise function - effectively the brains of the alignment!
// ==================================================================

131
void AlignmentTrackChi2::finalise() {
Simon Spannagel's avatar
Simon Spannagel committed
132
133
134
135
136

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

    // Make the fitting object
137
    TVirtualFitter* residualFitter = TVirtualFitter::Fitter(nullptr, 50);
138
    residualFitter->SetFCN(MinimiseTrackChi2);
Simon Spannagel's avatar
Simon Spannagel committed
139
140
141
142
143
144

    // Set the global parameters
    globalTracks = m_alignmenttracks;

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

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

152
153
154
155
156
157
158
159
    // 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;

160
161
    // 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.
162
    for(size_t iteration = 0; iteration < nIterations; iteration++) {
Simon Spannagel's avatar
Simon Spannagel committed
163

Simon Spannagel's avatar
Simon Spannagel committed
164
165
166
        LOG_PROGRESS(STATUS, "alignment_track") << "Alignment iteration " << (iteration + 1) << " of " << nIterations;

        int det = 0;
167
168
        for(auto& detector : get_detectors()) {
            string detectorID = detector->name();
169

Simon Spannagel's avatar
Simon Spannagel committed
170
            // Do not align the reference plane
171
            if(detector->isReference() || detector->isDUT()) {
Simon Spannagel's avatar
Simon Spannagel committed
172
                continue;
173
174
            }

Simon Spannagel's avatar
Simon Spannagel committed
175
            // Say that this is the detector we align
176
177
            globalDetector = detector;

Simon Spannagel's avatar
Simon Spannagel committed
178
179
            detNum = det;
            // Add the parameters to the fitter (z displacement not allowed to move!)
180
181
182
183
184
185
186
187
188
189
190
191
            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);
            }
192
            residualFitter->SetParameter(
193
                det * 6 + 2, (detectorID + "_displacementZ").c_str(), detector->displacement().Z(), 0, -10, 500);
Simon Spannagel's avatar
Simon Spannagel committed
194

195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
            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);
            }
210
211
212
            auto old_position = detector->displacement();
            auto old_orientation = detector->rotation();

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

216
217
218
219
220
221
222
223
            // 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);

224
            // Store corrections:
225
226
227
228
229
230
231
232
233
234
235
236
            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")));
237
238

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

Simon Spannagel's avatar
Simon Spannagel committed
242
243
            // Now that this device is fitted, set parameter errors to 0 so that they
            // are not fitted again
244
245
246
247
248
249
            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
250
251
252

            // Set the alignment parameters of this plane to be the optimised values
            // from the alignment
253
254
            detector->displacement(XYZPoint(displacementX, displacementY, displacementZ));
            detector->rotation(XYZVector(rotationX, rotationY, rotationZ));
255
            detector->update();
Simon Spannagel's avatar
Simon Spannagel committed
256
257
258
            det++;
        }
    }
259

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

Simon Spannagel's avatar
Simon Spannagel committed
262
    // Now list the new alignment parameters
263
    for(auto& detector : get_detectors()) {
Simon Spannagel's avatar
Simon Spannagel committed
264
        // Do not align the reference plane
265
        if(detector->isReference() || detector->isDUT()) {
Simon Spannagel's avatar
Simon Spannagel committed
266
267
            continue;
        }
Simon Spannagel's avatar
Simon Spannagel committed
268

269
        LOG(STATUS) << detector->name() << " new alignment: " << std::endl
270
271
                    << "T" << Units::display(detector->displacement(), {"mm", "um"}) << " R"
                    << Units::display(detector->rotation(), {"deg"});
272
273
274
275
276
277
278

        // 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()] =
279
            new TGraph(static_cast<int>(shiftsX[detector->name()].size()), &iterations[0], &shiftsX[detector->name()][0]);
280
281
        align_correction_shiftX[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_shiftX[detector->name()]->GetYaxis()->SetTitle("correction [#mum]");
282
283
284
285
        align_correction_shiftX[detector->name()]->Write(name.c_str());

        name = "alignment_correction_displacementY_" + detector->name();
        align_correction_shiftY[detector->name()] =
286
            new TGraph(static_cast<int>(shiftsY[detector->name()].size()), &iterations[0], &shiftsY[detector->name()][0]);
287
288
        align_correction_shiftY[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_shiftY[detector->name()]->GetYaxis()->SetTitle("correction [#mum]");
289
290
291
292
        align_correction_shiftY[detector->name()]->Write(name.c_str());

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

        name = "alignment_correction_rotationX_" + detector->name();
        align_correction_rotX[detector->name()] =
300
            new TGraph(static_cast<int>(rotX[detector->name()].size()), &iterations[0], &rotX[detector->name()][0]);
301
302
        align_correction_rotX[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_rotX[detector->name()]->GetYaxis()->SetTitle("correction [deg]");
303
304
305
306
        align_correction_rotX[detector->name()]->Write(name.c_str());

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

        name = "alignment_correction_rotationZ_" + detector->name();
        align_correction_rotZ[detector->name()] =
314
            new TGraph(static_cast<int>(rotZ[detector->name()].size()), &iterations[0], &rotZ[detector->name()][0]);
315
316
        align_correction_rotZ[detector->name()]->GetXaxis()->SetTitle("# iteration");
        align_correction_rotZ[detector->name()]->GetYaxis()->SetTitle("correction [deg]");
317
        align_correction_rotZ[detector->name()]->Write(name.c_str());
Daniel Hynds's avatar
Daniel Hynds committed
318
    }
319
}