CbmRoot
qa/CbmRichRecoQa.cxx
Go to the documentation of this file.
1 
2 #include "CbmRichRecoQa.h"
3 
4 #include "TCanvas.h"
5 #include "TClonesArray.h"
6 #include "TEllipse.h"
7 #include "TF1.h"
8 #include "TH1.h"
9 #include "TH1D.h"
10 #include "TMCProcess.h"
11 #include "TMarker.h"
12 #include "TStyle.h"
13 #include <TFile.h>
14 
15 #include "CbmDrawHist.h"
16 #include "CbmGlobalTrack.h"
17 #include "CbmMCTrack.h"
18 #include "CbmMatchRecoToMC.h"
19 #include "CbmRichDraw.h"
20 #include "CbmRichGeoManager.h"
21 #include "CbmRichHit.h"
22 #include "CbmRichPoint.h"
23 #include "CbmRichRing.h"
24 #include "CbmRichUtil.h"
25 #include "CbmTrackMatchNew.h"
26 #include "FairMCPoint.h"
27 #include "FairTrackParam.h"
28 #include "elid/CbmLitGlobalElectronId.h"
29 
30 #include "CbmDigiManager.h"
31 #include "CbmHistManager.h"
32 #include "CbmUtils.h"
33 
34 #include <boost/assign/list_of.hpp>
35 #include <iostream>
36 #include <sstream>
37 #include <string>
38 
39 using namespace std;
40 using boost::assign::list_of;
41 
42 CbmRichRecoQa::CbmRichRecoQa()
43  : FairTask("CbmRichRecoQa")
44  , fHM(nullptr)
45  , fEventNum(0)
46  , fOutputDir("")
47  , fMCTracks(nullptr)
48  , fRichPoints(nullptr)
49  , fRichHits(nullptr)
50  , fRichRings(nullptr)
51  , fRichRingMatches(nullptr)
52  , fGlobalTracks(nullptr)
53  , fStsTracks(nullptr)
54  , fStsTrackMatches(nullptr)
55  , fRichProjections(nullptr)
56  , fDigiMan(nullptr)
57  , fNofHitsInRingMap()
58  , fCanvas() {}
59 
60 
61 InitStatus CbmRichRecoQa::Init() {
62  cout << "CbmRichRecoQa::Init" << endl;
63  FairRootManager* ioman = FairRootManager::Instance();
64  if (nullptr == ioman) {
65  Fatal("CbmRichRecoQa::Init", "RootManager not instantised!");
66  }
67 
68  fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
69  if (nullptr == fMCTracks) { Fatal("CbmRichRecoQa::Init", "No MC Tracks!"); }
70 
71  fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
72  if (nullptr == fRichPoints) {
73  Fatal("CbmRichRecoQa::Init", "No Rich Points!");
74  }
75 
76  fDigiMan = CbmDigiManager::Instance();
77  fDigiMan->Init();
78 
79  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
80  if (nullptr == fRichHits) { Fatal("CbmRichRecoQa::Init", "No RichHits!"); }
81 
82  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
83  if (nullptr == fRichRings) { Fatal("CbmRichRecoQa::Init", "No RichRings!"); }
84 
85  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
86  if (nullptr == fRichRingMatches) {
87  Fatal("CbmRichRecoQa::Init", "No RichRingMatch array!");
88  }
89 
90  fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
91  if (nullptr == fGlobalTracks) {
92  Fatal("CbmRichRecoQa::Init", "No GlobalTrack array!");
93  }
94 
95  fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
96  if (nullptr == fStsTracks) {
97  Fatal("CbmRichRecoQa::Init", ": No StsTrack array!");
98  }
99 
100  fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
101  if (nullptr == fStsTrackMatches) {
102  Fatal("CbmRichRecoQa::Init", ": No StsTrackMatch array!");
103  }
104 
105  fRichProjections = (TClonesArray*) ioman->GetObject("RichProjection");
106  if (nullptr == fRichProjections) {
107  Fatal("CbmRichUrqmdTest::Init", "No fRichProjections array!");
108  }
109 
110  InitHistograms();
111 
112  // CbmLitGlobalElectronId::GetInstance();
113  cout << "CbmRichRecoQa::Init finished" << endl;
114 
115  return kSUCCESS;
116 }
117 
118 void CbmRichRecoQa::InitHistograms() {
119  fHM = new CbmHistManager();
120 
121  // 2D histogramm XY
122  double xMin = -120.;
123  double xMax = 120.;
124  int nBinsX = 30;
125  double yMin = -208;
126  double yMax = 208.;
127  int nBinsY = 52;
128 
129  // 1D histogram X or Y
130  int nBinsX1 = 60;
131  int xMin1 = -120;
132  int xMax1 = 120;
133  int nBinsY1 = 25;
134  int yMin1 = 100;
135  int yMax1 = 200;
136 
137  for (Int_t i = 0; i < 4; i++) {
138  string s;
139  if (i == 0) s = "Primel";
140  if (i == 1) s = "Pi";
141  if (i == 2) s = "PrimelPlus";
142  if (i == 3) s = "PrimelMinus";
143 
144  fHM->Create2<TH2D>("fhRingTrackDistVsMomTruematch" + s,
145  "fhRingTrackDistVsMomTruematch" + s
146  + ";P [GeV/c];Ring-track distance [cm];Yield (a.u.)",
147  20,
148  0.,
149  10.,
150  100,
151  0.,
152  5.);
153  fHM->Create2<TH2D>("fhRingTrackDistVsMomWrongmatch" + s,
154  "fhRingTrackDistVsMomWrongmatch" + s
155  + ";P [GeV/c];Ring-track distance [cm];Yield (a.u.)",
156  20,
157  0.,
158  10.,
159  100,
160  0.,
161  5.);
162 
163  fHM->Create2<TH2D>(
164  "fhRingTrackDistVsNofHitsTruematch" + s,
165  "fhRingTrackDistVsNofHitsTruematch" + s
166  + ";Nof hits in found ring;Ring-track distance [cm];Yield (a.u.)",
167  40,
168  -.5,
169  39.5,
170  100,
171  0.,
172  5.);
173  fHM->Create3<TH3D>("fhRingTrackDistVsXYTruematch" + s,
174  "fhRingTrackDistVsXYTruematch" + s
175  + ";X [cm];Y [cm];Ring-track distance [cm]",
176  nBinsX,
177  xMin,
178  xMax,
179  nBinsY,
180  yMin,
181  yMax,
182  100,
183  0.,
184  5.);
185  fHM->Create2<TH2D>("fhRingTrackDistVsXTruematch" + s,
186  "fhRingTrackDistVsXTruematch" + s
187  + ";X [cm];Ring-track distance [cm]",
188  nBinsX1,
189  xMin1,
190  xMax1,
191  100,
192  0.,
193  5.);
194  fHM->Create2<TH2D>("fhRingTrackDistVsYTruematch" + s,
195  "fhRingTrackDistVsYTruematch" + s
196  + ";Abs(Y) [cm];Ring-track distance [cm]",
197  nBinsY1,
198  yMin1,
199  yMax1,
200  100,
201  0.,
202  5.);
203  }
204 
205  // after electron identification
206  fHM->Create2<TH2D>("fhRingTrackDistVsMomTruematchElId",
207  "fhRingTrackDistVsMomTruematchElId;P [GeV/c];Ring-track "
208  "distance [cm];Yield (a.u.)",
209  20,
210  0.,
211  10.,
212  100,
213  0.,
214  5.);
215 
216  fHM->Create1<TH1D>("fhMismatchSource",
217  "fhMismatchSource;Global track category;% from MC",
218  13,
219  -0.5,
220  12.5);
221 
222  fHM->Create1<TH1D>("fhMismatchSourceMomMc",
223  "fhMismatchSourceMomMc;Momentum [GeV/c];Yield",
224  40,
225  0.,
226  10.);
227  fHM->Create1<TH1D>("fhMismatchSourceMomSts",
228  "fhMismatchSourceMomSts;Momentum [GeV/c];Yield",
229  40,
230  0.,
231  10.);
232  fHM->Create1<TH1D>("fhMismatchSourceMomStsAccRich",
233  "fhMismatchSourceMomStsAccRich;Momentum [GeV/c];Yield",
234  40,
235  0.,
236  10.);
237  fHM->Create1<TH1D>("fhMismatchSourceMomStsRich",
238  "fhMismatchSourceMomStsRich;Momentum [GeV/c];Yield",
239  40,
240  0.,
241  10.);
242  fHM->Create1<TH1D>("fhMismatchSourceMomStsRichTrue",
243  "fhMismatchSourceMomStsRichTrue;Momentum [GeV/c];Yield",
244  40,
245  0.,
246  10.);
247  fHM->Create1<TH1D>("fhMismatchSourceMomStsNoRich",
248  "fhMismatchSourceMomStsNoRich;Momentum [GeV/c];Yield",
249  40,
250  0.,
251  10.);
252  fHM->Create1<TH1D>("fhMismatchSourceMomStsNoRichRF",
253  "fhMismatchSourceMomStsNoRichRF;Momentum [GeV/c];Yield",
254  40,
255  0.,
256  10.);
257  fHM->Create1<TH1D>("fhMismatchSourceMomStsNoRichRM",
258  "fhMismatchSourceMomStsNoRichRM;Momentum [GeV/c];Yield",
259  40,
260  0.,
261  10.);
262  fHM->Create1<TH1D>("fhMismatchSourceMomStsNoRichNoRF",
263  "fhMismatchSourceMomStsNoRichNoRF;Momentum [GeV/c];Yield",
264  40,
265  0.,
266  10.);
267  fHM->Create1<TH1D>(
268  "fhMismatchSourceMomStsNoRichNoProj",
269  "fhMismatchSourceMomStsNoRichNoProj;Momentum [GeV/c];Yield",
270  40,
271  0.,
272  10.);
273  fHM->Create1<TH1D>("fhMismatchSourceMomStsRichWrong",
274  "fhMismatchSourceMomStsRichWrong;Momentum [GeV/c];Yield",
275  40,
276  0.,
277  10.);
278  fHM->Create1<TH1D>("fhMismatchSourceMomStsRichWrongRF",
279  "fhMismatchSourceMomStsRichWrongRF;Momentum [GeV/c];Yield",
280  40,
281  0.,
282  10.);
283  fHM->Create1<TH1D>("fhMismatchSourceMomStsRichWrongRM",
284  "fhMismatchSourceMomStsRichWrongRM;Momentum [GeV/c];Yield",
285  40,
286  0.,
287  10.);
288 }
289 
290 void CbmRichRecoQa::Exec(Option_t* /*option*/) {
291  fEventNum++;
292  cout << "CbmRichRecoQa, event No. " << fEventNum << endl;
293 
294  FillRichRingNofHits();
295  FillRingTrackDistance();
296  RingTrackMismatchSource();
297 }
298 
299 void CbmRichRecoQa::FillRichRingNofHits() {
300  fNofHitsInRingMap.clear();
301  Int_t nofRichHits = fRichHits->GetEntriesFast();
302  for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
303  CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
304  if (nullptr == hit) continue;
305 
306  vector<Int_t> motherIds = CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(
307  fDigiMan, hit, fRichPoints, fMCTracks);
308  for (UInt_t i = 0; i < motherIds.size(); i++) {
309  fNofHitsInRingMap[motherIds[i]]++;
310  }
311  }
312 }
313 
314 void CbmRichRecoQa::RingTrackMismatchSource() {
315  Int_t nofMcTracks = fMCTracks->GetEntriesFast();
316  for (Int_t iTrack = 0; iTrack < nofMcTracks; iTrack++) {
317  const CbmMCTrack* mcTrack =
318  static_cast<const CbmMCTrack*>(fMCTracks->At(iTrack));
319  if (mcTrack == nullptr) continue;
320  bool isEl = IsMcPrimaryElectron(mcTrack);
321  if (isEl) {
322  //MC
323  fHM->H1("fhMismatchSource")->Fill(0);
324  fHM->H1("fhMismatchSourceMomMc")->Fill(mcTrack->GetP());
325  }
326  }
327 
328 
329  Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
330  for (Int_t iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
331  const CbmGlobalTrack* globalTrack =
332  static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
333 
334  Int_t stsId = globalTrack->GetStsTrackIndex();
335  if (stsId < 0) continue;
336  const CbmTrackMatchNew* stsTrackMatch =
337  static_cast<const CbmTrackMatchNew*>(fStsTrackMatches->At(stsId));
338  if (stsTrackMatch == nullptr) continue;
339  int stsMcTrackId = stsTrackMatch->GetMatchedLink().GetIndex();
340 
341  CbmMCTrack* mctrack = static_cast<CbmMCTrack*>(fMCTracks->At(stsMcTrackId));
342  if (mctrack == nullptr) continue;
343  Double_t mom = mctrack->GetP();
344 
345  bool isEl = IsMcPrimaryElectron(mctrack);
346  if (!isEl) continue;
347 
348  //STS
349  fHM->H1("fhMismatchSource")->Fill(1);
350  fHM->H1("fhMismatchSourceMomSts")->Fill(mom);
351 
352  if (fNofHitsInRingMap[stsMcTrackId] >= 7) {
353  //Sts-AccRich
354  fHM->H1("fhMismatchSource")->Fill(2);
355  fHM->H1("fhMismatchSourceMomStsAccRich")->Fill(mctrack->GetP());
356  }
357 
358  Int_t richId = globalTrack->GetRichRingIndex();
359 
360  // No RICH segment
361  if (richId < 0) {
362  //STS-noRICH
363  fHM->H1("fhMismatchSource")->Fill(5);
364  fHM->H1("fhMismatchSourceMomStsNoRich")->Fill(mom);
365  bool ringFound = WasRingFound(stsMcTrackId);
366  bool ringMatched = WasRingMatched(stsMcTrackId);
367  bool hasProj = HasRichProjection(stsId);
368  if (ringFound) {
369  //STS-NoRich RF
370  fHM->H1("fhMismatchSource")->Fill(6);
371  fHM->H1("fhMismatchSourceMomStsNoRichRF")->Fill(mom);
372  } else {
373  //STS-NoRich NoRF
374  fHM->H1("fhMismatchSource")->Fill(8);
375  fHM->H1("fhMismatchSourceMomStsNoRichNoRF")->Fill(mom);
376  }
377 
378  if (ringMatched) {
379  //STS-NoRich RM
380  fHM->H1("fhMismatchSource")->Fill(7);
381  fHM->H1("fhMismatchSourceMomStsNoRichRM")->Fill(mom);
382  }
383 
384  if (!hasProj) {
385  //STS-NoRich NoProj
386  fHM->H1("fhMismatchSource")->Fill(9);
387  fHM->H1("fhMismatchSourceMomStsNoRichNoProj")->Fill(mom);
388  }
389 
390  continue;
391  }
392 
393  //STS-RICH
394  fHM->H1("fhMismatchSource")->Fill(3);
395  fHM->H1("fhMismatchSourceMomStsRich")->Fill(mom);
396 
397  const CbmTrackMatchNew* richRingMatch =
398  static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
399  if (richRingMatch == nullptr) continue;
400  int richMcTrackId = richRingMatch->GetMatchedLink().GetIndex();
401  const CbmRichRing* ring =
402  static_cast<const CbmRichRing*>(fRichRings->At(richId));
403  if (nullptr == ring) continue;
404 
405  if (stsMcTrackId == richMcTrackId) {
406  //STS-RICH true
407  fHM->H1("fhMismatchSource")->Fill(4);
408  fHM->H1("fhMismatchSourceMomStsRichTrue")->Fill(mom);
409  } else {
410  //STS-RICH wrong
411  fHM->H1("fhMismatchSource")->Fill(10);
412  fHM->H1("fhMismatchSourceMomStsRichWrong")->Fill(mom);
413  if (WasRingFound(stsMcTrackId)) {
414  //STS-RICH wrong RF
415  fHM->H1("fhMismatchSource")->Fill(11);
416  fHM->H1("fhMismatchSourceMomStsRichWrongRF")->Fill(mom);
417  }
418 
419  if (WasRingFound(stsMcTrackId)) {
420  //STS-RICH wrong RM
421  fHM->H1("fhMismatchSource")->Fill(12);
422  fHM->H1("fhMismatchSourceMomStsRichWrongRM")->Fill(mom);
423  }
424  }
425  }
426 }
427 
428 bool CbmRichRecoQa::WasRingFound(Int_t mcTrackId) {
429  Int_t nofRings = fRichRings->GetEntriesFast();
430  for (Int_t iR = 0; iR < nofRings; iR++) {
431  const CbmRichRing* ring =
432  static_cast<const CbmRichRing*>(fRichRings->At(iR));
433  if (ring == nullptr) continue;
434  const CbmTrackMatchNew* richRingMatch =
435  static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(iR));
436  if (richRingMatch == nullptr) continue;
437  int richMcTrackId = richRingMatch->GetMatchedLink().GetIndex();
438  if (richMcTrackId == mcTrackId) return true;
439  }
440  return false;
441 }
442 
443 bool CbmRichRecoQa::WasRingMatched(Int_t mcTrackId) {
444  Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
445  for (Int_t iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
446  const CbmGlobalTrack* globalTrack =
447  static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
448 
449  Int_t richId = globalTrack->GetRichRingIndex();
450  if (richId < 0) continue;
451  const CbmTrackMatchNew* richRingMatch =
452  static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
453  if (richRingMatch == nullptr) continue;
454  int richMcTrackId = richRingMatch->GetMatchedLink().GetIndex();
455  if (richMcTrackId == mcTrackId) { return true; }
456  }
457 
458  return false;
459 }
460 
461 
462 bool CbmRichRecoQa::HasRichProjection(Int_t stsTrackId) {
463  if (stsTrackId < 0) { return false; }
464  FairTrackParam* pTrack = (FairTrackParam*) fRichProjections->At(stsTrackId);
465  if (pTrack == nullptr) { return false; }
466 
467  if (pTrack->GetX() == 0. && pTrack->GetY() == 0.) {
468  return false;
469  } else {
470  return true;
471  }
472 }
473 
474 void CbmRichRecoQa::FillRingTrackDistance() {
475  Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
476  for (Int_t iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
477  const CbmGlobalTrack* globalTrack =
478  static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
479 
480  Int_t stsId = globalTrack->GetStsTrackIndex();
481  Int_t richId = globalTrack->GetRichRingIndex();
482  if (stsId < 0 || richId < 0) continue;
483 
484  const CbmTrackMatchNew* stsTrackMatch =
485  static_cast<const CbmTrackMatchNew*>(fStsTrackMatches->At(stsId));
486  if (stsTrackMatch == nullptr) continue;
487  int stsMcTrackId = stsTrackMatch->GetMatchedLink().GetIndex();
488 
489  const CbmTrackMatchNew* richRingMatch =
490  static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
491  if (richRingMatch == nullptr) continue;
492  int richMcTrackId = richRingMatch->GetMatchedLink().GetIndex();
493  const CbmRichRing* ring =
494  static_cast<const CbmRichRing*>(fRichRings->At(richId));
495  if (nullptr == ring) continue;
496 
497  double rtDistance = CbmRichUtil::GetRingTrackDistance(iTrack);
498  double xc = ring->GetCenterX();
499  double yc = ring->GetCenterY();
500  int nofHits = ring->GetNofHits();
501 
502  CbmMCTrack* mctrack = static_cast<CbmMCTrack*>(fMCTracks->At(stsMcTrackId));
503  if (mctrack == nullptr) continue;
504  double mom = mctrack->GetP();
505  int charge = mctrack->GetCharge();
506 
507  bool isEl = IsMcPrimaryElectron(mctrack);
508  bool isPi = IsMcPion(mctrack);
509 
510  if (!isEl && !isPi) continue;
511 
512  string s = "";
513 
514  for (int i = 0; i < 2; i++) {
515  if (i == 0) {
516  if (isEl) {
517  s = "Primel";
518  } else if (isPi) {
519  s = "Pi";
520  } else {
521  continue;
522  }
523  } else if (i == 1) {
524  if (isEl) {
525  if (charge > 0) {
526  s = "PrimelPlus";
527  } else if (charge < 0) {
528  s = "PrimelMinus";
529  } else {
530  continue;
531  }
532  } else {
533  continue;
534  }
535  }
536 
537  if (stsMcTrackId == richMcTrackId) {
538  fHM->H2("fhRingTrackDistVsMomTruematch" + s)->Fill(mom, rtDistance);
539  fHM->H3("fhRingTrackDistVsXYTruematch" + s)->Fill(xc, yc, rtDistance);
540  fHM->H2("fhRingTrackDistVsXTruematch" + s)->Fill(xc, rtDistance);
541  fHM->H2("fhRingTrackDistVsYTruematch" + s)->Fill(abs(yc), rtDistance);
542  fHM->H2("fhRingTrackDistVsNofHitsTruematch" + s)
543  ->Fill(nofHits, rtDistance);
544 
545  if (i == 0 && isEl) {
546  //if (CbmLitGlobalElectronId::GetInstance().IsRichElectron(iTrack, mom)){
547  fHM->H2("fhRingTrackDistVsMomTruematchElId")->Fill(mom, rtDistance);
548  //}
549  }
550 
551  } else {
552  fHM->H2("fhRingTrackDistVsMomWrongmatch" + s)->Fill(mom, rtDistance);
553  }
554  }
555  }
556 }
557 
558 
559 void CbmRichRecoQa::DrawHist() {
560  cout.precision(4);
561 
562 
563  SetDefaultDrawStyle();
564  //fHM->ScaleByPattern("fh_.*", 1./fEventNum);
565 
566  {
567  fHM->CreateCanvas("fh_ring_track_distance_truematch_comparison_primel",
568  "fh_ring_track_distance_truematch_comparison_primel",
569  800,
570  800);
571  TH1D* py_minus =
572  (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematchPrimelMinus")
573  ->ProjectionY("fhRingTrackDistVsMomTruematchPrimelMinus_py")
574  ->Clone());
575  py_minus->Scale(1. / py_minus->Integral());
576  TH1D* py_plus =
577  (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematchPrimelPlus")
578  ->ProjectionY("fhRingTrackDistVsMomTruematchPrimelPlus_py")
579  ->Clone());
580  py_plus->Scale(1. / py_plus->Integral());
581  DrawH1(
582  {py_minus, py_plus},
583  {"e_{prim}^{-} (" + this->GetMeanRmsOverflowString(py_minus, false) + ")",
584  "e_{prim}^{+} (" + this->GetMeanRmsOverflowString(py_plus, false) + ")"},
585  kLinear,
586  kLog,
587  true,
588  0.5,
589  0.85,
590  0.99,
591  0.99);
592  }
593 
594  {
595  fHM->CreateCanvas("fh_ring_track_distance_truematch_comparison_elpi",
596  "fh_ring_track_distance_truematch_comparison_elpi",
597  800,
598  800);
599  TH1D* py_el =
600  (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematchPrimel")
601  ->ProjectionY("fhRingTrackDistVsMomTruematchPrimel_py")
602  ->Clone());
603  py_el->Scale(1. / py_el->Integral());
604  TH1D* py_pi = (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematchPi")
605  ->ProjectionY("fhRingTrackDistVsMomTruematchPi_py")
606  ->Clone());
607  py_pi->Scale(1. / py_pi->Integral());
608  DrawH1(
609  {py_el, py_pi},
610  {"e_{prim}^{#pm} (" + this->GetMeanRmsOverflowString(py_el, false) + ")",
611  "#pi^{#pm} (" + this->GetMeanRmsOverflowString(py_pi, false) + ")"},
612  kLinear,
613  kLog,
614  true,
615  0.5,
616  0.85,
617  0.99,
618  0.99);
619  }
620 
621  {
622  gStyle->SetPaintTextFormat("4.1f");
623  fHM->CreateCanvas("fh_mismatch_source", "fh_mismatch_source", 1000, 800);
624  Double_t nofMcEl = fHM->H1("fhMismatchSource")->GetBinContent(1);
625  cout << "nofMcEl = " << nofMcEl << endl;
626  fHM->H1("fhMismatchSource")->Scale(100. / nofMcEl);
627  DrawH1(fHM->H1("fhMismatchSource"), kLinear, kLog, "hist text");
628  fHM->H1("fhMismatchSource")->SetMarkerSize(1.9);
629 
630  vector<string> labels = {"MC",
631  "STS",
632  "STS-Acc RICH",
633  "STS-RICH",
634  "STS-RICH true",
635  "STS-NoRICH",
636  "STS-NoRICH RF",
637  "STS-NoRICH RM",
638  "STS-NoRICH NoRF",
639  "STS-NoRICH NoPrj",
640  "STS-RICH wrong",
641  "STS-RICH wrong RF",
642  "STS-RICH wrong RM"};
643  for (unsigned int i = 0; i < labels.size(); i++) {
644  fHM->H1("fhMismatchSource")
645  ->GetXaxis()
646  ->SetBinLabel(i + 1, labels[i].c_str());
647  }
648  fHM->H1("fhMismatchSource")->GetXaxis()->SetLabelSize(0.03);
649  }
650 
651  {
652  fHM->CreateCanvas(
653  "fh_mismatch_source_mom", "fh_mismatch_source_mom", 1000, 800);
654  vector<string> labels = {"MC",
655  "STS",
656  "STS-Acc RICH",
657  "STS-RICH",
658  "STS-RICH true",
659  "STS-NoRICH",
660  "STS-NoRICH RF",
661  "STS-NoRICH RM",
662  "STS-NoRICH NoRF",
663  "STS-NoRICH NoPrj",
664  "STS-RICH wrong",
665  "STS-RICH wrong RF",
666  "STS-RICH wrong RM"};
667  vector<TH1*> hists = {fHM->H1("fhMismatchSourceMomMc"),
668  fHM->H1("fhMismatchSourceMomSts"),
669  fHM->H1("fhMismatchSourceMomStsAccRich"),
670  fHM->H1("fhMismatchSourceMomStsRich"),
671  fHM->H1("fhMismatchSourceMomStsRichTrue"),
672  fHM->H1("fhMismatchSourceMomStsNoRich"),
673  fHM->H1("fhMismatchSourceMomStsNoRichRF"),
674  fHM->H1("fhMismatchSourceMomStsNoRichRM"),
675  fHM->H1("fhMismatchSourceMomStsNoRichNoRF"),
676  fHM->H1("fhMismatchSourceMomStsNoRichNoProj"),
677  fHM->H1("fhMismatchSourceMomStsRichWrong"),
678  fHM->H1("fhMismatchSourceMomStsRichWrongRF"),
679  fHM->H1("fhMismatchSourceMomStsRichWrongRM")};
680 
681  DrawH1(hists, labels, kLinear, kLog, true, 0.8, 0.35, 0.99, 0.99);
682  fHM->H1("fhMismatchSourceMomMc")->SetMinimum(0.9);
683  }
684 
685  DrawRingTrackDistHistWithSuffix("Primel");
686  DrawRingTrackDistHistWithSuffix("PrimelPlus");
687  DrawRingTrackDistHistWithSuffix("PrimelMinus");
688  DrawRingTrackDistHistWithSuffix("Pi");
689 
690  // before and after electron identification
691  {
692  TCanvas* c = fHM->CreateCanvas("fh_ring_track_distance_truematch_elid",
693  "fh_ring_track_distance_truematch_elid",
694  1800,
695  600);
696  c->Divide(3, 1);
697  c->cd(1);
698  DrawH2WithProfile(
699  fHM->H2("fhRingTrackDistVsMomTruematchPrimel"), false, true);
700  c->cd(2);
701  DrawH2WithProfile(
702  fHM->H2("fhRingTrackDistVsMomTruematchElId"), false, true);
703  c->cd(3);
704  TH1D* py =
705  (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematchPrimel")
706  ->ProjectionY(
707  string("fhRingTrackDistVsMomTruematchPrimel_py2").c_str())
708  ->Clone());
709  TH1D* pyElId =
710  (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematchElId")
711  ->ProjectionY(
712  string("fhRingTrackDistVsMomTruematchElId_py").c_str())
713  ->Clone());
714  DrawH1(
715  {py, pyElId},
716  {string("before ElId (" + this->GetMeanRmsOverflowString(py) + ")"),
717  string("after ElId (" + this->GetMeanRmsOverflowString(pyElId) + ")")},
718  kLinear,
719  kLog,
720  true,
721  0.3,
722  0.75,
723  0.99,
724  0.99);
725  fHM->H2("fhRingTrackDistVsMomTruematchPrimel")
726  ->GetYaxis()
727  ->SetRangeUser(0., 2.);
728  fHM->H2("fhRingTrackDistVsMomTruematchElId")
729  ->GetYaxis()
730  ->SetRangeUser(0., 2.);
731  }
732 }
733 
734 string CbmRichRecoQa::GetMeanRmsOverflowString(TH1* h, Bool_t withOverflow) {
735  if (withOverflow) {
736  double overflow = h->GetBinContent(h->GetNbinsX() + 1);
737  return Cbm::NumberToString<Double_t>(h->GetMean(), 2) + " / "
738  + Cbm::NumberToString<Double_t>(h->GetRMS(), 2) + " / "
739  + Cbm::NumberToString<Double_t>(
740  100. * overflow / h->Integral(0, h->GetNbinsX() + 1), 2)
741  + "%";
742  } else {
743  return Cbm::NumberToString<Double_t>(h->GetMean(), 2) + " / "
744  + Cbm::NumberToString<Double_t>(h->GetRMS(), 2);
745  }
746 }
747 
748 
749 void CbmRichRecoQa::DrawRingTrackDistHistWithSuffix(const string& s) {
750  {
751  TCanvas* c = fHM->CreateCanvas("fh_ring_track_distance_truematch_" + s,
752  "fh_ring_track_distance_truematch_" + s,
753  1800,
754  600);
755  c->Divide(3, 1);
756  c->cd(1);
757  DrawH2WithProfile(
758  fHM->H2("fhRingTrackDistVsMomTruematch" + s), false, true);
759  c->cd(2);
760  DrawH2WithProfile(
761  fHM->H2("fhRingTrackDistVsNofHitsTruematch" + s), false, true);
762  c->cd(3);
763  TH1D* py =
764  (TH1D*) (fHM->H2("fhRingTrackDistVsMomTruematch" + s)
765  ->ProjectionY(("fhRingTrackDistVsMomTruematch_py" + s).c_str())
766  ->Clone());
767  DrawH1(py);
768  DrawTextOnPad(this->GetMeanRmsOverflowString(py), 0.2, 0.9, 0.8, 0.99);
769  gPad->SetLogy(true);
770 
771  fHM->H2("fhRingTrackDistVsMomTruematch" + s)
772  ->GetYaxis()
773  ->SetRangeUser(0., 2.);
774  fHM->H2("fhRingTrackDistVsNofHitsTruematch" + s)
775  ->GetYaxis()
776  ->SetRangeUser(0., 2.);
777  }
778 
779  {
780  fHM->CreateCanvas("fh_ring_track_distance_wrongmatch" + s,
781  "fh_ring_track_distance_wrongmatch" + s,
782  600,
783  600);
784  TH1D* py = (TH1D*) (fHM->H2("fhRingTrackDistVsMomWrongmatch" + s)
785  ->ProjectionY(
786  ("fhRingTrackDistVsMomWrongmatch_py" + s).c_str())
787  ->Clone());
788  DrawH1(py);
789  DrawTextOnPad(this->GetMeanRmsOverflowString(py), 0.2, 0.9, 0.8, 0.99);
790  gPad->SetLogy(true);
791  }
792 
793  {
794  TCanvas* c = fHM->CreateCanvas("fh_ring_track_distance_vs_xy_truematch" + s,
795  "fh_ring_track_distance_vs_xy_truematch" + s,
796  1800,
797  600);
798  c->Divide(3, 1);
799  c->cd(1);
800  DrawH3Profile(
801  fHM->H3("fhRingTrackDistVsXYTruematch" + s), true, false, 0., .4);
802  c->cd(2);
803  DrawH2WithProfile(fHM->H2("fhRingTrackDistVsXTruematch" + s), false, true);
804  fHM->H2("fhRingTrackDistVsXTruematch" + s)
805  ->GetYaxis()
806  ->SetRangeUser(0., 2.);
807  c->cd(3);
808  DrawH2WithProfile(fHM->H2("fhRingTrackDistVsYTruematch" + s), false, true);
809  fHM->H2("fhRingTrackDistVsYTruematch" + s)
810  ->GetYaxis()
811  ->SetRangeUser(0., 2.);
812  }
813 }
814 
815 bool CbmRichRecoQa::IsMcPrimaryElectron(const CbmMCTrack* mctrack) {
816  if (mctrack == nullptr) return false;
817  int pdg = TMath::Abs(mctrack->GetPdgCode());
818  if (mctrack->GetGeantProcessId() == kPPrimary && pdg == 11) return true;
819  return false;
820 }
821 
822 bool CbmRichRecoQa::IsMcPion(const CbmMCTrack* mctrack) {
823  if (mctrack == nullptr) return false;
824  int pdg = TMath::Abs(mctrack->GetPdgCode());
825  if (pdg == 211) return true;
826  return false;
827 }
828 
829 void CbmRichRecoQa::Finish() {
830  DrawHist();
831  fHM->SaveCanvasToImage(fOutputDir);
832 
833  TDirectory* oldir = gDirectory;
834  TFile* outFile = FairRootManager::Instance()->GetOutFile();
835  if (outFile != NULL) {
836  outFile->cd();
837  fHM->WriteToFile();
838  }
839  gDirectory->cd(oldir->GetPath());
840 }
841 
842 void CbmRichRecoQa::DrawFromFile(const string& fileName,
843  const string& outputDir) {
844  fOutputDir = outputDir;
845 
846  if (fHM != nullptr) delete fHM;
847 
848  fHM = new CbmHistManager();
849  TFile* file = new TFile(fileName.c_str());
850  fHM->ReadFromFile(file);
851 
852  DrawHist();
853 
854  fHM->SaveCanvasToImage(fOutputDir);
855 }
856 
857 ClassImp(CbmRichRecoQa)
CbmRichRecoQa::WasRingMatched
bool WasRingMatched(Int_t mcTrackId)
Definition: alignment/CbmRichRecoQa.cxx:594
CbmRichPoint.h
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
CbmRichRecoQa::fGlobalTracks
TClonesArray * fGlobalTracks
Definition: alignment/CbmRichRecoQa.h:148
CbmHistManager::WriteToFile
void WriteToFile()
Write all histograms to current opened file.
Definition: core/base/CbmHistManager.cxx:103
CbmRichRecoQa::fStsTracks
TClonesArray * fStsTracks
Definition: alignment/CbmRichRecoQa.h:149
CbmHistManager::Create3
void Create3(const std::string &name, const std::string &title, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY, Int_t nofBinsZ, Double_t minBinZ, Double_t maxBinZ)
Helper function for creation of 3-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:140
CbmRichRecoQa::fEventNum
Int_t fEventNum
Definition: alignment/CbmRichRecoQa.h:138
CbmRichRecoQa::fOutputDir
string fOutputDir
Definition: alignment/CbmRichRecoQa.h:140
CbmHistManager::ReadFromFile
void ReadFromFile(TFile *file)
Read histograms from file.
Definition: core/base/CbmHistManager.cxx:110
DrawH3Profile
TH2D * DrawH3Profile(TH3 *h, Bool_t drawMean, Bool_t doGaussFit, Double_t zMin, Double_t zMax)
Definition: CbmDrawHist.cxx:369
CbmRichRecoQa::IsMcPion
static Bool_t IsMcPion(const CbmMCTrack *mctrack)
Definition: alignment/CbmRichRecoQa.cxx:1200
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition: CbmDigiManager.cxx:71
CbmHistManager::Create2
void Create2(const std::string &name, const std::string &title, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY)
Helper function for creation of 2-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:104
CbmGlobalTrack::GetRichRingIndex
Int_t GetRichRingIndex() const
Definition: CbmGlobalTrack.h:41
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition: CbmHistManager.h:190
CbmGlobalTrack.h
CbmRichRing::GetNofHits
Int_t GetNofHits() const
Definition: CbmRichRing.h:40
CbmRichRing
Definition: CbmRichRing.h:17
CbmRichRecoQa
Definition: alignment/CbmRichRecoQa.h:21
CbmRichRecoQa::FillRichRingNofHits
void FillRichRingNofHits()
Fill map mcTrackId -> nof RICH hits.
Definition: alignment/CbmRichRecoQa.cxx:450
CbmRichRing.h
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichDraw.h
CbmHistManager.h
Histogram manager.
CbmMatchRecoToMC.h
FairTask for matching RECO data to MC.
CbmHistManager::H3
TH3 * H3(const std::string &name) const
Return pointer to TH3 histogram.
Definition: CbmHistManager.h:210
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition: CbmDigiManager.h:93
CbmRichRecoQa::fMCTracks
TClonesArray * fMCTracks
Definition: alignment/CbmRichRecoQa.h:142
CbmMCTrack::GetCharge
Double_t GetCharge() const
Charge of the associated particle.
Definition: CbmMCTrack.cxx:146
CbmRichGeoManager.h
CbmRichRecoQa::fRichRings
TClonesArray * fRichRings
Definition: alignment/CbmRichRecoQa.h:146
CbmRichUtil::GetRingTrackDistance
static Double_t GetRingTrackDistance(Int_t globalTrackId)
Definition: alignment/CbmRichUtil.h:20
DrawTextOnPad
void DrawTextOnPad(const string &text, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: CbmDrawHist.cxx:253
DrawH1
void DrawH1(TH1 *hist, HistScale logx, HistScale logy, const string &drawOpt, Int_t color, Int_t lineWidth, Int_t lineStyle, Int_t markerSize, Int_t markerStyle)
Definition: CbmDrawHist.cxx:49
CbmGlobalTrack::GetStsTrackIndex
Int_t GetStsTrackIndex() const
Definition: CbmGlobalTrack.h:38
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmRichRecoQa::DrawFromFile
void DrawFromFile(const string &fileName, const string &outputDir)
Draw histogram from file.
Definition: alignment/CbmRichRecoQa.cxx:1222
CbmRichRecoQa::fRichHits
TClonesArray * fRichHits
Definition: alignment/CbmRichRecoQa.h:145
CbmRichRecoQa.h
CbmRichRecoQa::WasRingFound
bool WasRingFound(Int_t mcTrackId)
Definition: alignment/CbmRichRecoQa.cxx:579
kLinear
@ kLinear
Definition: CbmDrawHist.h:79
CbmTrackMatchNew.h
CbmRichRecoQa::GetMeanRmsOverflowString
string GetMeanRmsOverflowString(TH1 *h, Bool_t withOverflow=true)
Return string with mean, RMS and overflow percent for input TH1.
Definition: alignment/CbmRichRecoQa.cxx:955
CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit
static std::vector< std::pair< Int_t, Int_t > > GetMcTrackMotherIdsForRichHit(CbmDigiManager *digiMan, const CbmRichHit *hit, CbmMCDataArray *richPoints, CbmMCDataArray *mcTracks, Int_t eventNumber)
Return McTrack Ids for RICH hit C++11 efficient way to return vector.
Definition: CbmMatchRecoToMC.cxx:821
CbmRichRecoQa::Finish
virtual void Finish()
Inherited from FairTask.
Definition: alignment/CbmRichRecoQa.cxx:1207
DrawH2WithProfile
void DrawH2WithProfile(TH2 *hist, Bool_t doGaussFit, Bool_t drawOnlyMean, const string &drawOpt2D, Int_t profileColor, Int_t profileLineWidth)
Definition: CbmDrawHist.cxx:296
CbmHistManager::Create1
void Create1(const std::string &name, const std::string &title, Int_t nofBins, Double_t minBin, Double_t maxBin)
Helper function for creation of 1-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:81
CbmRichRecoQa::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition: alignment/CbmRichRecoQa.h:150
CbmRichRecoQa::DrawHist
void DrawHist()
Draw histograms.
Definition: alignment/CbmRichRecoQa.cxx:743
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmMCTrack::GetGeantProcessId
UInt_t GetGeantProcessId() const
Definition: CbmMCTrack.h:69
CbmRichRecoQa::fHM
CbmHistManager * fHM
Definition: alignment/CbmRichRecoQa.h:136
CbmRichRecoQa::CbmRichRecoQa
CbmRichRecoQa()
Standard constructor.
Definition: alignment/CbmRichRecoQa.cxx:40
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition: CbmHistManager.h:170
CbmRichRecoQa::InitHistograms
void InitHistograms()
Initialize histograms.
Definition: alignment/CbmRichRecoQa.cxx:115
CbmRichRecoQa::FillRingTrackDistance
void FillRingTrackDistance()
Fill histogramms related to ring track distance.
Definition: alignment/CbmRichRecoQa.cxx:624
fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmTofAnaTestbeam.cxx:88
CbmUtils.h
CbmHistManager::CreateCanvas
TCanvas * CreateCanvas(const std::string &name, const std::string &title, Int_t width, Int_t height)
Create and draw TCanvas and store pointer to it.
Definition: core/base/CbmHistManager.cxx:267
CbmHistManager::SaveCanvasToImage
void SaveCanvasToImage(const std::string &outputDir, const std::string &options="png,eps")
Save all stored canvases to images.
Definition: core/base/CbmHistManager.cxx:276
CbmRichRecoQa::RingTrackMismatchSource
void RingTrackMismatchSource()
Fill histograms related to study of the source of ring-track mismatch.
Definition: alignment/CbmRichRecoQa.cxx:465
CbmRichRecoQa::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: alignment/CbmRichRecoQa.cxx:59
CbmGlobalTrack
Definition: CbmGlobalTrack.h:26
CbmRichRecoQa::IsMcPrimaryElectron
static Bool_t IsMcPrimaryElectron(const CbmMCTrack *mctrack)
Definition: alignment/CbmRichRecoQa.cxx:1193
CbmMCTrack.h
CbmDigiManager.h
CbmRichRecoQa::DrawRingTrackDistHistWithSuffix
void DrawRingTrackDistHistWithSuffix(const string &suffix)
Draw histograms related to ring-track distance for pions or electrons (+/-).
Definition: alignment/CbmRichRecoQa.cxx:969
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmRichRecoQa::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: alignment/CbmRichRecoQa.cxx:441
CbmRichRecoQa::fRichProjections
TClonesArray * fRichProjections
Definition: alignment/CbmRichRecoQa.h:151
CbmRichRing::GetCenterY
Float_t GetCenterY() const
Definition: CbmRichRing.h:81
CbmRichRecoQa::fRichPoints
TClonesArray * fRichPoints
Definition: alignment/CbmRichRecoQa.h:143
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
SetDefaultDrawStyle
void SetDefaultDrawStyle()
Definition: CbmDrawHist.cxx:33
CbmLitGlobalElectronId.h
CbmRichRecoQa::fDigiMan
CbmDigiManager * fDigiMan
Definition: qa/CbmRichRecoQa.h:148
CbmRichHit.h
CbmRichRecoQa::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: alignment/CbmRichRecoQa.h:147
CbmRichRecoQa::fNofHitsInRingMap
std::map< Int_t, Int_t > fNofHitsInRingMap
Definition: alignment/CbmRichRecoQa.h:154
CbmMCTrack::GetP
Double_t GetP() const
Definition: CbmMCTrack.h:100
kLog
@ kLog
Definition: CbmDrawHist.h:78
CbmRichRing::GetCenterX
Float_t GetCenterX() const
Definition: CbmRichRing.h:80
CbmRichHit
Definition: CbmRichHit.h:19
CbmRichRecoQa::HasRichProjection
bool HasRichProjection(Int_t stsTrackId)
Definition: alignment/CbmRichRecoQa.cxx:612