BmnRoot
Loading...
Searching...
No Matches
BmnGlobalAlignment.cxx
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#include <iconv.h>
14#include <TGraph.h>
15#include <TStyle.h>
16#include <TRandom2.h>
17#include "BmnGlobalAlignment.h"
18#include "CbmVertex.h"
19#include "BmnSiliconHit.h"
20
21Int_t BmnGlobalAlignment::fCurrentEvent = 0;
22
25
26BmnGlobalAlignment::~BmnGlobalAlignment() {
27 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
28 delete [] fixedGemElements[iStat];
29 delete [] fixedGemElements;
30
31 delete [] fDetectorSet;
32 delete fDetectorGEM;
33 delete fDetectorSI;
34 if (Labels) {
35 delete [] Labels;
36 system("cp millepede.res Millepede.res");
37 // Remove useless text files finalizing the code execution
38 system("rm millepede.*");
39 }
40}
41
42BmnGlobalAlignment::BmnGlobalAlignment(Int_t nEvents, TString inFileName, Int_t period, Int_t run, TString misAlignFile, Bool_t doTest)
43: fIsField(kFALSE),
44 fRunPeriod(period),
45 fRunId(run),
46 fGemAlignCorr(nullptr),
47 fGemHits(nullptr),
48 fGemTracks(nullptr),
49 fSilTracks(nullptr),
50 fGlobalTracks(nullptr),
51 fUseRealHitErrors(kFALSE),
52 fChi2MaxPerNDF(LDBL_MAX),
53 fMinHitsAccepted(3),
54 fTxMin(-LDBL_MAX),
55 fTxMax(LDBL_MAX),
56 fTyMin(-LDBL_MAX),
57 fTyMax(LDBL_MAX),
58 fIsExcludedTx(kFALSE),
59 fIsExcludedTy(kFALSE),
60 fTxLeft(0.),
61 fTxRight(0.),
62 fTyLeft(0.),
63 fTyRight(0.),
64 fNumOfIterations(500000),
65 fAccuracy(1e-3),
66 fPreSigma(1.),
67 fUseRegularization(kFALSE),
68 fHugecut(50.),
69 fEntries(10),
70 fOutlierdownweighting(0),
71 fDwfractioncut(0.0),
72 fNGL(0),
73 fNLC(4),
74 nDetectors(2),
75 fDebug(true),
76 Labels(nullptr),
77 fUseConstraints(kTRUE),
78 fIsDoTest(doTest),
79 fMisAlignFile(misAlignFile),
80 fBmnGemMisalign(nullptr)
81{
82 fNEvents = nEvents;
83 fNTracks = 0;
84 fRecoFileName = inFileName;
85
86 fBranchSiHits = "BmnSiliconHit";
87 fBranchGemHits = "BmnGemStripHit";
88
89 fBranchGemTracks = "BmnGemTrack";
90 fBranchSilTracks = "BmnSiliconTrack";
91 fBranchGlobalTracks = "BmnGlobalTrack";
92
93 fBranchGemAlignCorr = "BmnGemAlignCorrections";
94 fBranchSiAlignCorr = "BmnSiliconAlignCorrections";
95
96 fBranchGemResiduals = "BmnResiduals";
97 fBranchFairEventHeader = "EventHeader.";
98
99 CreateDetectorGeometries();
100 FillMaps();
101}
102
103InitStatus BmnGlobalAlignment::Init() {
104 cout << " BmnGlobalAlignment::Init() " << endl;
105 cout << "Use detectors: GEM - " << fDetectorSet[0] << " SILICON - " << fDetectorSet[1] << endl;
106
107 TChain* chain = new TChain("bmndata");
108 chain->Add(fRecoFileName.Data());
109 FairEventHeader* evHeader = NULL;
110 chain->SetBranchAddress(fBranchFairEventHeader.Data(), &evHeader);
111 chain->GetEntry(0);
112 delete chain;
113
114 Double_t fieldVolt = 0.;
115 UniRun* runInfo = NULL;
116 if (fRunId != 0) {
117 runInfo = UniRun::GetRun(fRunPeriod, fRunId);
118 if (!runInfo)
119 throw;
120 fieldVolt = *runInfo->GetFieldVoltage();
121 fIsField = (fieldVolt > 10.) ? kTRUE : kFALSE;
122 }
123
124 FairRootManager* ioman = FairRootManager::Instance();
125
126 fSiHits = (TClonesArray*) ioman->GetObject(fBranchSiHits.Data());
127 fGemHits = (TClonesArray*) ioman->GetObject(fBranchGemHits.Data());
128
129 fSilTracks = (TClonesArray*) ioman->GetObject(fBranchSilTracks.Data());
130 fGemTracks = (TClonesArray*) ioman->GetObject(fBranchGemTracks.Data());
131 fGlobalTracks = (TClonesArray*) ioman->GetObject(fBranchGlobalTracks.Data());
132
133 fFairEventHeader = (FairEventHeader*) ioman->GetObject(fBranchFairEventHeader.Data());
134
135 fGemAlignCorr = new TClonesArray(fBranchGemAlignCorr.Data());
136 fSiAlignCorr = new TClonesArray(fBranchSiAlignCorr.Data());
137
138 ioman->Register(fBranchGemAlignCorr.Data(), "GEM", fGemAlignCorr, kTRUE);
139 ioman->Register(fBranchSiAlignCorr.Data(), "SI", fSiAlignCorr, kTRUE);
140
141 fChain = ioman->GetInChain();
142
143 Char_t* geoFileName = (Char_t*) "current_geo_file.root";
144 Int_t res_code = UniRun::ReadGeometryFile(fRunPeriod, fRunId, geoFileName);
145 if (res_code != 0) {
146 cout << "Geometry file can't be read from the database" << endl;
147 exit(-1);
148 }
149 TGeoManager::Import(geoFileName);
150
151 if (fIsDoTest) {
152 TChain* ch = new TChain("bmndata");
153 ch->Add(fMisAlignFile.Data());
154 ch->SetBranchAddress(fBranchGemAlignCorr.Data(), &fBmnGemMisalign);
155 ch->GetEntry(0);
156 delete ch;
157 }
158 return kSUCCESS;
159}
160
161void BmnGlobalAlignment::Exec(Option_t* opt) {
162 fFairEventHeader->SetMCEntryNumber(fCurrentEvent);
163 fFairEventHeader->SetRunId(fRunId);
164 fCurrentEvent++;
165 if (fCurrentEvent % 1000 == 0)
166 cout << "Event# = " << fCurrentEvent << endl;
167
168 if (!fGlobalTracks)
169 return;
170
171 for (Int_t iGlobTrack = 0; iGlobTrack < fGlobalTracks->GetEntriesFast(); iGlobTrack++) {
172 BmnGlobalTrack* globTrack = (BmnGlobalTrack*) fGlobalTracks->UncheckedAt(iGlobTrack);
173 if (fDetectorSet[1] && globTrack->GetSilTrackIndex() == -1)
174 continue;
175
176 FairTrackParam* params = globTrack->GetParamFirst();
177
178 Double_t chi2 = globTrack->GetChi2();
179 Int_t ndf = globTrack->GetNDF();
180 Int_t nHits = globTrack->GetNHits();
181 Double_t Tx = params->GetTx();
182 Double_t Ty = params->GetTy();
183
184 // Use track constraints if necessary
185 if (Tx < fTxMin || Tx > fTxMax || Ty < fTyMin || Ty > fTyMax || nHits < fMinHitsAccepted || chi2 / ndf > fChi2MaxPerNDF)
186 continue;
187
188 // Exclude a range from the selected range of track. params (in order not to take into account tracks with almost zero values of track params.)
189 if (fIsExcludedTx && Tx > fTxLeft && Tx < fTxRight)
190 continue;
191 if (fIsExcludedTy && Ty > fTyLeft && Ty < fTyRight)
192 continue;
193
194 Int_t* idx = new Int_t[nDetectors];
195 idx[0] = globTrack->GetGemTrackIndex();
196 idx[1] = globTrack->GetSilTrackIndex();
197
198 vector <BmnMilleContainer*> GEM;
199 vector <BmnMilleContainer*> SILICON;
200
201 for (Int_t iDet = 0; iDet < nDetectors; iDet++) {
202 if (!fDetectorSet[iDet] || idx[iDet] == -1 || fIsField) //{
203 // delete [] idx;
204 continue;
205 //}
206 MilleNoFieldRuns(globTrack, idx[iDet], iDet, GEM, SILICON);
207 }
208
209 fCONTAINER[fNTracks] = pair <vector <BmnMilleContainer*>, vector < BmnMilleContainer*>> (SILICON, GEM);
210 fNTracks++;
211 delete [] idx;
212 }
213
214 if (fNEvents == fCurrentEvent) {
215 MakeBinFile();
216 MakeSteerFile();
217 Pede();
218 }
219}
220
221BmnMilleContainer* BmnGlobalAlignment::FillMilleContainer(BmnGlobalTrack* glTrack, BmnHit* hit) {
222 BmnMilleContainer* mille = new BmnMilleContainer();
223
224 mille->SetStation(hit->GetStation());
225 mille->SetModule(hit->GetModule());
226
227 vector <Double_t> locDerX;
228 vector <Double_t> locDerY;
229 locDerX.push_back(1.);
230 locDerX.push_back(hit->GetZ());
231 locDerX.push_back(0.);
232 locDerX.push_back(0.);
233
234 locDerY.push_back(0.);
235 locDerY.push_back(0.);
236 locDerY.push_back(1.);
237 locDerY.push_back(hit->GetZ());
238
239 vector <Double_t> globDerX;
240 vector <Double_t> globDerY;
241 globDerX.push_back(1.);
242 globDerX.push_back(0.);
243 globDerX.push_back(glTrack->GetParamFirst()->GetTx());
244
245 globDerY.push_back(0.);
246 globDerY.push_back(1.);
247 globDerY.push_back(glTrack->GetParamFirst()->GetTy());
248
249 mille->SetLocDers(locDerX, locDerY);
250 mille->SetGlobDers(globDerX, globDerY);
251
252 mille->SetMeasures(hit->GetX(), hit->GetY());
253 mille->SetDMeasures(fUseRealHitErrors ? 2. * hit->GetDx() : 1., fUseRealHitErrors ? hit->GetDy() : 1.);
254
255 return mille;
256}
257
258void BmnGlobalAlignment::MilleNoFieldRuns(BmnGlobalTrack* glTrack, Int_t idx, Int_t iDet, vector <BmnMilleContainer*>& GEM, vector <BmnMilleContainer*>& SILICON) {
259 // GEM
260 if (iDet == 0) {
261 BmnTrack* track = (BmnTrack*) fGemTracks->UncheckedAt(idx);
262 for (Int_t iHit = 0; iHit < track->GetNHits(); iHit++) {
263 BmnGemStripHit* hit = (BmnGemStripHit*) fGemHits->UncheckedAt(track->GetHitIndex(iHit));
264 GEM.push_back(FillMilleContainer(glTrack, hit));
265 }
266 }
267
268 // SILICON
269 if (iDet == 1) {
270 BmnTrack* track = (BmnTrack*) fSilTracks->UncheckedAt(idx);
271 for (Int_t iHit = 0; iHit < track->GetNHits(); iHit++) {
272 BmnSiliconHit* hit = (BmnSiliconHit*) fSiHits->UncheckedAt(track->GetHitIndex(iHit));
273 SILICON.push_back(FillMilleContainer(glTrack, hit));
274 }
275 }
276}
277
278void BmnGlobalAlignment::_Mille(Double_t* DerLc, Double_t* DerGl, BmnMille* Mille, Int_t iTrack) {
279 fITERATOR = next(fCONTAINER.begin(), iTrack);
280
281 TString* dets = new TString[nDetectors];
282 dets[0] = "GEM";
283 dets[1] = "SILICON";
284
285 for (Int_t iDet = 0; iDet < nDetectors; iDet++) {
286 vector <BmnMilleContainer*> cont;
287 if (dets[iDet].Contains("GEM"))
288 cont = fITERATOR->second.second;
289 else if (dets[iDet].Contains("SILICON"))
290 cont = fITERATOR->second.first;
291 else {
292 cout << "Wrong input data! Exiting ..." << endl;
293 throw;
294 }
295
296 const Int_t nLays = 2;
297 TString lays[nLays] = {"X", "Y"};
298
299 for (size_t iEle = 0; iEle < cont.size(); iEle++)
300 {
301 BmnMilleContainer* _cont = cont[iEle];
302 for (Int_t iLoc = 0; iLoc < fNLC; iLoc++)
303 DerLc[iLoc] = 0.;
304 for (Int_t iGlob = 0; iGlob < fNGL; iGlob++)
305 DerGl[iGlob] = 0.;
306
307 for (Int_t iLay = 0; iLay < nLays; iLay++) {
308 vector <Double_t> locDers = _cont->GetLocDers(lays[iLay]);
309 DerLc[0] = locDers[0];
310 DerLc[1] = locDers[1];
311 DerLc[2] = locDers[2];
312 DerLc[3] = locDers[3];
313
314 Int_t idx = dets[iDet].Contains("GEM") ? GemStatModLabel(_cont->GetStation(), _cont->GetModule()) :
315 SiliconStatModLabel(_cont->GetStation(), _cont->GetModule());
316 Int_t idxUp = idx - 1; // Upper label value
317 Int_t idxMed = idxUp - 1;
318 Int_t idxLow = idxMed - 1;
319
320 vector <Double_t> globDers = _cont->GetGlobDers(lays[iLay]);
321 DerGl[idxLow] = globDers[0];
322 DerGl[idxMed] = globDers[1];
323 DerGl[idxUp] = globDers[2];
324
325 Double_t meas = lays[iLay].Contains("X") ? _cont->GetMeasures().X() : _cont->GetMeasures().Y();
326 Double_t dmeas = lays[iLay].Contains("X") ? _cont->GetDMeasures().X() : _cont->GetDMeasures().Y();
327 Mille->mille(fNLC, DerLc, fNGL, DerGl, Labels, meas, dmeas);
328
329 if (fDebug) {
330 cout << dets[iDet] << endl;
331 cout << lays[iLay] << endl;
332 cout << _cont->GetStation() << " " << _cont->GetModule() << endl;
333 cout << "Loc. ders : ";
334 for (Int_t iLoc = 0; iLoc < fNLC; iLoc++)
335 cout << DerLc[iLoc] << " ";
336 cout << endl;
337 cout << "Glob. ders : ";
338 for (Int_t iGlob = 0; iGlob < fNGL; iGlob++)
339 cout << DerGl[iGlob] << " ";
340 cout << endl;
341 cout << "Meas = " << meas << " dMeas = " << dmeas << endl;
342 }
343 }
344 }
345 }
346 delete [] dets;
347}
348
349const Int_t BmnGlobalAlignment::MakeBinFile() {
350 const Int_t ngl_per_subdetector = 3; // x, y and z corrs to each det. subsyst. (GEM, SILICON at the moment)
351
352 // GEM
353 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
354 fNGL += ngl_per_subdetector * fDetectorGEM->GetGemStation(iStat)->GetNModules();
355
356 // SILICON
357 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
358 fNGL += ngl_per_subdetector * fDetectorSI->GetSiliconStation(iStat)->GetNModules();
359
360 // Array with labels
361 Labels = new Int_t[fNGL];
362 for (Int_t iEle = 0; iEle < fNGL; iEle++)
363 Labels[iEle] = 1 + iEle;
364
365 // Read collected tracks
366 // Arrays with loc. and glob. derivatives
367 Int_t nTracks = fCONTAINER.size();
368 Double_t DerLc[fNLC];
369 Double_t DerGl[fNGL];
370
371 BmnMille* Mille = new BmnMille("alignment.bin", kTRUE, kFALSE);
372 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
373 if (fDebug)
374 cout << "Track# " << iTrack << endl;
375 _Mille(DerLc, DerGl, Mille, iTrack);
376 Mille->end();
377 if (iTrack % 1000 == 0 && !fDebug)
378 cout << "Track processed# " << iTrack << endl;
379 }
380
381 delete Mille;
382 return 0;
383}
384
385void BmnGlobalAlignment::MakeSteerFile() {
386 // Double_t shiftX = 0., shiftY = 0., shiftZ = 0.;
387 // Double_t*** misAlign = new Double_t**[fDetectorGEM->GetNStations()];
388 // for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++) {
389 // misAlign[iStat] = new Double_t*[fDetectorGEM->GetGemStation(iStat)->GetNModules()];
390 // for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++) {
391 // misAlign[iStat][iMod] = new Double_t[3];
392 // for (Int_t iPar = 0; iPar < 3; iPar++)
393 // misAlign[iStat][iMod][iPar] = 0.;
394 // }
395 // }
396
397 // if (fIsDoTest && fBmnGemMisalign) {
398 // for (Int_t iCorr = 0; iCorr < fBmnGemMisalign->GetEntriesFast(); iCorr++) {
399 // BmnGemAlignCorrections* align = (BmnGemAlignCorrections*) fBmnGemMisalign->UncheckedAt(iCorr);
400 // Int_t iStat = align->GetStation();
401 // Int_t iMod = align->GetModule();
402 // TVector3 corrsXYZ = align->GetCorrections();
403 // misAlign[iStat][iMod][0] = corrsXYZ.X();
404 // shiftX += misAlign[iStat][iMod][0];
405 // misAlign[iStat][iMod][1] = corrsXYZ.Y();
406 // shiftY += misAlign[iStat][iMod][1];
407 // misAlign[iStat][iMod][2] = corrsXYZ.Z();
408 // shiftZ += misAlign[iStat][iMod][2];
409 // }
410 // }
411
412 FILE* steer = fopen("steer.txt", "w");
413 TString alignType = "alignment.bin";
414 fprintf(steer, "%s\n", alignType.Data());
415 fprintf(steer, "method inversion %d %f\n", fNumOfIterations, fAccuracy);
416 if (fUseRegularization)
417 fprintf(steer, "regularization 1.0\n");
418 fprintf(steer, "hugecut %G\n", fHugecut);
419 if (fChisqcut[0] * fChisqcut[1] != 0)
420 fprintf(steer, "chisqcut %G %G\n", fChisqcut[0], fChisqcut[1]);
421 fprintf(steer, "entries %d\n", fEntries);
422 fprintf(steer, "outlierdownweighting %d\n", fOutlierdownweighting);
423 fprintf(steer, "dwfractioncut %G\n", fDwfractioncut);
424 fprintf(steer, "Parameter\n");
425
426 Int_t parCounterGem = 0;
427 Int_t parCounterSi = 0;
428 const Int_t nParams = 3;
429
430 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
431 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
432 parCounterGem++;
433
434 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
435 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
436 parCounterSi++;
437
438 Int_t startIdx = 0;
439 for (Int_t iDet = 0; iDet < nDetectors; iDet++) {
440 if (iDet == 0) {// Process GEMs to mark fixed stations if exist
441 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++) {
442 for (Int_t iPar = 0; iPar < fDetectorGEM->GetGemStation(iStat)->GetNModules() * nParams; iPar++) {
443 if (!fDetectorSet[0])
444 fprintf(steer, "%d %G %G\n", Labels[startIdx + iPar], 0., -1.);
445 else
446 fprintf(steer, "%d %G %G\n", Labels[startIdx + iPar], 0., (fixedGemElements[iStat][iPar / nParams]) ? -1. : fPreSigma);
447 }
448 startIdx += fDetectorGEM->GetGemStation(iStat)->GetNModules() * nParams;
449 }
450 } else if (iDet == 1) { // Process SILICON to mark fixed modules if exist
451 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++) {
452 for (Int_t iPar = 0; iPar < fDetectorSI->GetSiliconStation(iStat)->GetNModules() * nParams; iPar++) {
453 if (!fDetectorSet[1])
454 fprintf(steer, "%d %G %G\n", Labels[startIdx + iPar], 0., -1.);
455 else
456 fprintf(steer, "%d %G %G\n", Labels[startIdx + iPar], 0., (fixedSiElements[iStat][iPar / nParams]) ? -1. : (fIsDoTest) ? -1. : fPreSigma);
457 }
458 startIdx += fDetectorSI->GetSiliconStation(iStat)->GetNModules() * nParams;
459 }
460 }
461 }
462
463 if (!fUseConstraints) {
464 fclose(steer);
465 return;
466 }
467
468 // Calculate center-of-gravity along Z-axis (GEM and SI)
469 map <pair <Int_t, Int_t>, Double_t> z_GEM;
470 map <pair <Int_t, Int_t>, Double_t> z_SI;
471
472 if (fDetectorSet[0]) {
473 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
474 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
475 if (!fixedGemElements[iStat][iMod])
476 z_GEM[pair <Int_t, Int_t> (iStat, iMod)] = fDetectorGEM->GetGemStation(iStat)->GetModule(iMod)->GetZPositionRegistered();
477 }
478
479 if (fDetectorSet[1]) {
480 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
481 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
482 if (!fixedSiElements[iStat][iMod])
483 z_SI[pair <Int_t, Int_t> (iStat, iMod)] = fDetectorSI->GetSiliconStation(iStat)->GetModule(iMod)->GetZPositionRegistered();
484 }
485
486 Double_t zSum = 0.;
487 for (auto it : z_GEM)
488 zSum += it.second;
489
490 for (auto it : z_SI)
491 zSum += it.second;
492
493 Double_t zC = zSum / (z_GEM.size() + z_SI.size());
494 cout << zC << " " << z_GEM.size() + z_SI.size() << endl;
495
496 // Calculate dZ = Zpos - Zc (GEM and SI)
497 map <pair <Int_t, Int_t>, Double_t> dZ_GEM;
498 map <pair <Int_t, Int_t>, Double_t> dZ_SI;
499
500 for (auto it : z_GEM)
501 dZ_GEM[it.first] = it.second - zC;
502
503 for (auto it : z_SI)
504 dZ_SI[it.first] = it.second - zC;
505
506 map <Int_t, Double_t> deltas; // Label --> Correct dZ
507
508 for (auto it : dZ_GEM) {
509 Int_t stat = it.first.first;
510 Int_t mod = it.first.second;
511 Int_t labZ = GemStatModLabel(stat, mod);
512 Int_t labY = labZ - 1;
513 Int_t labX = labY - 1;
514 deltas[labX] = it.second;
515 deltas[labY] = it.second;
516 deltas[labZ] = it.second + zC;
517 }
518
519 for (auto it : dZ_SI) {
520 Int_t stat = it.first.first;
521 Int_t mod = it.first.second;
522 Int_t labZ = SiliconStatModLabel(stat, mod);
523 Int_t labY = labZ - 1;
524 Int_t labX = labY - 1;
525 deltas[labX] = it.second;
526 deltas[labY] = it.second;
527 deltas[labZ] = it.second + zC;
528 }
529
530 if (!fDetectorSet[0] && !fDetectorSet[1])
531 return;
532
533 // Apply constraints ...
534 // Wi * a_Xi = 0, (iStep = 0), Wi = 1.
535 // Wi * a_Yi = 0, (iStep = 1), Wi = 1.
536 // Wi * a_Zi = 0, (iStep = 2), Wi = 1.
537 for (Int_t iRemain = 0; iRemain < nParams; iRemain++) {
538 fprintf(steer, "constraint 0.0\n");
539 if (fDetectorSet[0])
540 for (Int_t iPar = 0; iPar < parCounterGem * nParams; iPar++) {
541 Int_t stat = GetGemStatMod(Labels[iPar])[0];
542 Int_t mod = GetGemStatMod(Labels[iPar])[1];
543 if (Labels[iPar] % nParams == iRemain && !fixedGemElements[stat][mod])
544 fprintf(steer, "%d %G\n", Labels[iPar], 1.);
545 }
546
547 if (fDetectorSet[1])
548 for (Int_t iPar = parCounterGem * nParams; iPar < fNGL; iPar++) {
549 Int_t stat = GetSiliconStatMod(Labels[iPar])[0];
550 Int_t mod = GetSiliconStatMod(Labels[iPar])[1];
551 if (Labels[iPar] % nParams == iRemain && !fixedSiElements[stat][mod])
552 fprintf(steer, "%d %G\n", Labels[iPar], 1.);
553 }
554 }
555
556 // Wi * a_Xi = 0, (iStep = 0), Wi = delta_Zi
557 // Wi * a_Yi = 0, (iStep = 1), Wi = delta_Zi
558 // Wi * a_Zi = 0, (iStep = 2), Wi = Zi
559 for (Int_t iRemain = 0; iRemain < nParams; iRemain++) {
560 fprintf(steer, "constraint 0.0\n");
561 for (auto it : deltas) {
562 if (it.first <= parCounterGem * nParams) {
563 Int_t stat = GetGemStatMod(it.first)[0];
564 Int_t mod = GetGemStatMod(it.first)[1];
565 if (it.first % nParams == iRemain && !fixedGemElements[stat][mod])
566 fprintf(steer, "%d %G\n", it.first, it.second);
567 } else {
568 Int_t stat = GetSiliconStatMod(it.first)[0];
569 Int_t mod = GetSiliconStatMod(it.first)[1];
570 if (it.first % nParams == iRemain && !fixedSiElements[stat][mod])
571 fprintf(steer, "%d %G\n", it.first, it.second);
572 }
573 }
574 }
575 fclose(steer);
576}
577
578void BmnGlobalAlignment::Pede() {
579 system("pede steer.txt");
580 ifstream resFile("millepede.res", ios::in);
581 ReadPedeOutput(resFile);
582 resFile.close();
583}
584
585void BmnGlobalAlignment::ReadPedeOutput(ifstream& resFile) {
586 if (!resFile) {
587 cout << "BmnGlobalAlignment::ReadPedeOutput" << " No input file found!!" << endl;
588 throw;
589 }
590 resFile.ignore(numeric_limits<streamsize>::max(), '\n');
591
592 const Int_t nParams = 3;
593 Double_t* corrs = new Double_t[nParams];
594
595 // Read GEMs
596 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++) {
597 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++) {
598 ExtractCorrValues(resFile, corrs);
599 BmnGemAlignCorrections* gemCorrs = new((*fGemAlignCorr)[fGemAlignCorr->GetEntriesFast()]) BmnGemAlignCorrections();
600 gemCorrs->SetStation(iStat);
601 gemCorrs->SetModule(iMod);
602 gemCorrs->SetCorrections(corrs);
603 }
604 }
605
606 // Read SILICON
607 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++) {
608 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++) {
609 ExtractCorrValues(resFile, corrs);
610 BmnSiliconAlignCorrections* siCorrs = new((*fSiAlignCorr)[fSiAlignCorr->GetEntriesFast()]) BmnSiliconAlignCorrections();
611 siCorrs->SetStation(iStat);
612 siCorrs->SetModule(iMod);
613 siCorrs->SetCorrections(corrs);
614 }
615 }
616
617 delete [] corrs;
618}
619
620void BmnGlobalAlignment::ExtractCorrValues(ifstream& resFile, Double_t* corrs) {
621 const Int_t nParams = 3;
622 TString parValue = "", dummy = "";
623 string line;
624
625 for (Int_t iCorr = 0; iCorr < nParams; iCorr++)
626 corrs[iCorr] = 0.;
627 for (Int_t iLine = 0; iLine < nParams; iLine++) {
628 getline(resFile, line);
629 stringstream ss(line);
630 Int_t size = ss.str().length();
631 // 40 and 68 symbols are fixed in the Pede-output by a given format
632 if (size == 40)
633 ss >> dummy >> parValue >> dummy;
634 else if (size == 68)
635 ss >> dummy >> parValue >> dummy >> dummy >> dummy;
636 else
637 cout << "Unsupported format given!" << endl;
638
639 Int_t idx = (iLine % nParams == 0) ? 0 : (iLine % nParams == 1) ? 1 : 2;
640 corrs[idx] = -parValue.Atof();
641 }
642}
643
644void BmnGlobalAlignment::CreateDetectorGeometries() {
645 fDetectorSet = new Bool_t[nDetectors]();
646
647 // Choose appropriate configuration (BM@N or SRC)
648 Bool_t isSRC = (fRunId < 3589) ? kTRUE : kFALSE; // FIXME
649 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
650 TString confSi = "SiliconRun" + TString(isSRC ? "SRC" : "") + "Spring2018.xml";
651 TString confGem = "GemRun" + TString(isSRC ? "SRC" : "") + "Spring2018.xml";
652
654 TString gPathSiliconConfig = gPathConfig + "/parameters/silicon/XMLConfigs/";
655 fDetectorSI = new BmnSiliconStationSet(gPathSiliconConfig + confSi);
656
657 // Define fixed elements of SI-detector...
658 fixedSiElements = new Bool_t*[fDetectorSI->GetNStations()];
659 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++) {
660 fixedSiElements[iStat] = new Bool_t[fDetectorSI->GetSiliconStation(iStat)->GetNModules()];
661 }
662
663 for (Int_t iStat = 0; iStat < fDetectorSI->GetNStations(); iStat++)
664 for (Int_t iMod = 0; iMod < fDetectorSI->GetSiliconStation(iStat)->GetNModules(); iMod++)
665 fixedSiElements[iStat][iMod] = kFALSE;
666
668 TString gPathGemConfig = gPathConfig + "/parameters/gem/XMLConfigs/";
669 fDetectorGEM = new BmnGemStripStationSet(gPathGemConfig + confGem);
670
671 // Define fixed elements of GEM-detector...
672 fixedGemElements = new Bool_t*[fDetectorGEM->GetNStations()];
673 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++) {
674 fixedGemElements[iStat] = new Bool_t[fDetectorGEM->GetGemStation(iStat)->GetNModules()];
675 }
676
677 for (Int_t iStat = 0; iStat < fDetectorGEM->GetNStations(); iStat++)
678 for (Int_t iMod = 0; iMod < fDetectorGEM->GetGemStation(iStat)->GetNModules(); iMod++)
679 fixedGemElements[iStat][iMod] = kFALSE;
680}
681
BmnAlignCorrections::SetCorrections
void SetCorrections(Double_t x, Double_t y, Double_t z)
Definition BmnAlignCorrections.h:20
BmnGemStripModule::GetZPositionRegistered
Double_t GetZPositionRegistered()
Definition BmnGemStripModule.cxx:73
BmnGemStripStationSet::GetGemStation
BmnGemStripStation * GetGemStation(Int_t station_num)
Definition BmnGemStripStationSet.cxx:150
BmnGemStripStation::GetModule
BmnGemStripModule * GetModule(Int_t module_num)
Definition BmnGemStripStation.cxx:183
BmnGlobalAlignment::Init
virtual InitStatus Init()
Definition BmnGlobalAlignment.cxx:103
BmnGlobalAlignment::Exec
virtual void Exec(Option_t *opt)
Definition BmnGlobalAlignment.cxx:161
BmnGlobalTrack::GetGemTrackIndex
Int_t GetGemTrackIndex() const
Definition BmnGlobalTrack.h:36
BmnGlobalTrack::GetSilTrackIndex
Int_t GetSilTrackIndex() const
Definition BmnGlobalTrack.h:37
BmnHit::GetModule
Int_t GetModule()
Definition BmnHit.h:77
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:45
BmnMilleContainer::SetModule
void SetModule(Int_t mod)
Definition BmnMilleContainer.h:24
BmnMilleContainer::GetGlobDers
vector< Double_t > GetGlobDers(TString lay)
Definition BmnMilleContainer.h:46
BmnMilleContainer::SetMeasures
void SetMeasures(Double_t measX, Double_t measY)
Definition BmnMilleContainer.h:50
BmnMilleContainer::SetLocDers
void SetLocDers(vector< Double_t > vecX, vector< Double_t > vecY)
Definition BmnMilleContainer.h:32
BmnMilleContainer::GetLocDers
vector< Double_t > GetLocDers(TString lay)
Definition BmnMilleContainer.h:37
BmnMilleContainer::SetGlobDers
void SetGlobDers(vector< Double_t > vecX, vector< Double_t > vecY)
Definition BmnMilleContainer.h:41
BmnMilleContainer::SetDMeasures
void SetDMeasures(Double_t dMeasX, Double_t dMeasY)
Definition BmnMilleContainer.h:59
BmnMilleContainer::SetStation
void SetStation(Int_t stat)
Definition BmnMilleContainer.h:16
BmnMille::end
void end()
Write buffer (set of derivatives with same local parameters) to file.
Definition BmnMille.cxx:135
BmnMille::mille
void mille(Int_t NLC, const Double_t *derLc, Int_t NGL, const Double_t *derGl, const Int_t *label, Double_t rMeas, Double_t sigma)
Add measurement to buffer.
Definition BmnMille.cxx:44
BmnSiliconModule::GetZPositionRegistered
Double_t GetZPositionRegistered()
Definition BmnSiliconModule.cxx:53
BmnSiliconStationSet::GetSiliconStation
BmnSiliconStation * GetSiliconStation(Int_t station_num)
Definition BmnSiliconStationSet.cxx:89
BmnSiliconStation::GetModule
BmnSiliconModule * GetModule(Int_t module_num)
Definition BmnSiliconStation.cxx:95
BmnTrack::GetNDF
Int_t GetNDF() const
Definition BmnTrack.h:60
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:72
BmnTrack::GetChi2
Float_t GetChi2() const
Definition BmnTrack.h:56
BmnTrack::GetNHits
Int_t GetNHits() const
Definition BmnTrack.h:44
BmnTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition BmnTrack.h:48
UniRun::GetRun
static UniRun * GetRun(int period_number, int run_number)
get run from the database
Definition UniRun.cxx:177
UniRun::ReadGeometryFile
static int ReadGeometryFile(int period_number, int run_number, char *geo_file_path, bool usePrevGeometryIfMissing=false)
Definition UniRun.cxx:1105
UniRun::GetFieldVoltage
double * GetFieldVoltage()
get field voltage of the current run
Definition UniRun.h:125
run
-clang-format