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