13#ifndef BMNGLOBALALIGNMENT_H
14#define BMNGLOBALALIGNMENT_H 1
23#include <TClonesArray.h>
24#include <TFitResult.h>
25#include <TFitResultPtr.h>
26#include <TGeoManager.h>
30#include "FairRootManager.h"
31#include "FairEventHeader.h"
59 BmnGlobalAlignment(Int_t, TString, Int_t, Int_t, TString misAlignFile =
"misAlignment.root", Bool_t doTest = kFALSE);
64 virtual void Exec(Option_t* opt);
73 fDetectorSet[0] = gem;
78 fUseRealHitErrors = flag;
90 fUseRegularization = flag;
94 fOutlierdownweighting = n;
98 fDwfractioncut = fraction;
102 fChi2MaxPerNDF = val;
106 fMinHitsAccepted = val;
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) {
140 cout <<
"GEM" << endl;
141 const Int_t nModulMax = 2;
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];
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;
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;
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;
171 for (Int_t iStat = 0; iStat < fDetectorSI->
GetNStations(); iStat++)
172 si[iStat] =
new Bool_t[nModulMax];
182 for (Int_t iStat = 0; iStat < fDetectorSI->
GetNStations(); iStat++)
183 for (Int_t iMod = 0; iMod < fDetectorSI->
GetSiliconStation(iStat)->GetNModules(); iMod++)
185 fixedSiElements[iStat][iMod] = kTRUE;
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;
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;
199 for (Int_t iStat = 0; iStat < fDetectorSI->
GetNStations(); iStat++)
200 si[iStat] =
new Bool_t[nModulMax];
222 for (Int_t iStat = 0; iStat < fDetectorSI->
GetNStations(); iStat++)
223 for (Int_t iMod = 0; iMod < fDetectorSI->
GetSiliconStation(iStat)->GetNModules(); iMod++)
225 fixedSiElements[iStat][iMod] = kTRUE;
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);
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) {
242 cout <<
"GEM" << endl;
243 const Int_t nModulMax = 2;
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];
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;
291 fIsExcludedTx = kTRUE;
297 fIsExcludedTy = kTRUE;
303 fUseConstraints = flag;
307 void CreateDetectorGeometries();
308 const Int_t MakeBinFile();
309 void MakeSteerFile();
310 void MilleNoFieldRuns(
BmnGlobalTrack*, Int_t, Int_t, vector <BmnMilleContainer*>&, vector <BmnMilleContainer*>&);
312 void ReadPedeOutput(ifstream&);
313 void ExtractCorrValues(ifstream&, Double_t*);
316 inline Int_t GemStatModLabel(Int_t st, Int_t mod) {
318 return gemStatModLabel.find(pair <Int_t, Int_t> (st, mod))->second;
321 inline Int_t SiliconStatModLabel(Int_t st, Int_t mod) {
324 return silStatModLabel.find(pair <Int_t, Int_t> (st, mod))->second;
327 inline vector <Int_t> GetSiliconStatMod(Int_t label) {
328 while (label % 3 != 0)
332 for (
auto it : silStatModLabel) {
333 if (it.second != label)
335 out.push_back(it.first.first);
336 out.push_back(it.first.second);
342 inline vector <Int_t> GetGemStatMod(Int_t label) {
343 while (label % 3 != 0)
347 for (
auto it : gemStatModLabel) {
348 if (it.second != label)
350 out.push_back(it.first.first);
351 out.push_back(it.first.second);
359 Bool_t isSRC = (fRunId < 3589) ? kTRUE : kFALSE;
361 Int_t srcGEM[10][2] = {
374 Int_t bmnGEM[6][2] = {
383 const Int_t nGemStats = isSRC ? 10 : 6;
384 const Int_t nGemMods = 2;
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];
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];
397 const Int_t nSilStats = 3;
398 const Int_t nSilMods = 8;
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}
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}
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];
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];
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;
431 cout <<
"SILICON" << endl;
432 for (
auto it : silStatModLabel)
433 cout << it.first.first <<
", " << it.first.second <<
" --> " << it.second << endl;
436 void _Mille(Double_t*, Double_t*,
BmnMille*, Int_t);
438 static Int_t fCurrentEvent;
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;
452 TClonesArray* fGemAlignCorr;
453 TClonesArray* fSiAlignCorr;
455 TString fBranchGemAlignCorr;
456 TString fBranchSiAlignCorr;
458 Bool_t* fDetectorSet;
460 TString fBranchSiHits;
461 TString fBranchGemHits;
463 TString fBranchGemTracks;
464 TString fBranchSilTracks;
466 TString fBranchGlobalTracks;
467 TString fBranchGemResiduals;
468 TString fBranchFairEventHeader;
470 TClonesArray* fSiHits;
471 TClonesArray* fGemHits;
473 TClonesArray* fGemTracks;
474 TClonesArray* fSilTracks;
475 TClonesArray* fGlobalTracks;
477 TString fRecoFileName;
480 FairEventHeader* fFairEventHeader;
482 Bool_t fUseRealHitErrors;
485 Double_t fChi2MaxPerNDF;
486 Double_t fMinHitsAccepted;
492 Bool_t fIsExcludedTx;
493 Bool_t fIsExcludedTy;
501 UInt_t fNumOfIterations;
506 Bool_t fUseRegularization;
508 Double_t fChisqcut[2];
510 Int_t fOutlierdownweighting;
511 Double_t fDwfractioncut;
517 vector <Int_t> fFixedStats;
518 Bool_t** fixedGemElements;
519 Bool_t** fixedSiElements;
525 Bool_t fUseConstraints;
528 TString fMisAlignFile;
529 TClonesArray* fBmnGemMisalign;
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
BmnGemStripStation * GetGemStation(Int_t station_num)
virtual InitStatus Init()
void SetDwfractionCut(Double_t fraction)
void SetUseConstraints(Bool_t flag)
void SetExclusionRangeTy(Double_t min, Double_t max)
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)
virtual void Exec(Option_t *opt)
void SetDebug(Bool_t flag)
void SetChi2MaxPerNDF(Double_t val)
virtual ~BmnGlobalAlignment()
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)
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)
BmnGlobalAlignment(Int_t, TString, Int_t, Int_t, TString misAlignFile="misAlignment.root", Bool_t doTest=kFALSE)
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)