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