CbmRoot
CbmRichUrqmdTest.cxx
Go to the documentation of this file.
1 
8 #include "CbmRichUrqmdTest.h"
9 
10 #include "TCanvas.h"
11 #include "TClonesArray.h"
12 #include "TH1.h"
13 #include "TH1D.h"
14 
15 #include "CbmDigiManager.h"
16 #include "CbmDrawHist.h"
17 #include "CbmMCDataArray.h"
18 #include "CbmMCDataManager.h"
19 #include "CbmMCEventList.h"
20 #include "CbmMCTrack.h"
21 #include "CbmMatchRecoToMC.h"
22 #include "CbmRichDigi.h"
23 #include "CbmRichDraw.h"
24 #include "CbmRichGeoManager.h"
25 #include "CbmRichHit.h"
26 #include "CbmRichPoint.h"
27 #include "CbmRichRing.h"
28 #include "CbmTrackMatchNew.h"
29 #include "FairMCPoint.h"
30 #include "FairTrackParam.h"
31 #include "TStyle.h"
32 #include <TFile.h>
33 
34 #include "CbmHistManager.h"
35 #include "CbmUtils.h"
36 
37 #include <boost/assign/list_of.hpp>
38 #include <iostream>
39 #include <sstream>
40 #include <string>
41 
42 using namespace std;
43 using boost::assign::list_of;
44 
45 CbmRichUrqmdTest::CbmRichUrqmdTest()
46  : FairTask("CbmRichUrqmdTest")
47  , fHM(NULL)
48  , fOutputDir("")
49  , fRichHits(NULL)
50  , fRichRings(NULL)
51  , fRichPoints(NULL)
52  , fMcTracks(NULL)
53  , fRichRingMatches(NULL)
54  , fRichProjections(NULL)
55  , fDigiMan(nullptr)
56  , fEventList(NULL)
57  , fEventNum(0)
58  , fMinNofHits(7)
59  , fNofHitsInRingMap() {}
60 
61 CbmRichUrqmdTest::~CbmRichUrqmdTest() {}
62 
63 InitStatus CbmRichUrqmdTest::Init() {
64  cout << "CbmRichUrqmdTest::Init" << endl;
65  FairRootManager* ioman = FairRootManager::Instance();
66  if (NULL == ioman) {
67  Fatal("CbmRichUrqmdTest::Init", "RootManager not instantised!");
68  }
69 
70  CbmMCDataManager* mcManager =
71  (CbmMCDataManager*) ioman->GetObject("MCDataManager");
72  if (mcManager == nullptr)
73  LOG(fatal) << "CbmRichUrqmdTest::Init() NULL MCDataManager.";
74 
75  fMcTracks = mcManager->InitBranch("MCTrack");
76  if (NULL == fMcTracks) { LOG(fatal) << "CbmRichUrqmdTest::Init No MCTrack!"; }
77 
78  fRichPoints = mcManager->InitBranch("RichPoint");
79  if (NULL == fRichPoints) {
80  LOG(fatal) << "CbmRichUrqmdTest::Init No RichPoint!";
81  }
82 
83  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
84  if (NULL == fRichHits) { LOG(fatal) << "CbmRichUrqmdTest::Init No RichHit!"; }
85 
86  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
87  if (NULL == fRichRings) {
88  LOG(fatal) << "CbmRichUrqmdTest::Init No RichRing!";
89  }
90 
91  fDigiMan = CbmDigiManager::Instance();
92  fDigiMan->Init();
93 
94  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
95  if (NULL == fRichRingMatches) {
96  LOG(fatal) << "CbmRichUrqmdTest::Init No RichRingMatch!";
97  }
98 
99  fRichProjections = (TClonesArray*) ioman->GetObject("RichProjection");
100  if (NULL == fRichProjections) {
101  LOG(fatal) << "CbmRichUrqmdTest::Init No fRichProjection !";
102  }
103 
104  fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList.");
105  if (NULL == fEventList) {
106  LOG(fatal) << "CbmRichUrqmdTest::Init No MCEventList!";
107  }
108 
109  InitHistograms();
110 
111  return kSUCCESS;
112 }
113 
114 void CbmRichUrqmdTest::Exec(Option_t* /*option*/) {
115  fEventNum++;
116 
117  cout << "CbmRichUrqmdTest, event No. " << fEventNum << endl;
118 
119  FillRichRingNofHits();
120  NofRings();
121  NofHitsAndPoints();
122  NofProjections();
123  Vertex();
124  PmtXYSource();
125 }
126 
127 void CbmRichUrqmdTest::InitHistograms() {
128  fHM = new CbmHistManager();
129 
130  fHM->Create1<TH1D>(
131  "fh_vertex_z", "fh_vertex_z;z [cm];# vertices per event", 350, -1., 350);
132  fHM->Create1<TH1D>("fh_vertex_z_sts",
133  "fh_vertex_z_sts;z [cm];# vertices per event",
134  222,
135  -1.,
136  110.);
137  fHM->Create2<TH2D>("fh_vertex_xy",
138  "fh_vertex_xy;x [cm];y [cm];# vertices per event",
139  100,
140  -200.,
141  200.,
142  100,
143  -200.,
144  200.);
145  fHM->Create2<TH2D>("fh_vertex_zy",
146  "fh_vertex_zy;z [cm];y [cm];# vertices per event",
147  350,
148  -1.,
149  350,
150  100,
151  -200.,
152  200.);
153  fHM->Create2<TH2D>("fh_vertex_zx",
154  "fh_vertex_zx;z [cm];x [cm];# vertices per event",
155  350,
156  -1.,
157  350,
158  100,
159  -200.,
160  200.);
161  fHM->Create2<TH2D>("fh_vertex_xy_z100_180",
162  "fh_vertex_xy_z100_180;x [cm];y [cm];# vertices per event",
163  100,
164  -200.,
165  200.,
166  100,
167  -200.,
168  200.);
169  fHM->Create2<TH2D>("fh_vertex_xy_z180_370",
170  "fh_vertex_xy_z180_370;x [cm];y [cm];# vertices per event",
171  100,
172  -200.,
173  200.,
174  100,
175  -200.,
176  200.);
177  fHM->Create2<TH2D>("fh_vertex_xy_z180_230",
178  "fh_vertex_xy_z180_230;x [cm];y [cm];# vertices per event",
179  100,
180  -200.,
181  200.,
182  100,
183  -200.,
184  200.);
185 
186  fHM->Create2<TH2D>("fh_vertex_xy_z5",
187  "fh_vertex_xy_z5;x [cm];y [cm];# vertices per event",
188  100,
189  -100.,
190  100.,
191  100,
192  -100.,
193  100.);
194  fHM->Create2<TH2D>("fh_vertex_xy_z5_15",
195  "fh_vertex_xy_z5_15;x [cm];y [cm];# vertices per event",
196  100,
197  -100.,
198  100.,
199  100,
200  -100.,
201  100.);
202  fHM->Create2<TH2D>("fh_vertex_xy_z15_25",
203  "fh_vertex_xy_z15_25;x [cm];y [cm];# vertices per event",
204  100,
205  -100.,
206  100.,
207  100,
208  -100.,
209  100.);
210  fHM->Create2<TH2D>("fh_vertex_xy_z25_35",
211  "fh_vertex_xy_z25_35;x [cm];y [cm];# vertices per event",
212  100,
213  -100.,
214  100.,
215  100,
216  -100.,
217  100.);
218  fHM->Create2<TH2D>("fh_vertex_xy_z35_45",
219  "fh_vertex_xy_z35_45;x [cm];y [cm];# vertices per event",
220  100,
221  -100.,
222  100.,
223  100,
224  -100.,
225  100.);
226  fHM->Create2<TH2D>("fh_vertex_xy_z45_55",
227  "fh_vertex_xy_z45_55;x [cm];y [cm];# vertices per event",
228  100,
229  -100.,
230  100.,
231  100,
232  -100.,
233  100.);
234  fHM->Create2<TH2D>("fh_vertex_xy_z55_65",
235  "fh_vertex_xy_z55_65;x [cm];y [cm];# vertices per event",
236  100,
237  -100.,
238  100.,
239  100,
240  -100.,
241  100.);
242  fHM->Create2<TH2D>("fh_vertex_xy_z65_75",
243  "fh_vertex_xy_z65_75;x [cm];y [cm];# vertices per event",
244  100,
245  -100.,
246  100.,
247  100,
248  -100.,
249  100.);
250  fHM->Create2<TH2D>("fh_vertex_xy_z75_85",
251  "fh_vertex_xy_z75_85;x [cm];y [cm];# vertices per event",
252  100,
253  -100.,
254  100.,
255  100,
256  -100.,
257  100.);
258  fHM->Create2<TH2D>("fh_vertex_xy_z85_95",
259  "fh_vertex_xy_z85_95;x [cm];y [cm];# vertices per event",
260  100,
261  -100.,
262  100.,
263  100,
264  -100.,
265  100.);
266  fHM->Create2<TH2D>("fh_vertex_xy_z95_105",
267  "fh_vertex_xy_z95_105;x [cm];y [cm];# vertices per event",
268  100,
269  -100.,
270  100.,
271  100,
272  -100.,
273  100.);
274 
275  fHM->Create1<TH1D>("fh_nof_rings_1hit",
276  "fh_nof_rings_1hit;# detected particles/event;Yield",
277  250,
278  -.5,
279  249.5);
280  fHM->Create1<TH1D>("fh_nof_rings_7hits",
281  "fh_nof_rings_7hits;# detected particles/event;Yield",
282  250,
283  -.5,
284  249.5);
285  fHM->Create1<TH1D>("fh_nof_rings_prim_1hit",
286  "fh_nof_rings_prim_1hit;# detected particles/event;Yield",
287  50,
288  -.5,
289  69.5);
290  fHM->Create1<TH1D>("fh_nof_rings_prim_7hits",
291  "fh_nof_rings_prim_7hits;# detected particles/event;Yield",
292  50,
293  -.5,
294  69.5);
295  fHM->Create1<TH1D>(
296  "fh_nof_rings_target_1hit",
297  "fh_nof_rings_target_1hit;# detected particles/event;Yield",
298  60,
299  -.5,
300  79.5);
301  fHM->Create1<TH1D>(
302  "fh_nof_rings_target_7hits",
303  "fh_nof_rings_target_7hits;# detected particles/event;Yield",
304  60,
305  -.5,
306  79.5);
307 
308  fHM->Create1<TH1D>(
309  "fh_secel_mom", "fh_secel_mom;p [GeV/c];Number per event", 100, 0., 20);
310  fHM->Create1<TH1D>("fh_gamma_target_mom",
311  "fh_gamma_target_mom;p [GeV/c];Number per event",
312  100,
313  0.,
314  20);
315  fHM->Create1<TH1D>("fh_gamma_nontarget_mom",
316  "fh_gamma_nontarget_mom;p [GeV/c];Number per event",
317  100,
318  0.,
319  20);
320  fHM->Create1<TH1D>(
321  "fh_pi_mom", "fh_pi_mom;p [GeV/c];Number per event", 100, 0., 20);
322  fHM->Create1<TH1D>(
323  "fh_kaon_mom", "fh_kaon_mom;p [GeV/c];Number per event", 100, 0., 20);
324  fHM->Create1<TH1D>(
325  "fh_mu_mom", "fh_mu_mom;p [GeV/c];Number per event", 100, 0., 20);
326 
327  fHM->Create1<TH1D>("fh_nof_points_per_event",
328  "fh_nof_points_per_event;Particle;# MC points per event",
329  7,
330  .5,
331  7.5);
332  fHM->Create1<TH1D>("fh_nof_hits_per_event",
333  "fh_nof_hits_per_event;# hits per event;Yield",
334  100,
335  0,
336  4000);
337  fHM->Create1<TH1D>("fh_nof_hits_per_pmt",
338  "fh_nof_hits_per_pmt;# hits per PMT;% of total",
339  65,
340  -0.5,
341  64.5);
342 
343  vector<Double_t> xPmtBins = CbmRichDraw::GetPmtHistXbins();
344  vector<Double_t> yPmtBins = CbmRichDraw::GetPmtHistYbins();
345 
346  // before drawing must be normalized by 1/64
347  TH2D* fh_hitrate_xy = new TH2D("fh_hitrate_xy",
348  "fh_hitrate_xy;X [cm];Y [cm];# hits/pixel/s",
349  xPmtBins.size() - 1,
350  &xPmtBins[0],
351  yPmtBins.size() - 1,
352  &yPmtBins[0]);
353  fHM->Add("fh_hitrate_xy", fh_hitrate_xy);
354 
355  TH2D* fh_hits_xy = new TH2D("fh_hits_xy",
356  "fh_hits_xy;X [cm];Y [cm];# hits/PMT/event",
357  xPmtBins.size() - 1,
358  &xPmtBins[0],
359  yPmtBins.size() - 1,
360  &yPmtBins[0]);
361  fHM->Add("fh_hits_xy", fh_hits_xy);
362 
363  TH2D* fh_points_xy =
364  new TH2D("fh_points_xy",
365  "fh_points_xy;X [cm];Y [cm];# MC points/PMT/event",
366  xPmtBins.size() - 1,
367  &xPmtBins[0],
368  yPmtBins.size() - 1,
369  &yPmtBins[0]);
370  fHM->Add("fh_points_xy", fh_points_xy);
371 
372  TH2D* fh_points_xy_pions =
373  new TH2D("fh_points_xy_pions",
374  "fh_points_xy_pions;X [cm];Y [cm];# MC points/PMT/event",
375  xPmtBins.size() - 1,
376  &xPmtBins[0],
377  yPmtBins.size() - 1,
378  &yPmtBins[0]);
379  fHM->Add("fh_points_xy_pions", fh_points_xy_pions);
380 
381  TH2D* fh_points_xy_gamma_target =
382  new TH2D("fh_points_xy_gamma_target",
383  "fh_points_xy_gamma_target;X [cm];Y [cm];# MC points/PMT/event",
384  xPmtBins.size() - 1,
385  &xPmtBins[0],
386  yPmtBins.size() - 1,
387  &yPmtBins[0]);
388  fHM->Add("fh_points_xy_gamma_target", fh_points_xy_gamma_target);
389 
390  TH2D* fh_points_xy_gamma_nontarget =
391  new TH2D("fh_points_xy_gamma_nontarget",
392  "fh_points_xy_gamma_nontarget;X [cm];Y [cm];# MC points/PMT/event",
393  xPmtBins.size() - 1,
394  &xPmtBins[0],
395  yPmtBins.size() - 1,
396  &yPmtBins[0]);
397  fHM->Add("fh_points_xy_gamma_nontarget", fh_points_xy_gamma_nontarget);
398 
399  TH2D* fh_skipped_pmt_10_xy =
400  new TH2D("fh_skipped_pmt_10_xy",
401  "fh_skipped_pmt_10_xy;X [cm];Y [cm];# skipped PMTs (>10 hits) [%]",
402  xPmtBins.size() - 1,
403  &xPmtBins[0],
404  yPmtBins.size() - 1,
405  &yPmtBins[0]);
406  fHM->Add("fh_skipped_pmt_10_xy", fh_skipped_pmt_10_xy);
407 
408  TH2D* fh_skipped_pmt_20_xy =
409  new TH2D("fh_skipped_pmt_20_xy",
410  "fh_skipped_pmt_20_xy;X [cm];Y [cm];# skipped PMTs (>20 hits) [%]",
411  xPmtBins.size() - 1,
412  &xPmtBins[0],
413  yPmtBins.size() - 1,
414  &yPmtBins[0]);
415  fHM->Add("fh_skipped_pmt_20_xy", fh_skipped_pmt_20_xy);
416 
417  TH2D* fh_skipped_pmt_30_xy =
418  new TH2D("fh_skipped_pmt_30_xy",
419  "fh_skipped_pmt_30_xy;X [cm];Y [cm];# skipped PMTs (>30 hits) [%]",
420  xPmtBins.size() - 1,
421  &xPmtBins[0],
422  yPmtBins.size() - 1,
423  &yPmtBins[0]);
424  fHM->Add("fh_skipped_pmt_30_xy", fh_skipped_pmt_30_xy);
425 
426  fHM->Create1<TH1D>("fh_nof_proj_per_event",
427  "fh_nof_proj_per_event;# tracks per event;Yield",
428  50,
429  0,
430  1000);
431  fHM->Create2<TH2D>("fh_proj_xy",
432  "fh_proj_xy;X [cm];Y [cm];# tracks/cm^{2}/event",
433  240,
434  -120,
435  120,
436  420,
437  -210,
438  210);
439 }
440 
441 void CbmRichUrqmdTest::FillRichRingNofHits() {
442  fNofHitsInRingMap.clear();
443  Int_t nofRichHits = fRichHits->GetEntriesFast();
444  for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
445  CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
446  if (NULL == hit) continue;
447  Int_t digiIndex = hit->GetRefId();
448  if (digiIndex < 0) continue;
449  const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(digiIndex);
450  if (NULL == digi) continue;
451  const CbmMatch* digiMatch =
452  fDigiMan->GetMatch(ECbmModuleId::kRich, digiIndex);
453  if (NULL == digiMatch) continue;
454  Int_t eventId = digiMatch->GetMatchedLink().GetEntry();
455 
456  vector<pair<Int_t, Int_t>> motherIds =
457  CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(
458  fDigiMan, hit, fRichPoints, fMcTracks, eventId);
459  for (UInt_t i = 0; i < motherIds.size(); i++) {
460  fNofHitsInRingMap[motherIds[i]]++;
461  }
462  }
463 }
464 
465 void CbmRichUrqmdTest::NofRings() {
466  Int_t fileId = 0;
467  Int_t nofRings = fRichRings->GetEntriesFast();
468  int nRings1hit = 0, nRings7hits = 0;
469  int nRingsPrim1hit = 0, nRingsPrim7hits = 0;
470  int nRingsTarget1hit = 0, nRingsTarget7hits = 0;
471  for (Int_t iR = 0; iR < nofRings; iR++) {
472  const CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
473  if (NULL == ring) continue;
474  const CbmTrackMatchNew* ringMatch =
475  static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iR));
476  if (NULL == ringMatch) continue;
477 
478  Int_t mcEventId = ringMatch->GetMatchedLink().GetEntry();
479  Int_t mcTrackId = ringMatch->GetMatchedLink().GetIndex();
480  const CbmMCTrack* mcTrack =
481  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, mcEventId, mcTrackId));
482  if (NULL == mcTrack) continue;
483  Int_t motherId = mcTrack->GetMotherId();
484  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
485  double mom = mcTrack->GetP();
486  TVector3 vert;
487  mcTrack->GetStartVertex(vert);
488  double dZ = vert.Z();
489 
490  if (motherId == -1 && pdg == 11)
491  continue; // do not calculate embedded electrons
492 
493  int nofHits = ring->GetNofHits();
494  if (nofHits >= 1) nRings1hit++;
495  if (nofHits >= fMinNofHits) nRings7hits++;
496 
497  if (motherId == -1 && nofHits >= 1) nRingsPrim1hit++;
498  if (motherId == -1 && nofHits >= fMinNofHits) nRingsPrim7hits++;
499 
500  if (dZ < 0.1 && nofHits >= 1) nRingsTarget1hit++;
501  if (dZ < 0.1 && nofHits >= fMinNofHits) nRingsTarget7hits++;
502 
503  if (nofHits >= 1) {
504  if (motherId != -1) {
505  int motherPdg = -1;
506  const CbmMCTrack* mother =
507  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, mcEventId, motherId));
508  if (NULL != mother) motherPdg = mother->GetPdgCode();
509  if (motherId != -1 && pdg == 11 && motherPdg != 22)
510  fHM->H1("fh_secel_mom")->Fill(mom);
511 
512  if (motherId != -1 && pdg == 11 && motherPdg == 22) {
513  if (dZ < 0.1) {
514  fHM->H1("fh_gamma_target_mom")->Fill(mom);
515  } else {
516  fHM->H1("fh_gamma_nontarget_mom")->Fill(mom);
517  }
518  }
519  }
520  if (pdg == 211) fHM->H1("fh_pi_mom")->Fill(mom);
521  if (pdg == 321) fHM->H1("fh_kaon_mom")->Fill(mom);
522  if (pdg == 13) fHM->H1("fh_mu_mom")->Fill(mom);
523  }
524  }
525  fHM->H1("fh_nof_rings_1hit")->Fill(nRings1hit);
526  fHM->H1("fh_nof_rings_7hits")->Fill(nRings7hits);
527 
528  fHM->H1("fh_nof_rings_prim_1hit")->Fill(nRingsPrim1hit);
529  fHM->H1("fh_nof_rings_prim_7hits")->Fill(nRingsPrim7hits);
530 
531  fHM->H1("fh_nof_rings_target_1hit")->Fill(nRingsTarget1hit);
532  fHM->H1("fh_nof_rings_target_7hits")->Fill(nRingsTarget7hits);
533 }
534 
535 void CbmRichUrqmdTest::NofHitsAndPoints() {
536  int nofHits = fRichHits->GetEntriesFast();
537  fHM->H1("fh_nof_hits_per_event")->Fill(nofHits);
538  map<Int_t, Int_t> digisPerPmtMap;
539  for (int iH = 0; iH < nofHits; iH++) {
540  CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iH));
541  if (nullptr == hit) continue;
542  fHM->H2("fh_hits_xy")->Fill(hit->GetX(), hit->GetY());
543  fHM->H2("fh_hitrate_xy")->Fill(hit->GetX(), hit->GetY());
544 
545  Int_t digiInd = hit->GetRefId();
546  const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(digiInd);
547  if (nullptr == richDigi) continue;
548  CbmRichPixelData* pixelData =
549  CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(
550  richDigi->GetAddress());
551  if (nullptr == pixelData) continue;
552  Int_t pmtId = pixelData->fPmtId;
553  digisPerPmtMap[pmtId]++;
554  }
555 
556  for (auto const& it : digisPerPmtMap) {
557  Int_t pmtId = it.first;
558  Int_t nofDigis = it.second;
559  CbmRichPmtData* pmtData =
560  CbmRichDigiMapManager::GetInstance().GetPmtDataById(pmtId);
561  TVector3 inPos(pmtData->fX, pmtData->fY, pmtData->fZ);
562  TVector3 outPos;
563  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
564  if (nofDigis > 10)
565  fHM->H2("fh_skipped_pmt_10_xy")->Fill(outPos.X(), outPos.Y());
566  if (nofDigis > 20)
567  fHM->H2("fh_skipped_pmt_20_xy")->Fill(outPos.X(), outPos.Y());
568  if (nofDigis > 30)
569  fHM->H2("fh_skipped_pmt_30_xy")->Fill(outPos.X(), outPos.Y());
570  }
571 
572  vector<Int_t> allPmtIds = CbmRichDigiMapManager::GetInstance().GetPmtIds();
573  for (Int_t pmtId : allPmtIds) {
574  fHM->H1("fh_nof_hits_per_pmt")->Fill(digisPerPmtMap[pmtId]);
575  }
576 
577  Int_t nofEvents = fEventList->GetNofEvents();
578  for (Int_t iE = 0; iE < nofEvents; iE++) {
579  Int_t fileId = fEventList->GetFileIdByIndex(iE);
580  Int_t eventId = fEventList->GetEventIdByIndex(iE);
581  int nofPoints = fRichPoints->Size(fileId, eventId);
582  for (int i = 0; i < nofPoints; i++) {
583  const CbmRichPoint* point =
584  static_cast<CbmRichPoint*>(fRichPoints->Get(fileId, eventId, i));
585  if (NULL == point) continue;
586  fHM->H1("fh_nof_points_per_event")->Fill(1);
587 
588  Int_t mcPhotonTrackId = point->GetTrackID();
589  if (mcPhotonTrackId < 0) continue;
590  const CbmMCTrack* mcPhotonTrack = static_cast<CbmMCTrack*>(
591  fMcTracks->Get(fileId, eventId, mcPhotonTrackId));
592  if (NULL == mcPhotonTrack) continue;
593  Int_t motherPhotonId = mcPhotonTrack->GetMotherId();
594  if (motherPhotonId < 0) continue;
595  const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(
596  fMcTracks->Get(fileId, eventId, motherPhotonId));
597  if (NULL == mcTrack) continue;
598  Int_t motherId = mcTrack->GetMotherId();
599 
600  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
601  TVector3 vert;
602  mcTrack->GetStartVertex(vert);
603  double dZ = vert.Z();
604 
605  if (motherId == -1 && pdg == 11)
606  continue; // do not calculate embedded electrons
607 
608  if (motherId != -1) {
609  int motherPdg = -1;
610  const CbmMCTrack* mother =
611  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, motherId));
612  if (NULL != mother) motherPdg = mother->GetPdgCode();
613  if (motherId != -1 && pdg == 11 && motherPdg != 22)
614  fHM->H1("fh_nof_points_per_event")->Fill(2);
615 
616  if (motherId != -1 && pdg == 11 && motherPdg == 22) {
617  if (dZ < 0.1) {
618  fHM->H1("fh_nof_points_per_event")->Fill(3);
619  } else {
620  fHM->H1("fh_nof_points_per_event")->Fill(4);
621  }
622  }
623  }
624  if (pdg == 211) fHM->H1("fh_nof_points_per_event")->Fill(5);
625  if (pdg == 321) fHM->H1("fh_nof_points_per_event")->Fill(6);
626  if (pdg == 13) fHM->H1("fh_nof_points_per_event")->Fill(7);
627  }
628  }
629 }
630 
631 void CbmRichUrqmdTest::PmtXYSource() {
632  Int_t nofEvents = fEventList->GetNofEvents();
633  for (Int_t iE = 0; iE < nofEvents; iE++) {
634  Int_t fileId = fEventList->GetFileIdByIndex(iE);
635  Int_t eventId = fEventList->GetEventIdByIndex(iE);
636  int nofPoints = fRichPoints->Size(fileId, eventId);
637  for (int i = 0; i < nofPoints; i++) {
638  const CbmRichPoint* point =
639  static_cast<CbmRichPoint*>(fRichPoints->Get(fileId, eventId, i));
640  if (NULL == point) continue;
641 
642  Int_t iMCTrack = point->GetTrackID();
643  const CbmMCTrack* track =
644  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iMCTrack));
645  if (NULL == track) continue;
646 
647  Int_t iMother = track->GetMotherId();
648  if (iMother == -1) continue;
649 
650  const CbmMCTrack* track2 =
651  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iMother));
652  if (NULL == track2) continue;
653  int pdg = TMath::Abs(track2->GetPdgCode());
654  int motherId = track2->GetMotherId();
655  TVector3 inPos(point->GetX(), point->GetY(), point->GetZ());
656  TVector3 outPos;
657  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
658 
659  fHM->H2("fh_points_xy")->Fill(outPos.X(), outPos.Y());
660  if (motherId != -1) {
661  int motherPdg = -1;
662  const CbmMCTrack* mother =
663  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, motherId));
664  if (NULL != mother) motherPdg = mother->GetPdgCode();
665  TVector3 vert;
666  track2->GetStartVertex(vert);
667  if (motherId != -1 && pdg == 11 && motherPdg == 22) {
668  if (vert.Z() < 0.1) {
669  fHM->H2("fh_points_xy_gamma_target")->Fill(outPos.X(), outPos.Y());
670  } else {
671  fHM->H2("fh_points_xy_gamma_nontarget")
672  ->Fill(outPos.X(), outPos.Y());
673  }
674  }
675  }
676  if (pdg == 211)
677  fHM->H2("fh_points_xy_pions")->Fill(outPos.X(), outPos.Y());
678  }
679  }
680 }
681 
682 void CbmRichUrqmdTest::NofProjections() {
683  if (fRichProjections == NULL) return;
684  int nofProj = fRichProjections->GetEntriesFast();
685  int nofGoodProj = 0;
686  for (int i = 0; i < nofProj; i++) {
687  FairTrackParam* proj = (FairTrackParam*) fRichProjections->At(i);
688  if (NULL == proj) continue;
689  fHM->H2("fh_proj_xy")->Fill(proj->GetX(), proj->GetY());
690  if (proj->GetX() != 0 && proj->GetY() != 0) nofGoodProj++;
691  }
692  fHM->H1("fh_nof_proj_per_event")->Fill(nofGoodProj);
693 }
694 
695 void CbmRichUrqmdTest::Vertex() {
696  Int_t nofEvents = fEventList->GetNofEvents();
697  for (Int_t iE = 0; iE < nofEvents; iE++) {
698  Int_t fileId = fEventList->GetFileIdByIndex(iE);
699  Int_t eventId = fEventList->GetEventIdByIndex(iE);
700 
701  int nMcTracks = fMcTracks->Size(fileId, eventId);
702  for (int i = 0; i < nMcTracks; i++) {
703  //At least one hit in RICH
704  pair<Int_t, Int_t> val = std::make_pair(eventId, i);
705  if (fNofHitsInRingMap[val] < 1) continue;
706  const CbmMCTrack* mctrack =
707  static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, i));
708  TVector3 v;
709  mctrack->GetStartVertex(v);
710  fHM->H1("fh_vertex_z")->Fill(v.Z());
711  fHM->H1("fh_vertex_z_sts")->Fill(v.Z());
712  fHM->H2("fh_vertex_xy")->Fill(v.X(), v.Y());
713  fHM->H2("fh_vertex_zy")->Fill(v.Z(), v.Y());
714  fHM->H2("fh_vertex_zx")->Fill(v.Z(), v.X());
715  if (v.Z() >= 100 && v.Z() <= 180)
716  fHM->H2("fh_vertex_xy_z100_180")->Fill(v.X(), v.Y());
717  if (v.Z() >= 180 && v.Z() <= 370)
718  fHM->H2("fh_vertex_xy_z180_370")->Fill(v.X(), v.Y());
719  if (v.Z() >= 180 && v.Z() <= 230)
720  fHM->H2("fh_vertex_xy_z180_230")->Fill(v.X(), v.Y());
721 
722  if (v.Z() <= 5) fHM->H2("fh_vertex_xy_z5")->Fill(v.X(), v.Y());
723  if (v.Z() > 5 && v.Z() <= 15)
724  fHM->H2("fh_vertex_xy_z5_15")->Fill(v.X(), v.Y());
725  if (v.Z() > 15 && v.Z() <= 25)
726  fHM->H2("fh_vertex_xy_z15_25")->Fill(v.X(), v.Y());
727  if (v.Z() > 25 && v.Z() <= 35)
728  fHM->H2("fh_vertex_xy_z25_35")->Fill(v.X(), v.Y());
729  if (v.Z() > 35 && v.Z() <= 45)
730  fHM->H2("fh_vertex_xy_z35_45")->Fill(v.X(), v.Y());
731  if (v.Z() > 45 && v.Z() <= 55)
732  fHM->H2("fh_vertex_xy_z45_55")->Fill(v.X(), v.Y());
733  if (v.Z() > 55 && v.Z() <= 65)
734  fHM->H2("fh_vertex_xy_z55_65")->Fill(v.X(), v.Y());
735  if (v.Z() > 65 && v.Z() <= 75)
736  fHM->H2("fh_vertex_xy_z65_75")->Fill(v.X(), v.Y());
737  if (v.Z() > 75 && v.Z() <= 85)
738  fHM->H2("fh_vertex_xy_z75_85")->Fill(v.X(), v.Y());
739  if (v.Z() > 85 && v.Z() <= 95)
740  fHM->H2("fh_vertex_xy_z85_95")->Fill(v.X(), v.Y());
741  if (v.Z() > 95 && v.Z() <= 105)
742  fHM->H2("fh_vertex_xy_z95_105")->Fill(v.X(), v.Y());
743  }
744  }
745 }
746 
747 void CbmRichUrqmdTest::DrawHist() {
748  cout.precision(4);
749 
750  SetDefaultDrawStyle();
751 
752  {
753  fHM->H1("fh_vertex_z")->Scale(1. / fEventNum);
754  fHM->CreateCanvas("rich_urqmd_vertex_z", "rich_urqmd_vertex_z", 800, 800);
755  DrawH1(fHM->H1("fh_vertex_z"), kLinear, kLog, "hist");
756  }
757 
758  {
759  fHM->H1("fh_vertex_z_sts")->Scale(1. / fEventNum);
760  fHM->CreateCanvas(
761  "rich_urqmd_vertex_z_sts", "rich_urqmd_vertex_z_sts", 800, 800);
762  DrawH1(fHM->H1("fh_vertex_z_sts"), kLinear, kLog, "hist");
763  }
764 
765 
766  {
767  fHM->H2("fh_vertex_xy")->Scale(1. / fEventNum);
768  fHM->H2("fh_vertex_zy")->Scale(1. / fEventNum);
769  fHM->H2("fh_vertex_zx")->Scale(1. / fEventNum);
770  TCanvas* c = fHM->CreateCanvas(
771  "rich_urqmd_vertex_xyz", "rich_urqmd_vertex_xyz", 1500, 500);
772  c->Divide(3, 1);
773  c->cd(1);
774  DrawH2(fHM->H2("fh_vertex_xy"), kLinear, kLinear, kLog);
775  c->cd(2);
776  DrawH2(fHM->H2("fh_vertex_zy"), kLinear, kLinear, kLog);
777  c->cd(3);
778  DrawH2(fHM->H2("fh_vertex_zx"), kLinear, kLinear, kLog);
779  }
780 
781  {
782  gStyle->SetOptTitle(1);
783 
784  fHM->H2("fh_vertex_xy_z5")->Scale(1. / fEventNum);
785  fHM->H2("fh_vertex_xy_z5_15")->Scale(1. / fEventNum);
786  fHM->H2("fh_vertex_xy_z15_25")->Scale(1. / fEventNum);
787  fHM->H2("fh_vertex_xy_z25_35")->Scale(1. / fEventNum);
788  fHM->H2("fh_vertex_xy_z35_45")->Scale(1. / fEventNum);
789  fHM->H2("fh_vertex_xy_z45_55")->Scale(1. / fEventNum);
790  fHM->H2("fh_vertex_xy_z55_65")->Scale(1. / fEventNum);
791  fHM->H2("fh_vertex_xy_z65_75")->Scale(1. / fEventNum);
792  fHM->H2("fh_vertex_xy_z75_85")->Scale(1. / fEventNum);
793  fHM->H2("fh_vertex_xy_z85_95")->Scale(1. / fEventNum);
794  fHM->H2("fh_vertex_xy_z95_105")->Scale(1. / fEventNum);
795 
796 
797  TCanvas* c = fHM->CreateCanvas(
798  "rich_urqmd_vertex_sts_xyz", "rich_urqmd_vertex_sts_xyz", 1600, 1200);
799  c->Divide(4, 3);
800  c->cd(1);
801  DrawH2(fHM->H2("fh_vertex_xy_z5"), kLinear, kLinear, kLog);
802  DrawTextOnPad("Z < 5 cm", 0.3, 0.9, 0.7, 0.98);
803  c->cd(2);
804  DrawH2(fHM->H2("fh_vertex_xy_z5_15"), kLinear, kLinear, kLog);
805  DrawTextOnPad("5 < Z < 15 cm", 0.3, 0.9, 0.7, 0.98);
806  c->cd(3);
807  DrawH2(fHM->H2("fh_vertex_xy_z15_25"), kLinear, kLinear, kLog);
808  DrawTextOnPad("15 < Z < 25 cm", 0.3, 0.9, 0.7, 0.98);
809  c->cd(4);
810  DrawH2(fHM->H2("fh_vertex_xy_z25_35"), kLinear, kLinear, kLog);
811  DrawTextOnPad("25 < Z < 35 cm", 0.3, 0.9, 0.7, 0.98);
812  c->cd(5);
813  DrawH2(fHM->H2("fh_vertex_xy_z35_45"), kLinear, kLinear, kLog);
814  DrawTextOnPad("35 < Z < 45 cm", 0.3, 0.9, 0.7, 0.98);
815  c->cd(6);
816  DrawH2(fHM->H2("fh_vertex_xy_z45_55"), kLinear, kLinear, kLog);
817  DrawTextOnPad("45 < Z < 55 cm", 0.3, 0.9, 0.7, 0.98);
818  c->cd(7);
819  DrawH2(fHM->H2("fh_vertex_xy_z55_65"), kLinear, kLinear, kLog);
820  DrawTextOnPad("55 < Z < 65 cm", 0.3, 0.9, 0.7, 0.98);
821  c->cd(8);
822  DrawH2(fHM->H2("fh_vertex_xy_z65_75"), kLinear, kLinear, kLog);
823  DrawTextOnPad("65 < Z < 75 cm", 0.3, 0.9, 0.7, 0.98);
824  c->cd(9);
825  DrawH2(fHM->H2("fh_vertex_xy_z75_85"), kLinear, kLinear, kLog);
826  DrawTextOnPad("75 < Z < 85 cm", 0.3, 0.9, 0.7, 0.98);
827  c->cd(10);
828  DrawH2(fHM->H2("fh_vertex_xy_z85_95"), kLinear, kLinear, kLog);
829  DrawTextOnPad("85 < Z < 95 cm", 0.3, 0.9, 0.7, 0.98);
830  c->cd(11);
831  DrawH2(fHM->H2("fh_vertex_xy_z95_105"), kLinear, kLinear, kLog);
832  DrawTextOnPad("95 < Z < 105 cm", 0.3, 0.9, 0.7, 0.98);
833 
834  gStyle->SetOptTitle(0);
835  }
836 
837  {
838  fHM->H2("fh_vertex_xy_z100_180")->Scale(1. / fEventNum);
839  fHM->CreateCanvas("rich_urqmd_vertex_xy_z100_180",
840  "rich_urqmd_vertex_xy_z100_180",
841  800,
842  800);
843  DrawH2(fHM->H2("fh_vertex_xy_z100_180"));
844  }
845 
846  {
847  fHM->H2("fh_vertex_xy_z180_370")->Scale(1. / fEventNum);
848  fHM->CreateCanvas("rich_urqmd_vertex_xy_z180_370",
849  "rich_urqmd_vertex_xy_z180_370",
850  800,
851  800);
852  DrawH2(fHM->H2("fh_vertex_xy_z180_370"));
853  }
854 
855  {
856  fHM->H2("fh_vertex_xy_z180_230")->Scale(1. / fEventNum);
857  fHM->CreateCanvas("rich_urqmd_vertex_xy_z180_230",
858  "rich_urqmd_vertex_xy_z180_230",
859  800,
860  800);
861  DrawH2(fHM->H2("fh_vertex_xy_z180_230"));
862  }
863 
864  {
865  fHM->H1("fh_nof_points_per_event")->Scale(1. / fEventNum);
866  fHM->CreateCanvas("rich_urqmd_nof_points_per_event",
867  "rich_urqmd_nof_points_per_event",
868  800,
869  800);
870  //gStyle->SetPaintTextFormat("4.1f");
871  string labels[7] = {"all",
872  "e^{#pm}_{sec} other",
873  "e^{#pm}_{target} from #gamma",
874  "e^{#pm}_{not target} from #gamma",
875  "#pi^{#pm}",
876  "K^{#pm}",
877  "#mu^{#pm}"};
878  for (Int_t i = 1; i <= 7; i++) {
879  fHM->H1("fh_nof_points_per_event")
880  ->GetXaxis()
881  ->SetBinLabel(i, labels[i - 1].c_str());
882  }
883  fHM->H1("fh_nof_points_per_event")->GetXaxis()->SetLabelSize(0.05);
884  //fHM->H1("fh_nof_points_per_event")->SetMarkerSize(0.15);
885  DrawH1(fHM->H1("fh_nof_points_per_event"), kLinear, kLog, "hist text0");
886  }
887 
888  {
889  fHM->CreateCanvas("rich_urqmd_nof_rings", "rich_urqmd_nof_rings", 800, 800);
890  fHM->H1("fh_nof_rings_1hit")
891  ->Scale(1. / fHM->H1("fh_nof_rings_1hit")->Integral());
892  fHM->H1("fh_nof_rings_7hits")
893  ->Scale(1. / fHM->H1("fh_nof_rings_7hits")->Integral());
894  stringstream ss1, ss2;
895  ss1 << "At least 1 hit detected ("
896  << fHM->H1("fh_nof_rings_1hit")->GetMean() << ")";
897  ss2 << "At least 7 hits detected ("
898  << fHM->H1("fh_nof_rings_7hits")->GetMean() << ")";
899  DrawH1({fHM->H1("fh_nof_rings_1hit"), fHM->H1("fh_nof_rings_7hits")},
900  {ss1.str(), ss2.str()},
901  kLinear,
902  kLinear,
903  true,
904  0.3,
905  0.85,
906  0.99,
907  0.99,
908  "hist");
909  }
910 
911  {
912  fHM->CreateCanvas(
913  "rich_urqmd_nof_rings_prim", "rich_urqmd_nof_rings_prim", 800, 800);
914  fHM->H1("fh_nof_rings_prim_1hit")
915  ->Scale(1. / fHM->H1("fh_nof_rings_prim_1hit")->Integral());
916  fHM->H1("fh_nof_rings_prim_7hits")
917  ->Scale(1. / fHM->H1("fh_nof_rings_prim_7hits")->Integral());
918  stringstream ss1, ss2;
919  ss1 << "At least 1 hit detected ("
920  << fHM->H1("fh_nof_rings_prim_1hit")->GetMean() << ")";
921  ss2 << "At least 7 hits detected ("
922  << fHM->H1("fh_nof_rings_prim_7hits")->GetMean() << ")";
923  DrawH1(
924  {fHM->H1("fh_nof_rings_prim_1hit"), fHM->H1("fh_nof_rings_prim_7hits")},
925  {ss1.str(), ss2.str()},
926  kLinear,
927  kLinear,
928  true,
929  0.3,
930  0.85,
931  0.99,
932  0.99,
933  "hist");
934  }
935 
936  {
937  fHM->CreateCanvas(
938  "rich_urqmd_nof_rings_target", "rich_urqmd_nof_rings_target", 800, 800);
939  fHM->H1("fh_nof_rings_target_1hit")
940  ->Scale(1. / fHM->H1("fh_nof_rings_target_1hit")->Integral());
941  fHM->H1("fh_nof_rings_target_7hits")
942  ->Scale(1. / fHM->H1("fh_nof_rings_target_7hits")->Integral());
943  stringstream ss1, ss2;
944  ss1 << "At least 1 hit detected ("
945  << fHM->H1("fh_nof_rings_target_1hit")->GetMean() << ")";
946  ss2 << "At least 7 hits detected ("
947  << fHM->H1("fh_nof_rings_target_7hits")->GetMean() << ")";
948  DrawH1({fHM->H1("fh_nof_rings_target_1hit"),
949  fHM->H1("fh_nof_rings_target_7hits")},
950  {ss1.str(), ss2.str()},
951  kLinear,
952  kLinear,
953  true,
954  0.3,
955  0.85,
956  0.99,
957  0.99,
958  "hist");
959  }
960 
961  {
962  fHM->CreateCanvas(
963  "rich_urqmd_sources_mom", "rich_urqmd_sources_mom", 800, 800);
964  fHM->H1("fh_gamma_target_mom")->Scale(1. / fEventNum);
965  fHM->H1("fh_gamma_nontarget_mom")->Scale(1. / fEventNum);
966  fHM->H1("fh_secel_mom")->Scale(1. / fEventNum);
967  fHM->H1("fh_pi_mom")->Scale(1. / fEventNum);
968  fHM->H1("fh_kaon_mom")->Scale(1. / fEventNum);
969  fHM->H1("fh_mu_mom")->Scale(1. / fEventNum);
970  stringstream ss1, ss2, ss3, ss4, ss5, ss6;
971  ss1 << "e^{#pm}_{target} from #gamma ("
972  << fHM->H1("fh_gamma_target_mom")->GetEntries() / fEventNum << ")";
973  ss2 << "e^{#pm}_{not target} from #gamma ("
974  << fHM->H1("fh_gamma_nontarget_mom")->GetEntries() / fEventNum << ")";
975  ss3 << "e^{#pm}_{sec} other ("
976  << fHM->H1("fh_secel_mom")->GetEntries() / fEventNum << ")";
977  ss4 << "#pi^{#pm} (" << fHM->H1("fh_pi_mom")->GetEntries() / fEventNum
978  << ")";
979  ss5 << "K^{#pm} (" << fHM->H1("fh_kaon_mom")->GetEntries() / fEventNum
980  << ")";
981  ss6 << "#mu^{#pm} (" << fHM->H1("fh_mu_mom")->GetEntries() / fEventNum
982  << ")";
983  DrawH1(
984  {fHM->H1("fh_gamma_target_mom"),
985  fHM->H1("fh_gamma_nontarget_mom"),
986  fHM->H1("fh_secel_mom"),
987  fHM->H1("fh_pi_mom"),
988  fHM->H1("fh_kaon_mom"),
989  fHM->H1("fh_mu_mom")},
990  list_of(ss1.str())(ss2.str())(ss3.str())(ss4.str())(ss5.str())(ss6.str()),
991  kLinear,
992  kLog,
993  true,
994  0.5,
995  0.7,
996  0.99,
997  0.99,
998  "hist");
999  }
1000 
1001  {
1002  TCanvas* c =
1003  fHM->CreateCanvas("rich_urqmd_hits_xy", "rich_urqmd_hits_xy", 800, 800);
1004  TH2D* clone = (TH2D*) fHM->H2("fh_hits_xy")->Clone();
1005  clone->Scale(1. / (fEventNum));
1006  CbmRichDraw::DrawPmtH2(clone, c, true);
1007  }
1008 
1009  {
1010  TCanvas* c = fHM->CreateCanvas(
1011  "rich_urqmd_occupancy_xy", "rich_urqmd_occupancy_xy", 800, 800);
1012  TH2D* clone = (TH2D*) fHM->H2("fh_hits_xy")->Clone();
1013  clone->GetZaxis()->SetTitle("Occupancy:# hits/PMT/ev./64 [%]");
1014  clone->Scale(100. / (fEventNum * 64.));
1015  CbmRichDraw::DrawPmtH2(clone, c, true);
1016  }
1017 
1018  {
1019  TCanvas* c = fHM->CreateCanvas(
1020  "rich_urqmd_skipped_pmt_10_xy", "rich_urqmd_skipped_pmt_10_xy", 800, 800);
1021  fHM->H2("fh_skipped_pmt_10_xy")->Scale(100. / (fEventNum));
1022  CbmRichDraw::DrawPmtH2(fHM->H2("fh_skipped_pmt_10_xy"), c, true);
1023  }
1024 
1025  {
1026  TCanvas* c = fHM->CreateCanvas(
1027  "rich_urqmd_skipped_pmt_20_xy", "rich_urqmd_skipped_pmt_20_xy", 800, 800);
1028  fHM->H2("fh_skipped_pmt_20_xy")->Scale(100. / (fEventNum));
1029  CbmRichDraw::DrawPmtH2(fHM->H2("fh_skipped_pmt_20_xy"), c, true);
1030  }
1031 
1032  {
1033  TCanvas* c = fHM->CreateCanvas(
1034  "rich_urqmd_skipped_pmt_30_xy", "rich_urqmd_skipped_pmt_30_xy", 800, 800);
1035  fHM->H2("fh_skipped_pmt_30_xy")->Scale(100. / (fEventNum));
1036  CbmRichDraw::DrawPmtH2(fHM->H2("fh_skipped_pmt_30_xy"), c, true);
1037  }
1038 
1039  {
1040  fHM->CreateCanvas(
1041  "rich_urqmd_nof_hits_per_pmt", "rich_urqmd_nof_hits_per_pmt", 800, 800);
1042  fHM->H1("fh_nof_hits_per_pmt")
1043  ->Scale(100. / fHM->H1("fh_nof_hits_per_pmt")->Integral());
1044  DrawH1(fHM->H1("fh_nof_hits_per_pmt"), kLinear, kLog, "hist");
1045  }
1046 
1047  {
1048  TCanvas* c = fHM->CreateCanvas(
1049  "rich_urqmd_points_xy", "rich_urqmd_points_xy", 800, 800);
1050  fHM->H2("fh_points_xy")->Scale(1. / fEventNum);
1051  CbmRichDraw::DrawPmtH2(fHM->H2("fh_points_xy"), c, true);
1052  }
1053 
1054  {
1055  TCanvas* c = fHM->CreateCanvas(
1056  "rich_urqmd_points_xy_pions", "rich_urqmd_points_xy_pions", 800, 800);
1057  fHM->H2("fh_points_xy_pions")->Scale(1. / fEventNum);
1058  CbmRichDraw::DrawPmtH2(fHM->H2("fh_points_xy_pions"), c, true);
1059  }
1060 
1061  {
1062  TCanvas* c = fHM->CreateCanvas("rich_urqmd_points_xy_gamma_target",
1063  "rich_urqmd_points_xy_gamma_target",
1064  800,
1065  800);
1066  fHM->H2("fh_points_xy_gamma_target")->Scale(1. / fEventNum);
1067  CbmRichDraw::DrawPmtH2(fHM->H2("fh_points_xy_gamma_target"), c, true);
1068  }
1069 
1070  {
1071  TCanvas* c = fHM->CreateCanvas("rich_urqmd_points_xy_gamma_nontarget",
1072  "rich_urqmd_points_xy_gamma_nontarget",
1073  800,
1074  800);
1075  fHM->H2("fh_points_xy_gamma_nontarget")->Scale(1. / fEventNum);
1076  CbmRichDraw::DrawPmtH2(fHM->H2("fh_points_xy_gamma_nontarget"), c, true);
1077  }
1078 
1079  {
1080  fHM->CreateCanvas("rich_urqmd_nof_hits_per_event",
1081  "rich_urqmd_nof_hits_per_event",
1082  800,
1083  800);
1084  fHM->H1("fh_nof_hits_per_event")
1085  ->Scale(1. / fHM->H1("fh_nof_hits_per_event")->Integral());
1086  DrawH1andFitGauss(fHM->H1("fh_nof_hits_per_event"));
1087  cout << "Mean number of hits per event = "
1088  << fHM->H1("fh_nof_hits_per_event")->GetMean() << endl;
1089  }
1090 
1091  {
1092  TCanvas* c = fHM->CreateCanvas(
1093  "rich_urqmd_hitrate_xy", "rich_urqmd_hitrate_xy", 800, 800);
1094  fHM->H2("fh_hitrate_xy")->Scale(1e7 / (fEventNum * 64.));
1095  CbmRichDraw::DrawPmtH2(fHM->H2("fh_hitrate_xy"), c, true);
1096  //DrawH2(fHM->H2("fh_hitrate_xy"));
1097  }
1098 
1099  {
1100  TCanvas* c =
1101  fHM->CreateCanvas("rich_urqmd_proj_xy", "rich_urqmd_proj_xy", 800, 800);
1102  fHM->H2("fh_proj_xy")->Scale(1. / fEventNum);
1103  CbmRichDraw::DrawPmtH2(fHM->H2("fh_proj_xy"), c);
1104  }
1105 
1106  {
1107  fHM->CreateCanvas("rich_urqmd_nof_proj_per_event",
1108  "rich_urqmd_nof_proj_per_event",
1109  800,
1110  800);
1111  fHM->H1("fh_nof_proj_per_event")->Scale(1. / fEventNum);
1112  DrawH1andFitGauss(fHM->H1("fh_nof_proj_per_event"));
1113  cout << "Number of track projections per event = "
1114  << fHM->H1("fh_nof_proj_per_event")->GetMean() << endl;
1115  }
1116 }
1117 
1118 void CbmRichUrqmdTest::Finish() {
1119  DrawHist();
1120  fHM->SaveCanvasToImage(fOutputDir);
1121  TDirectory* oldir = gDirectory;
1122  TFile* outFile = FairRootManager::Instance()->GetOutFile();
1123  if (outFile != NULL) {
1124  outFile->cd();
1125  fHM->WriteToFile();
1126  }
1127  gDirectory->cd(oldir->GetPath());
1128 }
1129 
1130 ClassImp(CbmRichUrqmdTest)
CbmRichPoint.h
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmRichPixelData
Definition: CbmRichDetectorData.h:17
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition: CbmMCDataManager.cxx:137
CbmMCEventList::GetFileIdByIndex
Int_t GetFileIdByIndex(UInt_t index)
File number by index @value File number for event at given index in list.
Definition: CbmMCEventList.cxx:106
CbmMatch
Definition: CbmMatch.h:22
CbmRichPmtData::fY
Double_t fY
Definition: CbmRichDetectorData.h:42
CbmMCDataManager.h
CbmRichDigiMapManager::GetPmtIds
std::vector< Int_t > GetPmtIds()
Definition: CbmRichDigiMapManager.cxx:285
CbmHistManager::WriteToFile
void WriteToFile()
Write all histograms to current opened file.
Definition: core/base/CbmHistManager.cxx:103
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
CbmRichDigi
Definition: CbmRichDigi.h:25
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition: CbmMCDataArray.cxx:133
CbmRichUrqmdTest::NofHitsAndPoints
void NofHitsAndPoints()
Definition: CbmRichUrqmdTest.cxx:535
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
CbmRichUrqmdTest::fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmRichUrqmdTest.h:130
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichDigiMapManager::GetPmtDataById
CbmRichPmtData * GetPmtDataById(Int_t id)
Definition: CbmRichDigiMapManager.cxx:288
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmRichPmtData::fZ
Double_t fZ
Definition: CbmRichDetectorData.h:43
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition: CbmHistManager.h:190
CbmMCDataArray.h
CbmRichRing::GetNofHits
Int_t GetNofHits() const
Definition: CbmRichRing.h:40
CbmRichRing
Definition: CbmRichRing.h:17
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
CbmRichUrqmdTest::~CbmRichUrqmdTest
virtual ~CbmRichUrqmdTest()
Standard destructor.
Definition: CbmRichUrqmdTest.cxx:61
CbmRichUrqmdTest::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmRichUrqmdTest.h:128
CbmRichPmtData
Definition: CbmRichDetectorData.h:26
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.
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition: CbmDigiManager.h:93
CbmRichUrqmdTest::fRichPoints
CbmMCDataArray * fRichPoints
Definition: CbmRichUrqmdTest.h:126
CbmRichUrqmdTest::Vertex
void Vertex()
Definition: CbmRichUrqmdTest.cxx:695
CbmRichPixelData::fPmtId
Int_t fPmtId
Definition: CbmRichDetectorData.h:23
CbmRichGeoManager.h
CbmRichDigi.h
CbmRichDraw::GetPmtHistYbins
static std::vector< Double_t > GetPmtHistYbins()
Definition: CbmRichDraw.h:73
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
CbmRichUrqmdTest::NofRings
void NofRings()
Definition: CbmRichUrqmdTest.cxx:465
CbmRichUrqmdTest::fRichRings
TClonesArray * fRichRings
Definition: CbmRichUrqmdTest.h:125
CbmMCEventList::GetNofEvents
std::size_t GetNofEvents() const
Number of events in the list @value Number of events.
Definition: CbmMCEventList.h:90
CbmRichDraw::DrawPmtH2
static void DrawPmtH2(TH2 *h, TCanvas *c, Bool_t usePmtBins=false)
Definition: CbmRichDraw.h:21
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmRichUrqmdTest::fRichHits
TClonesArray * fRichHits
Definition: CbmRichUrqmdTest.h:124
CbmRichUrqmdTest.h
RICH geometry testing in Urqmd collisions.
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmRichUrqmdTest::CbmRichUrqmdTest
CbmRichUrqmdTest()
Standard constructor.
Definition: CbmRichUrqmdTest.cxx:45
CbmRichUrqmdTest::FillRichRingNofHits
void FillRichRingNofHits()
Definition: CbmRichUrqmdTest.cxx:441
CbmRichUrqmdTest::fEventList
CbmMCEventList * fEventList
Definition: CbmRichUrqmdTest.h:131
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmRichUrqmdTest::fHM
CbmHistManager * fHM
Definition: CbmRichUrqmdTest.h:120
kLinear
@ kLinear
Definition: CbmDrawHist.h:79
CbmTrackMatchNew.h
CbmRichUrqmdTest::fEventNum
Int_t fEventNum
Definition: CbmRichUrqmdTest.h:133
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
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:407
CbmRichPmtData::fX
Double_t fX
Definition: CbmRichDetectorData.h:41
ECbmModuleId::kRich
@ kRich
Ring-Imaging Cherenkov Detector.
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmMCEventList
Container class for MC events with number, file and start time.
Definition: CbmMCEventList.h:38
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition: CbmMCTrack.h:182
CbmRichUrqmdTest::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichUrqmdTest.cxx:114
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition: CbmHistManager.h:170
CbmRichUrqmdTest::fNofHitsInRingMap
std::map< pair< Int_t, Int_t >, Int_t > fNofHitsInRingMap
Definition: CbmRichUrqmdTest.h:138
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
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
CbmRichDigi::GetAddress
Int_t GetAddress() const
Definition: CbmRichDigi.h:37
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
CbmRichUrqmdTest::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichUrqmdTest.cxx:1118
CbmRichUrqmdTest
RICH geometry testing in Urqmd collisions.
Definition: CbmRichUrqmdTest.h:35
CbmRichUrqmdTest::NofProjections
void NofProjections()
Definition: CbmRichUrqmdTest.cxx:682
CbmMCTrack.h
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
CbmMCEventList::GetEventIdByIndex
Int_t GetEventIdByIndex(UInt_t index)
Event number by index @value Event number for event at given index in list.
Definition: CbmMCEventList.cxx:77
CbmRichUrqmdTest::fOutputDir
string fOutputDir
Definition: CbmRichUrqmdTest.h:122
CbmDigiManager.h
CbmRichUrqmdTest::InitHistograms
void InitHistograms()
Initialize histograms.
Definition: CbmRichUrqmdTest.cxx:127
CbmRichUrqmdTest::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichUrqmdTest.cxx:63
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmRichUrqmdTest::DrawHist
void DrawHist()
Draw histograms.
Definition: CbmRichUrqmdTest.cxx:747
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
DrawH2
void DrawH2(TH2 *hist, HistScale logx, HistScale logy, HistScale logz, const string &drawOpt)
Definition: CbmDrawHist.cxx:84
SetDefaultDrawStyle
void SetDefaultDrawStyle()
Definition: CbmDrawHist.cxx:33
CbmRichHit.h
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmRichDigiMapManager::GetInstance
static CbmRichDigiMapManager & GetInstance()
Definition: CbmRichDigiMapManager.h:29
CbmRichDraw::GetPmtHistXbins
static std::vector< Double_t > GetPmtHistXbins()
Definition: CbmRichDraw.h:69
CbmRichUrqmdTest::fMinNofHits
Int_t fMinNofHits
Definition: CbmRichUrqmdTest.h:135
CbmMCTrack::GetP
Double_t GetP() const
Definition: CbmMCTrack.h:100
kLog
@ kLog
Definition: CbmDrawHist.h:78
CbmRichUrqmdTest::fRichProjections
TClonesArray * fRichProjections
Definition: CbmRichUrqmdTest.h:129
CbmRichUrqmdTest::fMcTracks
CbmMCDataArray * fMcTracks
Definition: CbmRichUrqmdTest.h:127
DrawH1andFitGauss
void DrawH1andFitGauss(TH1 *hist, Bool_t drawResults, Bool_t doScale, Double_t userRangeMin, Double_t userRangeMax)
Definition: CbmDrawHist.cxx:266
CbmRichHit
Definition: CbmRichHit.h:19
CbmRichDigiMapManager::GetPixelDataByAddress
CbmRichPixelData * GetPixelDataByAddress(Int_t address)
Definition: CbmRichDigiMapManager.cxx:267
CbmRichUrqmdTest::PmtXYSource
void PmtXYSource()
Definition: CbmRichUrqmdTest.cxx:631
CbmHistManager::Add
void Add(const std::string &name, TNamed *object)
Add new named object to manager.
Definition: CbmHistManager.h:58