Project
Loading...
Searching...
No Matches
MchAlignment.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
12//-----------------------------------------------------------------------------
22//-----------------------------------------------------------------------------
23
24#include "MCHAlign/Alignment.h"
25#include "MCHAlign/MillePede2.h"
26#include "MCHAlign/MillePedeRecord.h"
27#include <iostream>
28
29#include "MCHTracking/Track.h"
30#include "MCHTracking/TrackParam.h"
31#include "MCHTracking/Cluster.h"
32#include "TGeoManager.h"
33
34// #include "DataFormatsMCH/ROFRecord.h"
35// #include "DataFormatsMCH/TrackMCH.h"
36// #include "DataFormatsMCH/Cluster.h"
37// #include "DataFormatsMCH/Digit.h"
38
39// #include "AliMUONGeometryTransformer.h"
40// #include "AliMUONGeometryModuleTransformer.h"
41// #include "MCHAlign/AliMUONGeometryDetElement.h"
42// #include "AliMUONGeometryBuilder.h"
43#include "MCHGeometryCreator/Geometry.h"
44#include "MCHGeometryTest/Helpers.h"
45#include "MCHGeometryTransformer/Transformations.h"
46#include "TGeoManager.h"
47
48// #include "Align/Millepede2Record.h" //to be replaced
49// #include "AliMpExMap.h"
50// #include "AliMpExMapIterator.h"
51
52#include "DetectorsCommonDataFormats/AlignParam.h"
53#include "Framework/Logger.h"
54
55#include <TMath.h>
56#include <TMatrixDSym.h>
57#include <TMatrixD.h>
58#include <TClonesArray.h>
59#include <TGraphErrors.h>
60#include <TObject.h>
61
62namespace o2
63{
64namespace mch
65{
66
67using namespace std;
68
69//_____________________________________________________________________
70// static variables
71const Int_t Alignment::fgNDetElemCh[Alignment::fgNCh] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26};
72const Int_t Alignment::fgSNDetElemCh[Alignment::fgNCh + 1] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156};
73
74// number of detector elements in each half-chamber
75const Int_t Alignment::fgNDetElemHalfCh[Alignment::fgNHalfCh] = {2, 2, 2, 2, 2, 2, 2, 2, 9, 9, 9, 9, 13, 13, 13, 13, 13, 13, 13, 13};
76
77// list of detector elements for each half chamber
78const Int_t Alignment::fgDetElemHalfCh[Alignment::fgNHalfCh][Alignment::fgNDetHalfChMax] =
79 {
80 {100, 103, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
81 {101, 102, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
82
83 {200, 203, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
84 {201, 202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
85
86 {300, 303, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
87 {301, 302, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
88
89 {400, 403, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
90 {401, 402, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
91
92 {500, 501, 502, 503, 504, 514, 515, 516, 517, 0, 0, 0, 0},
93 {505, 506, 507, 508, 509, 510, 511, 512, 513, 0, 0, 0, 0},
94
95 {600, 601, 602, 603, 604, 614, 615, 616, 617, 0, 0, 0, 0},
96 {605, 606, 607, 608, 609, 610, 611, 612, 613, 0, 0, 0, 0},
97
98 {700, 701, 702, 703, 704, 705, 706, 720, 721, 722, 723, 724, 725},
99 {707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719},
100
101 {800, 801, 802, 803, 804, 805, 806, 820, 821, 822, 823, 824, 825},
102 {807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819},
103
104 {900, 901, 902, 903, 904, 905, 906, 920, 921, 922, 923, 924, 925},
105 {907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919},
106
107 {1000, 1001, 1002, 1003, 1004, 1005, 1006, 1020, 1021, 1022, 1023, 1024, 1025},
108 {1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019}
109
110};
111
112//_____________________________________________________________________
114class Array
115{
116
117 public:
119 Array(void)
120 {
121 for (Int_t i = 0; i < Alignment::fNGlobal; ++i) {
122 values[i] = 0;
123 }
124 }
125
127 Double_t values[Alignment::fNGlobal];
128
129 private:
131 Array(const Array&);
132
134 Array& operator=(const Array&);
135};
136
137//________________________________________________________________________
138Double_t Square(Double_t x) { return x * x; }
139
140//_____________________________________________________________________
141Alignment::Alignment()
142 : TObject(),
143 fInitialized(kFALSE),
144 fRunNumber(0),
145 fBFieldOn(kFALSE),
146 fRefitStraightTracks(kFALSE),
147 fStartFac(256),
148 fResCutInitial(100),
149 fResCut(100),
150 fMillepede(0L), // to be modified
151 fCluster(0L),
152 fNStdDev(3),
153 fDetElemNumber(0),
154 fTrackRecord(),
155 fTransformCreator(),
156 fGeoCombiTransInverse(),
157 fDoEvaluation(kFALSE),
158 fTrackParamOrig(0),
159 fTrackParamNew(0),
160 fTFile(0),
161 fTTree(0)
162{
164 fSigma[0] = 1.5e-1;
165 fSigma[1] = 1.0e-2;
166
167 // default allowed variations
168 fAllowVar[0] = 0.5; // x
169 fAllowVar[1] = 0.5; // y
170 fAllowVar[2] = 0.01; // phi_z
171 fAllowVar[3] = 5; // z
172
173 // initialize millepede
174 fMillepede = new MillePede2();
175 // fMillepede = new o2::align::Mille("theMilleFile.txt"); // To be replaced by MillePede2
176
177 // initialize degrees of freedom
178 // by default all parameters are free
179 for (Int_t iPar = 0; iPar < fNGlobal; ++iPar) {
180 fGlobalParameterStatus[iPar] = kFreeParId;
181 }
182
183 // initialize local equations
184 for (int i = 0; i < fNLocal; ++i) {
185 fLocalDerivatives[i] = 0.0;
186 }
187
188 for (int i = 0; i < fNGlobal; ++i) {
189 fGlobalDerivatives[i] = 0.0;
190 }
191}
192
193//_____________________________________________________________________
194// Alignment::~Alignment()
195//{
196// /// destructor
197//}
198// Alignment::~Alignment() = default;
199//_____________________________________________________________________
200void Alignment::init(void)
201{
202
204
209 if (fInitialized) {
210 LOG(fatal) << "Millepede already initialized";
211 }
212
213 // assign proper groupID to free parameters
214 Int_t nGlobal = 0;
215 for (Int_t iPar = 0; iPar < fNGlobal; ++iPar) {
216
217 if (fGlobalParameterStatus[iPar] == kFixedParId) {
218 // fixed parameters are left unchanged
219 continue;
220
221 } else if (fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId) {
222
223 // free parameters or first element of group are assigned a new group id
224 fGlobalParameterStatus[iPar] = nGlobal++;
225 continue;
226
227 } else if (fGlobalParameterStatus[iPar] < kGroupBaseId) {
228
229 // get detector element id from status, get chamber parameter id
230 const Int_t iDeBase(kGroupBaseId - 1 - fGlobalParameterStatus[iPar]);
231 const Int_t iParBase = iPar % fgNParCh;
232
233 // check
234 if (iDeBase < 0 || iDeBase >= iPar / fgNParCh) {
235 LOG(fatal) << "Group for parameter index " << iPar << " has wrong base detector element: " << iDeBase;
236 }
237
238 // assign identical group id to current
239 fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase * fgNParCh + iParBase];
240 LOG(info) << "Parameter " << iPar << " grouped to detector " << iDeBase << " (" << GetParameterMaskString(1 << iParBase).Data() << ")";
241
242 } else
243 LOG(fatal) << "Unrecognized parameter status for index " << iPar << ": " << fGlobalParameterStatus[iPar];
244 }
245
246 LOG(info) << "Free Parameters: " << nGlobal << " out of " << fNGlobal;
247
248 // initialize millepede
249 // fMillepede->InitMille(fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial, fGlobalParameterStatus);
250 fMillepede->InitMille(fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial); // MillePede2 implementation
251
252 fInitialized = kTRUE;
253
254 // some debug output
255 for (Int_t iPar = 0; iPar < fgNParCh; ++iPar) {
256 LOG(info) << "fAllowVar[" << iPar << "]= " << fAllowVar[iPar];
257 }
258
259 // set allowed variations for all parameters
260 for (Int_t iDet = 0; iDet < fgNDetElem; ++iDet) {
261 for (Int_t iPar = 0; iPar < fgNParCh; ++iPar) {
262 fMillepede->SetParSigma(iDet * fgNParCh + iPar, fAllowVar[iPar]);
263 }
264 }
265
266 // Set iterations
267 if (fStartFac > 1) {
268 fMillepede->SetIterations(fStartFac);
269 }
270 // setup monitoring TFile
271 if (fDoEvaluation && fRefitStraightTracks) {
272 fTFile = new TFile("Alignment.root", "RECREATE");
273 fTTree = new TTree("TreeE", "Evaluation");
274
275 const Int_t kSplitlevel = 98;
276 const Int_t kBufsize = 32000;
277
278 fTrackParamOrig = new LocalTrackParam();
279 fTTree->Branch("fTrackParamOrig", "LocalTrackParam", &fTrackParamOrig, kBufsize, kSplitlevel);
280
281 fTrackParamNew = new LocalTrackParam();
282 fTTree->Branch("fTrackParamNew", "LocalTrackParam", &fTrackParamNew, kBufsize, kSplitlevel);
283 }
284}
285
286//_____________________________________________________
287void Alignment::terminate(void)
288{
289 LOG(info) << "Closing Evaluation TFile";
290 if (fTFile && fTTree) {
291 fTFile->cd();
292 fTTree->Write();
293 fTFile->Close();
294 }
295}
296
297//_____________________________________________________
298MillePedeRecord* Alignment::ProcessTrack(Track& track, Bool_t doAlignment, Double_t weight)
299{
301
306 // reset track records
307 fTrackRecord.Reset();
308 if (fMillepede->GetRecord()) {
309 fMillepede->GetRecord()->Reset();
310 }
311
312 // loop over clusters to get starting values
313 Bool_t first(kTRUE);
314 // if (!trackParam)
315 // continue;
316 for (auto itTrackParam(track.begin()); itTrackParam != track.end(); ++itTrackParam) {
317
318 // get cluster
319 const Cluster* Cluster = itTrackParam->getClusterPtr();
320 if (!cluster)
321 continue;
322
323 // for first valid cluster, save track position as "starting" values
324 if (first) {
325
326 first = kFALSE;
327 FillTrackParamData(&*itTrackParam);
328 fTrackPos0[0] = fTrackPos[0];
329 fTrackPos0[1] = fTrackPos[1];
330 fTrackPos0[2] = fTrackPos[2];
331 fTrackSlope0[0] = fTrackSlope[0];
332 fTrackSlope0[1] = fTrackSlope[1];
333
334 break;
335 }
336 }
337
338 // redo straight track fit
339 if (fRefitStraightTracks) {
340
341 // refit straight track
342 const LocalTrackParam trackParam(RefitStraightTrack(track, fTrackPos0[2]));
343
344 // fill evaluation tree
345 if (fTrackParamOrig) {
346 fTrackParamOrig->fTrackX = fTrackPos0[0];
347 fTrackParamOrig->fTrackY = fTrackPos0[1];
348 fTrackParamOrig->fTrackZ = fTrackPos0[2];
349 fTrackParamOrig->fTrackSlopeX = fTrackSlope[0];
350 fTrackParamOrig->fTrackSlopeY = fTrackSlope[1];
351 }
352
353 // new ones
354 if (fTrackParamNew) {
355 fTrackParamNew->fTrackX = trackParam.fTrackX;
356 fTrackParamNew->fTrackY = trackParam.fTrackY;
357 fTrackParamNew->fTrackZ = trackParam.fTrackZ;
358 fTrackParamNew->fTrackSlopeX = trackParam.fTrackSlopeX;
359 fTrackParamNew->fTrackSlopeY = trackParam.fTrackSlopeY;
360 }
361
362 if (fTTree)
363 fTTree->Fill();
364
365 /*
366 copy new parameters to stored ones for derivatives calculation
367 this is done only if BFieldOn is false, for which these parameters are used
368 */
369 if (!fBFieldOn) {
370 fTrackPos0[0] = trackParam.fTrackX;
371 fTrackPos0[1] = trackParam.fTrackY;
372 fTrackPos0[2] = trackParam.fTrackZ;
373 fTrackSlope[0] = trackParam.fTrackSlopeX;
374 fTrackSlope[1] = trackParam.fTrackSlopeY;
375 }
376 }
377
378 // second loop to perform alignment
379 for (auto itTrackParam(track.begin()); itTrackParam != track.end(); ++itTrackParam) {
380
381 // get track parameters
382 if (!&*itTrackParam)
383 continue;
384
385 // get cluster
386 const Cluster* cluster = itTrackParam->getClusterPtr();
387 if (!cluster)
388 continue;
389
390 // fill local variables for this position --> one measurement
391 FillDetElemData(cluster);
392 FillRecPointData(cluster);
393 FillTrackParamData(&*itTrackParam);
394
395 // 'inverse' (GlobalToLocal) rotation matrix
396 const Double_t* r(fGeoCombiTransInverse.GetRotationMatrix());
397
398 // calculate measurements
399 if (fBFieldOn) {
400
401 // use residuals (cluster - track) for measurement
402 fMeas[0] = r[0] * (fClustPos[0] - fTrackPos[0]) + r[1] * (fClustPos[1] - fTrackPos[1]);
403 fMeas[1] = r[3] * (fClustPos[0] - fTrackPos[0]) + r[4] * (fClustPos[1] - fTrackPos[1]);
404
405 } else {
406
407 // use cluster position for measurement
408 fMeas[0] = (r[0] * fClustPos[0] + r[1] * fClustPos[1]);
409 fMeas[1] = (r[3] * fClustPos[0] + r[4] * fClustPos[1]);
410 }
411
412 // Set local equations
413 LocalEquationX();
414 LocalEquationY();
415 }
416
417 // copy track record
418 fMillepede->SetRecordRun(fRunNumber);
419 fMillepede->SetRecordWeight(weight);
420 fTrackRecord = *fMillepede->GetRecord();
421
422 // save record data
423 if (doAlignment) {
424 fMillepede->SaveRecordData();
425 fMillepede->CloseDataRecStorage();
426 }
427
428 // return record
429 return &fTrackRecord;
430}
431
432//______________________________________________________________________________
433void Alignment::ProcessTrack(MillePedeRecord* trackRecord)
434{
435 LOG(fatal) << __PRETTY_FUNCTION__ << " is disabled";
436
438 if (!trackRecord)
439 return;
440
441 // // make sure record storage is initialized
442 if (!fMillepede->GetRecord()) {
443 fMillepede->InitDataRecStorage(kFalse);
444 }
445 // // copy content
446 *fMillepede->GetRecord() = *trackRecord;
447
448 // save record
449 fMillepede->SaveRecordData();
450 // write to local file
451 fMillepede->CloseDataRecStorage();
452
453 return;
454}
455
456//_____________________________________________________________________
457void Alignment::FixAll(UInt_t mask)
458{
460 LOG(info) << "Fixing " << GetParameterMaskString(mask).Data() << " for all detector elements";
461
462 // fix all stations
463 for (Int_t i = 0; i < fgNDetElem; ++i) {
464 if (mask & ParX)
465 FixParameter(i, 0);
466 if (mask & ParY)
467 FixParameter(i, 1);
468 if (mask & ParTZ)
469 FixParameter(i, 2);
470 if (mask & ParZ)
471 FixParameter(i, 3);
472 }
473}
474
475//_____________________________________________________________________
476void Alignment::FixChamber(Int_t iCh, UInt_t mask)
477{
479
480 // check boundaries
481 if (iCh < 1 || iCh > 10) {
482 LOG(fatal) << "Invalid chamber index " << iCh;
483 }
484
485 // get first and last element
486 const Int_t iDetElemFirst = fgSNDetElemCh[iCh - 1];
487 const Int_t iDetElemLast = fgSNDetElemCh[iCh];
488 for (Int_t i = iDetElemFirst; i < iDetElemLast; ++i) {
489
490 LOG(info) << "Fixing " << GetParameterMaskString(mask).Data() << " for detector element " << i;
491
492 if (mask & ParX)
493 FixParameter(i, 0);
494 if (mask & ParY)
495 FixParameter(i, 1);
496 if (mask & ParTZ)
497 FixParameter(i, 2);
498 if (mask & ParZ)
499 FixParameter(i, 3);
500 }
501}
502
503//_____________________________________________________________________
504void Alignment::FixDetElem(Int_t iDetElemId, UInt_t mask)
505{
507 const Int_t iDet(GetDetElemNumber(iDetElemId));
508 if (mask & ParX)
509 FixParameter(iDet, 0);
510 if (mask & ParY)
511 FixParameter(iDet, 1);
512 if (mask & ParTZ)
513 FixParameter(iDet, 2);
514 if (mask & ParZ)
515 FixParameter(iDet, 3);
516}
517
518//_____________________________________________________________________
519void Alignment::FixHalfSpectrometer(const Bool_t* lChOnOff, UInt_t sidesMask, UInt_t mask)
520{
521
523 for (Int_t i = 0; i < fgNDetElem; ++i) {
524
525 // get chamber matching detector
526 const Int_t iCh(GetChamberId(i));
527 if (!lChOnOff[iCh - 1])
528 continue;
529
530 // get detector element in chamber
531 Int_t lDetElemNumber = i - fgSNDetElemCh[iCh - 1];
532
533 // skip detector if its side is off
534 // stations 1 and 2
535 if (iCh >= 1 && iCh <= 4) {
536 if (lDetElemNumber == 0 && !(sidesMask & SideTopRight))
537 continue;
538 if (lDetElemNumber == 1 && !(sidesMask & SideTopLeft))
539 continue;
540 if (lDetElemNumber == 2 && !(sidesMask & SideBottomLeft))
541 continue;
542 if (lDetElemNumber == 3 && !(sidesMask & SideBottomRight))
543 continue;
544 }
545
546 // station 3
547 if (iCh >= 5 && iCh <= 6) {
548 if (lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(sidesMask & SideTopRight))
549 continue;
550 if (lDetElemNumber >= 5 && lDetElemNumber <= 10 && !(sidesMask & SideTopLeft))
551 continue;
552 if (lDetElemNumber >= 11 && lDetElemNumber <= 13 && !(sidesMask & SideBottomLeft))
553 continue;
554 if (lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(sidesMask & SideBottomRight))
555 continue;
556 }
557
558 // stations 4 and 5
559 if (iCh >= 7 && iCh <= 10) {
560 if (lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(sidesMask & SideTopRight))
561 continue;
562 if (lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(sidesMask & SideTopLeft))
563 continue;
564 if (lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(sidesMask & SideBottomLeft))
565 continue;
566 if (lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(sidesMask & SideBottomRight))
567 continue;
568 }
569
570 // detector is accepted, fix it
571 FixDetElem(i, mask);
572 }
573}
574
575//______________________________________________________________________
576void Alignment::FixParameter(Int_t iPar)
577{
578
580 if (fInitialized) {
581 LOG(fatal) << "Millepede already initialized";
582 }
583
584 fGlobalParameterStatus[iPar] = kFixedParId;
585}
586
587//_____________________________________________________________________
588void Alignment::ReleaseChamber(Int_t iCh, UInt_t mask)
589{
591
592 // check boundaries
593 if (iCh < 1 || iCh > 10) {
594 LOG(fatal) << "Invalid chamber index " << iCh;
595 }
596
597 // get first and last element
598 const Int_t iDetElemFirst = fgSNDetElemCh[iCh - 1];
599 const Int_t iDetElemLast = fgSNDetElemCh[iCh];
600 for (Int_t i = iDetElemFirst; i < iDetElemLast; ++i) {
601
602 LOG(info) << "Releasing " << GetParameterMaskString(mask).Data() << " for detector element " << i;
603
604 if (mask & ParX)
605 ReleaseParameter(i, 0);
606 if (mask & ParY)
607 ReleaseParameter(i, 1);
608 if (mask & ParTZ)
609 ReleaseParameter(i, 2);
610 if (mask & ParZ)
611 ReleaseParameter(i, 3);
612 }
613}
614
615//_____________________________________________________________________
616void Alignment::ReleaseDetElem(Int_t iDetElemId, UInt_t mask)
617{
619 const Int_t iDet(GetDetElemNumber(iDetElemId));
620 if (mask & ParX)
621 ReleaseParameter(iDet, 0);
622 if (mask & ParY)
623 ReleaseParameter(iDet, 1);
624 if (mask & ParTZ)
625 ReleaseParameter(iDet, 2);
626 if (mask & ParZ)
627 ReleaseParameter(iDet, 3);
628}
629
630//______________________________________________________________________
631void Alignment::ReleaseParameter(Int_t iPar)
632{
633
635 if (fInitialized) {
636 LOG(fatal) << "Millepede already initialized";
637 }
638
639 fGlobalParameterStatus[iPar] = kFreeParId;
640}
641
642//_____________________________________________________________________
643void Alignment::GroupChamber(Int_t iCh, UInt_t mask)
644{
646 if (iCh < 1 || iCh > fgNCh) {
647 LOG(fatal) << "Invalid chamber index " << iCh;
648 }
649
650 const Int_t detElemMin = 100 * iCh;
651 const Int_t detElemMax = 100 * iCh + fgNDetElemCh[iCh] - 1;
652 GroupDetElems(detElemMin, detElemMax, mask);
653}
654
655//_____________________________________________________________________
656void Alignment::GroupHalfChamber(Int_t iCh, Int_t iHalf, UInt_t mask)
657{
659 if (iCh < 1 || iCh > fgNCh) {
660 LOG(fatal) << "Invalid chamber index " << iCh;
661 }
662
663 if (iHalf < 0 || iHalf > 1) {
664 LOG(fatal) << "Invalid half chamber index " << iHalf;
665 }
666
667 const Int_t iHalfCh = 2 * (iCh - 1) + iHalf;
668 GroupDetElems(&fgDetElemHalfCh[iHalfCh][0], fgNDetElemHalfCh[iHalfCh], mask);
669}
670
671//_____________________________________________________________________
672void Alignment::GroupDetElems(Int_t detElemMin, Int_t detElemMax, UInt_t mask)
673{
675 // check number of detector elements
676 const Int_t nDetElem = detElemMax - detElemMin + 1;
677 if (nDetElem < 2) {
678 LOG(fatal) << "Requested group of DEs " << detElemMin << "-" << detElemMax << " contains less than 2 DE's";
679 }
680
681 // create list
682 Int_t* detElemList = new int[nDetElem];
683 for (Int_t i = 0; i < nDetElem; ++i) {
684 detElemList[i] = detElemMin + i;
685 }
686
687 // group
688 GroupDetElems(detElemList, nDetElem, mask);
689 delete[] detElemList;
690}
691
692//_____________________________________________________________________
693void Alignment::GroupDetElems(const Int_t* detElemList, Int_t nDetElem, UInt_t mask)
694{
696 if (fInitialized) {
697 LOG(fatal) << "Millepede already initialized";
698 }
699
700 const Int_t iDeBase(GetDetElemNumber(detElemList[0]));
701 for (Int_t i = 0; i < nDetElem; ++i) {
702 const Int_t iDeCurrent(GetDetElemNumber(detElemList[i]));
703 if (mask & ParX)
704 fGlobalParameterStatus[iDeCurrent * fgNParCh + 0] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1);
705 if (mask & ParY)
706 fGlobalParameterStatus[iDeCurrent * fgNParCh + 1] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1);
707 if (mask & ParTZ)
708 fGlobalParameterStatus[iDeCurrent * fgNParCh + 2] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1);
709 if (mask & ParZ)
710 fGlobalParameterStatus[iDeCurrent * fgNParCh + 3] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1);
711
712 if (i == 0)
713 LOG(info) << "Creating new group for detector " << detElemList[i] << " and variable " << GetParameterMaskString(mask).Data();
714 else
715 LOG(info) << "Adding detector element " << detElemList[i] << " to current group";
716 }
717}
718
719//______________________________________________________________________
720void Alignment::SetChamberNonLinear(Int_t iCh, UInt_t mask)
721{
723 const Int_t iDetElemFirst = fgSNDetElemCh[iCh - 1];
724 const Int_t iDetElemLast = fgSNDetElemCh[iCh];
725 for (Int_t i = iDetElemFirst; i < iDetElemLast; ++i) {
726
727 if (mask & ParX)
728 SetParameterNonLinear(i, 0);
729 if (mask & ParY)
730 SetParameterNonLinear(i, 1);
731 if (mask & ParTZ)
732 SetParameterNonLinear(i, 2);
733 if (mask & ParZ)
734 SetParameterNonLinear(i, 3);
735 }
736}
737
738//_____________________________________________________________________
739void Alignment::SetDetElemNonLinear(Int_t iDetElemId, UInt_t mask)
740{
742 const Int_t iDet(GetDetElemNumber(iDetElemId));
743 if (mask & ParX)
744 SetParameterNonLinear(iDet, 0);
745 if (mask & ParY)
746 SetParameterNonLinear(iDet, 1);
747 if (mask & ParTZ)
748 SetParameterNonLinear(iDet, 2);
749 if (mask & ParZ)
750 SetParameterNonLinear(iDet, 3);
751}
752
753//______________________________________________________________________
754void Alignment::SetParameterNonLinear(Int_t iPar)
755{
757 if (!fInitialized) {
758 LOG(fatal) << "Millepede not initialized";
759 }
760
761 fMillepede->SetNonLinear(iPar);
762 LOG(info) << "Parameter " << iPar << " set to non linear ";
763}
764
765//______________________________________________________________________
766void Alignment::AddConstraints(const Bool_t* lChOnOff, UInt_t mask)
767{
769
770 Array fConstraintX;
771 Array fConstraintY;
772 Array fConstraintTZ;
773 Array fConstraintZ;
774
775 for (Int_t i = 0; i < fgNDetElem; ++i) {
776
777 // get chamber matching detector
778 const Int_t iCh(GetChamberId(i));
779 if (lChOnOff[iCh - 1]) {
780
781 if (mask & ParX)
782 fConstraintX.values[i * fgNParCh + 0] = 1.0;
783 if (mask & ParY)
784 fConstraintY.values[i * fgNParCh + 1] = 1.0;
785 if (mask & ParTZ)
786 fConstraintTZ.values[i * fgNParCh + 2] = 1.0;
787 if (mask & ParZ)
788 fConstraintTZ.values[i * fgNParCh + 3] = 1.0;
789 }
790 }
791
792 if (mask & ParX)
793 AddConstraint(fConstraintX.values, 0.0);
794 if (mask & ParY)
795 AddConstraint(fConstraintY.values, 0.0);
796 if (mask & ParTZ)
797 AddConstraint(fConstraintTZ.values, 0.0);
798 if (mask & ParZ)
799 AddConstraint(fConstraintZ.values, 0.0);
800}
801
802//______________________________________________________________________
803void Alignment::AddConstraints(const Bool_t* lChOnOff, const Bool_t* lVarXYT, UInt_t sidesMask)
804{
805 /*
806 questions:
807 - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ?
808 - why is weight ignored for ConstrainT and ConstrainB
809 - why is there no constrain on z
810 */
811
813 Double_t lMeanY = 0.;
814 Double_t lSigmaY = 0.;
815 Double_t lMeanZ = 0.;
816 Double_t lSigmaZ = 0.;
817 Int_t lNDetElem = 0;
818
819 for (Int_t i = 0; i < fgNDetElem; ++i) {
820
821 // get chamber matching detector
822 const Int_t iCh(GetChamberId(i));
823
824 // skip detector if chamber is off
825 if (lChOnOff[iCh - 1])
826 continue;
827
828 // get detector element id from detector element number
829 const Int_t lDetElemNumber = i - fgSNDetElemCh[iCh - 1];
830 const Int_t lDetElemId = iCh * 100 + lDetElemNumber;
831
832 // skip detector if its side is off
833 // stations 1 and 2
834 if (iCh >= 1 && iCh <= 4) {
835 if (lDetElemNumber == 0 && !(sidesMask & SideTopRight))
836 continue;
837 if (lDetElemNumber == 1 && !(sidesMask & SideTopLeft))
838 continue;
839 if (lDetElemNumber == 2 && !(sidesMask & SideBottomLeft))
840 continue;
841 if (lDetElemNumber == 3 && !(sidesMask & SideBottomRight))
842 continue;
843 }
844
845 // station 3
846 if (iCh >= 5 && iCh <= 6) {
847 if (lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(sidesMask & SideTopRight))
848 continue;
849 if (lDetElemNumber >= 5 && lDetElemNumber <= 10 && !(sidesMask & SideTopLeft))
850 continue;
851 if (lDetElemNumber >= 11 && lDetElemNumber <= 13 && !(sidesMask & SideBottomLeft))
852 continue;
853 if (lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(sidesMask & SideBottomRight))
854 continue;
855 }
856
857 // stations 4 and 5
858 if (iCh >= 7 && iCh <= 10) {
859 if (lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(sidesMask & SideTopRight))
860 continue;
861 if (lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(sidesMask & SideTopLeft))
862 continue;
863 if (lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(sidesMask & SideBottomLeft))
864 continue;
865 if (lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(sidesMask & SideBottomRight))
866 continue;
867 }
868
869 // get global x, y and z position
870 Double_t lDetElemGloX = 0.;
871 Double_t lDetElemGloY = 0.;
872 Double_t lDetElemGloZ = 0.;
873
874 auto fTransform = fTransformCreator(lDetElemId);
875 o2::math_utils::Point3D<double> SlatPos{0.0, 0.0, 0.0};
876 o2::math_utils::Point3D<double> GlobalPos;
877
878 fTransform.LocalToMaster(SlatPos, GlobalPos);
879 lDetElemGloX = GlobalPos.x();
880 lDetElemGloY = GlobalPos.y();
881 lDetElemGloZ = GlobalPos.z();
882 // fTransform->Local2Global(lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ);
883
884 // increment mean Y, mean Z, sigmas and number of accepted detectors
885 lMeanY += lDetElemGloY;
886 lSigmaY += lDetElemGloY * lDetElemGloY;
887 lMeanZ += lDetElemGloZ;
888 lSigmaZ += lDetElemGloZ * lDetElemGloZ;
889 lNDetElem++;
890 }
891
892 // calculate mean values
893 lMeanY /= lNDetElem;
894 lSigmaY /= lNDetElem;
895 lSigmaY = TMath::Sqrt(lSigmaY - lMeanY * lMeanY);
896 lMeanZ /= lNDetElem;
897 lSigmaZ /= lNDetElem;
898 lSigmaZ = TMath::Sqrt(lSigmaZ - lMeanZ * lMeanZ);
899 LOG(info) << "Used " << lNDetElem << " DetElem, MeanZ= " << lMeanZ << ", SigmaZ= " << lSigmaZ;
900
901 // create all possible arrays
902 Array fConstraintX[4]; // Array for constraint equation X
903 Array fConstraintY[4]; // Array for constraint equation Y
904 Array fConstraintP[4]; // Array for constraint equation P
905 Array fConstraintXZ[4]; // Array for constraint equation X vs Z
906 Array fConstraintYZ[4]; // Array for constraint equation Y vs Z
907 Array fConstraintPZ[4]; // Array for constraint equation P vs Z
908
909 // do we really need these ?
910 Array fConstraintXY[4]; // Array for constraint equation X vs Y
911 Array fConstraintYY[4]; // Array for constraint equation Y vs Y
912 Array fConstraintPY[4]; // Array for constraint equation P vs Y
913
914 // fill Bool_t sides array based on masks, for convenience
915 Bool_t lDetTLBR[4];
916 lDetTLBR[0] = sidesMask & SideTop;
917 lDetTLBR[1] = sidesMask & SideLeft;
918 lDetTLBR[2] = sidesMask & SideBottom;
919 lDetTLBR[3] = sidesMask & SideRight;
920
921 for (Int_t i = 0; i < fgNDetElem; ++i) {
922
923 // get chamber matching detector
924 const Int_t iCh(GetChamberId(i));
925
926 // skip detector if chamber is off
927 if (!lChOnOff[iCh - 1])
928 continue;
929
930 // get detector element id from detector element number
931 const Int_t lDetElemNumber = i - fgSNDetElemCh[iCh - 1];
932 const Int_t lDetElemId = iCh * 100 + lDetElemNumber;
933
934 // get global x, y and z position
935 Double_t lDetElemGloX = 0.;
936 Double_t lDetElemGloY = 0.;
937 Double_t lDetElemGloZ = 0.;
938
939 auto fTransform = fTransformCreator(lDetElemId);
940 o2::math_utils::Point3D<double> SlatPos{0.0, 0.0, 0.0};
941 o2::math_utils::Point3D<double> GlobalPos;
942
943 fTransform.LocalToMaster(SlatPos, GlobalPos);
944 lDetElemGloX = GlobalPos.x();
945 lDetElemGloY = GlobalPos.y();
946 lDetElemGloZ = GlobalPos.z();
947 // fTransform->Local2Global(lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ);
948
949 // loop over sides
950 for (Int_t iSide = 0; iSide < 4; iSide++) {
951
952 // skip if side is not selected
953 if (!lDetTLBR[iSide])
954 continue;
955
956 // skip detector if it is not in the selected side
957 // stations 1 and 2
958 if (iCh >= 1 && iCh <= 4) {
959 if (lDetElemNumber == 0 && !(iSide == 0 || iSide == 3))
960 continue; // top-right
961 if (lDetElemNumber == 1 && !(iSide == 0 || iSide == 1))
962 continue; // top-left
963 if (lDetElemNumber == 2 && !(iSide == 2 || iSide == 1))
964 continue; // bottom-left
965 if (lDetElemNumber == 3 && !(iSide == 2 || iSide == 3))
966 continue; // bottom-right
967 }
968
969 // station 3
970 if (iCh >= 5 && iCh <= 6) {
971 if (lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3))
972 continue; // top-right
973 if (lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1))
974 continue; // top-left
975 if (lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1))
976 continue; // bottom-left
977 if (lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3))
978 continue; // bottom-right
979 }
980
981 // stations 4 and 5
982 if (iCh >= 7 && iCh <= 10) {
983 if (lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3))
984 continue; // top-right
985 if (lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1))
986 continue; // top-left
987 if (lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1))
988 continue; // bottom-left
989 if (lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3))
990 continue; // bottom-right
991 }
992
993 // constrain x
994 if (lVarXYT[0])
995 fConstraintX[iSide].values[i * fgNParCh + 0] = 1;
996
997 // constrain y
998 if (lVarXYT[1])
999 fConstraintY[iSide].values[i * fgNParCh + 1] = 1;
1000
1001 // constrain phi (rotation around z)
1002 if (lVarXYT[2])
1003 fConstraintP[iSide].values[i * fgNParCh + 2] = 1;
1004
1005 // x-z shearing
1006 if (lVarXYT[3])
1007 fConstraintXZ[iSide].values[i * fgNParCh + 0] = (lDetElemGloZ - lMeanZ) / lSigmaZ;
1008
1009 // y-z shearing
1010 if (lVarXYT[4])
1011 fConstraintYZ[iSide].values[i * fgNParCh + 1] = (lDetElemGloZ - lMeanZ) / lSigmaZ;
1012
1013 // phi-z shearing
1014 if (lVarXYT[5])
1015 fConstraintPZ[iSide].values[i * fgNParCh + 2] = (lDetElemGloZ - lMeanZ) / lSigmaZ;
1016
1017 // x-y shearing
1018 if (lVarXYT[6])
1019 fConstraintXY[iSide].values[i * fgNParCh + 0] = (lDetElemGloY - lMeanY) / lSigmaY;
1020
1021 // y-y shearing
1022 if (lVarXYT[7])
1023 fConstraintYY[iSide].values[i * fgNParCh + 1] = (lDetElemGloY - lMeanY) / lSigmaY;
1024
1025 // phi-y shearing
1026 if (lVarXYT[8])
1027 fConstraintPY[iSide].values[i * fgNParCh + 2] = (lDetElemGloY - lMeanY) / lSigmaY;
1028 }
1029 }
1030
1031 // pass constraints to millepede
1032 for (Int_t iSide = 0; iSide < 4; iSide++) {
1033 // skip if side is not selected
1034 if (!lDetTLBR[iSide])
1035 continue;
1036
1037 if (lVarXYT[0])
1038 AddConstraint(fConstraintX[iSide].values, 0.0);
1039 if (lVarXYT[1])
1040 AddConstraint(fConstraintY[iSide].values, 0.0);
1041 if (lVarXYT[2])
1042 AddConstraint(fConstraintP[iSide].values, 0.0);
1043 if (lVarXYT[3])
1044 AddConstraint(fConstraintXZ[iSide].values, 0.0);
1045 if (lVarXYT[4])
1046 AddConstraint(fConstraintYZ[iSide].values, 0.0);
1047 if (lVarXYT[5])
1048 AddConstraint(fConstraintPZ[iSide].values, 0.0);
1049 if (lVarXYT[6])
1050 AddConstraint(fConstraintXY[iSide].values, 0.0);
1051 if (lVarXYT[7])
1052 AddConstraint(fConstraintYY[iSide].values, 0.0);
1053 if (lVarXYT[8])
1054 AddConstraint(fConstraintPY[iSide].values, 0.0);
1055 }
1056}
1057
1058//______________________________________________________________________
1059void Alignment::InitGlobalParameters(Double_t* par)
1060{
1062 if (!fInitialized) {
1063 LOG(fatal) << "Millepede is not initialized";
1064 }
1065
1066 fMillepede->SetGlobalParameters(par);
1067}
1068
1069//______________________________________________________________________
1070void Alignment::SetAllowedVariation(Int_t iPar, Double_t value)
1071{
1073 // check initialization
1074 if (fInitialized) {
1075 LOG(fatal) << "Millepede already initialized";
1076 }
1077
1078 // check initialization
1079 if (!(iPar >= 0 && iPar < fgNParCh)) {
1080 LOG(fatal) << "Invalid index: " << iPar;
1081 }
1082
1083 fAllowVar[iPar] = value;
1084}
1085
1086//______________________________________________________________________
1087void Alignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
1088{
1089
1091 fSigma[0] = sigmaX;
1092 fSigma[1] = sigmaY;
1093
1094 // print
1095 for (Int_t i = 0; i < 2; ++i) {
1096 LOG(info) << "fSigma[" << i << "] =" << fSigma[i];
1097 }
1098}
1099
1100//_____________________________________________________
1101void Alignment::GlobalFit(Double_t* parameters, Double_t* errors, Double_t* pulls)
1102{
1103
1105 fMillepede->GlobalFit(parameters, errors, pulls);
1106
1107 LOG(info) << "Done fitting global parameters";
1108 for (int iDet = 0; iDet < fgNDetElem; ++iDet) {
1109 LOG(info) << iDet << " " << parameters[iDet * fgNParCh + 0] << " " << parameters[iDet * fgNParCh + 1] << " " << parameters[iDet * fgNParCh + 3] << " " << parameters[iDet * fgNParCh + 2];
1110 }
1111}
1112
1113//_____________________________________________________
1114void Alignment::PrintGlobalParameters() const
1115{
1116 fMillepede->PrintGlobalParameters();
1117}
1118
1119//_____________________________________________________
1120Double_t Alignment::GetParError(Int_t iPar) const
1121{
1122 return fMillepede->GetParError(iPar);
1123}
1124
1125// //______________________________________________________________________
1126// AliMUONGeometryTransformer* Alignment::ReAlign(
1127// const AliMUONGeometryTransformer* transformer,
1128// const double* misAlignments, Bool_t)
1129// {
1130
1131// /// Returns a new AliMUONGeometryTransformer with the found misalignments
1132// /// applied.
1133
1134// // Takes the internal geometry module transformers, copies them
1135// // and gets the Detection Elements from them.
1136// // Takes misalignment parameters and applies these
1137// // to the local transform of the Detection Element
1138// // Obtains the global transform by multiplying the module transformer
1139// // transformation with the local transformation
1140// // Applies the global transform to a new detection element
1141// // Adds the new detection element to a new module transformer
1142// // Adds the new module transformer to a new geometry transformer
1143// // Returns the new geometry transformer
1144
1145// Double_t lModuleMisAlignment[fgNParCh] = {0};
1146// Double_t lDetElemMisAlignment[fgNParCh] = {0};
1147// const TClonesArray* oldMisAlignArray(transformer->GetMisAlignmentData());
1148
1149// AliMUONGeometryTransformer* newGeometryTransformer = new AliMUONGeometryTransformer();
1150// for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); ++iMt) {
1151
1152// // module transformers
1153// const AliMUONGeometryModuleTransformer* kModuleTransformer = transformer->GetModuleTransformer(iMt, kTRUE);
1154
1155// AliMUONGeometryModuleTransformer* newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt);
1156// newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1157
1158// // get transformation
1159// TGeoHMatrix deltaModuleTransform(DeltaTransform(lModuleMisAlignment));
1160
1161// // update module
1162// TGeoHMatrix moduleTransform(*kModuleTransformer->GetTransformation());
1163// TGeoHMatrix newModuleTransform(AliMUONGeometryBuilder::Multiply(deltaModuleTransform, moduleTransform));
1164// newModuleTransformer->SetTransformation(newModuleTransform);
1165
1166// // Get matching old alignment and update current matrix accordingly
1167// if (oldMisAlignArray) {
1168
1169// const AliAlignObjMatrix* oldAlignObj(0);
1170// const Int_t moduleId(kModuleTransformer->GetModuleId());
1171// const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, moduleId);
1172// for (Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos) {
1173
1174// const AliAlignObjMatrix* localAlignObj(dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At(pos)));
1175// if (localAlignObj && localAlignObj->GetVolUID() == volId) {
1176// oldAlignObj = localAlignObj;
1177// break;
1178// }
1179// }
1180
1181// // multiply
1182// if (oldAlignObj) {
1183
1184// TGeoHMatrix oldMatrix;
1185// oldAlignObj->GetMatrix(oldMatrix);
1186// deltaModuleTransform.Multiply(&oldMatrix);
1187// }
1188// }
1189
1190// // Create module mis alignment matrix
1191// newGeometryTransformer->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
1192
1193// AliMpExMap* detElements = kModuleTransformer->GetDetElementStore();
1194
1195// TIter next(detElements->CreateIterator());
1196// AliMUONGeometryDetElement* detElement;
1197// Int_t iDe(-1);
1198// while ((detElement = static_cast<AliMUONGeometryDetElement*>(next()))) {
1199// ++iDe;
1200// // make a new detection element
1201// AliMUONGeometryDetElement* newDetElement = new AliMUONGeometryDetElement(detElement->GetId(), detElement->GetVolumePath());
1202// TString lDetElemName(detElement->GetDEName());
1203// lDetElemName.ReplaceAll("DE", "");
1204
1205// // store detector element id and number
1206// const Int_t iDetElemId = lDetElemName.Atoi();
1207// if (DetElemIsValid(iDetElemId)) {
1208
1209// const Int_t iDetElemNumber(GetDetElemNumber(iDetElemId));
1210
1211// for (int i = 0; i < fgNParCh; ++i) {
1212// lDetElemMisAlignment[i] = 0.0;
1213// if (iMt < fgNTrkMod) {
1214// lDetElemMisAlignment[i] = misAlignments[iDetElemNumber * fgNParCh + i];
1215// }
1216// }
1217
1218// // get transformation
1219// TGeoHMatrix deltaGlobalTransform(DeltaTransform(lDetElemMisAlignment));
1220
1221// // update module
1222// TGeoHMatrix globalTransform(*detElement->GetGlobalTransformation());
1223// TGeoHMatrix newGlobalTransform(AliMUONGeometryBuilder::Multiply(deltaGlobalTransform, globalTransform));
1224// newDetElement->SetGlobalTransformation(newGlobalTransform);
1225// newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
1226
1227// // Get matching old alignment and update current matrix accordingly
1228// if (oldMisAlignArray) {
1229
1230// const AliAlignObjMatrix* oldAlignObj(0);
1231// const int detElemId(detElement->GetId());
1232// const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId);
1233// for (Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos) {
1234
1235// const AliAlignObjMatrix* localAlignObj(dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At(pos)));
1236// if (localAlignObj && localAlignObj->GetVolUID() == volId) {
1237// oldAlignObj = localAlignObj;
1238// break;
1239// }
1240// }
1241
1242// // multiply
1243// if (oldAlignObj) {
1244
1245// TGeoHMatrix oldMatrix;
1246// oldAlignObj->GetMatrix(oldMatrix);
1247// deltaGlobalTransform.Multiply(&oldMatrix);
1248// }
1249// }
1250
1251// // Create misalignment matrix
1252// newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
1253
1254// } else {
1255
1256// // "invalid" detector elements come from MTR and are left unchanged
1257// Aliinfo(Form("Keeping detElement %i unchanged", iDetElemId));
1258
1259// // update module
1260// TGeoHMatrix globalTransform(*detElement->GetGlobalTransformation());
1261// newDetElement->SetGlobalTransformation(globalTransform);
1262// newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
1263
1264// // Get matching old alignment and update current matrix accordingly
1265// if (oldMisAlignArray) {
1266
1267// const AliAlignObjMatrix* oldAlignObj(0);
1268// const int detElemId(detElement->GetId());
1269// const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId);
1270// for (Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos) {
1271
1272// const AliAlignObjMatrix* localAlignObj(dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At(pos)));
1273// if (localAlignObj && localAlignObj->GetVolUID() == volId) {
1274// oldAlignObj = localAlignObj;
1275// break;
1276// }
1277// }
1278
1279// // multiply
1280// if (oldAlignObj) {
1281
1282// TGeoHMatrix oldMatrix;
1283// oldAlignObj->GetMatrix(oldMatrix);
1284// newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), oldMatrix);
1285// }
1286// }
1287// }
1288// }
1289
1290// newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1291// }
1292
1293// return newGeometryTransformer;
1294// }
1295
1296//______________________________________________________________________
1297void Alignment::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY)
1298{
1299
1302 TMatrixDSym mChCorrMatrix(6);
1303 mChCorrMatrix[0][0] = chResX * chResX;
1304 mChCorrMatrix[1][1] = chResY * chResY;
1305
1306 TMatrixDSym mDECorrMatrix(6);
1307 mDECorrMatrix[0][0] = deResX * deResX;
1308 mDECorrMatrix[1][1] = deResY * deResY;
1309
1310 o2::detectors::AlignParam* alignMat = 0x0;
1311
1312 for (Int_t chId = 0; chId <= 9; ++chId) {
1313
1314 // skip chamber if selection is valid, and does not match
1315 if (rChId > 0 && chId + 1 != rChId)
1316 continue;
1317
1318 TString chName1;
1319 TString chName2;
1320 if (chId < 4) {
1321
1322 chName1 = Form("GM%d", chId);
1323 chName2 = Form("GM%d", chId);
1324
1325 } else {
1326
1327 chName1 = Form("GM%d", 4 + (chId - 4) * 2);
1328 chName2 = Form("GM%d", 4 + (chId - 4) * 2 + 1);
1329 }
1330
1331 for (int i = 0; i < misAlignArray->GetEntries(); ++i) {
1332
1333 alignMat = (o2::detectors::AlignParam*)misAlignArray->At(i);
1334 TString volName(alignMat->getSymName());
1335 if ((volName.Contains(chName1) &&
1336 ((volName.Last('/') == volName.Index(chName1) + chName1.Length()) ||
1337 (volName.Length() == volName.Index(chName1) + chName1.Length()))) ||
1338 (volName.Contains(chName2) &&
1339 ((volName.Last('/') == volName.Index(chName2) + chName2.Length()) ||
1340 (volName.Length() == volName.Index(chName2) + chName2.Length())))) {
1341
1342 volName.Remove(0, volName.Last('/') + 1);
1343 // if (volName.Contains("GM")){
1344 // alignMat->SetCorrMatrix(mChCorrMatrix);
1345 // }else if (volName.Contains("DE")){
1346 // alignMat->SetCorrMatrix(mDECorrMatrix);
1347 // }
1348 }
1349 }
1350 }
1351}
1352
1353//_____________________________________________________
1354LocalTrackParam Alignment::RefitStraightTrack(Track& track, Double_t z0) const
1355{
1356
1357 // initialize matrices
1358 TMatrixD AtGASum(4, 4);
1359 AtGASum.Zero();
1360
1361 TMatrixD AtGMSum(4, 1);
1362 AtGMSum.Zero();
1363
1364 // loop over clusters
1365 for (auto itTrackParam(track.begin()); itTrackParam != track.end(); ++itTrackParam) {
1366
1367 // get track parameters
1368 if (!&*itTrackParam)
1369 continue;
1370
1371 // get cluster
1372 const Cluster* cluster = itTrackParam->getClusterPtr();
1373 if (!cluster)
1374 continue;
1375
1376 // projection matrix
1377 TMatrixD A(2, 4);
1378 A.Zero();
1379 A(0, 0) = 1;
1380 A(0, 2) = (cluster->getZ() - z0);
1381 A(1, 1) = 1;
1382 A(1, 3) = (cluster->getZ() - z0);
1383
1384 TMatrixD At(TMatrixD::kTransposed, A);
1385
1386 // gain matrix
1387 TMatrixD G(2, 2);
1388 G.Zero();
1389 G(0, 0) = 1.0 / Square(cluster->getEx());
1390 G(1, 1) = 1.0 / Square(cluster->getEy());
1391
1392 const TMatrixD AtG(At, TMatrixD::kMult, G);
1393 const TMatrixD AtGA(AtG, TMatrixD::kMult, A);
1394 AtGASum += AtGA;
1395
1396 // measurement
1397 TMatrixD M(2, 1);
1398 M(0, 0) = cluster->getX();
1399 M(1, 0) = cluster->getY();
1400 const TMatrixD AtGM(AtG, TMatrixD::kMult, M);
1401 AtGMSum += AtGM;
1402 }
1403
1404 // perform inversion
1405 TMatrixD AtGASumInv(TMatrixD::kInverted, AtGASum);
1406 TMatrixD X(AtGASumInv, TMatrixD::kMult, AtGMSum);
1407
1408 // // TODO: compare with initial track parameters
1409 // Aliinfo( Form( "x: %.3f vs %.3f", fTrackPos0[0], X(0,0) ) );
1410 // Aliinfo( Form( "y: %.3f vs %.3f", fTrackPos0[1], X(1,0) ) );
1411 // Aliinfo( Form( "dxdz: %.6g vs %.6g", fTrackSlope0[0], X(2,0) ) );
1412 // Aliinfo( Form( "dydz: %.6g vs %.6g\n", fTrackSlope0[1], X(3,0) ) );
1413
1414 // fill output parameters
1415 LocalTrackParam out;
1416 out.fTrackX = X(0, 0);
1417 out.fTrackY = X(1, 0);
1418 out.fTrackZ = z0;
1419 out.fTrackSlopeX = X(2, 0);
1420 out.fTrackSlopeY = X(3, 0);
1421
1422 return out;
1423}
1424
1425//_____________________________________________________
1426void Alignment::FillDetElemData(const Cluster* cluster)
1427{
1428 // LOG(fatal) << __PRETTY_FUNCTION__ << " is disabled";
1429 LOG(info) << __PRETTY_FUNCTION__ << " is enabled";
1430
1432 // get detector element number from Alice ID
1433 const Int_t detElemId = cluster->getDEId();
1434 fDetElemNumber = GetDetElemNumber(detElemId);
1435
1436 // get detector element
1437 // const AliMUONGeometryDetElement detElement(detElemId);
1438 auto fTransform = fTransformCreator(detElemId);
1439 /*
1440 get the global transformation matrix and store its inverse, in order to manually perform
1441 the global to Local transformations needed to calculate the derivatives
1442 */
1443 // fTransform = fTransform.Inverse();
1444 // fTransform.GetTransformMatrix(fGeoCombiTransInverse);
1445}
1446
1447//______________________________________________________________________
1448void Alignment::FillRecPointData(const Cluster* cluster)
1449{
1450
1452 fClustPos[0] = cluster->getX();
1453 fClustPos[1] = cluster->getY();
1454 fClustPos[2] = cluster->getZ();
1455}
1456
1457//______________________________________________________________________
1458void Alignment::FillTrackParamData(const TrackParam* trackParam)
1459{
1460
1462 fTrackPos[0] = trackParam->getNonBendingCoor();
1463 fTrackPos[1] = trackParam->getBendingCoor();
1464 fTrackPos[2] = trackParam->getZ();
1465 fTrackSlope[0] = trackParam->getNonBendingSlope();
1466 fTrackSlope[1] = trackParam->getBendingSlope();
1467}
1468
1469//______________________________________________________________________
1470void Alignment::LocalEquationX(void)
1471{
1473
1474 // 'inverse' (GlobalToLocal) rotation matrix
1475 const Double_t* r(fGeoCombiTransInverse.GetRotationMatrix());
1476
1477 // local derivatives
1478 SetLocalDerivative(0, r[0]);
1479 SetLocalDerivative(1, r[0] * (fTrackPos[2] - fTrackPos0[2]));
1480 SetLocalDerivative(2, r[1]);
1481 SetLocalDerivative(3, r[1] * (fTrackPos[2] - fTrackPos0[2]));
1482
1483 // global derivatives
1484 /*
1485 alignment parameters are
1486 0: delta_x
1487 1: delta_y
1488 2: delta_phiz
1489 3: delta_z
1490 */
1491
1492 SetGlobalDerivative(fDetElemNumber * fgNParCh + 0, -r[0]);
1493 SetGlobalDerivative(fDetElemNumber * fgNParCh + 1, -r[1]);
1494
1495 if (fBFieldOn) {
1496
1497 // use local position for derivatives vs 'delta_phi_z'
1498 SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[1] * fTrackPos[0] + r[0] * fTrackPos[1]);
1499
1500 // use local slopes for derivatives vs 'delta_z'
1501 SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[0] * fTrackSlope[0] + r[1] * fTrackSlope[1]);
1502
1503 } else {
1504
1505 // local copy of extrapolated track positions
1506 const Double_t trackPosX = fTrackPos0[0] + fTrackSlope0[0] * (fTrackPos[2] - fTrackPos0[2]);
1507 const Double_t trackPosY = fTrackPos0[1] + fTrackSlope0[1] * (fTrackPos[2] - fTrackPos0[2]);
1508
1509 // use properly extrapolated position for derivatives vs 'delta_phi_z'
1510 SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[1] * trackPosX + r[0] * trackPosY);
1511
1512 // use slopes at origin for derivatives vs 'delta_z'
1513 SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[0] * fTrackSlope0[0] + r[1] * fTrackSlope0[1]);
1514 }
1515
1516 // store local equation
1517 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]);
1518}
1519
1520//______________________________________________________________________
1521void Alignment::LocalEquationY(void)
1522{
1524
1525 // 'inverse' (GlobalToLocal) rotation matrix
1526 const Double_t* r(fGeoCombiTransInverse.GetRotationMatrix());
1527
1528 // store local derivatives
1529 SetLocalDerivative(0, r[3]);
1530 SetLocalDerivative(1, r[3] * (fTrackPos[2] - fTrackPos0[2]));
1531 SetLocalDerivative(2, r[4]);
1532 SetLocalDerivative(3, r[4] * (fTrackPos[2] - fTrackPos0[2]));
1533
1534 // set global derivatives
1535 SetGlobalDerivative(fDetElemNumber * fgNParCh + 0, -r[3]);
1536 SetGlobalDerivative(fDetElemNumber * fgNParCh + 1, -r[4]);
1537
1538 if (fBFieldOn) {
1539
1540 // use local position for derivatives vs 'delta_phi'
1541 SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[4] * fTrackPos[0] + r[3] * fTrackPos[1]);
1542
1543 // use local slopes for derivatives vs 'delta_z'
1544 SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[3] * fTrackSlope[0] + r[4] * fTrackSlope[1]);
1545
1546 } else {
1547
1548 // local copy of extrapolated track positions
1549 const Double_t trackPosX = fTrackPos0[0] + fTrackSlope0[0] * (fTrackPos[2] - fTrackPos0[2]);
1550 const Double_t trackPosY = fTrackPos0[1] + fTrackSlope0[1] * (fTrackPos[2] - fTrackPos0[2]);
1551
1552 // use properly extrapolated position for derivatives vs 'delta_phi'
1553 SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[4] * trackPosX + r[3] * trackPosY);
1554
1555 // use slopes at origin for derivatives vs 'delta_z'
1556 SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[3] * fTrackSlope0[0] + r[4] * fTrackSlope0[1]);
1557 }
1558
1559 // store local equation
1560 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]);
1561}
1562
1563//_________________________________________________________________________
1564TGeoCombiTrans Alignment::DeltaTransform(const double* lMisAlignment) const
1565{
1567
1568 // translation
1569 const TGeoTranslation deltaTrans(lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]);
1570
1571 // rotation
1572 TGeoRotation deltaRot;
1573 deltaRot.RotateZ(lMisAlignment[2] * 180. / TMath::Pi());
1574
1575 // combined rotation and translation.
1576 return TGeoCombiTrans(deltaTrans, deltaRot);
1577}
1578
1579//______________________________________________________________________
1580void Alignment::AddConstraint(Double_t* par, Double_t value)
1581{
1583 if (!fInitialized) {
1584 LOG(fatal) << "Millepede is not initialized";
1585 }
1586
1587 fMillepede->SetGlobalConstraint(par, value);
1588}
1589
1590//______________________________________________________________________
1591Bool_t Alignment::DetElemIsValid(Int_t iDetElemId) const
1592{
1594 const Int_t iCh = iDetElemId / 100;
1595 const Int_t iDet = iDetElemId % 100;
1596 return (iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh - 1]);
1597}
1598
1599//______________________________________________________________________
1600Int_t Alignment::GetDetElemNumber(Int_t iDetElemId) const
1601{
1603 // get chamber and element number in chamber
1604 const Int_t iCh = iDetElemId / 100;
1605 const Int_t iDet = iDetElemId % 100;
1606
1607 // make sure detector index is valid
1608 if (!(iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh - 1])) {
1609 LOG(fatal) << "Invalid detector element id: " << iDetElemId;
1610 }
1611
1612 // add number of detectors up to this chamber
1613 return iDet + fgSNDetElemCh[iCh - 1];
1614}
1615
1616//______________________________________________________________________
1617Int_t Alignment::GetChamberId(Int_t iDetElemNumber) const
1618{
1620 Int_t iCh(0);
1621 for (iCh = 0; iCh < fgNCh; iCh++) {
1622 if (iDetElemNumber < fgSNDetElemCh[iCh])
1623 break;
1624 }
1625
1626 return iCh;
1627}
1628
1629//______________________________________________________________________
1630TString Alignment::GetParameterMaskString(UInt_t mask) const
1631{
1632 TString out;
1633 if (mask & ParX)
1634 out += "X";
1635 if (mask & ParY)
1636 out += "Y";
1637 if (mask & ParZ)
1638 out += "Z";
1639 if (mask & ParTZ)
1640 out += "T";
1641 return out;
1642}
1643
1644//______________________________________________________________________
1645TString Alignment::GetSidesMaskString(UInt_t mask) const
1646{
1647 TString out;
1648 if (mask & SideTop)
1649 out += "T";
1650 if (mask & SideLeft)
1651 out += "L";
1652 if (mask & SideBottom)
1653 out += "B";
1654 if (mask & SideRight)
1655 out += "R";
1656 return out;
1657}
1658
1659} // namespace mch
1660} // namespace o2
AlignParam.h
Definition of the base alignment parameters class.
fgSNDetElemCh
const int fgSNDetElemCh[fgNCh+1]
Definition AlignRecordSpec.cxx:61
fgNCh
const int fgNCh
Definition AlignRecordSpec.cxx:59
fgNDetElemCh
const int fgNDetElemCh[fgNCh]
Definition AlignRecordSpec.cxx:60
Track.h
Definition of the MCH track for internal use.
i
int32_t i
Definition GPUCommonAlgorithm.h:443
Geometry.h
Interface for MCH geometry creation.
TrackParam.h
Definition of the MCH track parameters for internal use.
A
Definition A.h:16
o2::detectors::AlignParam::getSymName
const std::string & getSymName() const
return symbolic name of the volume
Definition AlignParam.h:45
o2::hmpid::Cluster
HMPID cluster implementation.
Definition Cluster.h:27
o2::mch::Array
self initialized array, used for adding constraints
Definition MchAlignment.cxx:115
o2::mch::Array::Array
Array(void)
contructor
Definition MchAlignment.cxx:119
x
GLint GLenum GLint x
Definition glcorearb.h:403
first
GLint first
Definition glcorearb.h:399
weight
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
value
GLsizei const GLfloat * value
Definition glcorearb.h:819
values
GLenum GLsizei GLsizei GLint * values
Definition glcorearb.h:1576
r
GLboolean r
Definition glcorearb.h:1233
mask
GLint GLuint mask
Definition glcorearb.h:291
o2::mch::dcs::detElemId
int detElemId(Chamber chamber, Side side, int number)
Definition DetectionElement.cxx:63
o2::mch::Square
Double_t Square(Double_t x)
Definition MchAlignment.cxx:138
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition BitstreamReader.h:24
std
Defining DataPointCompositeObject explicitly as copiable.
Definition CcdbObjectInfo.h:121
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"