CbmRoot
CbmRichRecoTbQa.cxx
Go to the documentation of this file.
1 #include "CbmRichRecoTbQa.h"
2 #include "TCanvas.h"
3 #include "TClonesArray.h"
4 #include "TEllipse.h"
5 #include "TF1.h"
6 #include "TH1.h"
7 #include "TH1D.h"
8 #include "TMCProcess.h"
9 #include "TMarker.h"
10 #include "TStyle.h"
11 #include <TFile.h>
12 
13 #include "CbmDigiManager.h"
14 #include "CbmDrawHist.h"
15 #include "CbmGlobalTrack.h"
16 #include "CbmMCDataArray.h"
17 #include "CbmMCDataManager.h"
18 #include "CbmMCTrack.h"
19 #include "CbmMatchRecoToMC.h"
20 #include "CbmRichDigi.h"
21 #include "CbmRichDraw.h"
22 #include "CbmRichGeoManager.h"
23 #include "CbmRichHit.h"
24 #include "CbmRichPoint.h"
25 #include "CbmRichRing.h"
26 #include "CbmRichUtil.h"
27 #include "CbmStsPoint.h"
28 #include "CbmTrackMatchNew.h"
29 #include "CbmUtils.h"
30 #include "FairMCPoint.h"
31 #include "FairTrackParam.h"
32 #include "elid/CbmLitGlobalElectronId.h"
33 
34 #include "FairEventHeader.h"
35 #include "FairLogger.h"
36 #include "FairRunAna.h"
37 #include "FairRunSim.h"
38 
39 #include "CbmHistManager.h"
40 #include "CbmUtils.h"
41 
42 #include <iostream>
43 #include <sstream>
44 #include <string>
45 
46 #include "CbmMCEventList.h"
47 
48 using namespace std;
49 
50 CbmRichRecoTbQa::CbmRichRecoTbQa()
51  : FairTask("CbmRichRecoTbQa")
52  , fHM(nullptr)
53  , fTimeSliceNum(0.)
54  , fNofLogEvents(150)
55  , fOutputDir("")
56  , fMCTracks(nullptr)
57  , fRichPoints(nullptr)
58  , fStsPoints(nullptr)
59  , fDigiMan(nullptr)
60  , fRichHits(nullptr)
61  , fRichRings(nullptr)
62  , fRichRingMatches(nullptr)
63  , fEventList(nullptr)
64  , fRecRings() {}
65 
66 
67 InitStatus CbmRichRecoTbQa::Init() {
68  cout << "CbmRichRecoTbQa::Init" << endl;
69  FairRootManager* ioman = FairRootManager::Instance();
70  if (nullptr == ioman) {
71  Fatal("CbmRichRecoTbQa::Init", "RootManager not instantised!");
72  }
73 
74  CbmMCDataManager* mcManager =
75  (CbmMCDataManager*) ioman->GetObject("MCDataManager");
76  if (mcManager == nullptr)
77  LOG(fatal)
78  << "CbmRichRecoTbQa::ReadAndCreateDataBranches() NULL MCDataManager.";
79 
80  fMCTracks = mcManager->InitBranch("MCTrack");
81  fRichPoints = mcManager->InitBranch("RichPoint");
82  fStsPoints = mcManager->InitBranch("StsPoint");
83 
84  fDigiMan = CbmDigiManager::Instance();
85  fDigiMan->Init();
86 
87  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
88  if (nullptr == fRichHits) { Fatal("CbmRichRecoTbQa::Init", "No Rich hits!"); }
89 
90  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
91  if (nullptr == fRichRings) {
92  Fatal("CbmRichRecoTbQa::Init", "No Rich rings!");
93  }
94 
95  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
96  if (nullptr == fRichRingMatches) {
97  Fatal("CbmRichRecoTbQa::Init", "No RichRingMatch!");
98  }
99 
100  fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList.");
101  if (nullptr == fEventList) {
102  Fatal("CbmRichRecoTbQa::Init", "No MCEventList!");
103  }
104 
105  InitHistograms();
106 
107  return kSUCCESS;
108 }
109 
110 void CbmRichRecoTbQa::InitHistograms() {
111  fHM = new CbmHistManager();
112 
113  Int_t nBinsLog = 20000;
114  Double_t minLog = 0.;
115  Double_t maxLog = 10000.;
116 
117 
118  fHM->Create1<TH1D>("fhNofRichDigisPerTS",
119  "fhNofRichDigisPerTS;RICH digis per time slice;Yield",
120  100,
121  -1,
122  -1.);
123  fHM->Create1<TH1D>("fhNofRichHitsPerTS",
124  "fhNofRichHitsPerTS;RICH hits per time slice;Yield",
125  100,
126  -1.,
127  -1.);
128  fHM->Create1<TH1D>("fhNofRichRingsPerTS",
129  "fhNofRichRingsPerTS;RICH rings per time slice;Yield",
130  100,
131  -1.,
132  -1.);
133 
134  fHM->Create1<TH1D>("fhRichDigiTimeLog",
135  "fhRichDigiTimeLog;RICH digi time [ns];Yield",
136  nBinsLog,
137  minLog,
138  maxLog);
139 
140  fHM->Create1<TH1D>("fhRichRingTimeLog",
141  "fhRichRingTimeLog;RICH ring time [ns];Yield",
142  nBinsLog,
143  minLog,
144  maxLog);
145 
146  fHM->Create1<TH1D>("fhRichPointTime",
147  "fhRichPointTime;RICH MC point time [ns];Yield",
148  100,
149  0.,
150  100.);
151  fHM->Create1<TH1D>("fhRichPointTimeChPhoton",
152  "fhRichPointTimeChPhoton;RICH MC point time [ns];Yield",
153  100,
154  0.,
155  100.);
156  fHM->Create1<TH1D>("fhRichPointTimeNotChPhoton",
157  "fhRichPointTimeNotChPhoton;RICH MC point time [ns];Yield",
158  100,
159  0.,
160  100.);
161  fHM->Create1<TH1D>("fhRichPointTimePrimEl",
162  "fhRichPointTimePrimEl;RICH MC point time [ns];Yield",
163  100,
164  0.,
165  100.);
166  fHM->Create1<TH1D>("fhRichPointTimeSecEl",
167  "fhRichPointTimeSecEl;RICH MC point time [ns];Yield",
168  100,
169  0.,
170  100.);
171  fHM->Create1<TH1D>("fhRichPointTimeOther",
172  "fhRichPointTimeOther;RICH MC point time [ns];Yield",
173  100,
174  0.,
175  100.);
176  fHM->Create1<TH1D>("fhRichPointTimePion",
177  "fhRichPointTimePion;RICH MC point time [ns];Yield",
178  100,
179  0.,
180  100.);
181 
182  fHM->Create1<TH1D>("fhStsPointTime",
183  "fhStsPointTime;STS MC point time [ns];Yield",
184  400,
185  0.,
186  400.);
187 
188  // -1 for noise hits
189  for (int iEv = -1; iEv < fNofLogEvents; iEv++) {
190  string richDigiStr = "fhRichDigiTimeLog_" + to_string(iEv);
191  fHM->Create1<TH1D>(
192  richDigiStr, richDigiStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
193 
194  string richRingStr = "fhRichRingTimeLog_" + to_string(iEv);
195  fHM->Create1<TH1D>(
196  richRingStr, richRingStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
197 
198  string richStr = "fhRichPointTimeLog_" + to_string(iEv);
199  fHM->Create1<TH1D>(
200  richStr, richStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
201 
202  string stsStr = "fhStsPointTimeLog_" + to_string(iEv);
203  fHM->Create1<TH1D>(
204  stsStr, stsStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
205  }
206 
207  fHM->Create1<TH1D>("fhTimeBetweenEvents",
208  "fhTimeBetweenEvents;Time between events [ns];Yield",
209  100,
210  0.,
211  500.);
212 
213 
214  // efficiency histograms
215  fHM->Create1<TH1D>("fhMomElAcc", "fhMomElAcc;P [GeV/c];Yield", 10, 0., 10.);
216  fHM->Create1<TH1D>("fhMomElRec", "fhMomElRec;P [GeV/c];Yield", 10, 0., 10.);
217  fHM->Create1<TH1D>(
218  "fhMomRefElAcc", "fhMomRefElAcc;P [GeV/c];Yield", 10, 0., 10.);
219  fHM->Create1<TH1D>(
220  "fhMomRefElRec", "fhMomRefElRec;P [GeV/c];Yield", 10, 0., 10.);
221 
222 
223  fHM->Create1<TH1D>(
224  "fhNofHitsElAcc", "fhNofHitsElAcc;Nof hits in ring;Yield", 20, -.5, 39.5);
225  fHM->Create1<TH1D>(
226  "fhNofHitsElRec", "fhNofHitsElRec;Nof hits in ring;Yield", 20, -.5, 39.5);
227  fHM->Create1<TH1D>("fhNofHitsRefElAcc",
228  "fhNofHitsRefElAcc;Nof hits in ring;Yield",
229  20,
230  -.5,
231  39.5);
232  fHM->Create1<TH1D>("fhNofHitsRefElRec",
233  "fhNofHitsRefElRec;Nof hits in ring;Yield",
234  20,
235  -.5,
236  39.5);
237 
238  fHM->Create1<TH1D>(
239  "fhEventMultElAcc", "fhEventMultElAcc;Nof MC tracks;Yield", 10, 0., 2000);
240  fHM->Create1<TH1D>(
241  "fhEventMultElRec", "fhEventMultElRec;Nof MC tracks;Yield", 10, 0, 2000);
242 }
243 
244 void CbmRichRecoTbQa::Exec(Option_t* /*option*/) {
245  fTimeSliceNum++;
246  int nofMcEvents = fEventList->GetNofEvents();
247 
248  cout << "CbmRichRecoTbQa, time slice:" << fTimeSliceNum
249  << " nofMcEvents:" << nofMcEvents << endl;
250 
251  Process();
252 
253  RingRecoEfficiency();
254 }
255 
256 void CbmRichRecoTbQa::Process() {
257  Int_t nofRichDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
258  Int_t nofRichHits = fRichHits->GetEntries();
259  Int_t nofRichRings = fRichRings->GetEntries();
260  fHM->H1("fhNofRichDigisPerTS")->Fill(nofRichDigis);
261  fHM->H1("fhNofRichHitsPerTS")->Fill(nofRichHits);
262  fHM->H1("fhNofRichRingsPerTS")->Fill(nofRichRings);
263 
264  LOG(debug) << "nofRichDigis:" << nofRichDigis;
265  for (Int_t iDigi = 0; iDigi < nofRichDigis; iDigi++) {
266  const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(iDigi);
267  const CbmMatch* digiMatch = fDigiMan->GetMatch(ECbmModuleId::kRich, iDigi);
268  fHM->H1("fhRichDigiTimeLog")->Fill(digi->GetTime());
269 
270  Int_t eventNum = digiMatch->GetMatchedLink().GetEntry();
271  Int_t index = digiMatch->GetMatchedLink().GetIndex();
272 
273  if (eventNum < 0 || index < 0) {
274  fHM->H1("fhRichDigiTimeLog_-1")->Fill(digi->GetTime());
275  } else {
276  if (eventNum < fNofLogEvents) {
277  string hName = "fhRichDigiTimeLog_" + to_string(eventNum);
278  fHM->H1(hName)->Fill(digi->GetTime());
279  }
280  }
281  }
282 
283  for (Int_t iR = 0; iR < nofRichRings; iR++) {
284  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
285  CbmTrackMatchNew* ringMatch =
286  static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iR));
287  if (ring == nullptr || ringMatch == nullptr) continue;
288 
289  fHM->H1("fhRichRingTimeLog")->Fill(ring->GetTime());
290 
291  Int_t eventNum = ringMatch->GetMatchedLink().GetEntry();
292  Int_t index = ringMatch->GetMatchedLink().GetIndex();
293 
294  if (eventNum < 0 || index < 0) {
295  fHM->H1("fhRichRingTimeLog_-1")->Fill(ring->GetTime());
296  } else {
297  if (eventNum < fNofLogEvents) {
298  string hName = "fhRichRingTimeLog_" + to_string(eventNum);
299  fHM->H1(hName)->Fill(ring->GetTime());
300  }
301  }
302  }
303 }
304 
305 void CbmRichRecoTbQa::ProcessMc() {
306  int fileId = 0;
307 
308  // int nMCEvents = fEventList->GetNofEvents();
309  for (int iEv = 0; fEventList->GetEventTime(iEv, fileId) > 0.; iEv++) {
310  double dT = (iEv == 0) ? fEventList->GetEventTime(iEv, fileId)
311  : fEventList->GetEventTime(iEv, fileId)
312  - fEventList->GetEventTime(iEv - 1, fileId);
313  fHM->H1("fhTimeBetweenEvents")->Fill(dT);
314  }
315 
316  for (Int_t iEv = 0; fRichPoints->Size(fileId, iEv) >= 0; iEv++) {
317  Int_t nofRichPoints = fRichPoints->Size(fileId, iEv);
318  for (Int_t iP = 0; iP < nofRichPoints; iP++) {
319  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->Get(fileId, iEv, iP);
320  fHM->H1("fhRichPointTime")->Fill(point->GetTime());
321  if (IsCherenkovPhoton(point, fileId, iEv)) {
322  fHM->H1("fhRichPointTimeChPhoton")->Fill(point->GetTime());
323  } else {
324  fHM->H1("fhRichPointTimeNotChPhoton")->Fill(point->GetTime());
325  }
326 
327  if (IsCherenkovPhotonFromPrimaryElectron(point, fileId, iEv)) {
328  fHM->H1("fhRichPointTimePrimEl")->Fill(point->GetTime());
329  } else if (IsCherenkovPhotonFromPion(point, fileId, iEv)) {
330  fHM->H1("fhRichPointTimePion")->Fill(point->GetTime());
331  }
332  if (IsCherenkovPhotonFromSecondaryElectron(point, fileId, iEv)) {
333  fHM->H1("fhRichPointTimeSecEl")->Fill(point->GetTime());
334  } else {
335  fHM->H1("fhRichPointTimeOther")->Fill(point->GetTime());
336  }
337 
338  if (iEv < fNofLogEvents) {
339  string richStr = "fhRichPointTimeLog_" + to_string(iEv);
340  double eventTime = fEventList->GetEventTime(iEv, fileId);
341  fHM->H1(richStr)->Fill(eventTime + point->GetTime());
342  }
343  }
344  }
345 
346 
347  for (Int_t iEv = 0; fStsPoints->Size(fileId, iEv) >= 0; iEv++) {
348  Int_t nofStsPoints = fStsPoints->Size(fileId, iEv);
349  for (Int_t j = 0; j < nofStsPoints; j++) {
350  const CbmStsPoint* point = (CbmStsPoint*) fStsPoints->Get(fileId, iEv, j);
351  fHM->H1("fhStsPointTime")->Fill(point->GetTime());
352  if (iEv < fNofLogEvents) {
353  string stsStr = "fhStsPointTimeLog_" + to_string(iEv);
354  double eventTime = fEventList->GetEventTime(iEv, fileId);
355  fHM->H1(stsStr)->Fill(eventTime + point->GetTime());
356  }
357  }
358  }
359 }
360 
361 Int_t CbmRichRecoTbQa::GetNofPrimaryMcTracks(Int_t iEv) {
362  Int_t fileId = 0;
363  Int_t counter = 0;
364  Int_t nofMcTracks = fMCTracks->Size(fileId, iEv);
365  for (Int_t j = 0; j < nofMcTracks; j++) {
366  const CbmMCTrack* track =
367  static_cast<CbmMCTrack*>(fMCTracks->Get(fileId, iEv, j));
368  if (track->GetGeantProcessId() == kPPrimary) counter++;
369  }
370  return counter;
371 }
372 
373 void CbmRichRecoTbQa::RingRecoEfficiency() {
374  map<pair<Int_t, Int_t>, Int_t> nofHitsInRing;
375  Int_t nofRichHits = fRichHits->GetEntriesFast();
376  for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
377  const CbmRichHit* hit = static_cast<const CbmRichHit*>(fRichHits->At(iHit));
378  if (nullptr == hit) continue;
379  vector<pair<Int_t, Int_t>> motherIds =
380  CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(
381  fDigiMan, hit, fRichPoints, fMCTracks, -1);
382  for (UInt_t i = 0; i < motherIds.size(); i++) {
383  nofHitsInRing[motherIds[i]]++;
384  }
385  }
386 
387  Int_t fileId = 0;
388  for (Int_t iEv = 0; fMCTracks->Size(fileId, iEv) >= 0; iEv++) {
389  Int_t nofMcTracks = fMCTracks->Size(fileId, iEv);
390  Int_t nofPrimMcTracks = GetNofPrimaryMcTracks(iEv);
391  for (Int_t j = 0; j < nofMcTracks; j++) {
392  const CbmMCTrack* track =
393  static_cast<CbmMCTrack*>(fMCTracks->Get(fileId, iEv, j));
394  pair<Int_t, Int_t> val = std::make_pair(iEv, j);
395  Int_t nofRichHitsInRing = nofHitsInRing[val];
396  Double_t mom = track->GetP();
397  if (nofRichHitsInRing >= 7 && IsMcPrimaryElectron(track)) {
398  fHM->H1("fhMomElAcc")->Fill(mom);
399  fHM->H1("fhNofHitsElAcc")->Fill(nofRichHitsInRing);
400  fHM->H1("fhEventMultElAcc")->Fill(nofPrimMcTracks);
401  if (nofRichHitsInRing >= 15) {
402  fHM->H1("fhMomRefElAcc")->Fill(mom);
403  fHM->H1("fhNofHitsRefElAcc")->Fill(nofRichHitsInRing);
404  }
405  }
406  }
407  }
408 
409  Int_t nofRings = fRichRings->GetEntriesFast();
410  for (Int_t iRing = 0; iRing < nofRings; iRing++) {
411  const CbmTrackMatchNew* richRingMatch =
412  static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(iRing));
413  Bool_t isRichOk = richRingMatch->GetTrueOverAllHitsRatio() >= 0.6;
414 
415  const CbmMCTrack* track = static_cast<CbmMCTrack*>(
416  fMCTracks->Get(fileId,
417  richRingMatch->GetMatchedLink().GetEntry(),
418  richRingMatch->GetMatchedLink().GetIndex()));
419  Int_t nofPrimMcTracks =
420  GetNofPrimaryMcTracks(richRingMatch->GetMatchedLink().GetEntry());
421  pair<Int_t, Int_t> val =
422  std::make_pair(richRingMatch->GetMatchedLink().GetEntry(),
423  richRingMatch->GetMatchedLink().GetIndex());
424  Int_t nofRichHitsInRing = nofHitsInRing[val];
425  Double_t mom = track->GetP();
426  Bool_t isClone =
427  (std::find(fRecRings.begin(), fRecRings.end(), val) != fRecRings.end());
428  if (nofRichHitsInRing >= 7 && isRichOk && IsMcPrimaryElectron(track)
429  && !isClone) {
430  fHM->H1("fhMomElRec")->Fill(mom);
431  fHM->H1("fhNofHitsElRec")->Fill(nofRichHitsInRing);
432  fHM->H1("fhEventMultElRec")->Fill(nofPrimMcTracks);
433  fRecRings.push_back(val);
434  if (nofRichHitsInRing >= 15) {
435  fHM->H1("fhMomRefElRec")->Fill(mom);
436  fHM->H1("fhNofHitsRefElRec")->Fill(nofRichHitsInRing);
437  }
438  }
439  }
440 }
441 
442 
443 Bool_t CbmRichRecoTbQa::IsCherenkovPhoton(const CbmRichPoint* point,
444  Int_t fileId,
445  Int_t eventId) {
446  if (point == nullptr) return false;
447  Int_t trackId = point->GetTrackID();
448  if (trackId < 0) return false;
449  CbmMCTrack* p = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
450  return (p != nullptr && TMath::Abs(p->GetPdgCode()) == 50000050);
451 }
452 
453 Bool_t
454 CbmRichRecoTbQa::IsCherenkovPhotonFromPrimaryElectron(const CbmRichPoint* point,
455  Int_t fileId,
456  Int_t eventId) {
457  if (point == nullptr) return false;
458  Int_t trackId = point->GetTrackID();
459  if (trackId < 0) return false;
460  CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
461  if (mcTrack == nullptr || TMath::Abs(mcTrack->GetPdgCode()) != 50000050)
462  return false;
463 
464  Int_t motherId = mcTrack->GetMotherId();
465  if (motherId < 0) return false;
466  CbmMCTrack* mcMotherTrack =
467  (CbmMCTrack*) fMCTracks->Get(fileId, eventId, motherId);
468  return IsMcPrimaryElectron(mcMotherTrack);
469 }
470 
471 Bool_t CbmRichRecoTbQa::IsCherenkovPhotonFromSecondaryElectron(
472  const CbmRichPoint* point,
473  Int_t fileId,
474  Int_t eventId) {
475  if (point == nullptr) return false;
476  Int_t trackId = point->GetTrackID();
477  if (trackId < 0) return false;
478  CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
479  if (mcTrack == nullptr || TMath::Abs(mcTrack->GetPdgCode()) != 50000050)
480  return false;
481 
482  Int_t motherId = mcTrack->GetMotherId();
483  if (motherId < 0) return false;
484  CbmMCTrack* mcMotherTrack =
485  (CbmMCTrack*) fMCTracks->Get(fileId, eventId, motherId);
486  if (mcMotherTrack == nullptr) return false;
487  int pdg = TMath::Abs(mcMotherTrack->GetPdgCode());
488  if (mcMotherTrack->GetGeantProcessId() != kPPrimary && pdg == 11) return true;
489 
490  return false;
491 }
492 
493 Bool_t CbmRichRecoTbQa::IsMcPrimaryElectron(const CbmMCTrack* mctrack) {
494  if (mctrack == nullptr) return false;
495  int pdg = TMath::Abs(mctrack->GetPdgCode());
496  if (mctrack->GetGeantProcessId() == kPPrimary && pdg == 11) return true;
497  return false;
498 }
499 
500 Bool_t CbmRichRecoTbQa::IsMcPion(const CbmMCTrack* mctrack) {
501  if (mctrack == nullptr) return false;
502  int pdg = TMath::Abs(mctrack->GetPdgCode());
503  if (pdg == 211) return true;
504  return false;
505 }
506 
507 Bool_t CbmRichRecoTbQa::IsCherenkovPhotonFromPion(const CbmRichPoint* point,
508  Int_t fileId,
509  Int_t eventId) {
510  if (point == nullptr) return false;
511  Int_t trackId = point->GetTrackID();
512  if (trackId < 0) return false;
513  CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
514  if (mcTrack == nullptr || TMath::Abs(mcTrack->GetPdgCode()) != 50000050)
515  return false;
516 
517  Int_t motherId = mcTrack->GetMotherId();
518  if (motherId < 0) return false;
519  CbmMCTrack* mcMotherTrack =
520  (CbmMCTrack*) fMCTracks->Get(fileId, eventId, motherId);
521  return IsMcPion(mcMotherTrack);
522 }
523 
524 void CbmRichRecoTbQa::DrawHist() {
525  SetDefaultDrawStyle();
526 
527 
528  {
529  fHM->CreateCanvas(
530  "rich_tb_fhTimeBetweenEvents", "rich_tb_fhTimeBetweenEvents", 800, 800);
531  DrawH1(fHM->H1("fhTimeBetweenEvents"), kLinear, kLog);
532  Double_t min = 0.;
533  Double_t max = 10.;
534  Double_t partIntegral =
535  fHM->H1("fhTimeBetweenEvents")
536  ->Integral(fHM->H1("fhTimeBetweenEvents")->FindBin(min),
537  fHM->H1("fhTimeBetweenEvents")->FindBin(max));
538  Double_t allIntegral = fHM->H1("fhTimeBetweenEvents")->Integral();
539  Double_t ratio = 100. * partIntegral / allIntegral;
540  cout << ratio << "% of all time between events are in range (" << min
541  << " ," << max << ") ns" << endl;
542  }
543 
544 
545  {
546  fHM->CreateCanvas(
547  "rich_tb_fhRichPointTime", "rich_tb_fhRichPointTime", 800, 800);
548  DrawH1(fHM->H1("fhRichPointTime"), kLinear, kLog);
549  }
550 
551  {
552  TCanvas* c = fHM->CreateCanvas("rich_tb_fhRichPointTimeSources",
553  "rich_tb_fhRichPointTimeSources",
554  1600,
555  800);
556  c->Divide(2, 1);
557  c->cd(1);
558  DrawH1({fHM->H1("fhRichPointTimeChPhoton"),
559  fHM->H1("fhRichPointTimeNotChPhoton")},
560  {"Ch. Photons", "Not Ch. Photons"},
561  kLinear,
562  kLog);
563 
564  c->cd(2);
565  DrawH1({fHM->H1("fhRichPointTimeSecEl"),
566  fHM->H1("fhRichPointTimePrimEl"),
567  fHM->H1("fhRichPointTimePion"),
568  fHM->H1("fhRichPointTimeOther")},
569  {"e^{#pm}_{sec}", "e^{#pm}_{prim}", "#pi^{#pm}", "Other"},
570  kLinear,
571  kLog);
572  }
573 
574 
575  {
576  fHM->CreateCanvas(
577  "rich_tb_fhStsPointTime", "rich_tb_fhStsPointTime", 800, 800);
578  DrawH1(fHM->H1("fhStsPointTime"), kLinear, kLog);
579  }
580 
581  {
582  fHM->CreateCanvas(
583  "rich_tb_fhRichPointTimeLog", "rich_tb_fhRichPointTimeLog", 1600, 600);
584  DrawTimeLog("fhRichPointTimeLog", fNofLogEvents);
585  }
586 
587  {
588  fHM->CreateCanvas(
589  "rich_tb_fhStsPointTimeLog", "rich_tb_fhStsPointTimeLog", 1600, 600);
590  DrawTimeLog("fhStsPointTimeLog", fNofLogEvents);
591  }
592 
593  {
594  fHM->CreateCanvas("rich_tb_reco_fhRichDigiTimeLogAll",
595  "rich_tb_reco_fhRichDigiTimeLogAll",
596  1600,
597  600);
598  DrawH1(fHM->H1("fhRichDigiTimeLog"), kLinear, kLog);
599  gPad->SetLeftMargin(0.07);
600  gPad->SetRightMargin(0.05);
601  }
602 
603  {
604  fHM->CreateCanvas("rich_tb_reco_fhRichDigiTimeLog",
605  "rich_tb_reco_fhRichDigiTimeLog",
606  1600,
607  600);
608  DrawTimeLog("fhRichDigiTimeLog", fNofLogEvents, true);
609  }
610 
611  {
612  fHM->CreateCanvas("rich_tb_reco_fhRichRingTimeLogAll",
613  "rich_tb_reco_fhRichRingTimeLogAll",
614  1600,
615  600);
616  DrawH1(fHM->H1("fhRichRingTimeLog"), kLinear, kLog);
617  gPad->SetLeftMargin(0.07);
618  gPad->SetRightMargin(0.05);
619  }
620 
621  {
622  fHM->CreateCanvas("rich_tb_reco_timelog_hits_rings",
623  "rich_tb_reco_timelog_hits_rings",
624  1600,
625  600);
626  DrawH1({fHM->H1("fhRichDigiTimeLog"), fHM->H1("fhRichRingTimeLog")},
627  {"Hits", "Rings"});
628  gPad->SetLeftMargin(0.07);
629  gPad->SetRightMargin(0.05);
630  }
631 
632  {
633  fHM->CreateCanvas("rich_tb_reco_fhRichRingTimeLog",
634  "rich_tb_reco_fhRichRingTimeLog",
635  1600,
636  600);
637  DrawTimeLog("fhRichRingTimeLog", fNofLogEvents, true);
638  }
639 
640  {
641  TCanvas* c = fHM->CreateCanvas("rich_tb_reco_fhRichObectsPerTimeSlice",
642  "rich_tb_reco_fhRichObectsPerTimeSlice",
643  1500,
644  500);
645  c->Divide(3, 1);
646  c->cd(1);
647  DrawH1(fHM->H1("fhNofRichDigisPerTS"), kLinear, kLog);
648  c->cd(2);
649  DrawH1(fHM->H1("fhNofRichHitsPerTS"), kLinear, kLog);
650  c->cd(3);
651  DrawH1(fHM->H1("fhNofRichRingsPerTS"), kLinear, kLog);
652  }
653 
654  {
655  TCanvas* c = fHM->CreateCanvas(
656  "rich_tb_reco_acc_rec", "rich_tb_reco_acc_rec", 1600, 600);
657  c->Divide(2, 1);
658  c->cd(1);
659  Int_t nofAccMom = fHM->H1("fhMomElAcc")->GetEntries();
660  Int_t nofRecMom = fHM->H1("fhMomElRec")->GetEntries();
661  Int_t nofRefAccMom = fHM->H1("fhMomRefElAcc")->GetEntries();
662  Int_t nofRefRecMom = fHM->H1("fhMomRefElRec")->GetEntries();
663  Double_t effMom = 100. * (Double_t) nofRecMom / (Double_t) nofAccMom;
664  Double_t effRefMom =
665  100. * (Double_t) nofRefRecMom / (Double_t) nofRefAccMom;
666  DrawH1({fHM->H1("fhMomElAcc"), fHM->H1("fhMomElRec")},
667  {"Acc (" + to_string(nofAccMom) + ")",
668  "Rec (" + to_string(nofRecMom) + ")"});
669  c->cd(2);
670  Int_t nofAccNofHits = fHM->H1("fhNofHitsElAcc")->GetEntries();
671  Int_t nofRecNofHits = fHM->H1("fhNofHitsElRec")->GetEntries();
672  Int_t nofRefAccNofHits = fHM->H1("fhNofHitsRefElAcc")->GetEntries();
673  Int_t nofRefRecNofHits = fHM->H1("fhNofHitsRefElRec")->GetEntries();
674  Double_t effNofHits =
675  100. * (Double_t) nofRecNofHits / (Double_t) nofAccNofHits;
676  Double_t effRefNofHits =
677  100. * (Double_t) nofRefRecNofHits / (Double_t) nofRefAccNofHits;
678  DrawH1({fHM->H1("fhNofHitsElAcc"), fHM->H1("fhNofHitsElRec")},
679  {"Acc (" + to_string(nofAccNofHits) + ")",
680  "Rec" + to_string(nofRecNofHits) + ")"});
681 
682  fHM->CreateCanvas("rich_tb_reco_eff_mom", "rich_tb_reco_eff_mom", 800, 800);
683  TH1D* hEffMom = Cbm::DivideH1(fHM->H1("fhMomElRec"), fHM->H1("fhMomElAcc"));
684  TH1D* hRefEffMom =
685  Cbm::DivideH1(fHM->H1("fhMomRefElRec"), fHM->H1("fhMomRefElAcc"));
686  DrawH1({hEffMom, hRefEffMom},
687  {"Electrons (" + Cbm::NumberToString(effMom, 2) + "%)",
688  "Reference electrons (" + Cbm::NumberToString(effRefMom, 2) + "%)"},
689  kLinear,
690  kLinear,
691  true,
692  0.5,
693  0.85,
694  0.99,
695  0.99);
696  hEffMom->SetMinimum(0.);
697  hEffMom->SetMaximum(105.);
698 
699  fHM->CreateCanvas(
700  "rich_tb_reco_eff_nofHits", "rich_tb_reco_eff_nofHits", 800, 800);
701  TH1D* hEffHits =
702  Cbm::DivideH1(fHM->H1("fhNofHitsElRec"), fHM->H1("fhNofHitsElAcc"));
703  TH1D* hRefEffHits =
704  Cbm::DivideH1(fHM->H1("fhNofHitsRefElRec"), fHM->H1("fhNofHitsRefElAcc"));
705  DrawH1(
706  {hEffHits, hRefEffHits},
707  {"Electron (" + Cbm::NumberToString(effNofHits, 2) + "%)",
708  "ElectronReference (" + Cbm::NumberToString(effRefNofHits, 2) + "%)"},
709  kLinear,
710  kLinear,
711  true,
712  0.5,
713  0.9,
714  0.99,
715  0.99);
716  hEffHits->SetMinimum(0.);
717  hEffHits->SetMaximum(105.);
718 
719  fHM->CreateCanvas(
720  "rich_tb_reco_eff_eventMult", "rich_tb_reco_eff_eventMult", 800, 800);
721  TH1D* hEffEventMult =
722  Cbm::DivideH1(fHM->H1("fhEventMultElRec"), fHM->H1("fhEventMultElAcc"));
723  DrawH1({hEffEventMult},
724  {"Electron (" + Cbm::NumberToString(effNofHits, 2) + "%)"},
725  kLinear,
726  kLinear,
727  true,
728  0.5,
729  0.9,
730  0.99,
731  0.99);
732  hEffHits->SetMinimum(0.);
733  hEffHits->SetMaximum(105.);
734  }
735 }
736 
737 void CbmRichRecoTbQa::DrawTimeLog(const string& hMainName,
738  Int_t nofLogEvents,
739  bool withNoise) {
740  Double_t max = std::numeric_limits<Double_t>::min();
741  Int_t startInd = (withNoise) ? -1 : 0;
742  for (int iEv = startInd; iEv < nofLogEvents; iEv++) {
743  string hName = hMainName + "_" + to_string(iEv);
744  string option =
745  ((iEv == 0 && !withNoise) || (iEv == -1 && withNoise)) ? "" : "same";
746  int ind = iEv % 6;
747  int color = (ind == 0)
748  ? kBlue
749  : (ind == 1)
750  ? kRed + 1
751  : (ind == 2)
752  ? kGreen + 3
753  : (ind == 3) ? kMagenta + 2
754  : (ind == 4) ? kYellow + 2 : kOrange + 8;
755  if (ind == -1) color = kAzure + 6;
756  max = std::max(max, fHM->H1(hName)->GetMaximum());
757  DrawH1(fHM->H1(hName), kLinear, kLog, option, color);
758  }
759  fHM->H1(hMainName + "_" + to_string(startInd))->SetMaximum(max * 1.10);
760  gPad->SetLeftMargin(0.07);
761  gPad->SetRightMargin(0.05);
762 }
763 
764 
765 void CbmRichRecoTbQa::Finish() {
766  ProcessMc();
767  DrawHist();
768  fHM->SaveCanvasToImage(fOutputDir);
769 
770  TDirectory* oldir = gDirectory;
771  TFile* outFile = FairRootManager::Instance()->GetOutFile();
772  if (outFile != NULL) {
773  outFile->cd();
774  fHM->WriteToFile();
775  }
776  gDirectory->cd(oldir->GetPath());
777 }
778 
779 
780 ClassImp(CbmRichRecoTbQa)
CbmRichPoint.h
CbmRichRing::GetTime
Double_t GetTime() const
Definition: CbmRichRing.h:102
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition: CbmMCDataManager.cxx:137
CbmMatch
Definition: CbmMatch.h:22
CbmMCDataManager.h
CbmRichRecoTbQa::InitHistograms
void InitHistograms()
Initialize histograms.
Definition: CbmRichRecoTbQa.cxx:110
CbmTrackMatchNew::GetTrueOverAllHitsRatio
Double_t GetTrueOverAllHitsRatio() const
Definition: CbmTrackMatchNew.h:35
CbmHistManager::WriteToFile
void WriteToFile()
Write all histograms to current opened file.
Definition: core/base/CbmHistManager.cxx:103
CbmRichDigi
Definition: CbmRichDigi.h:25
CbmRichRecoTbQa::fHM
CbmHistManager * fHM
Definition: CbmRichRecoTbQa.h:109
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition: CbmMCDataArray.cxx:133
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition: CbmDigiManager.cxx:71
CbmRichRecoTbQa::DrawHist
void DrawHist()
Definition: CbmRichRecoTbQa.cxx:524
CbmRichRecoTbQa::fMCTracks
CbmMCDataArray * fMCTracks
Definition: CbmRichRecoTbQa.h:116
CbmRichRecoTbQa::IsMcPion
Bool_t IsMcPion(const CbmMCTrack *mctrack)
Definition: CbmRichRecoTbQa.cxx:500
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRecoTbQa::GetNofPrimaryMcTracks
Int_t GetNofPrimaryMcTracks(Int_t iEv)
Definition: CbmRichRecoTbQa.cxx:361
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmRichRecoTbQa::Process
void Process()
Definition: CbmRichRecoTbQa.cxx:256
CbmRichRecoTbQa::fRichHits
TClonesArray * fRichHits
Definition: CbmRichRecoTbQa.h:120
CbmMCDataArray.h
CbmGlobalTrack.h
CbmDigiManager::GetNofDigis
static Int_t GetNofDigis(ECbmModuleId systemId)
Definition: CbmDigiManager.cxx:62
CbmRichRing
Definition: CbmRichRing.h:17
CbmRichRecoTbQa::RingRecoEfficiency
void RingRecoEfficiency()
Definition: CbmRichRecoTbQa.cxx:373
CbmRichRing.h
CbmStsPoint
Definition: CbmStsPoint.h:27
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichDraw.h
CbmRichRecoTbQa::fNofLogEvents
Int_t fNofLogEvents
Definition: CbmRichRecoTbQa.h:113
CbmRichRecoTbQa::DrawTimeLog
void DrawTimeLog(const string &hMainName, Int_t nofLogEvents, bool withNoise=false)
Definition: CbmRichRecoTbQa.cxx:737
CbmHistManager.h
Histogram manager.
CbmMatchRecoToMC.h
FairTask for matching RECO data to MC.
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition: CbmDigiManager.h:93
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:33
CbmRichGeoManager.h
Cbm::DivideH1
TH1D * DivideH1(TH1 *h1, TH1 *h2, const string &histName, double scale, const string &titleYaxis)
Definition: CbmUtils.cxx:70
CbmRichDigi.h
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
CbmMCEventList::GetNofEvents
std::size_t GetNofEvents() const
Number of events in the list @value Number of events.
Definition: CbmMCEventList.h:90
CbmRichRecoTbQa::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichRecoTbQa.cxx:244
CbmRichRecoTbQa::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichRecoTbQa.cxx:67
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmRichRecoTbQa::ProcessMc
void ProcessMc()
Definition: CbmRichRecoTbQa.cxx:305
CbmRichDigi::GetTime
Double_t GetTime() const
Definition: CbmRichDigi.h:61
CbmRichRecoTbQa::fTimeSliceNum
Int_t fTimeSliceNum
Definition: CbmRichRecoTbQa.h:111
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmRichRecoTbQa::IsCherenkovPhotonFromSecondaryElectron
Bool_t IsCherenkovPhotonFromSecondaryElectron(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition: CbmRichRecoTbQa.cxx:471
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
kLinear
@ kLinear
Definition: CbmDrawHist.h:79
CbmTrackMatchNew.h
CbmMCEventList::GetEventTime
Double_t GetEventTime(UInt_t event, UInt_t file)
Event start time.
Definition: CbmMCEventList.cxx:86
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
CbmDigiManager::GetMatch
const CbmMatch * GetMatch(ECbmModuleId systemId, UInt_t index) const
Get a match object.
Definition: CbmDigiManager.cxx:54
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
ECbmModuleId::kRich
@ kRich
Ring-Imaging Cherenkov Detector.
CbmRichRecoTbQa::fEventList
CbmMCEventList * fEventList
Definition: CbmRichRecoTbQa.h:123
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmMCTrack::GetGeantProcessId
UInt_t GetGeantProcessId() const
Definition: CbmMCTrack.h:69
CbmMCEventList
Container class for MC events with number, file and start time.
Definition: CbmMCEventList.h:38
CbmRichRecoTbQa::fRichRings
TClonesArray * fRichRings
Definition: CbmRichRecoTbQa.h:121
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition: CbmHistManager.h:170
CbmRichRecoTbQa::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichRecoTbQa.cxx:765
fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmTofAnaTestbeam.cxx:88
CbmRichRecoTbQa::fStsPoints
CbmMCDataArray * fStsPoints
Definition: CbmRichRecoTbQa.h:118
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
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
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
CbmMCEventList.h
counter
int counter
Definition: CbmMvdSensorDigiToHitTask.cxx:72
CbmStsPoint.h
CbmMCTrack.h
CbmRichRecoTbQa::IsCherenkovPhotonFromPion
Bool_t IsCherenkovPhotonFromPion(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition: CbmRichRecoTbQa.cxx:507
CbmDigiManager.h
CbmRichRecoTbQa::fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmRichRecoTbQa.h:119
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmRichRecoTbQa
Definition: CbmRichRecoTbQa.h:19
CbmRichRecoTbQa::CbmRichRecoTbQa
CbmRichRecoTbQa()
Standard constructor.
Definition: CbmRichRecoTbQa.cxx:50
CbmRichRecoTbQa::IsCherenkovPhoton
Bool_t IsCherenkovPhoton(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition: CbmRichRecoTbQa.cxx:443
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
SetDefaultDrawStyle
void SetDefaultDrawStyle()
Definition: CbmDrawHist.cxx:33
CbmLitGlobalElectronId.h
CbmRichHit.h
CbmRichRecoTbQa::fRecRings
vector< pair< Int_t, Int_t > > fRecRings
Definition: CbmRichRecoTbQa.h:125
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmRichRecoTbQa::IsMcPrimaryElectron
Bool_t IsMcPrimaryElectron(const CbmMCTrack *mctrack)
Definition: CbmRichRecoTbQa.cxx:493
CbmRichRecoTbQa::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmRichRecoTbQa.h:122
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
Cbm::NumberToString
std::string NumberToString(const T &value, int precision=1)
Definition: CbmUtils.h:23
CbmRichRecoTbQa.h
CbmRichRecoTbQa::IsCherenkovPhotonFromPrimaryElectron
Bool_t IsCherenkovPhotonFromPrimaryElectron(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition: CbmRichRecoTbQa.cxx:454
CbmMCTrack::GetP
Double_t GetP() const
Definition: CbmMCTrack.h:100
kLog
@ kLog
Definition: CbmDrawHist.h:78
CbmRichRecoTbQa::fRichPoints
CbmMCDataArray * fRichPoints
Definition: CbmRichRecoTbQa.h:117
CbmRichHit
Definition: CbmRichHit.h:19
CbmRichRecoTbQa::fOutputDir
string fOutputDir
Definition: CbmRichRecoTbQa.h:114