BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnTrackingQaOffline.cxx
Go to the documentation of this file.
1#include "BmnTrackingQaOffline.h"
2
3#include "GeoSetup.h"
4
5BmnTrackingQaOffline::BmnTrackingQaOffline()
6 : fPeriod(7)
7 , histoMan(nullptr)
8{
9 // Defining trigger map ...
10 fTriggerMap["Beam Trigger + BD ( > 0) & Si (> 2)"] = "BT+BD0+FD2";
11 fTriggerMap["Beam Trigger + BD (>3) & Si(>3)"] = "BT+BD3+FD3";
12 fTriggerMap["Beam Trigger + BD > 1 & Si > 3"] = "BT+BD1+FD3";
13 fTriggerMap["Beam Trigger + BD >5"] = "BT+BD5";
14 fTriggerMap["Beam Trigger + BD(>1) and Si(>2)"] = "BT+BD1+FD2";
15 fTriggerMap["Beam Trigger + BD(>=2) and FD(>=2)"] = "BT+BD1+FD1";
16 fTriggerMap["Beam Trigger + BD(>=2) and FD(>=3)"] = "BT+BD1+FD2";
17 fTriggerMap["Beam Trigger + BD(>=2) and FD(>=4)"] = "BT+BD1+FD3";
18 fTriggerMap["Beam Trigger + BD(>=3)"] = "BT+BD2";
19 fTriggerMap["Beam Trigger + BD(>=3) and FD(>=3)"] = "BT+BD2+FD2";
20 fTriggerMap["Beam Trigger + BD(>=3) and FD(>=4)"] = "BT+BD2+FD3";
21 fTriggerMap["Beam Trigger + BD>3"] = "BT+BD3";
22 fTriggerMap["Beam Trigger + Si (> 4)"] = "BT+FD4";
23 fTriggerMap["Beam Trigger + Si >3"] = "BT+FD3";
24 fTriggerMap["Beam Trigger + Si >= 3"] = "BT+FD2";
25 fTriggerMap["Beam Trigger + Veto Off + BD(>=4)"] = "BT+BD3+VetoOff";
26 fTriggerMap["Beam Trigger (BC1+BC2+T0+Veto)"] = "BT";
27 fTriggerMap["Beam BC1 + BC2"] = "BC1+BC2";
28 fTriggerMap["Beam BC1"] = "BC1";
29}
30
31BmnTrackingQaOffline::BmnTrackingQaOffline(TString setup, Int_t period, TString ext)
32 : fPeriod(period)
33 , fSetup(setup)
34 , start(-1)
35 , finish(-1)
36 , isEventAnal(kFALSE)
37 , isTrackAnal(kFALSE)
38 , isMatchAnal(kFALSE)
39 , isEoS(kFALSE)
40 , fEosHost("root://ncm.jinr.ru/")
41 , histoMan(nullptr)
42 , fPath(nullptr)
43 , hDst(nullptr)
44 , fGlobTracks(nullptr)
45 , fTof400Hits(nullptr)
46 , fTof700Hits(nullptr)
47 , fDchHits(nullptr)
48 , fCscHits(nullptr)
49 , fVertices(nullptr)
50 , fAna(nullptr)
51 , fKalman(nullptr)
52 , fMagField(nullptr)
53{
54 fAna = new FairRunAna();
55 GeoSetup gSetup;
56 string setupTag = gSetup.GetLastApprovedSetupTag(TString::Format("run%d", fPeriod).Data());
57 Int_t res_code = gSetup.LoadSetupToGeoManager(setupTag);
58 if (res_code != 0) {
59 cout << "Geometry '" << setupTag << "' can't be loaded from the Geometry Database" << endl;
60 exit(-1);
61 }
62
63 // Defining start and finish ...
64 if (period == 7 && setup.Contains("BM@N")) {
65 start = 3589;
66 finish = 5186;
67 } else if (period == 7 && setup.Contains("SRC")) {
68 start = 2041;
69 finish = 3589;
70 }
71
72 // Defining trigger map ...
73 fTriggerMap["Beam Trigger + BD ( > 0) & Si (> 2)"] = "BT+BD0+FD2";
74 fTriggerMap["Beam Trigger + BD (>3) & Si(>3)"] = "BT+BD3+FD3";
75 fTriggerMap["Beam Trigger + BD > 1 & Si > 3"] = "BT+BD1+FD3";
76 fTriggerMap["Beam Trigger + BD >5"] = "BT+BD5";
77 fTriggerMap["Beam Trigger + BD(>1) and Si(>2)"] = "BT+BD1+FD2";
78 fTriggerMap["Beam Trigger + BD(>=2) and FD(>=2)"] = "BT+BD1+FD1";
79 fTriggerMap["Beam Trigger + BD(>=2) and FD(>=3)"] = "BT+BD1+FD2";
80 fTriggerMap["Beam Trigger + BD(>=2) and FD(>=4)"] = "BT+BD1+FD3";
81 fTriggerMap["Beam Trigger + BD(>=3)"] = "BT+BD2";
82 fTriggerMap["Beam Trigger + BD(>=3) and FD(>=3)"] = "BT+BD2+FD2";
83 fTriggerMap["Beam Trigger + BD(>=3) and FD(>=4)"] = "BT+BD2+FD3";
84 fTriggerMap["Beam Trigger + BD>3"] = "BT+BD3";
85 fTriggerMap["Beam Trigger + Si (> 4)"] = "BT+FD4";
86 fTriggerMap["Beam Trigger + Si >3"] = "BT+FD3";
87 fTriggerMap["Beam Trigger + Si >= 3"] = "BT+FD2";
88 fTriggerMap["Beam Trigger + Veto Off + BD(>=4)"] = "BT+BD3+VetoOff";
89 fTriggerMap["Beam Trigger (BC1+BC2+T0+Veto)"] = "BT";
90 fTriggerMap["Beam BC1 + BC2"] = "BC1+BC2";
91 fTriggerMap["Beam BC1"] = "BC1";
92
93 // Calling histo manager ...
94 histoMan = new BmnTrackingQaOfflineDraw();
95 if (ext.Contains("root"))
96 histoMan->usedOutExtension = ext;
97}
98
99TString BmnTrackingQaOffline::GetTrigger(Int_t run)
100{
101 TObjArray* runRecord = ElogRecord::GetRecords(fPeriod, run);
102
103 TString trigger = "";
104
105 TIter it(runRecord);
106 ElogRecord* curRecord;
107 while ((curRecord = (ElogRecord*)it())) {
108 if (curRecord->GetTriggerId())
109 trigger = ElogTrigger::GetTrigger(*(curRecord->GetTriggerId()))->GetTriggerInfo();
110 }
111
112 auto trigIter = fTriggerMap.find(trigger);
113
114 if (trigIter != fTriggerMap.end())
115 trigger = trigIter->second;
116
117 return trigger;
118}
119
120pair<Int_t, Int_t> BmnTrackingQaOffline::FindBinXY(TString condition)
121{
122
123 pair<Int_t, Int_t> binXY = make_pair(-1, -1);
124
125 // Simple matching
126 if (condition == "tof400")
127 binXY = make_pair(1, 1);
128 else if (condition == "tof700")
129 binXY = make_pair(2, 2);
130 else if (condition == "csc")
131 binXY = make_pair(3, 3);
132 else if (condition == "dch")
133 binXY = make_pair(4, 4);
134
135 // Matching OneToOne
136 else if (condition == "tof400+tof700")
137 binXY = make_pair(1, 2);
138 else if (condition == "tof400+csc")
139 binXY = make_pair(1, 3);
140 else if (condition == "tof400+dch")
141 binXY = make_pair(1, 4);
142
143 else if (condition == "tof700+csc")
144 binXY = make_pair(2, 3);
145 else if (condition == "tof700+dch")
146 binXY = make_pair(2, 4);
147 else if (condition == "csc+dch")
148 binXY = make_pair(3, 4);
149
150 // Matching TwoToOne
151 else if (condition == "tof400+tof700+csc")
152 binXY = make_pair(1, 1);
153 else if (condition == "tof400+tof700+dch")
154 binXY = make_pair(1, 2);
155 else if (condition == "tof700+csc+dch")
156 binXY = make_pair(2, 2);
157 else if (condition == "tof400+csc+dch")
158 binXY = make_pair(3, 2);
159
160 // Matching ThreeToOne
161 else if (condition == "tof400+tof700+csc+dch")
162 binXY = make_pair(1, 1);
163
164 return binXY;
165}
166
167void BmnTrackingQaOffline::DoAnalisys(Bool_t anal1, Bool_t anal2, Bool_t anal3)
168{
169 if (fPath.Contains("eos"))
170 isEoS = kTRUE;
171
172 isEventAnal = anal1;
173 isTrackAnal = anal2;
174 isMatchAnal = anal3;
175
176 histoMan->fBeams = fBeams;
177 histoMan->fTargets = fTargets;
178 histoMan->fTriggers = fTriggers;
179
180 histoMan->PrepareHistos();
181
182 TH1F***** h = histoMan->GetH1();
183 TH2F* h2 = histoMan->GetH2();
184 TH2F* h3 = histoMan->GetH3();
185 TH2F***** h4 = histoMan->GetH4();
186 TH1F***** h5 = histoMan->GetH5();
187
188 Long_t nTotalGlobTrack[fBeams.size()][fTargets.size()]
189 [fTriggers.size()]; // Total number of glob. tracks for all files to be processed
190 for (size_t iBeam = 0; iBeam < fBeams.size(); iBeam++)
191 for (size_t iTarget = 0; iTarget < fTargets.size(); iTarget++)
192 for (size_t iTrigger = 0; iTrigger < fTriggers.size(); iTrigger++)
193 nTotalGlobTrack[iBeam][iTarget][iTrigger] = 0;
194
195 if (!anal1 && !anal2 && !anal3)
196 return;
197
198 for (Int_t iDst = start; iDst < finish; iDst++) {
199 UniRun* pCurrentRun = UniRun::GetRun(fPeriod, iDst);
200
201 // Check presense of the current run in DB ...
202 if (!pCurrentRun)
203 continue;
204
205 TString targ = pCurrentRun->GetTargetParticle();
206
207 // Skipping runs w/o target ...
208 if (targ == "")
209 continue;
210
211 pair<Bool_t, Int_t> isBeamInAnal = make_pair(kFALSE, -1);
212 pair<Bool_t, Int_t> isTargetInAnal = make_pair(kFALSE, -1);
213 pair<Bool_t, Int_t> isTriggerInAnal = make_pair(kFALSE, -1);
214
215 TString target = targ;
216 TString beam = pCurrentRun->GetBeamParticle();
217 TString trigger = GetTrigger(iDst);
218
219 // Check beam, target and trigger ...
220
221 for (size_t iEle = 0; iEle < fBeams.size(); iEle++) {
222 if (fBeams[iEle] == beam) {
223 isBeamInAnal.first = kTRUE;
224 isBeamInAnal.second = iEle;
225 break;
226 }
227 }
228
229 for (size_t iEle = 0; iEle < fTargets.size(); iEle++) {
230 if (fTargets[iEle] == target) {
231 isTargetInAnal.first = kTRUE;
232 isTargetInAnal.second = iEle;
233 break;
234 }
235 }
236
237 for (size_t iEle = 0; iEle < fTriggers.size(); iEle++) {
238 if (fTriggers[iEle] == trigger) {
239 isTriggerInAnal.first = kTRUE;
240 isTriggerInAnal.second = iEle;
241 break;
242 }
243 }
244
245 if (!isBeamInAnal.first || !isTargetInAnal.first || !isTriggerInAnal.first)
246 continue;
247
248 // Trying to open dst file ...
249 TChain* dst = new TChain("bmndata");
250
251 // Trying to get dst file from EoS ...
252 if (isEoS) {
253 TString path = fPath + "/" + TString::Format("bmndst_%d.root", iDst);
254 TString command = "xrdcp -DIRedirCntTimeout 30 " + fEosHost + path + " .";
255
256 gSystem->Exec(command.Data());
257 TString dstFile = TString::Format("bmndst_%d.root", iDst);
258 dst->Add(dstFile.Data());
259 } else {
260 TString path = fPath + "/" + TString::Format("bmndst_%d.root", iDst);
261 dst->Add(path.Data());
262 }
263
264 if (!dst->GetFile()) {
265 delete dst;
266 continue;
267 }
268
269 if (isMatchAnal) {
270 // Getting mag. field ...
271 fMagField = new BmnNewFieldMap("field_sp41v4_ascii_Extrap.root");
272 fMagField->SetScale(*pCurrentRun->GetFieldVoltage() / 55.87);
273 fAna->SetField(fMagField);
274 fMagField->Init();
275 fIsField = kTRUE; // FIXME
276
277 fKalman = new BmnKalmanFilter();
278 }
279
280 if (isEventAnal)
281 h2->Fill(fTriggers[isTriggerInAnal.second], fTargets[isTargetInAnal.second], dst->GetEntries() / 1e6);
282
283 dst->SetBranchAddress("DstEventHeader.", &hDst);
284 dst->SetBranchAddress("BmnGlobalTrack", &fGlobTracks);
285 dst->SetBranchAddress("BmnVertex", &fVertices);
286 dst->SetBranchAddress("BmnTof400Hit", &fTof400Hits);
287 dst->SetBranchAddress("BmnTof700Hit", &fTof700Hits);
288 dst->SetBranchAddress("BmnDchHit", &fDchHits);
289 dst->SetBranchAddress("BmnCSCHit", &fCscHits);
290
291 // Loop over events ...
292 UInt_t nEventsWithMoreThanOneTrack = 0;
293
294 if (isEventAnal) {
295 for (Int_t iEvent = 0; iEvent < dst->GetEntries(); iEvent++) {
296
297 dst->GetEntry(iEvent);
298
299 if (iEvent % 10000 == 0)
300 cout << "Processing file# " << iDst << " Event# " << iEvent << endl;
301
302 Int_t totMult = fGlobTracks->GetEntriesFast();
303 if (totMult > 1)
304 nEventsWithMoreThanOneTrack++;
305 }
306
307 h3->Fill(fTriggers[isTriggerInAnal.second], fTargets[isTargetInAnal.second],
308 nEventsWithMoreThanOneTrack / 1e6);
309 }
310
311 if (isTrackAnal)
312 for (Int_t iEvent = 0; iEvent < dst->GetEntries(); iEvent++) {
313
314 dst->GetEntry(iEvent);
315
316 if (iEvent % 10000 == 0)
317 cout << "Processing file# " << iDst << " Event# " << iEvent << endl;
318
319 CbmVertex* vertex = (CbmVertex*)fVertices->UncheckedAt(0);
320 h[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][0]->Fill(vertex->GetZ());
321
322 Int_t totMult = fGlobTracks->GetEntriesFast();
323 h[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][1]->Fill(totMult);
324
325 if (totMult > 1)
326 nEventsWithMoreThanOneTrack++;
327
328 Int_t nPos = 0;
329 Int_t nNeg = 0;
330
331 // Loop over glob. tracks ...
332 for (Int_t iTrack = 0; iTrack < fGlobTracks->GetEntriesFast(); iTrack++) {
333 BmnGlobalTrack* track = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(iTrack);
334
335 if (track->GetParamLast()->GetQp() > 0.)
336 nPos++;
337 else
338 nNeg++;
339
340 h[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][4]->Fill(track->GetP());
341 }
342
343 if (nPos != 0)
344 h[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][2]->Fill(nPos);
345
346 if (nNeg != 0)
347 h[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][3]->Fill(nNeg);
348 }
349
350 if (isMatchAnal) {
351 // Matching ...
352 TH2F* H0 = h4[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][0];
353 TH2F* H1 = h4[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][1];
354 TH2F* H2 = h4[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][2];
355
356 // Beta ...
357 TH2F* beta400 = h4[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][3];
358 TH2F* beta700 = h4[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second][4];
359
360 // Residuals ...
361 TH1F** residuals = h5[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second];
362
363 for (Int_t iEvent = 0; iEvent < dst->GetEntries(); iEvent++) {
364
365 dst->GetEntry(iEvent);
366
367 if (iEvent % 10000 == 0)
368 cout << "Processing file# " << iDst << " Event# " << iEvent << endl;
369
370 // Loop over glob. tracks ...
371 nTotalGlobTrack[isBeamInAnal.second][isTargetInAnal.second][isTriggerInAnal.second] +=
372 fGlobTracks
373 ->GetEntriesFast(); // Increasing num. of glob. tracks by each file event to be processed
374 for (Int_t iTrack = 0; iTrack < fGlobTracks->GetEntriesFast(); iTrack++) {
375 BmnGlobalTrack* track = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(iTrack);
376
377 Double_t b400 = -100; // track->GetBeta(1);
378 Double_t b700 = -100; // track->GetBeta(2);
379
380 if (b400 < 1.5)
381 beta400->Fill(TMath::Abs(track->GetP()), b400);
382 if (b700 < 1.5)
383 beta700->Fill(TMath::Abs(track->GetP()), b700);
384
385 Int_t mTof1 = track->GetTof1HitIndex();
386 Int_t mTof2 = track->GetTof2HitIndex();
387
388 Int_t mCsc = -1; // track->GetCscHitIndex(0);
389
390 Int_t mDch1 = -1; // track->GetDch1TrackIndex();
391 Int_t mDch2 = -1; // track->GetDch2TrackIndex();
392 Int_t mDch = -1; // track->GetDchTrackIndex();
393
394 BmnHit* tof400Hit = nullptr;
395 BmnHit* tof700Hit = nullptr;
396 BmnHit* dchHit = nullptr;
397 BmnHit* cscHit = nullptr;
398
399 // Simple matching ...
400 if (mTof1 != -1) {
401 tof400Hit = (BmnHit*)fTof400Hits->UncheckedAt(mTof1);
402 Double_t alreadyFilled =
403 H0->GetBinContent(FindBinXY("tof400").first, FindBinXY("tof400").second);
404 H0->SetBinContent(FindBinXY("tof400").first, FindBinXY("tof400").second, alreadyFilled + 1.);
405 }
406
407 if (mTof2 != -1) {
408 tof700Hit = (BmnHit*)fTof700Hits->UncheckedAt(mTof2);
409 Double_t alreadyFilled =
410 H0->GetBinContent(FindBinXY("tof700").first, FindBinXY("tof700").second);
411 H0->SetBinContent(FindBinXY("tof700").first, FindBinXY("tof700").second, alreadyFilled + 1.);
412 }
413
414 if (mCsc != -1) {
415 cscHit = (BmnHit*)fCscHits->UncheckedAt(mCsc);
416 Double_t alreadyFilled = H0->GetBinContent(FindBinXY("csc").first, FindBinXY("csc").second);
417 H0->SetBinContent(FindBinXY("csc").first, FindBinXY("csc").second, alreadyFilled + 1.);
418 }
419
420 if (mDch1 != -1 || mDch2 != -1 || mDch != -1) {
421 dchHit = (BmnHit*)fDchHits->UncheckedAt(mDch);
422 Double_t alreadyFilled = H0->GetBinContent(FindBinXY("dch").first, FindBinXY("dch").second);
423 H0->SetBinContent(FindBinXY("dch").first, FindBinXY("dch").second, alreadyFilled + 1.);
424 }
425
426 // OneTwoOne matching ...
427 if (mTof1 != -1 && mTof2 != -1) {
428 Double_t alreadyFilled =
429 H0->GetBinContent(FindBinXY("tof400+tof700").first, FindBinXY("tof400+tof700").second);
430 H0->SetBinContent(FindBinXY("tof400+tof700").first, FindBinXY("tof400+tof700").second,
431 alreadyFilled + 1.);
432 }
433
434 if (mTof1 != -1 && mCsc != -1) {
435 Double_t alreadyFilled =
436 H0->GetBinContent(FindBinXY("tof400+csc").first, FindBinXY("tof400+csc").second);
437 H0->SetBinContent(FindBinXY("tof400+csc").first, FindBinXY("tof400+csc").second,
438 alreadyFilled + 1.);
439 }
440
441 if (mTof1 != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
442 Double_t alreadyFilled =
443 H0->GetBinContent(FindBinXY("tof400+dch").first, FindBinXY("tof400+dch").second);
444 H0->SetBinContent(FindBinXY("tof400+dch").first, FindBinXY("tof400+dch").second,
445 alreadyFilled + 1.);
446 }
447
448 if (mTof2 != -1 && mCsc != -1) {
449 Double_t alreadyFilled =
450 H0->GetBinContent(FindBinXY("tof700+csc").first, FindBinXY("tof700+csc").second);
451 H0->SetBinContent(FindBinXY("tof700+csc").first, FindBinXY("tof700+csc").second,
452 alreadyFilled + 1.);
453 }
454
455 if (mTof2 != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
456 Double_t alreadyFilled =
457 H0->GetBinContent(FindBinXY("tof700+dch").first, FindBinXY("tof700+dch").second);
458 H0->SetBinContent(FindBinXY("tof700+dch").first, FindBinXY("tof700+dch").second,
459 alreadyFilled + 1.);
460 }
461
462 if (mCsc != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
463 Double_t alreadyFilled =
464 H0->GetBinContent(FindBinXY("csc+dch").first, FindBinXY("csc+dch").second);
465 H0->SetBinContent(FindBinXY("csc+dch").first, FindBinXY("csc+dch").second, alreadyFilled + 1.);
466 }
467
468 // TwoToOne matching ...
469 if (mTof1 != -1 && mTof2 != -1 && mCsc != -1) {
470 Double_t alreadyFilled = H1->GetBinContent(FindBinXY("tof400+tof700+csc").first,
471 FindBinXY("tof400+tof700+csc").second);
472 H1->SetBinContent(FindBinXY("tof400+tof700+csc").first, FindBinXY("tof400+tof700+csc").second,
473 alreadyFilled + 1.);
474 }
475
476 if (mTof1 != -1 && mTof2 != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
477 Double_t alreadyFilled = H1->GetBinContent(FindBinXY("tof400+tof700+dch").first,
478 FindBinXY("tof400+tof700+dch").second);
479 H1->SetBinContent(FindBinXY("tof400+tof700+dch").first, FindBinXY("tof400+tof700+dch").second,
480 alreadyFilled + 1.);
481 }
482
483 if (mTof2 != -1 && mCsc != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
484 Double_t alreadyFilled =
485 H1->GetBinContent(FindBinXY("tof700+csc+dch").first, FindBinXY("tof700+csc+dch").second);
486 H1->SetBinContent(FindBinXY("tof700+csc+dch").first, FindBinXY("tof700+csc+dch").second,
487 alreadyFilled + 1.);
488 }
489
490 if (mTof1 != -1 && mCsc != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
491 Double_t alreadyFilled =
492 H1->GetBinContent(FindBinXY("tof400+csc+dch").first, FindBinXY("tof400+csc+dch").second);
493 H1->SetBinContent(FindBinXY("tof400+csc+dch").first, FindBinXY("tof400+csc+dch").second,
494 alreadyFilled + 1.);
495 }
496
497 // Matching ThreeToOne
498 if (mTof1 != -1 && mTof2 != -1 && mCsc != -1 && (mDch != -1 || mDch1 != -1 || mDch2 != -1)) {
499 Double_t alreadyFilled = H2->GetBinContent(FindBinXY("tof400+tof700+csc+dch").first,
500 FindBinXY("tof400+tof700+csc+dch").second);
501 H2->SetBinContent(FindBinXY("tof400+tof700+csc+dch").first,
502 FindBinXY("tof400+tof700+csc+dch").second, alreadyFilled + 1.);
503 }
504
505 // Pushing vector with det. hits ...
506 vector<BmnHit*> outerTrackerHits{tof400Hit, tof700Hit, cscHit, dchHit};
507 map<Int_t, pair<Double_t, Double_t>> detIdxResXY = GetResiduals(track, outerTrackerHits);
508
509 for (auto detector : detIdxResXY) {
510 Int_t idx = detector.first;
511 pair<Double_t, Double_t> resid = detector.second;
512
513 residuals[2 * idx]->Fill(resid.first);
514 residuals[2 * idx + 1]->Fill(resid.second);
515 }
516 }
517 }
518 }
519
520 delete dst;
521
522 if (isEoS)
523 gSystem->Exec(Form("rm bmndst_%d.root", iDst));
524
525 if (isMatchAnal) {
526 delete fKalman;
527 delete fMagField;
528 }
529 }
530
531 // Normalizing in the end ...
532 for (size_t iBeam = 0; iBeam < fBeams.size(); iBeam++)
533 for (size_t iTarget = 0; iTarget < fTargets.size(); iTarget++)
534 for (size_t iTrigger = 0; iTrigger < fTriggers.size(); iTrigger++) {
535 if (nTotalGlobTrack[iBeam][iTarget][iTrigger] == 0)
536 continue;
537
538 TH2F* H0 = h4[iBeam][iTarget][iTrigger][0];
539 TH2F* H1 = h4[iBeam][iTarget][iTrigger][1];
540 TH2F* H2 = h4[iBeam][iTarget][iTrigger][2];
541
542 vector<TH2F*> H{H0, H1, H2};
543
544 for (TH2F* histo : H) {
545 for (Int_t iBinX = 1; iBinX < histo->GetNbinsX() + 1; iBinX++)
546 for (Int_t iBinY = 1; iBinY < histo->GetNbinsY() + 1; iBinY++) {
547 Double_t contentBefore = histo->GetBinContent(iBinX, iBinY);
548
549 Double_t contentAfter = contentBefore / nTotalGlobTrack[iBeam][iTarget][iTrigger];
550 histo->SetBinContent(iBinX, iBinY, 100. * contentAfter);
551 }
552 }
553 }
554}
555
556map<Int_t, pair<Double_t, Double_t>> BmnTrackingQaOffline::GetResiduals(BmnGlobalTrack* track, vector<BmnHit*> hits)
557{
558 // hits --> tof400, tof700, csc, dch
559
560 map<Int_t, pair<Double_t, Double_t>> detResidXY;
561 pair<Double_t, Double_t> residXY = make_pair(-1000., -1000.);
562
563 for (size_t iDet = 0; iDet < hits.size(); iDet++)
564 detResidXY[iDet] = residXY;
565
566 Double_t Qp = track->GetParamLast()->GetQp();
567 FairTrackParam par(*(track->GetParamFirst()));
568
569 // Loop over detectors ...
570 for (size_t iDet = 0; iDet < hits.size(); iDet++) {
571 if (!hits[iDet])
572 continue;
573
574 BmnHit* hit = hits[iDet];
575
576 if (fKalman->TGeoTrackPropagate(&par, hit->GetZ(), (Qp > 0) ? 2212 : -211, nullptr, nullptr, fIsField)
577 == kBMNERROR)
578 continue;
579
580 Double_t xRes = par.GetX() - hit->GetX();
581 Double_t yRes = par.GetY() - hit->GetY();
582
583 residXY = make_pair(xRes, yRes);
584
585 detResidXY[iDet] = residXY;
586 }
587
588 return detResidXY;
589}
Bool_t
bool Bool_t
Int_t
int Int_t
Long_t
long Long_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Form
char * Form(const char *fmt,...)
gSystem
R__EXTERN TSystem * gSystem
kBMNERROR
@ kBMNERROR
Error.
Definition BmnEnums.h:95
BmnGlobalTrack
BmnGlobalTrack.h.
Definition BmnGlobalTrack.h:14
BmnGlobalTrack::GetTof2HitIndex
Int_t GetTof2HitIndex() const
Definition BmnGlobalTrack.h:37
BmnGlobalTrack::GetTof1HitIndex
Int_t GetTof1HitIndex() const
Definition BmnGlobalTrack.h:36
BmnNewFieldMap
Field Map does not have any symmetries w.r.t coordinate axises.
Definition BmnNewFieldMap.h:16
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:63
BmnTrack::GetP
Double_t GetP()
Definition BmnTrack.h:67
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:65
BmnTrackingQaOffline::fTriggers
vector< TString > fTriggers
Definition BmnTrackingQaOffline.h:78
BmnTrackingQaOffline::FindBinXY
pair< Int_t, Int_t > FindBinXY(TString)
Definition BmnTrackingQaOffline.cxx:120
BmnTrackingQaOffline::fBeams
std::vector< TString > fBeams
Definition BmnTrackingQaOffline.h:76
BmnTrackingQaOffline::fTargets
vector< TString > fTargets
Definition BmnTrackingQaOffline.h:77
BmnTrackingQaOffline::fKalman
BmnKalmanFilter * fKalman
Definition BmnTrackingQaOffline.h:92
BmnTrackingQaOffline::histoMan
BmnTrackingQaOfflineDraw * histoMan
Definition BmnTrackingQaOffline.h:80
BmnTrackingQaOffline::GetResiduals
map< Int_t, pair< Double_t, Double_t > > GetResiduals(BmnGlobalTrack *, vector< BmnHit * >)
Definition BmnTrackingQaOffline.cxx:556
BmnTrackingQaOffline::fTriggerMap
map< TString, TString > fTriggerMap
Definition BmnTrackingQaOffline.h:74
BmnTrackingQaOffline::DoAnalisys
void DoAnalisys(Bool_t evAnal=kTRUE, Bool_t trackAnal=kTRUE, Bool_t matchAnal=kFALSE)
Definition BmnTrackingQaOffline.cxx:167
CbmVertex
CbmVertex.h.
Definition CbmVertex.h:21
CbmVertex::GetZ
Double_t GetZ() const
Definition CbmVertex.h:60
ElogRecord
ElogRecord.h Class for the table: record_.
Definition ElogRecord.h:20
ElogRecord::GetTriggerId
int * GetTriggerId()
get trigger id of the current record
Definition ElogRecord.h:145
ElogRecord::GetRecords
static TObjArray * GetRecords(int period_number, int run_number, int findNeighbourRun=0)
get array of ElogRecord-s for a given or a neighbour run (the same period) from the database if findN...
Definition ElogRecord.cxx:1215
ElogTrigger::GetTrigger
static ElogTrigger * GetTrigger(int trigger_id)
get trigger from the database
Definition ElogTrigger.cxx:94
ElogTrigger::GetTriggerInfo
TString GetTriggerInfo()
get trigger info of the current trigger
Definition ElogTrigger.h:61
GeoSetup
GeoSetup class.
Definition GeoSetup.h:40
GeoSetup::LoadSetupToGeoManager
int LoadSetupToGeoManager(string_view setup_tag)
Load selected geometry setup from the Local Geometry Database and create gGeoManager instance.
Definition GeoSetup.cxx:512
GeoSetup::GetLastApprovedSetupTag
string GetLastApprovedSetupTag(string_view tag_prefix)
Get tag for the last approved setup corresponding to the pattern.
Definition GeoSetup.cxx:122
TChain::SetBranchAddress
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
TChain::Add
virtual Int_t Add(const char *name, Long64_t nentries=TTree::kMaxEntries)
TChain::GetEntries
Long64_t GetEntries() const override
TChain::GetFile
TFile * GetFile() const
TChain::GetEntry
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH2::Fill
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
TH2::SetBinContent
void SetBinContent(Int_t bin, Double_t content) override
TH2::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
TString::Data
const char * Data() const
TString::Format
static TString Format(const char *fmt,...)
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
UniRun
UniRun.h Class for the table: run_.
Definition UniRun.h:21
UniRun::GetRun
static UniRun * GetRun(int period_number, int run_number)
get run from the database
Definition UniRun.cxx:161
UniRun::GetBeamParticle
TString GetBeamParticle()
get beam particle of the current run
Definition UniRun.h:109
UniRun::GetFieldVoltage
double * GetFieldVoltage()
get field voltage of the current run
Definition UniRun.h:121
UniRun::GetTargetParticle
TString GetTargetParticle()
get target particle of the current run
Definition UniRun.h:111
TMath::Abs
Double_t Abs(Double_t d)