BmnRoot
Loading...
Searching...
No Matches
BmnGlobalAlignment.h
Go to the documentation of this file.
1// @(#)bmnroot/alignment:$Id$
2// Author: Pavel Batyuk <pavel.batyuk@jinr.ru> 2017-03-31
3
5// //
6// BmnGlobalAlignment //
7// //
8// Alignment of tracking detectors. //
9// //
10// Uses Volker Blobel and Claus Kleinwort Millepede II //
11// //
13#ifndef BMNGLOBALALIGNMENT_H
14#define BMNGLOBALALIGNMENT_H 1
15
16#include <vector>
17#include <fstream>
18#include <sstream>
19
20#include <TSystem.h>
21#include <TMath.h>
22#include <TNamed.h>
23#include <TClonesArray.h>
24#include <TFitResult.h>
25#include <TFitResultPtr.h>
26#include <TGeoManager.h>
27#include <TString.h>
28
29#include "FairTask.h"
30#include "FairRootManager.h"
31#include "FairEventHeader.h"
32
33#include <BmnGlobalTrack.h>
34#include "BmnGemStripHit.h"
35#include "BmnSiliconHit.h"
36#include "BmnMille.h"
37#include "BmnMilleContainer.h"
40
41#include <BmnSiliconStationSet.h>
42#include <BmnGemStripStationSet.h>
43
44#include <UniDetectorParameter.h>
45#include <UniRun.h>
46
47#include <BmnFieldMap.h>
48#include <BmnNewFieldMap.h>
49#include <FairField.h>
50
51using namespace std;
52using namespace TMath;
53
54class BmnGlobalAlignment : public FairTask {
55public:
56
59 BmnGlobalAlignment(Int_t, TString, Int_t, Int_t, TString misAlignFile = "misAlignment.root", Bool_t doTest = kFALSE);
60 virtual ~BmnGlobalAlignment();
61
62 virtual InitStatus Init();
63
64 virtual void Exec(Option_t* opt);
65
66 virtual void Finish();
67
68 void SetDoTest(Bool_t flag) {
69 fIsDoTest = flag;
70 }
71
72 void SetDetectors(Bool_t gem, Bool_t si) {
73 fDetectorSet[0] = gem;
74 fDetectorSet[1] = si;
75 }
76
77 void SetUseRealHitErrors(Bool_t flag) {
78 fUseRealHitErrors = flag;
79 }
80
81 void SetPreSigma(Double_t presigma) {
82 fPreSigma = presigma;
83 }
84
85 void SetAccuracy(Double_t accuracy) {
86 fAccuracy = accuracy;
87 }
88
89 void SetUseRegularization(Bool_t flag) {
90 fUseRegularization = flag;
91 }
92
94 fOutlierdownweighting = n;
95 }
96
97 void SetDwfractionCut(Double_t fraction) {
98 fDwfractioncut = fraction;
99 }
100
101 void SetChi2MaxPerNDF(Double_t val) {
102 fChi2MaxPerNDF = val;
103 }
104
105 void SetMinHitsAccepted(Int_t val) {
106 fMinHitsAccepted = val;
107 }
108
109 void SetTxMinMax(Double_t min, Double_t max) {
110 fTxMin = min;
111 fTxMax = max;
112 }
113
114 void SetTyMinMax(Double_t min, Double_t max) {
115 fTyMin = min;
116 fTyMax = max;
117 }
118
119 void SetDebug(Bool_t flag) {
120 fDebug = flag;
121 }
122
123 void SetHugecut(Double_t val) {
124 fHugecut = val;
125 }
126
127 void SetChisqcut(Double_t val1, Double_t val2) {
128 fChisqcut[0] = val1;
129 fChisqcut[1] = val2;
130 }
131
132 void SetEntriesPerParam(Int_t entries) {
133 fEntries = entries;
134 }
135
136 void SetGemFixedRun6(Bool_t st0_0, Bool_t st1_0, Bool_t st2_0,
137 Bool_t st3_0, Bool_t st3_1, Bool_t st4_0,
138 Bool_t st4_1, Bool_t st5_0, Bool_t st5_1) {
139
140 cout << "GEM" << endl;
141 const Int_t nModulMax = 2; // To be fixed
142 Bool_t** gems = new Bool_t*[fDetectorGEM->GetNStations()];
143 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
144 gems[iStat] = new Bool_t[nModulMax];
145 gems[0][0] = st0_0;
146 gems[1][0] = st1_0;
147 gems[2][0] = st2_0;
148 gems[3][0] = st3_0;
149 gems[3][1] = st3_1;
150 gems[4][0] = st4_0;
151 gems[4][1] = st4_1;
152 gems[5][0] = st5_0;
153 gems[5][1] = st5_1;
154
155 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
156 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
157 if (gems[iStat][iMod])
158 fixedGemElements[iStat][iMod] = kTRUE;
159
160 if (fDebug)
161 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
162 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
163 cout << "Stat = " << iStat << " Mod = " << iMod << " isFixed (true / false) " << fixedGemElements[iStat][iMod] << endl;
164 }
165
166 void SetSiFixedRun6(Bool_t st0_0, Bool_t st0_1, Bool_t st0_2, Bool_t st0_3,
167 Bool_t st0_4, Bool_t st0_5, Bool_t st0_6, Bool_t st0_7) {
168 cout << "SI" << endl;
169 const Int_t nModulMax = 8; // FIXME
170 Bool_t** si = new Bool_t*[fDetectorSI->GetNStations()];
171 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
172 si[iStat] = new Bool_t[nModulMax];
173 si[0][0] = st0_0;
174 si[0][1] = st0_1;
175 si[0][2] = st0_2;
176 si[0][3] = st0_3;
177 si[0][4] = st0_4;
178 si[0][5] = st0_5;
179 si[0][6] = st0_6;
180 si[0][7] = st0_7;
181
182 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
183 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
184 if (si[iStat][iMod])
185 fixedSiElements[iStat][iMod] = kTRUE;
186
187 if (fDebug)
188 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
189 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
190 cout << "Stat = " << iStat << " Mod = " << iMod << " isFixed (true / false) " << fixedSiElements[iStat][iMod] << endl;
191 }
192
193 void SetSiFixedRun7(Bool_t st0_0, Bool_t st0_1, Bool_t st0_2, Bool_t st0_3,
194 Bool_t st1_0, Bool_t st1_1,
195 Bool_t st2_0, Bool_t st2_1, Bool_t st2_2, Bool_t st2_3, Bool_t st2_4, Bool_t st2_5, Bool_t st2_6, Bool_t st2_7) {
196 cout << "SI" << endl;
197 const Int_t nModulMax = 8; // FIXME
198 Bool_t** si = new Bool_t*[fDetectorSI->GetNStations()];
199 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
200 si[iStat] = new Bool_t[nModulMax];
201
202 // Stat #0
203 si[0][0] = st0_0;
204 si[0][1] = st0_1;
205 si[0][2] = st0_2;
206 si[0][3] = st0_3;
207
208 // Stat #1
209 si[1][0] = st1_0;
210 si[1][1] = st1_1;
211
212 // Stat #2
213 si[2][0] = st2_0;
214 si[2][1] = st2_1;
215 si[2][2] = st2_2;
216 si[2][3] = st2_3;
217 si[2][4] = st2_4;
218 si[2][5] = st2_5;
219 si[2][6] = st2_6;
220 si[2][7] = st2_7;
221
222 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
223 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
224 if (si[iStat][iMod])
225 fixedSiElements[iStat][iMod] = kTRUE;
226 }
227
228 void SetSiConfigSRC(Bool_t st0_0, Bool_t st0_1, Bool_t st0_2, Bool_t st0_3,
229 Bool_t st1_0, Bool_t st1_1,
230 Bool_t st2_0, Bool_t st2_1, Bool_t st2_2, Bool_t st2_3, Bool_t st2_4, Bool_t st2_5, Bool_t st2_6, Bool_t st2_7) {
231 SetSiFixedRun7(st0_0, st0_1, st0_2, st0_3, st1_0, st1_1, st2_0, st2_1, st2_2, st2_3, st2_4, st2_5, st2_6, st2_7);
232 }
233
234 void SetGemConfigSRC(Bool_t st0_0, Bool_t st1_0, Bool_t st2_0, Bool_t st3_0,
235 Bool_t st4_0, Bool_t st4_1,
236 Bool_t st5_0, Bool_t st5_1,
237 Bool_t st6_0, Bool_t st6_1,
238 Bool_t st7_0, Bool_t st7_1,
239 Bool_t st8_0, Bool_t st8_1,
240 Bool_t st9_0, Bool_t st9_1) {
241
242 cout << "GEM" << endl;
243 const Int_t nModulMax = 2; // FIXME
244 Bool_t** gem = new Bool_t*[fDetectorGEM->GetNStations()];
245 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
246 gem[iStat] = new Bool_t[nModulMax];
247
248 // Stat #0
249 gem[0][0] = st0_0;
250
251 // Stat #1
252 gem[1][0] = st1_0;
253
254 // Stat #2
255 gem[2][0] = st2_0;
256
257 // Stat #3
258 gem[3][0] = st3_0;
259
260 // Stat #4
261 gem[4][0] = st4_0;
262 gem[4][1] = st4_1;
263
264 // Stat #5
265 gem[5][0] = st5_0;
266 gem[5][1] = st5_1;
267
268 // Stat #6
269 gem[6][0] = st6_0;
270 gem[6][1] = st6_1;
271
272 // Stat #7
273 gem[7][0] = st7_0;
274 gem[7][1] = st7_1;
275
276 // Stat #8
277 gem[8][0] = st8_0;
278 gem[8][1] = st8_1;
279
280 // Stat #9
281 gem[9][0] = st9_0;
282 gem[9][1] = st9_1;
283
284 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
285 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
286 if (gem[iStat][iMod])
287 fixedGemElements[iStat][iMod] = kTRUE;
288 }
289
290 void SetExclusionRangeTx(Double_t min, Double_t max) {
291 fIsExcludedTx = kTRUE;
292 fTxLeft = min;
293 fTxRight = max;
294 }
295
296 void SetExclusionRangeTy(Double_t min, Double_t max) {
297 fIsExcludedTy = kTRUE;
298 fTyLeft = min;
299 fTyRight = max;
300 }
301
302 void SetUseConstraints(Bool_t flag) {
303 fUseConstraints = flag;
304 }
305
306private:
307 void CreateDetectorGeometries();
308 const Int_t MakeBinFile();
309 void MakeSteerFile();
310 void MilleNoFieldRuns(BmnGlobalTrack*, Int_t, Int_t, vector <BmnMilleContainer*>&, vector <BmnMilleContainer*>&);
311 void Pede();
312 void ReadPedeOutput(ifstream&);
313 void ExtractCorrValues(ifstream&, Double_t*);
314 BmnMilleContainer* FillMilleContainer(BmnGlobalTrack*, BmnHit*);
315
316 inline Int_t GemStatModLabel(Int_t st, Int_t mod) {
317 // return 3 * (2 * st + mod + 1); // FIXME!!!
318 return gemStatModLabel.find(pair <Int_t, Int_t> (st, mod))->second;
319 }
320
321 inline Int_t SiliconStatModLabel(Int_t st, Int_t mod) {
322 // Int_t coeff[3] = {0, 6, 3}; FIXME!!
323 // return GemStatModLabel(st, mod) + coeff[st] * st;
324 return silStatModLabel.find(pair <Int_t, Int_t> (st, mod))->second;
325 }
326
327 inline vector <Int_t> GetSiliconStatMod(Int_t label) {
328 while (label % 3 != 0)
329 label++;
330
331 vector <Int_t> out;
332 for (auto it : silStatModLabel) {
333 if (it.second != label)
334 continue;
335 out.push_back(it.first.first);
336 out.push_back(it.first.second);
337 break;
338 }
339 return out;
340 }
341
342 inline vector <Int_t> GetGemStatMod(Int_t label) {
343 while (label % 3 != 0)
344 label++;
345
346 vector <Int_t> out;
347 for (auto it : gemStatModLabel) {
348 if (it.second != label)
349 continue;
350 out.push_back(it.first.first);
351 out.push_back(it.first.second);
352 break;
353 }
354 return out;
355 }
356
357 void FillMaps() {
358 // Choose appropriate configuration (BM@N or SRC)
359 Bool_t isSRC = (fRunId < 3589) ? kTRUE : kFALSE;
360
361 Int_t srcGEM[10][2] = {
362 {3, -1},
363 {6, -1},
364 {9, -1},
365 {12, -1},
366 {15, 18},
367 {21, 24},
368 {27, 30},
369 {33, 36},
370 {39, 42},
371 {45, 48}
372 };
373
374 Int_t bmnGEM[6][2] = {
375 {3, 6},
376 {9, 12},
377 {15, 18},
378 {21, 24},
379 {27, 30},
380 {33, 36}
381 };
382
383 const Int_t nGemStats = isSRC ? 10 : 6;
384 const Int_t nGemMods = 2;
385
386 Int_t** gem = new Int_t*[nGemStats];
387 for (Int_t iStat = 0; iStat < nGemStats; iStat++) {
388 gem[iStat] = new Int_t[nGemMods];
389 for (Int_t iMod = 0; iMod < nGemMods; iMod++)
390 gem[iStat][iMod] = isSRC ? srcGEM[iStat][iMod] : bmnGEM[iStat][iMod];
391 }
392
393 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
394 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
395 gemStatModLabel[pair <Int_t, Int_t> (iStat, iMod)] = gem[iStat][iMod];
396
397 const Int_t nSilStats = 3; // FIXME!!
398 const Int_t nSilMods = 8;
399
400 Int_t bmnSIL[nSilStats][nSilMods] = {
401 {39, 42, 45, 48, -1, -1, -1, -1},
402 {51, 54, -1, -1, -1, -1, -1, -1},
403 {57, 60, 63, 66, 69, 72, 75, 78}
404 };
405
406 Int_t srcSIL[nSilStats][nSilMods] = {
407 {51, 54, 57, 60, -1, -1, -1, -1},
408 {63, 66, -1, -1, -1, -1, -1, -1},
409 {69, 72, 75, 78, 81, 84, 87, 90}
410 };
411
412 Int_t** sil = new Int_t*[nSilStats];
413 for (Int_t iStat = 0; iStat < nSilStats; iStat++) {
414 sil[iStat] = new Int_t[nSilMods];
415 for (Int_t iMod = 0; iMod < nSilMods; iMod++)
416 sil[iStat][iMod] = isSRC ? srcSIL[iStat][iMod] : bmnSIL[iStat][iMod];
417 }
418
419 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
420 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
421 silStatModLabel[pair <Int_t, Int_t> (iStat, iMod)] = sil[iStat][iMod];
422
423 // if (!fDebug)
424 // return;
425
426 cout << "GEM" << endl;
427 cout << gemStatModLabel.size() << endl;
428 for (auto it : gemStatModLabel)
429 cout << it.first.first << ", " << it.first.second << " --> " << it.second << endl;
430
431 cout << "SILICON" << endl;
432 for (auto it : silStatModLabel)
433 cout << it.first.first << ", " << it.first.second << " --> " << it.second << endl;
434 }
435
436 void _Mille(Double_t*, Double_t*, BmnMille*, Int_t);
437
438 static Int_t fCurrentEvent;
439 Int_t fNEvents;
440 Int_t fNTracks;
441 map <Int_t, pair <vector <BmnMilleContainer*>, vector < BmnMilleContainer*>>> fCONTAINER;
442 map <Int_t, pair <vector <BmnMilleContainer*>, vector < BmnMilleContainer*>>>::iterator fITERATOR;
443 map <pair <Int_t, Int_t>, Int_t> gemStatModLabel;
444 map <pair <Int_t, Int_t>, Int_t> silStatModLabel;
445 Bool_t fIsField;
446 Int_t fRunPeriod;
447 Int_t fRunId;
448
449 BmnGemStripStationSet* fDetectorGEM; // GEM-geometry
450 BmnSiliconStationSet* fDetectorSI; // SI-geometry
451
452 TClonesArray* fGemAlignCorr;
453 TClonesArray* fSiAlignCorr;
454
455 TString fBranchGemAlignCorr;
456 TString fBranchSiAlignCorr;
457
458 Bool_t* fDetectorSet;
459
460 TString fBranchSiHits;
461 TString fBranchGemHits;
462
463 TString fBranchGemTracks;
464 TString fBranchSilTracks;
465
466 TString fBranchGlobalTracks;
467 TString fBranchGemResiduals;
468 TString fBranchFairEventHeader;
469
470 TClonesArray* fSiHits;
471 TClonesArray* fGemHits;
472
473 TClonesArray* fGemTracks;
474 TClonesArray* fSilTracks;
475 TClonesArray* fGlobalTracks;
476
477 TString fRecoFileName;
478 TChain* fChain;
479
480 FairEventHeader* fFairEventHeader;
481
482 Bool_t fUseRealHitErrors; // errors are taken from hit finder algorithm
483
484 // Restrictions on track params
485 Double_t fChi2MaxPerNDF;
486 Double_t fMinHitsAccepted;
487 Double_t fTxMin;
488 Double_t fTxMax;
489 Double_t fTyMin;
490 Double_t fTyMax;
491 // Range to be exluded in case of "beam-target" alignment...
492 Bool_t fIsExcludedTx; // do exclusion (true) or not, manipulated by SetExclusionRangeTx
493 Bool_t fIsExcludedTy; // do exclusion (true) or not, manipulated by SetExclusionRangeTy
494 Double_t fTxLeft;
495 Double_t fTxRight;
496 Double_t fTyLeft;
497 Double_t fTyRight;
498
499 // Int_t nSelectedTracks;
500
501 UInt_t fNumOfIterations;
502 Double_t fAccuracy;
503 Double_t fPreSigma;
504
505 // Millepede params
506 Bool_t fUseRegularization; // use regularization or not
507 Double_t fHugecut; // cut factor in iteration 0
508 Double_t fChisqcut[2]; // cut factor in iterations 1 and 2
509 Int_t fEntries; // lower limit on number of entries/parameter
510 Int_t fOutlierdownweighting; // number of internal iterations (> 1)
511 Double_t fDwfractioncut; // reject all cases with a down-weight fraction >= val
512
513 Int_t fNGL;
514 Int_t fNLC;
515 Int_t nDetectors;
516
517 vector <Int_t> fFixedStats;
518 Bool_t** fixedGemElements;
519 Bool_t** fixedSiElements;
520
521 Bool_t fDebug;
522 Int_t* Labels; //array containing a fixed param. number for each detector.
523
524 // Use constraints or not
525 Bool_t fUseConstraints;
526
527 Bool_t fIsDoTest;
528 TString fMisAlignFile;
529 TClonesArray* fBmnGemMisalign;
530
531 ClassDef(BmnGlobalAlignment, 1)
532};
533
534#endif
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:31
BmnGemStripStation * GetGemStation(Int_t station_num)
void SetDwfractionCut(Double_t fraction)
void SetUseConstraints(Bool_t flag)
void SetExclusionRangeTy(Double_t min, Double_t max)
virtual InitStatus Init()
void SetOutlierDownweighting(Int_t n)
void SetGemConfigSRC(Bool_t st0_0, Bool_t st1_0, Bool_t st2_0, Bool_t st3_0, Bool_t st4_0, Bool_t st4_1, Bool_t st5_0, Bool_t st5_1, Bool_t st6_0, Bool_t st6_1, Bool_t st7_0, Bool_t st7_1, Bool_t st8_0, Bool_t st8_1, Bool_t st9_0, Bool_t st9_1)
void SetDebug(Bool_t flag)
void SetChi2MaxPerNDF(Double_t val)
void SetSiFixedRun7(Bool_t st0_0, Bool_t st0_1, Bool_t st0_2, Bool_t st0_3, Bool_t st1_0, Bool_t st1_1, Bool_t st2_0, Bool_t st2_1, Bool_t st2_2, Bool_t st2_3, Bool_t st2_4, Bool_t st2_5, Bool_t st2_6, Bool_t st2_7)
virtual void Exec(Option_t *opt)
void SetGemFixedRun6(Bool_t st0_0, Bool_t st1_0, Bool_t st2_0, Bool_t st3_0, Bool_t st3_1, Bool_t st4_0, Bool_t st4_1, Bool_t st5_0, Bool_t st5_1)
void SetSiConfigSRC(Bool_t st0_0, Bool_t st0_1, Bool_t st0_2, Bool_t st0_3, Bool_t st1_0, Bool_t st1_1, Bool_t st2_0, Bool_t st2_1, Bool_t st2_2, Bool_t st2_3, Bool_t st2_4, Bool_t st2_5, Bool_t st2_6, Bool_t st2_7)
void SetExclusionRangeTx(Double_t min, Double_t max)
void SetEntriesPerParam(Int_t entries)
void SetUseRegularization(Bool_t flag)
void SetTxMinMax(Double_t min, Double_t max)
void SetTyMinMax(Double_t min, Double_t max)
void SetDoTest(Bool_t flag)
void SetChisqcut(Double_t val1, Double_t val2)
void SetDetectors(Bool_t gem, Bool_t si)
void SetSiFixedRun6(Bool_t st0_0, Bool_t st0_1, Bool_t st0_2, Bool_t st0_3, Bool_t st0_4, Bool_t st0_5, Bool_t st0_6, Bool_t st0_7)
void SetAccuracy(Double_t accuracy)
void SetUseRealHitErrors(Bool_t flag)
void SetPreSigma(Double_t presigma)
void SetHugecut(Double_t val)
void SetMinHitsAccepted(Int_t val)
BmnSiliconStation * GetSiliconStation(Int_t station_num)
STL namespace.