BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnLambdaEmbeddingQa.cxx
Go to the documentation of this file.
1#include "BmnLambdaEmbeddingQa.h"
2
3#include "BmnSiliconPoint.h"
4#include "GeoSetup.h"
5#include "TObjString.h"
6
7BmnLambdaEmbeddingQa::BmnLambdaEmbeddingQa()
8 : fRunAna(nullptr)
9 , geoms(nullptr)
10 , fDigi(nullptr)
11 , fDst(nullptr)
12 , fEmb(nullptr)
13 , fPath(nullptr)
14 , fMagField(nullptr)
15 , fMon(nullptr)
16 , hDst(nullptr)
17 , hDigi(nullptr)
18 , fGemDigits(nullptr)
19 , fSiliconDigits(nullptr)
20 , fGemHits(nullptr)
21 , fSiliconHits(nullptr)
22 , fSiliconTracks(nullptr)
23 , fGlobTracks(nullptr)
24 , fMcTracks(nullptr)
25 , fSiliconPoints(nullptr)
26 , fGemPoints(nullptr)
27 , hXzMC(nullptr)
28 , hXzReco(nullptr)
29 , hYzMC(nullptr)
30 , hYzReco(nullptr)
31 , hGemXYProfiles(nullptr)
32 , hSiliconXYProfiles(nullptr)
33 , hXzRecoFromTracks(nullptr)
34 , hProtonMomentaEmb(nullptr)
35 , hPionMomentaEmb(nullptr)
36 , hProtonMomentaReco(nullptr)
37 , hPionMomentaReco(nullptr)
38 , hProtonNhitsEmb(nullptr)
39 , hPionNhitsEmb(nullptr)
40 , hProtonNhitsReco(nullptr)
41 , hPionNhitsReco(nullptr)
42 , hProtonNhitsRecoAll(nullptr)
43 , hPionNhitsRecoAll(nullptr)
44 , pEffProton(nullptr)
45 , pEffPion(nullptr)
46 , hStripChannel(nullptr)
47 , fPrefix("")
48 , histoMan(nullptr)
49{
50 // Get inn. tracker geometry ...
51 geoms = new BmnInnerTrackerGeometryDraw();
52}
53
54BmnLambdaEmbeddingQa::BmnLambdaEmbeddingQa(TString list)
55 : fRunAna(nullptr)
56 , geoms(nullptr)
57 , fDigi(nullptr)
58 , fDst(nullptr)
59 , fEmb(nullptr)
60 , fPath(nullptr)
61 , fMagField(nullptr)
62 , fMon(nullptr)
63 , hDst(nullptr)
64 , hDigi(nullptr)
65 , fGemDigits(nullptr)
66 , fSiliconDigits(nullptr)
67 , fGemHits(nullptr)
68 , fSiliconHits(nullptr)
69 , fSiliconTracks(nullptr)
70 , fGlobTracks(nullptr)
71 , fMcTracks(nullptr)
72 , fSiliconPoints(nullptr)
73 , fGemPoints(nullptr)
74 , hXzMC(nullptr)
75 , hXzReco(nullptr)
76 , hYzMC(nullptr)
77 , hYzReco(nullptr)
78 , hGemXYProfiles(nullptr)
79 , hSiliconXYProfiles(nullptr)
80 , hXzRecoFromTracks(nullptr)
81 , hProtonMomentaEmb(nullptr)
82 , hPionMomentaEmb(nullptr)
83 , hProtonMomentaReco(nullptr)
84 , hPionMomentaReco(nullptr)
85 , hProtonNhitsEmb(nullptr)
86 , hPionNhitsEmb(nullptr)
87 , hProtonNhitsReco(nullptr)
88 , hPionNhitsReco(nullptr)
89 , hProtonNhitsRecoAll(nullptr)
90 , hPionNhitsRecoAll(nullptr)
91 , pEffProton(nullptr)
92 , pEffPion(nullptr)
93 , hStripChannel(nullptr)
94 , fPrefix("")
95 , histoMan(nullptr)
96{
97 // Get inn. tracker geometry ...
98 geoms = new BmnInnerTrackerGeometryDraw();
99
100 // Get histo manager ...
101 histoMan = new BmnLambdaEmbeddingDrawHistos();
102
103 nInputs = -1;
104 drawFoundTracks = kFALSE;
105
106 // Strip lays to be shown by default ...
107 aLaysToShow[0] = 0;
108 aLaysToShow[1] = 1;
109
110 // Signal windows by default ...
111 fSignalWindow[GEM][0] = 0;
112 fSignalWindow[GEM][1] = LDBL_MAX;
113 fSignalWindow[SILICON][0] = 0;
114 fSignalWindow[SILICON][1] = LDBL_MAX;
115
116 // 1. Process list of passed files ...
117 ifstream f(list.Data());
118 Int_t nLines = 0;
119 for (string line; getline(f, line);)
120 nLines++;
121 f.close();
122
123 cout << "List file contains " << nLines << " inputs" << endl;
124
125 if (nLines != 0) {
126 fDigi = new TString[nLines];
127 fDst = new TString[nLines];
128 fEmb = new TString[nLines];
129 fPath = new TString[nLines];
130 }
131
132 nLines = 0;
133 ifstream g(list.Data());
134 while (!g.eof()) {
135 g >> fDigi[nLines] >> fDst[nLines] >> fEmb[nLines] >> fPath[nLines];
136 nLines++;
137 }
138 g.close();
139
140 nInputs = nLines;
141
142 for (Int_t iLine = 0; iLine < nLines; iLine++)
143 cout << fDigi[iLine] << " " << fDst[iLine] << " " << fEmb[iLine] << " " << fPath[iLine] << endl;
144}
145
146void BmnLambdaEmbeddingQa::DoEffAnalysis()
147{
148 histoMan->SetActive(1);
149
150 BmnGemStripStationSet* GEM = geoms->GetGemGeometry();
151 BmnSiliconStationSet* SILICON = geoms->GetSiliconGeometry();
152
153 TProfile2D* pEmbSilHitEff = histoMan->GetSilHitEff();
154 TProfile2D** pEmbGemHitEff = histoMan->GetGemHitEff();
155
156 TProfile2D** pEffSilStatXY = histoMan->GetSilHitEff2D();
157 TProfile2D** pEffGemStatXY = histoMan->GetGemHitEff2D();
158
159 TH1F*** hEtaLambda = histoMan->GetEtaLambda();
160 TH2F** hNhitsEmbReco = histoMan->GetNHitsEmbReco();
161
162 enum
163 {
164 proton,
165 pion
166 };
167
168 if (nInputs != -1)
169 for (Int_t iInput = 0; iInput < nInputs; iInput++) {
170 Double_t etaMin = 0.;
171 Double_t etaMax = 0.;
172
173 // Trying to get pseudorapidity range of lambda ...
174 TObjArray* tx = fPath[iInput].Tokenize("/");
175 for (Int_t i = 0; i < tx->GetEntries(); i++) {
176 TString str = ((TObjString*)(tx->At(i)))->String();
177
178 if (!str.Contains("lambda"))
179 continue;
180
181 TObjArray* ty = str.Tokenize("_");
182
183 etaMin = ((TObjString*)(ty->At(2)))->String().Atof();
184 etaMax = ((TObjString*)(ty->At(3)))->String().Atof();
185
186 break;
187 }
188
189 TChain* emb = new TChain("bmndata");
190 emb->Add(fEmb[iInput].Data());
191
192 emb->SetBranchAddress("EmbeddingMonitor.", &fMon);
193
194 TChain* dst = new TChain("bmndata");
195 dst->Add(fDst[iInput].Data());
196
197 dst->SetBranchAddress("DstEventHeader.", &hDst);
198 dst->SetBranchAddress("BmnGemStripHit", &fGemHits);
199 dst->SetBranchAddress("BmnSiliconHit", &fSiliconHits);
200 dst->SetBranchAddress("BmnGlobalTrack", &fGlobTracks);
201
202 const Int_t nCutBins = 4; // 0, 20, 40, 50 %
203 Double_t cutsNHitsFrac[nCutBins] = {0., 0.2, 0.4, 0.5};
204 Int_t nRecProtons[nCutBins];
205 Int_t nRecPions[nCutBins];
206
207 for (Int_t iBin = 0; iBin < nCutBins; iBin++) {
208 nRecProtons[iBin] = 0;
209 nRecPions[iBin] = 0;
210 }
211
212 Int_t nEmbeddedProtons = 0;
213 Int_t nEmbeddedPions = 0;
214
215 // Loop over events ...
216 Int_t counter = 0;
217 for (Int_t iEmb = 0; iEmb < emb->GetEntries(); iEmb++) {
218 emb->GetEntry(iEmb);
219
220 if (iEmb % 1000 == 0)
221 cout << "Event# " << iEmb << endl;
222
223 if (!fMon->IsEmbedded())
224 continue;
225
226 nEmbeddedProtons++;
227 nEmbeddedPions++;
228
229 UInt_t idEmb = fMon->GetEventId();
230
231 for (UInt_t iDst = idEmb - 1; iDst < idEmb + 1; iDst++) {
232
233 dst->GetEntry(iDst);
234
235 if (hDst->GetEventId() != idEmb)
236 continue;
237
238 counter++;
239
240 // Get event info ...
241 Int_t store = fMon->GetStoreVertexEvent()[0];
242 Int_t vertex = fMon->GetStoreVertexEvent()[1];
243 Int_t event = fMon->GetStoreVertexEvent()[2];
244
245 // Open lambda file ...
246 TChain* lambda = new TChain("bmndata");
247 lambda->Add(Form("%s/lambda%d_vertex%d.root", fPath[iInput].Data(), store, vertex));
248
249 lambda->SetBranchAddress("StsPoint", &fGemPoints);
250 lambda->SetBranchAddress("SiliconPoint", &fSiliconPoints);
251
252 lambda->GetEntry(event);
253
254 map<Int_t, vector<CbmStsPoint>> gemsMc;
255 for (Int_t iStat = 0; iStat < GEM->GetNStations(); iStat++) {
256 vector<CbmStsPoint> tmp;
257 gemsMc[iStat] = tmp;
258 }
259
260 map<Int_t, vector<BmnSiliconPoint>> siliconsMc;
261 for (Int_t iStat = 0; iStat < SILICON->GetNStations(); iStat++) {
262 vector<BmnSiliconPoint> tmp;
263 siliconsMc[iStat] = tmp;
264 }
265
266 // Get MC point info (GEM) ...
267 for (Int_t iPoint = 0; iPoint < fGemPoints->GetEntriesFast(); iPoint++) {
268 CbmStsPoint* point = (CbmStsPoint*)fGemPoints->UncheckedAt(iPoint);
269
270 gemsMc.find(point->GetStation())->second.push_back(*point);
271 }
272
273 // Get MC point info (SILICON) ...
274 for (Int_t iPoint = 0; iPoint < fSiliconPoints->GetEntriesFast(); iPoint++) {
275 BmnSiliconPoint* point = (BmnSiliconPoint*)fSiliconPoints->UncheckedAt(iPoint);
276
277 siliconsMc.find(point->GetStation())->second.push_back(*point);
278 }
279
280 map<Int_t, vector<BmnHit>> gemsReco;
281 map<Int_t, vector<BmnHit>> siliconsReco;
282
283 for (Int_t iStat = 0; iStat < GEM->GetNStations(); iStat++) {
284 vector<BmnHit> tmp;
285 gemsReco[iStat] = tmp;
286 }
287
288 for (Int_t iStat = 0; iStat < SILICON->GetNStations(); iStat++) {
289 vector<BmnHit> tmp;
290 siliconsReco[iStat] = tmp;
291 }
292
293 // Get GEM hit info ...
294 if (fGemHits)
295 for (Int_t iHit = 0; iHit < fGemHits->GetEntriesFast(); iHit++) {
296 BmnHit* hit = (BmnHit*)fGemHits->UncheckedAt(iHit);
297
298 gemsReco.find(hit->GetStation())->second.push_back(*hit);
299 }
300
301 // Get SILICON hit info ...
302 if (fSiliconHits)
303 for (Int_t iHit = 0; iHit < fSiliconHits->GetEntriesFast(); iHit++) {
304 BmnHit* hit = (BmnHit*)fSiliconHits->UncheckedAt(iHit);
305
306 siliconsReco.find(hit->GetStation())->second.push_back(*hit);
307 }
308
309 // Loop over siliconMc and siliconReco maps to find embedded MC-points (in digits) and reconstructed
310 // hits ...
311 for (auto mc : siliconsMc)
312 for (auto reco : siliconsReco) {
313
314 // Stations should be the same ...
315 if (mc.first != reco.first)
316 continue;
317
318 vector<BmnSiliconPoint> points = mc.second;
319 vector<BmnHit> hits = reco.second;
320
321 for (auto mcPoint : points) {
322 Bool_t isEmbedded = kFALSE;
323
324 vector<Double_t> distances;
325
326 if (hits.size() != 0)
327 for (auto recoHit : hits) {
328 Double_t xMC = mcPoint.GetX();
329 Double_t xReco = recoHit.GetX();
330
331 // Calculate distance between a hit and a point in X-direction ...
332 distances.push_back(Abs(xMC - xReco));
333 }
334
335 if (distances.size() != 0) {
336 Double_t min = *min_element(distances.begin(), distances.end());
337 if (min < 0.3)
338 isEmbedded = kTRUE;
339 }
340
341 pEmbSilHitEff->Fill(mcPoint.GetStation(), mcPoint.GetModule(), isEmbedded ? 1. : 0.);
342 pEffSilStatXY[mcPoint.GetStation()]->Fill(mcPoint.GetX(), mcPoint.GetY(),
343 isEmbedded ? 1. : 0.);
344 }
345 }
346
347 // Loop over gemMc and gemReco maps to find embedded and reconstructed hits ...
348 for (auto mc : gemsMc)
349 for (auto reco : gemsReco) {
350
351 // Stations should be the same ...
352 if (mc.first != reco.first)
353 continue;
354
355 vector<CbmStsPoint> points = mc.second;
356 vector<BmnHit> hits = reco.second;
357
358 if (hits.size() == 0)
359 continue;
360
361 for (auto mcPoint : points) {
362 Bool_t isEmbedded = kFALSE;
363
364 // Taking into account correct Z-readout position for GEM ...
365 BmnGemStripModule* mod =
366 GEM->GetStation(mcPoint.GetStation())->GetModule(mcPoint.GetModule());
367 Double_t driftCenterShift = 0.;
368 if (mod->GetElectronDriftDirection() == ForwardZAxisEDrift)
369 driftCenterShift = 0.15;
370 else
371 driftCenterShift = 0.75;
372
373 vector<Double_t> distances;
374
375 if (hits.size() != 0)
376 for (auto recoHit : hits) {
377 Double_t xMC = mcPoint.GetX(mcPoint.GetZ() + driftCenterShift);
378 Double_t xReco = recoHit.GetX();
379
380 // Calculate distance between a hit and a point in X-direction ...
381 distances.push_back(Abs(xMC - xReco));
382 }
383
384 if (distances.size() != 0) {
385 Double_t min = *min_element(distances.begin(), distances.end());
386 if (min < 0.5)
387 isEmbedded = kTRUE;
388 }
389
390 // Defining a zone on the stat, module we are considering ...
391 const Int_t nLayers = 4;
392 Bool_t isLayerActive[nLayers] = {kFALSE, kFALSE, kFALSE, kFALSE};
393
394 Bool_t isHotZone = kFALSE;
395 Bool_t isBigZone = kFALSE;
396
397 for (Int_t iLayer = 0; iLayer < GEM->GetStation(mcPoint.GetStation())
398 ->GetModule(mcPoint.GetModule())
399 ->GetNStripLayers();
400 iLayer++)
401 {
402 BmnGemStripLayer layer = GEM->GetStation(mcPoint.GetStation())
403 ->GetModule(mcPoint.GetModule())
404 ->GetStripLayer(iLayer);
405 isLayerActive[iLayer] = layer.IsPointInsideStripLayer(
406 (-1.) * mcPoint.GetX(mcPoint.GetZ() + driftCenterShift),
407 mcPoint.GetY(mcPoint.GetZ() + driftCenterShift));
408 }
409
410 if (isLayerActive[0] && isLayerActive[1])
411 isBigZone = kTRUE;
412 else if (isLayerActive[2] && isLayerActive[3])
413 isHotZone = kTRUE;
414
415 pEffGemStatXY[mcPoint.GetStation()]->Fill(
416 mcPoint.GetX(mcPoint.GetZ() + driftCenterShift),
417 mcPoint.GetY(mcPoint.GetZ() + driftCenterShift), isEmbedded ? 1. : 0.);
418
419 Int_t idx = isBigZone ? 1 : isHotZone ? 0 : -1;
420
421 if (idx != -1)
422 pEmbGemHitEff[idx]->Fill(mcPoint.GetStation(), mcPoint.GetModule(),
423 isEmbedded ? 1. : 0.);
424 }
425 }
426
427 delete lambda;
428
429 vector<BmnGlobalTrack*> protonCands;
430 vector<BmnGlobalTrack*> pionCands;
431
432 for (Int_t iGlobTrack = 0; iGlobTrack < fGlobTracks->GetEntriesFast(); iGlobTrack++) {
433 BmnGlobalTrack* glTrack = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(iGlobTrack);
434
435 if (glTrack->GetParamFirst()->GetQp() > 0.) {
436 protonCands.push_back(glTrack);
437 } else
438 pionCands.push_back(glTrack);
439 }
440
441 // Soft check (just momRes)
442 for (BmnGlobalTrack* cand : protonCands) {
443 if (!wasReconstructed(cand, 0.05))
444 continue;
445
446 hNhitsEmbReco[proton]->Fill(cand->GetNHits(), fMon->GetNHitsProton());
447
448 break;
449 }
450
451 // Soft check (just momRes)
452 for (BmnGlobalTrack* cand : pionCands) {
453 if (!wasReconstructed(cand, 0.05))
454 continue;
455
456 hNhitsEmbReco[pion]->Fill(cand->GetNHits(), fMon->GetNHitsPion());
457
458 break;
459 }
460
461 // Checking selected candidates once more ...
462 for (Int_t iCut = 0; iCut < nCutBins; iCut++) {
463 /* Protons */
464 // Strong check (momRes + nHits)
465 for (BmnGlobalTrack* cand : protonCands) {
466 if (!wasReconstructed(cand, cutsNHitsFrac[iCut], 0.05))
467 continue;
468
469 nRecProtons[iCut]++;
470 break;
471 }
472
473 /* Pions */
474 // Strong check (momRes + nHits)
475 for (BmnGlobalTrack* cand : pionCands) {
476 if (!wasReconstructed(cand, cutsNHitsFrac[iCut], 0.05))
477 continue;
478
479 nRecPions[iCut]++;
480 break;
481 }
482 }
483
484 break;
485 }
486 }
487
488 for (Int_t iCut = 0; iCut < nCutBins; iCut++) {
489 hEtaLambda[proton][iCut]->SetBinContent(hEtaLambda[proton][iCut]->FindBin(.5 * (etaMax + etaMin)),
490 1. * nRecProtons[iCut] / nEmbeddedProtons);
491 hEtaLambda[pion][iCut]->SetBinContent(hEtaLambda[pion][iCut]->FindBin(.5 * (etaMax + etaMin)),
492 1. * nRecPions[iCut] / nEmbeddedPions);
493 }
494
495 cout << "Events with embedded particles# " << counter << endl;
496 delete emb;
497 delete dst;
498 }
499}
500
501Bool_t BmnLambdaEmbeddingQa::wasReconstructed(BmnGlobalTrack* glTrack,
502 Double_t lostOrWrongConnectedHitsCut,
503 Double_t momResCut)
504{
505 Double_t pReco = 1. / glTrack->GetParamFirst()->GetQp();
506 Int_t nHitsReco = glTrack->GetNHits();
507
508 Int_t sign = (pReco > 0.) ? +1. : -1.;
509
510 Double_t pMc = (pReco > 0.) ? fMon->GetProtonP() : fMon->GetPionP();
511 Int_t nHitsMc = (pReco > 0.) ? fMon->GetNHitsProton() : fMon->GetNHitsPion();
512
513 Double_t momRes = TMath::Abs(pReco * sign - pMc) / pMc;
514
515 Double_t lostOrWrongConnectedHits = 1. * TMath::Abs(nHitsReco - nHitsMc) / nHitsMc;
516
517 if (momRes < momResCut && lostOrWrongConnectedHits <= lostOrWrongConnectedHitsCut)
518 return kTRUE;
519 else
520 return kFALSE;
521}
522
523Bool_t BmnLambdaEmbeddingQa::wasReconstructed(BmnGlobalTrack* glTrack, Double_t momResCut)
524{
525 Double_t pReco = 1. / glTrack->GetParamFirst()->GetQp();
526
527 Int_t sign = (pReco > 0.) ? +1. : -1.;
528
529 Double_t pMc = (pReco > 0.) ? fMon->GetProtonP() : fMon->GetPionP();
530
531 Double_t momRes = TMath::Abs(pReco * sign - pMc) / pMc;
532
533 if (momRes < momResCut)
534 return kTRUE;
535 else
536 return kFALSE;
537}
538
539void BmnLambdaEmbeddingQa::DoDigiAnalysis()
540{
541 histoMan->SetActive(0);
542
543 TH1F***** hGemStripInfo = histoMan->GetGemStripInfo();
544 TH1F***** hSiliconStripInfo = histoMan->GetSiliconStripInfo();
545
546 if (nInputs != -1)
547 for (Int_t iInput = 0; iInput < nInputs; iInput++) {
548
549 TChain* emb = new TChain("bmndata");
550 emb->Add(fEmb[iInput].Data());
551 emb->SetBranchAddress("EmbeddingMonitor.", &fMon);
552
553 TChain* digi = new TChain("bmndata");
554 digi->Add(fDigi[iInput].Data());
555 digi->SetBranchAddress("BmnEventHeader.", &hDigi);
556 digi->SetBranchAddress("GEM", &fGemDigits);
557 digi->SetBranchAddress("SILICON", &fSiliconDigits);
558
559 for (Int_t iEmb = 0; iEmb < emb->GetEntries(); iEmb++) {
560 emb->GetEntry(iEmb);
561
562 if (iEmb % 1000 == 0)
563 cout << "Event# " << iEmb << endl;
564
565 if (!fMon->IsEmbedded())
566 continue;
567
568 Int_t store = fMon->GetStoreVertexEvent()[0];
569 Int_t vertex = fMon->GetStoreVertexEvent()[1];
570 Int_t event = fMon->GetStoreVertexEvent()[2];
571
572 TChain* lambda = new TChain("bmndata");
573 lambda->Add(Form("%s/lambda%d_vertex%d.root", fPath[iInput].Data(), store, vertex));
574
575 lambda->SetBranchAddress("BmnGemStripDigit", &fGemDigits);
576 lambda->SetBranchAddress("BmnSiliconDigit", &fSiliconDigits);
577
578 lambda->GetEntry(event);
579
580 // Looping over GEM-digits ...
581 for (Int_t iDig = 0; iDig < fGemDigits->GetEntriesFast(); iDig++) {
582 BmnStripDigit* dig = (BmnStripDigit*)fGemDigits->UncheckedAt(iDig);
583
584 Int_t stat = dig->GetStation();
585 Int_t mod = dig->GetModule();
586 Int_t lay = dig->GetStripLayer();
587 Int_t strip = dig->GetStripNumber();
588
589 hGemStripInfo[0][stat][mod][lay]->Fill(strip);
590 }
591
592 // Looping over SILICON-digits ...
593 for (Int_t iDig = 0; iDig < fSiliconDigits->GetEntriesFast(); iDig++) {
594 BmnStripDigit* dig = (BmnStripDigit*)fSiliconDigits->UncheckedAt(iDig);
595
596 Int_t stat = dig->GetStation();
597 Int_t mod = dig->GetModule();
598 Int_t lay = dig->GetStripLayer();
599 Int_t strip = dig->GetStripNumber();
600
601 hSiliconStripInfo[0][stat][mod][lay]->Fill(strip);
602 }
603
604 delete lambda;
605
606 // Analyzing events with embedded digits ...
607 UInt_t idEmb = fMon->GetEventId();
608
609 for (UInt_t iEvent = idEmb - 1; iEvent < idEmb + 1; iEvent++) {
610 digi->GetEntry(iEvent);
611
612 for (Int_t iDigi = 0; iDigi < fGemDigits->GetEntriesFast(); iDigi++) {
613 BmnStripDigit* dig = (BmnStripDigit*)fGemDigits->UncheckedAt(iDigi);
614
615 Int_t stat = dig->GetStation();
616 Int_t mod = dig->GetModule();
617 Int_t lay = dig->GetStripLayer();
618 Int_t strip = dig->GetStripNumber();
619 Double_t signal = dig->GetStripSignal();
620
621 // Skipping digits not satisfied the signal window chosen ...
622 if (signal < fSignalWindow[GEM][0] || signal > fSignalWindow[GEM][1])
623 continue;
624
625 hGemStripInfo[1][stat][mod][lay]->Fill(strip);
626 }
627
628 for (Int_t iDigi = 0; iDigi < fSiliconDigits->GetEntriesFast(); iDigi++) {
629 BmnStripDigit* dig = (BmnStripDigit*)fSiliconDigits->UncheckedAt(iDigi);
630
631 Int_t stat = dig->GetStation();
632 Int_t mod = dig->GetModule();
633 Int_t lay = dig->GetStripLayer();
634 Int_t strip = dig->GetStripNumber();
635 Double_t signal = dig->GetStripSignal();
636
637 // Skipping digits not satisfied the signal window chosen ...
638 if (signal < fSignalWindow[SILICON][0] || signal > fSignalWindow[SILICON][1])
639 continue;
640
641 hSiliconStripInfo[1][stat][mod][lay]->Fill(strip);
642 }
643 }
644 }
645
646 delete emb;
647 delete digi;
648 }
649}
650
651void BmnLambdaEmbeddingQa::DoInnerTrackerRecoEfficiency()
652{
653 hProtonMomentaEmb = new TH1F("pMomentaEmb", "pMomentaEmb", 100, 1., 7.);
654 hPionMomentaEmb = new TH1F("pionMomentaEmb", "pionMomentaEmb", 100, 0., 2.);
655
656 hProtonMomentaReco = new TH1F("pMomentaReco", "pMomentaReco", 100, 1., 7.);
657 hPionMomentaReco = new TH1F("pionMomentaReco", "pionMomentaReco", 100, 0., 2.);
658
659 hProtonNhitsEmb = new TH1F("pNhitsEmb", "pNhitsEmb", 10, 0., 10.);
660 hPionNhitsEmb = new TH1F("pionNhitsEmb", "pionNhitsEmb", 10, 0., 10.);
661
662 hProtonNhitsReco = new TH1F("pNhitsReco", "pNhitsReco", 10, 0., 10.);
663 hPionNhitsReco = new TH1F("pionNhitsReco", "pionNhitsReco", 10, 0., 10.);
664
665 hProtonNhitsRecoAll = new TH1F("pNhitsRecoAll", "pNhitsRecoAll", 10, 0., 10.);
666 hPionNhitsRecoAll = new TH1F("pionNhitsRecoAll", "pionNhitsRecoAll", 10, 0., 10.);
667
668 const Int_t nHitsCut = 4;
669 const Double_t deltaPtoP = 0.05;
670
671 pEffProton = new TProfile("effProton vs. p",
672 Form("effProton vs. p, nHits > %d, Delta p / p < %G", nHitsCut, deltaPtoP), 100, 1., 7.);
673 pEffPion = new TProfile("effPion vs. p", Form("effPion vs. p, nHits > %d, Delta p / p < %G", nHitsCut, deltaPtoP),
674 100, 0., 2.);
675
676 if (nInputs != -1)
677 for (Int_t iInput = 0; iInput < nInputs; iInput++) {
678 TChain* emb = new TChain("bmndata");
679 emb->Add(fEmb[iInput].Data());
680
681 TChain* dst = new TChain("bmndata");
682 dst->Add(fDst[iInput].Data());
683
684 emb->SetBranchAddress("EmbeddingMonitor.", &fMon);
685 dst->SetBranchAddress("DstEventHeader.", &hDst);
686 dst->SetBranchAddress("BmnGlobalTrack", &fGlobTracks);
687 dst->SetBranchAddress("BmnGemStripHit", &fGemHits);
688 dst->SetBranchAddress("BmnSiliconHit", &fSiliconHits);
689
690 for (Int_t iEmb = 0; iEmb < emb->GetEntries(); iEmb++) {
691 emb->GetEntry(iEmb);
692
693 if (iEmb % 1000 == 0)
694 cout << "Event# " << iEmb << endl;
695
696 if (!fMon->IsEmbedded())
697 continue;
698
699 UInt_t idEmb = fMon->GetEventId();
700
701 for (UInt_t iDst = idEmb - 1; iDst < idEmb + 1; iDst++) {
702
703 dst->GetEntry(iDst);
704
705 if (hDst->GetEventId() != idEmb)
706 continue;
707
708 // Get event info ...
709 Int_t store = fMon->GetStoreVertexEvent()[0];
710 Int_t vertex = fMon->GetStoreVertexEvent()[1];
711 Int_t event = fMon->GetStoreVertexEvent()[2];
712
713 // Open lambda file ...
714 TChain* lambda = new TChain("bmndata");
715 lambda->Add(Form("%s/lambda%d_vertex%d.root", fPath[iInput].Data(), store, vertex));
716
717 lambda->SetBranchAddress("StsPoint", &fGemPoints);
718 lambda->SetBranchAddress("SiliconPoint", &fSiliconPoints);
719 lambda->SetBranchAddress("MCTrack", &fMcTracks);
720
721 lambda->GetEntry(event);
722
723 // [0] --> GEM, [1] --> SILICON
724 Int_t nRecoHitsProton[2] = {0, 0};
725 Int_t nRecoHitsPion[2] = {0, 0};
726
727 // Loop over GEM-points and reconstructed hits
728 for (Int_t iPoint = 0; iPoint < fGemPoints->GetEntriesFast(); iPoint++) {
729 CbmStsPoint* point = (CbmStsPoint*)fGemPoints->UncheckedAt(iPoint);
730
731 Bool_t isProton = kFALSE;
732 Bool_t isPion = kFALSE;
733
734 for (Int_t iHit = 0; iHit < fGemHits->GetEntriesFast(); iHit++) {
735 BmnHit* hit = (BmnHit*)fGemHits->UncheckedAt(iHit);
736
737 if (hit->GetStation() != point->GetStation())
738 continue;
739
740 Double_t dist = Abs(hit->GetX() - point->GetXIn());
741
742 if (dist < 0.5) {
743 // Probably this hit corresponds to the considering MC-point ...
744 // Here we should decide whether we consider proton or pion right now ...
745 CbmMCTrack* track = (CbmMCTrack*)fMcTracks->UncheckedAt(point->GetTrackID());
746 if (track->GetPdgCode() > 0)
747 isProton = kTRUE;
748 else
749 isPion = kTRUE;
750
751 break;
752 }
753 }
754
755 if (isProton)
756 nRecoHitsProton[0]++;
757 else if (isPion)
758 nRecoHitsPion[0]++;
759 }
760
761 // Loop over SILICON-points and reconstructed hits
762 for (Int_t iPoint = 0; iPoint < fSiliconPoints->GetEntriesFast(); iPoint++) {
763 BmnSiliconPoint* point = (BmnSiliconPoint*)fSiliconPoints->UncheckedAt(iPoint);
764
765 Bool_t isProton = kFALSE;
766 Bool_t isPion = kFALSE;
767
768 for (Int_t iHit = 0; iHit < fSiliconHits->GetEntriesFast(); iHit++) {
769 BmnHit* hit = (BmnHit*)fSiliconHits->UncheckedAt(iHit);
770
771 if (hit->GetStation() != point->GetStation())
772 continue;
773
774 Double_t dist = Abs(hit->GetX() - point->GetX());
775
776 if (dist < 0.3) {
777 // Probably this hit corresponds to the considering MC-point ...
778 // Here we should decide whether we consider proton or pion right now ...
779 CbmMCTrack* track = (CbmMCTrack*)fMcTracks->UncheckedAt(point->GetTrackID());
780 if (track->GetPdgCode() > 0)
781 isProton = kTRUE;
782 else
783 isPion = kTRUE;
784
785 break;
786 }
787 }
788
789 if (isProton)
790 nRecoHitsProton[1]++;
791 else if (isPion)
792 nRecoHitsPion[1]++;
793 }
794
795 hProtonMomentaEmb->Fill(fMon->GetProtonP());
796 hPionMomentaEmb->Fill(fMon->GetPionP());
797
798 hProtonNhitsEmb->Fill(fMon->GetNHitsProton());
799 hPionNhitsEmb->Fill(fMon->GetNHitsPion());
800
801 hProtonNhitsRecoAll->Fill(nRecoHitsProton[0] + nRecoHitsProton[1]);
802 hPionNhitsRecoAll->Fill(nRecoHitsPion[0] + nRecoHitsPion[1]);
803
804 // Loop over reconstructed glob. tracks ...
805 for (Int_t iGlobTrack = 0; iGlobTrack < fGlobTracks->GetEntriesFast(); iGlobTrack++) {
806 BmnGlobalTrack* glTrack = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(iGlobTrack);
807
808 Double_t p = glTrack->GetP();
809 Int_t nHits = glTrack->GetNHits();
810
811 if (p > 0.) {
812 hProtonMomentaReco->Fill(p);
813 hProtonNhitsReco->Fill(nHits);
814 } else {
815 hPionMomentaReco->Fill(-p);
816 hPionNhitsReco->Fill(nHits);
817 }
818
819 // Efficiency of track reconstruction ...
820 // Looking for a reconstr. track to be a proton ...
821 if (p > 0. && Abs((p - fMon->GetProtonP()) / p) < deltaPtoP && nHits > nHitsCut)
822 pEffProton->Fill(p, (1. * nHits) / fMon->GetNHitsProton());
823
824 if (p < 0. && Abs((-p - fMon->GetPionP()) / -p) < deltaPtoP && nHits > nHitsCut)
825 pEffPion->Fill(-p, (1. * nHits) / fMon->GetNHitsPion());
826 }
827 delete lambda;
828 }
829 }
830 delete emb;
831 delete dst;
832 }
833
834 // DrawHistos3();
835}
836
837void BmnLambdaEmbeddingQa::DoDrawEventsWithEmbeddedSignals()
838{
839 if (nInputs != 1) {
840 cout << endl;
841 cout << "*****************************************************" << endl;
842 cout << "*.list file contains more than one line! Exiting ...." << endl;
843 cout << "*****************************************************" << endl;
844 throw;
845 }
846
847 BmnGemStripStationSet* GEM = geoms->GetGemGeometry();
848 BmnSiliconStationSet* SILICON = geoms->GetSiliconGeometry();
849
850 TBox*** gemModBoxes = nullptr;
851 TBox**** gemLayBoxes = nullptr;
852 TPolyLine**** gemDeadZones = nullptr;
853 geoms->GetGemBorders(gemModBoxes, gemLayBoxes, gemDeadZones);
854
855 TBox*** siliconModBoxes = nullptr;
856 TBox**** siliconLayBoxes = nullptr;
857 TPolyLine**** siliconDeadZones = nullptr;
858 geoms->GetSiliconBorders(siliconModBoxes, siliconLayBoxes, siliconDeadZones);
859
860 TChain* emb = new TChain("bmndata");
861 emb->Add(fEmb[0].Data());
862
863 TChain* dst = new TChain("bmndata");
864 dst->Add(fDst[0].Data());
865
866 emb->SetBranchAddress("EmbeddingMonitor.", &fMon);
867 dst->SetBranchAddress("DstEventHeader.", &hDst);
868 dst->SetBranchAddress("BmnGemStripHit", &fGemHits);
869 dst->SetBranchAddress("BmnSiliconHit", &fSiliconHits);
870 dst->SetBranchAddress("BmnSiliconTrack", &fSiliconTracks);
871 dst->SetBranchAddress("BmnGlobalTrack", &fGlobTracks);
872
873 hXzMC = new TH2F("", "", 500, 0., 200., 500, -80., 20.);
874 hXzReco = new TH2F("", "", 500, 0., 200., 500, -80., 20.);
875 hYzMC = new TH2F("", "", 500, 0., 0., 500, 0., 0.);
876 hYzReco = new TH2F("", "", 500, 0., 0., 500, 0., 0.);
877
878 hGemXYProfiles = new TH2F*[GEM->GetNStations()];
879 for (Int_t iStat = 0; iStat < GEM->GetNStations(); iStat++) {
880 hGemXYProfiles[iStat] = new TH2F(Form("GEM: X vs. Y, stat# %d", iStat), Form("GEM: X vs. Y, stat# %d", iStat),
881 500, -90., +90., 500, -10., +40.);
882 hGemXYProfiles[iStat]->GetXaxis()->SetLabelSize(0.1);
883 hGemXYProfiles[iStat]->GetXaxis()->SetLabelSize(0.1);
884 hGemXYProfiles[iStat]->GetYaxis()->SetLabelSize(0.1);
885 }
886
887 hSiliconXYProfiles = new TH2F*[SILICON->GetNStations()];
888 for (Int_t iStat = 0; iStat < SILICON->GetNStations(); iStat++) {
889 Double_t xLow = (iStat == 0) ? -8. : (iStat == 1) ? -8. : -13.;
890 Double_t xUp = -xLow;
891
892 Double_t yLow = (iStat == 0) ? -12. : (iStat == 1) ? -13. : -20.;
893 Double_t yUp = (iStat == 0) ? +3. : (iStat == 1) ? +4. : +10.;
894
895 hSiliconXYProfiles[iStat] = new TH2F(Form("SILICON: X vs. Y, stat# %d", iStat),
896 Form("SILICON: X vs. Y, stat# %d", iStat), 500, xLow, xUp, 500, yLow, yUp);
897 hSiliconXYProfiles[iStat]->GetXaxis()->SetLabelSize(0.1);
898 hSiliconXYProfiles[iStat]->GetXaxis()->SetLabelSize(0.1);
899 hSiliconXYProfiles[iStat]->GetYaxis()->SetLabelSize(0.1);
900 }
901
902 if (drawFoundTracks) {
903 fRunAna = new FairRunAna();
904
905 // Read current geometry from the Geometry Database
906 dst->GetEntry(0);
907 GeoSetup gSetup;
908 string setupTag = gSetup.GetLastApprovedSetupTag("run7");
909 Int_t res_code = gSetup.LoadSetupToGeoManager(setupTag);
910 if (res_code != 0) {
911 cout << "Geometry '" << setupTag << "' can't be loaded from the Geometry Database" << endl;
912 exit(-1);
913 }
914
915 // Setting appropriate mag. field
916 fMagField = new BmnNewFieldMap("field_sp41v4_ascii_Extrap.root");
917 UniRun* pCurrentRun = UniRun::GetRun(7, hDst->GetRunId());
918 fMagField->SetScale(*pCurrentRun->GetFieldVoltage() / 55.7);
919 fRunAna->SetField(fMagField);
920 fMagField->Init();
921
922 hXzRecoFromTracks = new TH2F("", "", 500, 0., 200., 500, 0., 0.);
923 }
924
925 if (drawFoundTracks) {
926 for (Int_t iStat = 0; iStat < SILICON->GetNStations(); iStat++)
927 fSilGemZ[iStat] = SILICON->GetStation(iStat)->GetZPosition();
928
929 for (Int_t iStat = 0; iStat < GEM->GetNStations(); iStat++)
930 fSilGemZ[3 + iStat] = GEM->GetStation(iStat)->GetZPosition();
931 }
932
933 // Loop over events ...
934 for (Int_t iDst = 0; iDst < dst->GetEntries(); iDst++) {
935 dst->GetEntry(iDst);
936
937 hXzMC->Reset();
938 hXzReco->Reset();
939 hYzMC->Reset();
940 hYzReco->Reset();
941
942 if (drawFoundTracks)
943 hXzRecoFromTracks->Reset();
944
945 UInt_t idDst = hDst->GetEventId();
946
947 for (Int_t iEmb = 0; iEmb < emb->GetEntries(); iEmb++) {
948 emb->GetEntry(iEmb);
949
950 UInt_t idEmb = fMon->GetEventId();
951
952 if (idDst != idEmb)
953 continue;
954
955 UInt_t id = idDst;
956
957 if (!fMon->IsEmbedded())
958 continue;
959
960 hXzMC->SetTitle(Form(" XZ SI-GEM profile, id = %d", id));
961 hXzReco->SetTitle(Form("XZ SI-GEM profile, id = %d", id));
962 hYzMC->SetTitle(Form(" YZ SI-GEM profile, id = %d", id));
963 hYzReco->SetTitle(Form("YZ SI-GEM profile, id = %d", id));
964
965 if (drawFoundTracks)
966 DrawFoundTracks();
967
968 // Get MC-event info from file ...
969 Int_t store = fMon->GetStoreVertexEvent()[0];
970 Int_t vertex = fMon->GetStoreVertexEvent()[1];
971 Int_t event = fMon->GetStoreVertexEvent()[2];
972
973 TChain* lambda = new TChain("bmndata");
974 lambda->Add(Form("%s/lambda%d_vertex%d.root", fPath[0].Data(), store, vertex));
975
976 lambda->SetBranchAddress("StsPoint", &fGemPoints);
977 lambda->SetBranchAddress("SiliconPoint", &fSiliconPoints);
978
979 lambda->GetEntry(event);
980
981 for (Int_t iStat = 0; iStat < GEM->GetNStations(); iStat++) {
982 vector<pair<Double_t, Double_t>> tmp;
983 gems[0][iStat] = tmp;
984 }
985
986 for (Int_t iStat = 0; iStat < SILICON->GetNStations(); iStat++) {
987 vector<pair<Double_t, Double_t>> tmp;
988 silicons[0][iStat] = tmp;
989 }
990
991 // Get MC point info (GEM) ...
992 for (Int_t iPoint = 0; iPoint < fGemPoints->GetEntriesFast(); iPoint++) {
993 CbmStsPoint* point = (CbmStsPoint*)fGemPoints->UncheckedAt(iPoint);
994 Double_t x = point->GetXIn();
995 Double_t y = point->GetYIn();
996 Double_t z = point->GetZIn();
997 hXzMC->Fill(z, x);
998 hYzMC->Fill(z, y);
999
1000 gems[0].find(point->GetStation())->second.push_back(make_pair(x, y));
1001 }
1002
1003 // Get MC point info (SILICON) ...
1004 for (Int_t iPoint = 0; iPoint < fSiliconPoints->GetEntriesFast(); iPoint++) {
1005 BmnSiliconPoint* point = (BmnSiliconPoint*)fSiliconPoints->UncheckedAt(iPoint);
1006 Double_t x = point->GetX();
1007 Double_t y = point->GetY();
1008 Double_t z = point->GetZ();
1009 hXzMC->Fill(z, x);
1010 hYzMC->Fill(z, y);
1011
1012 silicons[0].find(point->GetStation())->second.push_back(make_pair(x, y));
1013 }
1014
1015 for (Int_t iStat = 0; iStat < GEM->GetNStations(); iStat++) {
1016 vector<pair<Double_t, Double_t>> tmp;
1017 gems[1][iStat] = tmp;
1018 }
1019
1020 for (Int_t iStat = 0; iStat < SILICON->GetNStations(); iStat++) {
1021 vector<pair<Double_t, Double_t>> tmp;
1022 silicons[1][iStat] = tmp;
1023 }
1024
1025 // Get GEM hit info ...
1026 if (fGemHits)
1027 for (Int_t iHit = 0; iHit < fGemHits->GetEntriesFast(); iHit++) {
1028 BmnGemStripHit* hit = (BmnGemStripHit*)fGemHits->UncheckedAt(iHit);
1029
1030 Double_t x = hit->GetX();
1031 Double_t y = hit->GetY();
1032 Double_t z = hit->GetZ();
1033 hXzReco->Fill(z, x);
1034 hYzReco->Fill(z, y);
1035
1036 gems[1].find(hit->GetStation())->second.push_back(make_pair(x, y));
1037 }
1038
1039 // Get SILICON hit info ...
1040 if (fSiliconHits)
1041 for (Int_t iHit = 0; iHit < fSiliconHits->GetEntriesFast(); iHit++) {
1042 BmnHit* hit = (BmnHit*)fSiliconHits->UncheckedAt(iHit);
1043 Double_t x = hit->GetX();
1044 Double_t y = hit->GetY();
1045 Double_t z = hit->GetZ();
1046 hXzReco->Fill(z, x);
1047 hYzReco->Fill(z, y);
1048
1049 silicons[1].find(hit->GetStation())->second.push_back(make_pair(x, y));
1050 }
1051
1052 delete lambda;
1053
1054 DrawHistos1();
1055 }
1056 }
1057}
1058
1059BmnLambdaEmbeddingQa::~BmnLambdaEmbeddingQa()
1060{
1061 if (fRunAna)
1062 delete fRunAna;
1063
1064 if (fMagField)
1065 delete fMagField;
1066
1067 if (geoms)
1068 delete geoms;
1069
1070 if (hXzMC)
1071 delete hXzMC;
1072 if (hXzReco)
1073 delete hXzReco;
1074 if (hYzMC)
1075 delete hYzMC;
1076 if (hYzReco)
1077 delete hYzReco;
1078 if (hGemXYProfiles)
1079 delete hGemXYProfiles;
1080 if (hSiliconXYProfiles)
1081 delete hSiliconXYProfiles;
1082
1083 if (hProtonMomentaEmb)
1084 delete hProtonMomentaEmb;
1085 if (hPionMomentaEmb)
1086 delete hPionMomentaEmb;
1087 if (hProtonMomentaReco)
1088 delete hProtonMomentaReco;
1089 if (hPionMomentaReco)
1090 delete hPionMomentaReco;
1091 if (hProtonNhitsEmb)
1092 delete hProtonNhitsEmb;
1093 if (hPionNhitsEmb)
1094 delete hPionNhitsEmb;
1095 if (hProtonNhitsReco)
1096 delete hProtonNhitsReco;
1097 if (hPionNhitsReco)
1098 delete hPionNhitsReco;
1099
1100 if (pEffProton)
1101 delete pEffProton;
1102 if (pEffPion)
1103 delete pEffPion;
1104
1105 if (histoMan) {
1106 histoMan->ProcessHistos();
1107 delete histoMan;
1108 }
1109}
1110
1111void BmnLambdaEmbeddingQa::DrawFoundTracks()
1112{
1113 for (Int_t iGlobTrack = 0; iGlobTrack < fGlobTracks->GetEntriesFast(); iGlobTrack++) {
1114 BmnGlobalTrack* glTrack = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(iGlobTrack);
1115
1116 BmnKalmanFilter* kalman = new BmnKalmanFilter();
1117
1118 FairTrackParam* first = glTrack->GetParamFirst();
1119 FairTrackParam* last_par = glTrack->GetParamLast();
1120
1121 Double_t zStart = first->GetZ();
1122 hXzRecoFromTracks->Fill(first->GetZ(), first->GetX());
1123
1124 map<Int_t, Int_t> statMod0;
1125 if (zStart < geoms->GetGemGeometry()->GetStation(0)->GetZPosition()) {
1126 BmnTrack* silTrack = (BmnTrack*)fSiliconTracks->UncheckedAt(glTrack->GetSilTrackIndex());
1127
1128 for (Int_t iHit = 0; iHit < silTrack->GetNHits(); iHit++) {
1129 BmnHit* silHit = (BmnHit*)fSiliconHits->UncheckedAt(silTrack->GetHitIndex(iHit));
1130 statMod0[silHit->GetStation()] = silHit->GetModule();
1131 }
1132 }
1133
1134 Double_t zFinish = last_par->GetZ();
1135
1136 Int_t pdg = (first->GetQp() > 0.) ? 2212 : -211;
1137
1138 // Try to locate a detector where the track has first point
1139 // in order to go forward by propagator ...
1140 Int_t idx = -1;
1141 if (statMod0.size() == 0)
1142 for (auto it : fSilGemZ) {
1143 if (TMath::Abs(it.second - zStart) < 1.) {
1144 idx = it.first;
1145 break;
1146 }
1147 }
1148 else {
1149 for (auto it : statMod0) {
1150 idx = it.first;
1151 break;
1152 }
1153 }
1154
1155 for (Int_t iStep = idx + 1; iStep < 9; iStep++) {
1156 Double_t zStat = -1000.;
1157 if (iStep < 3) {
1158 for (auto it : statMod0) {
1159 if (it.first != iStep)
1160 continue;
1161
1162 zStat = geoms->GetSiliconGeometry()
1163 ->GetStation(it.first)
1164 ->GetModule(it.second)
1165 ->GetZPositionRegistered();
1166 }
1167 } else
1168 zStat = fSilGemZ.find(iStep)->second;
1169
1170 if (zStat > zFinish)
1171 break;
1172
1173 kalman->TGeoTrackPropagate(first, zStat, pdg, nullptr, nullptr, kTRUE);
1174 hXzRecoFromTracks->Fill(first->GetZ(), first->GetX());
1175 }
1176
1177 delete kalman;
1178 }
1179}
1180
1181void BmnLambdaEmbeddingQa::DrawHistos1()
1182{
1183 gStyle->SetOptStat(0);
1184 TBox*** gemModBoxes = nullptr;
1185 TBox**** gemLayBoxes = nullptr;
1186 TPolyLine**** gemDeadZones = nullptr;
1187 geoms->GetGemBorders(gemModBoxes, gemLayBoxes, gemDeadZones);
1188
1189 TBox*** siliconModBoxes = nullptr;
1190 TBox**** siliconLayBoxes = nullptr;
1191 TPolyLine**** siliconDeadZones = nullptr;
1192 geoms->GetSiliconBorders(siliconModBoxes, siliconLayBoxes, siliconDeadZones);
1193
1194 TCanvas* c = new TCanvas("c1", "c1", 1200, 1200);
1195 gStyle->SetTitleFontSize(0.1);
1196 gStyle->SetTitleY(1.015);
1197 c->Divide(1, 2);
1198 TVirtualPad* pad = c->cd(1);
1199 pad->Divide(2, 1);
1200 pad->cd(1);
1201 hXzMC->Draw("P*");
1202 hXzMC->SetMarkerColor(kRed);
1203 hXzMC->SetMarkerSize(1.5);
1204 hXzMC->SetMarkerStyle(20);
1205 if (hXzReco->GetEntries() > 0.) {
1206 hXzReco->Draw("sameP*");
1207 hXzReco->SetMarkerColor(kBlue);
1208 hXzReco->SetMarkerSize(1.);
1209 hXzReco->SetMarkerStyle(19);
1210 }
1211
1212 if (hXzRecoFromTracks && hXzRecoFromTracks->GetEntries() > 0.) {
1213 hXzRecoFromTracks->Draw("sameP*");
1214 hXzRecoFromTracks->SetMarkerColor(kBlack);
1215 hXzRecoFromTracks->SetMarkerSize(2.5);
1216 hXzRecoFromTracks->SetMarkerStyle(kOpenStar);
1217 }
1218
1219 pad->cd(2);
1220 hYzMC->Draw("P*");
1221 hYzMC->SetMarkerColor(kRed);
1222 hYzMC->SetMarkerSize(1.5);
1223 hYzMC->SetMarkerStyle(20);
1224 if (hYzReco->GetEntries() > 0.) {
1225 hYzReco->Draw("sameP*");
1226 hYzReco->SetMarkerColor(kBlue);
1227 hYzReco->SetMarkerSize(1.);
1228 hYzReco->SetMarkerStyle(19);
1229 }
1230
1231 pad = c->cd(2);
1232 pad->Divide(3, 3);
1233
1234 for (Int_t iStat = 0; iStat < geoms->GetGemGeometry()->GetNStations(); iStat++) {
1235 pad->cd(iStat + 1);
1236
1237 hGemXYProfiles[iStat]->Draw();
1238
1239 Int_t nHitsPerStat = gems[0].find(iStat)->second.size();
1240 for (Int_t iHit = 0; iHit < nHitsPerStat; iHit++) {
1241 TMarker* marker =
1242 new TMarker(gems[0].find(iStat)->second[iHit].first, gems[0].find(iStat)->second[iHit].second, 19);
1243 marker->SetMarkerColor(kRed);
1244 marker->SetMarkerSize(1.5);
1245 marker->Draw();
1246 }
1247
1248 nHitsPerStat = gems[1].find(iStat)->second.size();
1249 for (Int_t iHit = 0; iHit < nHitsPerStat; iHit++) {
1250 TMarker* marker =
1251 new TMarker(gems[1].find(iStat)->second[iHit].first, gems[1].find(iStat)->second[iHit].second, 19);
1252 marker->SetMarkerColor(kBlue);
1253 marker->SetMarkerSize(.8);
1254 marker->Draw();
1255 }
1256
1257 for (Int_t iMod = 0; iMod < geoms->GetGemGeometry()->GetStation(iStat)->GetNModules(); iMod++) {
1258 gemModBoxes[iStat][iMod]->Draw("l");
1259
1260 for (Int_t iLayer = 0;
1261 iLayer < geoms->GetGemGeometry()->GetStation(iStat)->GetModule(iMod)->GetNStripLayers(); iLayer++)
1262 {
1263 gemLayBoxes[iStat][iMod][iLayer]->Draw("l");
1264 gemDeadZones[iStat][iMod][iLayer]->Draw("l");
1265 }
1266 }
1267 }
1268
1269 for (Int_t iStat = 0; iStat < geoms->GetSiliconGeometry()->GetNStations(); iStat++) {
1270 pad->cd(iStat + geoms->GetGemGeometry()->GetNStations() + 1);
1271
1272 hSiliconXYProfiles[iStat]->Draw();
1273
1274 Int_t nHitsPerStat = silicons[0].find(iStat)->second.size();
1275 for (Int_t iHit = 0; iHit < nHitsPerStat; iHit++) {
1276 TMarker* marker = new TMarker(silicons[0].find(iStat)->second[iHit].first,
1277 silicons[0].find(iStat)->second[iHit].second, 19);
1278 marker->SetMarkerColor(kRed);
1279 marker->SetMarkerSize(1.5);
1280 marker->Draw();
1281 }
1282
1283 nHitsPerStat = silicons[1].find(iStat)->second.size();
1284 for (Int_t iHit = 0; iHit < nHitsPerStat; iHit++) {
1285 TMarker* marker = new TMarker(silicons[1].find(iStat)->second[iHit].first,
1286 silicons[1].find(iStat)->second[iHit].second, 19);
1287 marker->SetMarkerColor(kBlue);
1288 marker->SetMarkerSize(.8);
1289 marker->Draw();
1290 }
1291
1292 for (Int_t iMod = 0; iMod < geoms->GetSiliconGeometry()->GetStation(iStat)->GetNModules(); iMod++) {
1293 siliconModBoxes[iStat][iMod]->Draw("l");
1294
1295 for (Int_t iLayer = 0;
1296 iLayer < geoms->GetSiliconGeometry()->GetStation(iStat)->GetModule(iMod)->GetNStripLayers(); iLayer++)
1297 {
1298 siliconLayBoxes[iStat][iMod][iLayer]->Draw("l");
1299
1300 if (iStat == 2 && iMod == 7 && (iLayer == 0 || iLayer == 2))
1301 siliconDeadZones[iStat][iMod][iLayer]->Draw("l");
1302 }
1303 }
1304 }
1305 c->SaveAs("tmp.pdf");
1306
1307 delete c;
1308
1309 getchar();
1310}
1311
1312void BmnLambdaEmbeddingQa::DrawHistos6(TClonesArray* map1, TClonesArray* map2)
1313{
1314 // map1 - main part
1315 // map2 - common ADC
1316
1317 const Int_t nZones = 2; // main zone and common ADC
1318
1319 hStripChannel = new TH2F***[nZones];
1320 TBox**** boundBoxesStrips = new TBox***[nZones];
1321 TBox**** boundBoxesChannels = new TBox***[nZones];
1322
1323 // Creating necessary histograms and help boxes on histograms...
1324 vector<Int_t> naturalGemOrder{11, 10, 5, 6, 8, 9};
1325
1326 for (Int_t iZone = 0; iZone < nZones; iZone++) {
1327 const Int_t nStats = geoms->GetGemGeometry()->GetNStations();
1328 boundBoxesStrips[iZone] = new TBox**[nStats];
1329 boundBoxesChannels[iZone] = new TBox**[nStats];
1330 hStripChannel[iZone] = new TH2F**[nStats];
1331
1332 for (Int_t iStat = 0; iStat < nStats; iStat++) {
1333
1334 const Int_t nMaps = 2; // Right, Left
1335 hStripChannel[iZone][iStat] = new TH2F*[nMaps];
1336 boundBoxesStrips[iZone][iStat] = new TBox*[nMaps];
1337 boundBoxesChannels[iZone][iStat] = new TBox*[nMaps];
1338
1339 for (Int_t iMap = 0; iMap < nMaps; iMap++) {
1340
1341 TString mapping = (iMap == 0) ? "Left Part" : "Right Part";
1342
1343 hStripChannel[iZone][iStat][iMap] =
1344 new TH2F(Form("Zone# %d Stat# %d Map# %d", iZone, iStat, iMap),
1345 Form("Stat# %d (GEM %d), %s", iStat, naturalGemOrder[iStat], mapping.Data()), 2060, 0.,
1346 2060., 1200, 0., 1200.);
1347 hStripChannel[iZone][iStat][iMap]->GetXaxis()->SetLabelSize(0.1);
1348 hStripChannel[iZone][iStat][iMap]->GetYaxis()->SetLabelSize(0.1);
1349 hStripChannel[iZone][iStat][iMap]->GetXaxis()->SetTitle("Channel#");
1350 hStripChannel[iZone][iStat][iMap]->GetXaxis()->SetTitleOffset(-.28);
1351 hStripChannel[iZone][iStat][iMap]->GetXaxis()->SetTitleSize(0.15);
1352 hStripChannel[iZone][iStat][iMap]->GetXaxis()->CenterTitle();
1353
1354 hStripChannel[iZone][iStat][iMap]->GetYaxis()->SetTitle("Strip#");
1355 hStripChannel[iZone][iStat][iMap]->GetYaxis()->SetTitleOffset(-.15);
1356 hStripChannel[iZone][iStat][iMap]->GetYaxis()->SetTitleSize(0.15);
1357 hStripChannel[iZone][iStat][iMap]->GetYaxis()->CenterTitle();
1358
1359 Int_t low = (iMap == 0) ? 985 : 1024;
1360 Int_t up = (iMap == 0) ? 1080 : 1129;
1361
1362 boundBoxesStrips[iZone][iStat][iMap] = new TBox(0., low, 2048., up);
1363 boundBoxesStrips[iZone][iStat][iMap]->SetFillColorAlpha(TColor::GetColor("#cccccc"), 0.3);
1364 boundBoxesStrips[iZone][iStat][iMap]->SetFillStyle(3001);
1365
1366 boundBoxesChannels[iZone][iStat][iMap] = new TBox(0., 0., 0., 1200.);
1367 boundBoxesChannels[iZone][iStat][iMap]->SetFillColorAlpha(TColor::GetColor("#cccccc"), 0.3);
1368 boundBoxesChannels[iZone][iStat][iMap]->SetFillStyle(3001);
1369 }
1370 }
1371 }
1372
1373 // Filling histograms ...
1374
1375 // Main part of big zone [0] ...
1376 for (Int_t iInfo = 0; iInfo < map1->GetEntriesFast(); iInfo++) {
1377 MappingInfo* info = (MappingInfo*)map1->UncheckedAt(iInfo);
1378
1379 Int_t stat = info->station;
1380 TString mapping = info->mapFile;
1381
1382 Int_t ch = info->channel;
1383 Int_t str = info->strip;
1384
1385 Int_t mapIdx = (mapping.Contains("Left")) ? 0 : 1;
1386
1387 hStripChannel[0][stat][mapIdx]->Fill(ch, str);
1388 }
1389
1390 // Common ADC [1] ...
1391 for (Int_t iInfo = 0; iInfo < map2->GetEntriesFast(); iInfo++) {
1392 MappingInfo* info = (MappingInfo*)map2->UncheckedAt(iInfo);
1393
1394 Int_t stat = info->station;
1395 TString mapping = info->mapFile;
1396
1397 Int_t mapIdx = (mapping.Contains("Left")) ? 0 : 1;
1398
1399 map<Int_t, Int_t> stripGlobChan = info->stripChan;
1400
1401 if ((stat == 0 || stat == 1) && mapIdx == 1)
1402 continue;
1403
1404 for (auto it : stripGlobChan)
1405 hStripChannel[1][stat][mapIdx]->Fill(it.second, it.first);
1406
1407 boundBoxesChannels[1][stat][mapIdx]->SetX1(info->channels.first);
1408 boundBoxesChannels[1][stat][mapIdx]->SetX2(info->channels.second);
1409 }
1410
1411 TCanvas* c = new TCanvas("c1", "c1", 1200, 800);
1412 c->Divide(2, 3);
1413
1414 gStyle->SetOptStat(0);
1415 gStyle->SetTitleFontSize(0.1);
1416 gStyle->SetTitleY(1.0);
1417
1418 const Int_t nPads = 6;
1419 Int_t statsC[nPads] = {0, 0, 3, 3, 5, 5}; // LeftToRight
1420 Int_t statsD[nPads] = {1, 1, 2, 2, 4, 4}; // RightToLeft
1421
1422 for (Int_t iPad = 1; iPad < nPads + 1; iPad++) {
1423 TVirtualPad* vPadC = c->cd(iPad);
1424 vPadC->Divide(1, 2, 0.01, 0.01);
1425
1426 vPadC->cd(1);
1427 Int_t mapping = 0;
1428
1429 Int_t stat = (iPad % 2 == 1) ? statsC[iPad - 1] : statsD[iPad - 1];
1430
1431 hStripChannel[0][stat][mapping]->Draw("");
1432 hStripChannel[1][stat][mapping]->Draw("same");
1433 boundBoxesStrips[1][stat][mapping]->Draw("l");
1434 boundBoxesChannels[1][stat][mapping]->Draw("l");
1435
1436 vPadC->cd(2);
1437 mapping = 1;
1438
1439 hStripChannel[0][stat][mapping]->Draw("");
1440 hStripChannel[1][stat][mapping]->Draw("same");
1441 boundBoxesStrips[1][stat][mapping]->Draw("l");
1442 boundBoxesChannels[1][stat][mapping]->Draw("l");
1443 }
1444
1445 c->SaveAs("strChan.pdf");
1446
1447 delete c;
1448}
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
i
int i
Definition P4_F32vec4.h:22
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
Bool_t
bool Bool_t
Int_t
int Int_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
kRed
kRed
kBlack
kBlack
kBlue
kBlue
kOpenStar
kOpenStar
Form
char * Form(const char *fmt,...)
gStyle
R__EXTERN TStyle * gStyle
BmnGemStripLayer::IsPointInsideStripLayer
Bool_t IsPointInsideStripLayer(Double_t x, Double_t y)
Definition BmnGemStripLayer.cxx:254
BmnGemStripModule::GetElectronDriftDirection
ElectronDriftDirectionInModule GetElectronDriftDirection()
Definition BmnGemStripModule.h:61
BmnGlobalTrack
BmnGlobalTrack.h.
Definition BmnGlobalTrack.h:14
BmnHit::GetModule
Int_t GetModule()
Definition BmnHit.h:36
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:35
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:438
BmnLambdaEmbeddingQa::wasReconstructed
Bool_t wasReconstructed(BmnGlobalTrack *, Double_t fracLostOrWrongConnectedHits, Double_t momResCut)
Definition BmnLambdaEmbeddingQa.cxx:501
BmnLambdaEmbeddingQa::histoMan
BmnLambdaEmbeddingDrawHistos * histoMan
Definition BmnLambdaEmbeddingQa.h:178
BmnLambdaEmbeddingQa::geoms
BmnInnerTrackerGeometryDraw * geoms
Definition BmnLambdaEmbeddingQa.h:96
BmnLambdaEmbeddingQa::fMon
BmnLambdaEmbeddingMonitor * fMon
Definition BmnLambdaEmbeddingQa.h:106
BmnLambdaEmbeddingQa::DrawHistos6
void DrawHistos6(TClonesArray *, TClonesArray *)
Definition BmnLambdaEmbeddingQa.cxx:1312
BmnLambdaEmbeddingQa::fSilGemZ
map< Int_t, Double_t > fSilGemZ
Definition BmnLambdaEmbeddingQa.h:120
BmnLambdaEmbeddingQa::gems
map< Int_t, vector< pair< Double_t, Double_t > > > gems[2]
Definition BmnLambdaEmbeddingQa.h:122
BmnLambdaEmbeddingQa::silicons
map< Int_t, vector< pair< Double_t, Double_t > > > silicons[2]
Definition BmnLambdaEmbeddingQa.h:123
BmnNewFieldMap
Field Map does not have any symmetries w.r.t coordinate axises.
Definition BmnNewFieldMap.h:16
BmnStripDigit
Class BmnStripDigit Base class for strip detector digits.
Definition BmnStripDigit.h:14
BmnStripDigit::GetStripNumber
Int_t GetStripNumber()
Definition BmnStripDigit.h:29
BmnStripDigit::GetStripLayer
Int_t GetStripLayer()
Definition BmnStripDigit.h:28
BmnStripDigit::GetStation
Int_t GetStation()
Definition BmnStripDigit.h:26
BmnStripDigit::GetStripSignal
Double_t GetStripSignal()
Definition BmnStripDigit.h:30
BmnStripDigit::GetModule
Int_t GetModule()
Definition BmnStripDigit.h:27
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:63
BmnTrack::GetNHits
Int_t GetNHits() const
Accessors.
Definition BmnTrack.h:49
BmnTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition BmnTrack.h:51
BmnTrack::GetP
Double_t GetP()
Definition BmnTrack.h:67
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:65
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Accessors.
Definition CbmMCTrack.h:56
CbmStsPoint::GetXIn
Double_t GetXIn() const
Accessors.
Definition CbmStsPoint.h:64
CbmStsPoint::GetStation
Int_t GetStation() const
Definition CbmStsPoint.h:76
CbmStsPoint::GetZIn
Double_t GetZIn() const
Definition CbmStsPoint.h:66
CbmStsPoint::GetYIn
Double_t GetYIn() const
Definition CbmStsPoint.h:65
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
MappingInfo::mapFile
TString mapFile
Definition BmnLambdaMisc.h:72
MappingInfo::channels
pair< Int_t, Int_t > channels
Definition BmnLambdaMisc.h:73
MappingInfo::stripChan
map< Int_t, Int_t > stripChan
Definition BmnLambdaMisc.h:80
TAttFill::SetFillColorAlpha
virtual void SetFillColorAlpha(Color_t fcolor, Float_t falpha)
TAttFill::SetFillStyle
virtual void SetFillStyle(Style_t fstyle)
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
TBox::Draw
void Draw(Option_t *option="") override
TBox::SetX1
virtual void SetX1(Double_t x1)
TBox::SetX2
virtual void SetX2(Double_t x2)
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
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::GetEntry
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
TColor::GetColor
static Int_t GetColor(const char *hexcolor)
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
TH2::Fill
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
TMarker::Draw
void Draw(Option_t *option="") override
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::GetEntries
Int_t GetEntries() const override
TObjArray::At
TObject * At(Int_t idx) const override
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
TObject::TClonesArray
friend friend class TClonesArray
TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
TPad::SaveAs
void SaveAs(const char *filename="", Option_t *option="") const override
TPolyLine::Draw
void Draw(Option_t *option="") override
TProfile2D::Fill
Int_t Fill(const char *namex, const char *namey, Double_t z) override
TString::Data
const char * Data() const
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
TVirtualPad::cd
virtual TVirtualPad * cd(Int_t subpadnumber=0)=0
TVirtualPad::Divide
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)=0
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::GetFieldVoltage
double * GetFieldVoltage()
get field voltage of the current run
Definition UniRun.h:121
ForwardZAxisEDrift
@ ForwardZAxisEDrift
Definition BmnGemStripModule.h:20
TMath::Abs
Double_t Abs(Double_t d)