CbmRoot
CbmRichMCbmQa.cxx
Go to the documentation of this file.
1 #include "CbmRichMCbmQa.h"
2 
3 #include "TCanvas.h"
4 #include "TClonesArray.h"
5 #include "TEllipse.h"
6 #include "TF1.h"
7 #include "TGeoBBox.h"
8 #include "TGeoManager.h"
9 #include "TGeoNode.h"
10 #include "TH1.h"
11 #include "TH1D.h"
12 #include "TMarker.h"
13 #include "TStyle.h"
14 #include <TFile.h>
15 
16 
17 #include "CbmDrawHist.h"
18 #include "CbmGlobalTrack.h"
19 #include "CbmMCTrack.h"
20 #include "CbmMatchRecoToMC.h"
21 #include "CbmRichDraw.h"
22 #include "CbmRichGeoManager.h"
23 #include "CbmRichHit.h"
24 #include "CbmRichPoint.h"
25 #include "CbmRichRing.h"
26 #include "CbmTofHit.h"
27 #include "CbmTofPoint.h"
28 #include "CbmTrackMatchNew.h"
29 #include "CbmTrdTrack.h"
30 #include "FairMCPoint.h"
31 #include "FairTrackParam.h"
32 
33 #include "CbmHistManager.h"
34 #include "CbmUtils.h"
35 
36 #include <boost/assign/list_of.hpp>
37 #include <iostream>
38 #include <sstream>
39 #include <string>
40 
41 using namespace std;
42 using boost::assign::list_of;
43 
44 CbmRichMCbmQa::CbmRichMCbmQa()
45  : FairTask("CbmRichMCbmQa")
46  , fHM(NULL)
47  , fEventNum(0)
48  , fOutputDir("")
49  , fMCTracks(NULL)
50  , fRichPoints(NULL)
51  , fRichDigis(NULL)
52  , fRichHits(NULL)
53  , fRichRings(NULL)
54  , fRichRingMatches(NULL)
55  , fRefPlanePoints(NULL)
56  , fGlobalTracks(NULL)
57  , fTrdTracks(NULL)
58  , fTofHits(NULL)
59  , fTofPoints(NULL)
60  , fTofHitMatches(NULL) {}
61 
62 InitStatus CbmRichMCbmQa::Init() {
63  cout << "CbmRichMCbmQa::Init" << endl;
64 
65  FairRootManager* ioman = FairRootManager::Instance();
66  if (NULL == ioman) {
67  Fatal("CbmRichMCbmQa::Init", "RootManager not instantised!");
68  }
69 
70  fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
71  if (NULL == fMCTracks) { Fatal("CbmRichMCbmQa::Init", "No MC Tracks!"); }
72 
73  fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
74  if (NULL == fRichPoints) { Fatal("CbmRichMCbmQa::Init", "No Rich Points!"); }
75 
76  fRichDigis = (TClonesArray*) ioman->GetObject("RichDigi");
77  if (NULL == fRichDigis) { Fatal("CbmRichMCbmQa::Init", "No Rich Digis!"); }
78 
79  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
80  if (NULL == fRichHits) { Fatal("CbmRichMCbmQa::Init", "No RichHits!"); }
81 
82  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
83  if (NULL == fRichRings) { Fatal("CbmRichMCbmQa::Init", "No RichRings!"); }
84 
85  fTofHits = (TClonesArray*) ioman->GetObject("TofHit");
86  if (NULL == fTofHits) { Fatal("CbmRichMCbmQa::Init", "No TofHits!"); }
87 
88  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
89  if (NULL == fRichRingMatches) {
90  Fatal("CbmRichMCbmQa::Init", "No RichRingMatch array!");
91  }
92 
93  fTofHitMatches = (TClonesArray*) ioman->GetObject("TofHitMatch");
94  if (NULL == fTofHitMatches) {
95  Fatal("CbmRichMCbmQa::Init", "No TofHitMatch array!");
96  }
97 
98  fTofPoints = (TClonesArray*) ioman->GetObject("TofPoint");
99  if (NULL == fTofPoints) { Fatal("CbmRichMCbmQa::Init", "No Tof Points!"); }
100 
101  fRefPlanePoints = (TClonesArray*) ioman->GetObject("RefPlanePoint");
102  if (NULL == fRefPlanePoints) {
103  Fatal("CbmRichMCbmQa::Init", "No RefPlanePoints!");
104  }
105 
106  // fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
107  // if ( NULL == fGlobalTracks) { Fatal("CbmRichMCbmQa::Init","No Global Tracks!");}
108 
109 
110  InitHistograms();
111 
112  return kSUCCESS;
113 }
114 /*
115 vector<Double_t> CbmRichMCbmQa::GetHistBins(Bool_t isX)
116 {
117  set<Double_t> setBins;
118  TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
119  TGeoNode* curNode;
120  TString pixelNameStr("pmt_pixel");
121  geoIterator.Reset();
122  Double_t halfPixelSize = 0.;
123  while (curNode == geoIterator()) {
124  TString nodeName(curNode->GetName());
125  TString nodePath;
126  if (TString(curNode->GetVolume()->GetName()).Contains(pixelNameStr)) {
127  geoIterator.GetPath(nodePath);
128  const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
129  const Double_t* curNodeTr = curMatrix->GetTranslation();
130  const TGeoBBox* shape = (const TGeoBBox*)(curNode->GetVolume()->GetShape());
131  Double_t lowEdge, upEdge;
132  if (isX) {
133  halfPixelSize = shape->GetDX();
134  lowEdge = curNodeTr[0] - halfPixelSize;
135  upEdge = curNodeTr[0] + halfPixelSize;
136  } else {
137  halfPixelSize = shape->GetDY();
138  lowEdge = curNodeTr[1] - halfPixelSize;
139  upEdge = curNodeTr[1] + halfPixelSize;
140  }
141  setBins.insert(lowEdge);
142  setBins.insert(upEdge);
143  }
144  }
145 
146  vector<double> bins(setBins.begin(), setBins.end());
147  sort (bins.begin(), bins.end());
148 
149  return bins;
150 }
151  */
152 
153 void CbmRichMCbmQa::InitHistograms() {
154  cout << "CbmRichMCbmQa::InitHistograms" << endl;
155 
156  fHM = new CbmHistManager();
157 
158 
159  //RICH
160  fHM->Create1<TH1D>(
161  "fh_nof_rich_points",
162  "fh_nof_rich_points;Nof RICH Points per event;Yield (a.u.)",
163  1000,
164  -.5,
165  999.5);
166  fHM->Create1<TH1D>("fh_nof_rich_hits",
167  "fh_nof_rich_hits;Nof RICH hits per event;Yield (a.u.)",
168  100,
169  -.5,
170  99.5);
171  fHM->Create1<TH1D>("fh_nof_rich_rings",
172  "fh_nof_rich_rings;Nof RICH rings;# per event",
173  20,
174  -0.5,
175  19.5);
176  fHM->Create1<TH1D>("fh_rich_ring_radius",
177  "fh_rich_ring_radius;Ring radius [cm];Yield (a.u.)",
178  200,
179  1.,
180  8.);
181 
182 
183  // vector<Double_t> xbins = GetHistBins(true);
184  // vector<Double_t> ybins = GetHistBins(false);
185  // fHM->Add("fh_rich_hits_xy", new TH2D("fh_rich_hits_xy", "fh_rich_hits_xy;X [cm];Y [cm];Yield (a.u.)", xbins.size() - 1, &xbins[0], ybins.size() - 1, &ybins[0]));
186  // cout << "binsize: " << xbins.size() << "bin Ende Anfang: " << &xbins[0] << endl;
187 
188  const Double_t xMin = -35.5;
189  const Double_t xMax = 35.5;
190  const Double_t yMin = -40.5;
191  const Double_t yMax = 40.5;
192 
193 
194  // RICH hists
195  fHM->Create2<TH2D>("fh_rich_hits_xy",
196  "fh_rich_hits_xy;X [cm];Y [cm];Yield",
197  71,
198  xMin,
199  xMax,
200  81,
201  yMin,
202  yMax);
203  fHM->Create2<TH2D>("fh_rich_points_xy",
204  "fh_rich_points_xy;X [cm];Y [cm];Yield (a.u.)",
205  250,
206  xMin,
207  xMax,
208  250,
209  yMin,
210  yMax);
211 
212  fHM->Create1<TH1D>(
213  "fh_hits_per_ring", "fh_hits_per_ring;Hits per Ring;Yield", 18, 2.5, 20.5);
214  fHM->Create1<TH1D>("fh_dR", "fh_dR;dR [cm];Yield (a.u.)", 80, -0.8, 0.8);
215  fHM->Create1<TH1D>(
216  "fh_photon_energy", "fh_photon_energy;Momentum [eV]", 100, 0., 10.);
217  fHM->Create2<TH2D>("fh_radius_momentum",
218  "fh_radius_momentum; Momentum [GeV];Radius [cm]; Yield",
219  500,
220  0.,
221  5.,
222  500,
223  0.,
224  5.);
225  fHM->Create2<TH2D>("fh_ring_center_xy",
226  "fh_ring_center_xy;X [cm];Y [cm];Yield (a.u.)",
227  100,
228  xMin,
229  xMax,
230  100,
231  yMin,
232  yMax);
233  fHM->Create2<TH2D>("fh_nonphoton_pmt_points_xy",
234  "fh_nonphoton_pmt_points_xy;X [cm]; Y [cm];Yield",
235  250.,
236  xMin,
237  xMax,
238  250,
239  yMin,
240  yMax);
241  fHM->Create2<TH2D>("fh_radius_beta",
242  "fh_radius_beta; beta; Radius [cm]; Yield",
243  400,
244  0.7,
245  1.1,
246  500,
247  0.,
248  5.);
249  fHM->Create1<TH1D>(
250  "fh_beta_dis_all", "fh_beta_dis_all;beta;Yield; ", 300, 0.7, 1.);
251  fHM->Create1<TH1D>(
252  "fh_beta_dis_ring", "fh_beta_dis_ring;beta;Yield;", 300, 0.7, 1.);
253 
254  // fHM->Create1<TH1D>("fh_nof_trd_hits","fh_nof_trd_hits; Nof Trd Hits; Yield(a.u.)", 5., -0.5, 4.5);
255 
256  // fHM->Create1<TH1D>("fh_nof_sts_hits","fh_nof_sts_hits; Nof Sts Hits; Yield(a.u.)", 5., -0.5, 4.5);
257 
258  // TOF hists
259  fHM->Create2<TH2D>("fh_eloss_tof",
260  "fh_eloss_tof; time of flight [ns]; energy loss; Yield",
261  500,
262  0.,
263  5.,
264  500,
265  0.,
266  5.);
267  fHM->Create2<TH2D>("fh_radius_mass2",
268  "fh_radius_mass2; mass^2 [GeV^2]; Radius [cm]; Yield",
269  450,
270  -0.3,
271  1.2,
272  500,
273  0.,
274  5.);
275 
276  fHM->Create1<TH1D>(
277  "number_of_events", "number_of_events; something", 1, 0.5, 1.5);
278 }
279 
280 void CbmRichMCbmQa::Exec(Option_t* /*option*/) {
281  fEventNum++;
282  fHM->H1("number_of_events")->Fill(1);
283 
284  cout << "CbmRichMCbmQa, event No. " << fEventNum << endl;
285 
286 
287  // Int_t nofMCTracks = fMCTracks->GetEntriesFast();
288  Int_t nofRichPoints = fRichPoints->GetEntriesFast();
289  // Int_t nofRichDigis = fRichDigis->GetEntriesFast();
290  Int_t nofRichHits = fRichHits->GetEntriesFast();
291  Int_t nofRichRings = fRichRings->GetEntriesFast();
292  // Int_t nofRefPlanePoints = fRefPlanePoints->GetEntriesFast();
293  Int_t nofTofHits = fTofHits->GetEntriesFast();
294  // Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
295 
296  fHM->H1("fh_nof_rich_points")->Fill(nofRichPoints);
297  fHM->H1("fh_nof_rich_hits")->Fill(nofRichHits);
298  fHM->H1("fh_nof_rich_rings")->Fill(nofRichRings);
299 
300  for (int i = 0; i < nofRichHits; i++) {
301  CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(i));
302  if (hit == NULL) continue;
303  fHM->H2("fh_rich_hits_xy")->Fill(hit->GetX(), hit->GetY());
304  }
305 
306 
307  //--------------------TofHits properties -----------------------------------
308 
309  for (int iT = 0; iT < nofTofHits; iT++) {
310  CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(iT));
311  if (NULL == tofHit) continue;
312  CbmTrackMatchNew* tofHitMatch = (CbmTrackMatchNew*) fTofHitMatches->At(iT);
313  if (NULL == tofHitMatch) continue;
314  Int_t tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
315  const CbmTofPoint* tofPoint =
316  static_cast<const CbmTofPoint*>(fTofPoints->At(tofPointIndex));
317  if (NULL == tofPoint) continue;
318  Int_t mcTrackIdTofHit = tofPoint->GetTrackID();
319  if (mcTrackIdTofHit < 0) continue;
320  CbmMCTrack* mcTrackTof = (CbmMCTrack*) fMCTracks->At(mcTrackIdTofHit);
321  if (NULL == mcTrackTof) continue;
322 
323 
324  TVector3 vec;
325  mcTrackTof->GetMomentum(vec);
326  Double_t momTotal = sqrt(vec.Px() * vec.Px() + vec.Py() * vec.Py()
327  + vec.Pz() * vec.Pz()); // GeV
328  Double_t time = tofHit->GetTime();
329  Double_t timect = 0.2998 * time; //time in ns, timect in m
330  Double_t trackLength = tofHit->GetR() / 100;
331  Double_t beta = trackLength / timect;
332  Double_t mass2 = TMath::Power(momTotal, 2.)
333  * (TMath::Power(1 / beta, 2) - 1); //m² = p²*((1/beta)²-1)
334 
335 
336  for (int i = 0; i < nofRichPoints; i++) {
337  CbmRichPoint* point = static_cast<CbmRichPoint*>(fRichPoints->At(i));
338  if (point == NULL) continue;
339 
340  Int_t mcTrackIdPoint = point->GetTrackID();
341  if (mcTrackIdPoint < 0) continue;
342  CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcTrackIdPoint);
343  if (mcTrack == NULL) continue;
344 
345  //if ( mcTrack->GetPdgCode() != 50000050) continue; // select only Cherenkov photons
346 
347  if (mcTrack->GetPdgCode() != 50000050) { // fill in non photons
348  if (mcTrackIdPoint == mcTrackIdTofHit) {
349  fHM->H1("fh_beta_dis_all")->Fill(beta);
350  break;
351  }
352  }
353 
354  Int_t motherId =
355  mcTrack->GetMotherId(); //fill in mothers of Cherenkov photons
356  CbmMCTrack* mcTrackMother = (CbmMCTrack*) fMCTracks->At(motherId);
357  if (mcTrackMother == NULL) continue;
358 
359  if (motherId == mcTrackIdTofHit) {
360  fHM->H1("fh_beta_dis_all")->Fill(beta);
361  // one entry per tof hit only
362  break;
363  }
364  }
365 
366  for (int iRing = 0; iRing < nofRichRings; iRing++) {
367  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iRing));
368  if (NULL == ring) continue;
369  CbmTrackMatchNew* ringMatch =
370  (CbmTrackMatchNew*) fRichRingMatches->At(iRing);
371  if (NULL == ringMatch) continue;
372  Int_t mcTrackIdRing = ringMatch->GetMatchedLink().GetIndex();
373  if (mcTrackIdRing < 0) continue;
374 
375  CbmMCTrack* mcTrackRing = (CbmMCTrack*) fMCTracks->At(mcTrackIdRing);
376  if (mcTrackRing == NULL) continue;
377  Double_t radius = ring->GetRadius();
378 
379  if (mcTrackIdTofHit == mcTrackIdRing) {
380  fHM->H2("fh_radius_mass2")->Fill(mass2, radius);
381  fHM->H2("fh_radius_beta")->Fill(beta, radius);
382  fHM->H1("fh_beta_dis_ring")->Fill(beta);
383  fHM->H2("fh_radius_momentum")->Fill(momTotal, radius);
384  // one entry per tof hit only
385  break;
386  }
387  }
388  }
389 
390 
391  /*
392  for(int i = 0; i < nofRichPoints; i++) {
393  CbmRichPoint* point= static_cast<CbmRichPoint*>(fRichPoints->At(i));
394  if (point == NULL) continue;
395 
396  Int_t mcTrackIdPoint = point->GetTrackID();
397 
398  if (mcTrackIdPoint == NULL) continue;
399  CbmMCTrack* mcTrack = (CbmMCTrack*)fMCTracks->At(mcTrackIdPoint);
400  if (mcTrack == NULL) continue;
401 
402 
403  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
404  // cout << "pdg: " << pdg << endl;
405  fHM->H2("fh_rich_points_xy")->Fill(point->GetX(), point->GetY());
406 
407 
408  for(int iT = 0; iT < nofTofHits; iT++){
409  CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(iT));
410  if(NULL == tofHit) continue;
411  CbmTrackMatchNew* tofHitMatch = (CbmTrackMatchNew*) fTofHitMatches->At(iT);
412  if(NULL == tofHitMatch) continue;
413  Int_t tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
414  const CbmTofPoint* tofPoint = static_cast<const CbmTofPoint*>(fTofPoints->At(tofPointIndex));
415  if(NULL == tofPoint) continue;
416  Int_t mcTrackIdTofHit = tofPoint->GetTrackID();
417  if(mcTrackIdTofHit < 0) continue;
418  CbmMCTrack* mcTrackTof = (CbmMCTrack*)fMCTracks->At(mcTrackIdTofHit);
419  if (NULL == mcTrackTof) continue;
420 
421 
422  Double_t time = tofHit->GetTime();
423  Double_t timect = 0.2998*time; //time in ns, timect in m
424  Double_t trackLength = tofHit->GetR()/100;
425  Double_t beta = trackLength/timect;
426 
427 
428  if(mcTrackIdPoint == mcTrackIdTofHit){
429  fHM->H1("fh_beta_dis_all")->Fill(beta);
430  }
431  }
432 
433  for(int iRing = 0; iRing < nofRichRings; iRing++){
434  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iRing));
435  if (NULL == ring) continue;
436  CbmTrackMatchNew* ringMatch = (CbmTrackMatchNew*) fRichRingMatches->At(iRing);
437  if (NULL == ringMatch) continue;
438  Int_t mcTrackIdRing = ringMatch->GetMatchedLink().GetIndex();
439  if (mcTrackIdRing < 0) continue;
440 
441 
442  CbmMCTrack* mcTrackRing = (CbmMCTrack*)fMCTracks->At(mcTrackIdRing);
443  if(mcTrackRing == NULL) continue;
444 
445  if(mcTrackIdTofHit == mcTrackIdRing){
446  // cout << "mcTrackIdTofHit: " << mcTrackIdTofHit << endl;
447  // cout << "mcTrackIdRing: " << mcTrackIdRing << endl;
448  Int_t pdg = TMath::Abs(mcTrackRing->GetPdgCode());
449  if(pdg != 50000050){
450  fHM->H2("fh_radius_mass2")->Fill(mass2, radius);
451  fHM->H2("fh_radius_beta")->Fill(beta, radius);
452  fHM->H1("fh_beta_dis_ring")->Fill(beta);
453  fHM->H2("fh_radius_momentum")->Fill(momTotal, radius);
454 
455  }
456  }
457  }
458  }
459 
460 */
461 
462 
463  //--------------------TofHits properties -----------------------------------
464  /*
465  for(int iRing = 0; iRing < nofRichRings; iRing++){
466  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iRing));
467  if (NULL == ring) continue;
468  CbmTrackMatchNew* ringMatch = (CbmTrackMatchNew*) fRichRingMatches->At(iRing);
469  if (NULL == ringMatch) continue;
470  Int_t mcTrackIdRing = ringMatch->GetMatchedLink().GetIndex();
471  if (mcTrackIdRing < 0) continue;
472 
473 
474  CbmMCTrack* mcTrackRing = (CbmMCTrack*)fMCTracks->At(mcTrackIdRing);
475  if(mcTrackRing == NULL) continue;
476  Double_t radius = ring->GetRadius();
477 
478  Int_t motherId = mcTrackRing->GetMotherId();
479  // Int_t pdg = TMath::Abs(mcTrackRing->GetPdgCode());
480  // cout << "pdg" << pdg << endl;
481 
482 
483  int nofHits = ring->GetNofHits();
484  TVector3 vec1;
485  mcTrackRing->GetMomentum(vec1);
486  Double_t momTotal1 = sqrt(vec1.Px()*vec1.Px() + vec1.Py()*vec1.Py() + vec1.Pz()*vec1.Pz()); // GeV
487 
488  Double_t cX = ring->GetCenterX();
489  Double_t cY = ring->GetCenterY();
490  fHM->H2("fh_ring_center_xy")->Fill(cX,cY);
491  fHM->H1("fh_hits_per_ring")->Fill(nofHits);
492  fHM->H1("fh_rich_ring_radius")->Fill(radius);
493  //fHM->H2("fh_radius_momentum")->Fill(momTotal1, radius);
494  cout << "richRingId: " << mcTrackIdRing << endl;
495 
496  for (int iH = 0; iH < nofHits; iH++) {
497  Int_t hitInd = ring->GetHit(iH);
498  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
499  if (NULL == hit) continue;
500  Double_t hitX = hit->GetX();
501  Double_t hitY = hit->GetY();
502 
503  Double_t dR = radius - TMath::Sqrt( (cX - hitX)*( cX - hitX) + (cY - hitY)*(cY - hitY) );
504  fHM->H1("fh_dR")->Fill(dR);
505  }
506 
507 
508  for(int iT = 0; iT < nofTofHits; iT++){
509 
510  CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(iT));
511  if(NULL == tofHit) continue;
512  CbmTrackMatchNew* tofHitMatch = (CbmTrackMatchNew*) fTofHitMatches->At(iT);
513  if(NULL == tofHitMatch) continue;
514  Int_t tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
515  const CbmTofPoint* tofPoint = static_cast<const CbmTofPoint*>(fTofPoints->At(tofPointIndex));
516  if(NULL == tofPoint) continue;
517  Int_t mcTrackIdTofHit = tofPoint->GetTrackID();
518  if(mcTrackIdTofHit < 0) continue;
519  CbmMCTrack* mcTrackTof = (CbmMCTrack*)fMCTracks->At(mcTrackIdTofHit);
520  if (NULL == mcTrackTof) continue;
521 
522  TVector3 vec;
523  mcTrackTof->GetMomentum(vec);
524  Double_t momTotal = sqrt(vec.Px()*vec.Px() + vec.Py()*vec.Py() + vec.Pz()*vec.Pz()); // GeV
525  Double_t time = tofHit->GetTime();
526  Double_t timect = 0.2998*time; //time in ns, timect in m
527  Double_t trackLength = tofHit->GetR()/100;
528  Double_t beta = trackLength/timect;
529 
530  Double_t mass2 = TMath::Power(momTotal, 2.) * (TMath::Power(1/beta, 2) - 1); //m² = p²*((1/beta)²-1)
531 
532 
533  if(mcTrackIdTofHit == mcTrackIdRing){
534  //cout << "mcTrackIdTofHit: " << mcTrackIdTofHit << endl;
535  //cout << "mcTrackIdRing: " << mcTrackIdRing << endl;
536  Int_t pdg = TMath::Abs(mcTrackRing->GetPdgCode());
537  if(pdg != 50000050){
538  // fHM->H2("fh_radius_mass2")->Fill(mass2, radius);
539  // fHM->H2("fh_radius_beta")->Fill(beta, radius);
540  // fHM->H1("fh_beta_dis_ring")->Fill(beta);
541  // fHM->H2("fh_radius_momentum")->Fill(momTotal, radius);
542 
543  }
544  }
545  }
546  }
547  */
548 }
549 
550 /*
551  if (ring) {
552  DrawEvent();
553  cout << "DrawEvent() called" << endl;
554  }
555  */
556 void CbmRichMCbmQa::DrawHist() {
557  cout.precision(4);
558 
559  SetDefaultDrawStyle();
560  int numberOfEvents = fHM->H1("number_of_events")->GetEntries();
561  // fHM->ScaleByPattern("fh_.*", 1./fEventNum);
562  fHM->ScaleByPattern("fh_.*", 1. / numberOfEvents);
563 
564  {
565  fHM->CreateCanvas("number_of_events", "number_of_events", 600, 600);
566  DrawH1(fHM->H1("number_of_events"));
567  }
568 
569  {
570  fHM->CreateCanvas(
571  "richsp_radius_momentum", "richsp_radius_momentum", 1200, 600);
572  DrawH2(fHM->H2("fh_radius_momentum"));
573  }
574 
575  {
576  fHM->CreateCanvas("richsp_radius_mass2", "richsp_radius_mass2", 1200, 600);
577  DrawH2(fHM->H2("fh_radius_mass2"));
578  }
579 
580  {
581  fHM->CreateCanvas("richsp_radius_beta", "richsp_radius_beta", 1200, 600);
582  DrawH2(fHM->H2("fh_radius_beta"));
583  }
584 
585  {
586  fHM->CreateCanvas("richsp_beta_dis_all", "richsp_beta_dis_all", 1200, 600);
587  DrawH1(fHM->H1("fh_beta_dis_all"));
588  }
589 
590  {
591  fHM->CreateCanvas(
592  "richsp_beta_dis_ring", "richsp_beta_dis_ring", 1200, 600);
593  DrawH1(fHM->H1("fh_beta_dis_ring"));
594  }
595 
596  {
597  fHM->CreateCanvas("fh_nof_rich_rings", "fh_nof_rich_rings", 1200, 600);
598  DrawH1(fHM->H1("fh_nof_rich_rings"));
599  }
600 
601 
602  {
603  fHM->CreateCanvas("richsp_eloss_tof", "richsp_eloss_tof", 1200, 600);
604  DrawH2(fHM->H2("fh_eloss_tof"));
605  }
606 
607  {
608  TCanvas* c = fHM->CreateCanvas(
609  "richsp_nof_rich_hits_points", "richsp_nof_rich_hits_points", 1800, 600);
610  c->Divide(3, 1);
611  c->cd(1);
612  DrawH1(fHM->H1("fh_nof_rich_points"));
613  c->cd(2);
614  DrawH1(fHM->H1("fh_nof_rich_hits"));
615  c->cd(3);
616  DrawH1(fHM->H1("fh_hits_per_ring"));
617  }
618 
619  {
620  TCanvas* c = fHM->CreateCanvas(
621  "richsp_rich_points_hits_xy", "richsp_rich_points_hits_xy", 1200, 600);
622  c->Divide(2, 1);
623  c->cd(1);
624  DrawH2(fHM->H2("fh_rich_points_xy"));
625  c->cd(2);
626  DrawH2(fHM->H2("fh_rich_hits_xy"));
627  }
628 
629  {
630  fHM->CreateCanvas(
631  "richsp_rich_ring_radius", "richsp_rich_ring_radius", 600, 600);
632  DrawH1(fHM->H1("fh_rich_ring_radius"));
633  fHM->H1("fh_rich_ring_radius")->SetTitle("Ring Radius");
634  }
635 
636 
637  {
638  fHM->CreateCanvas("richsp_dR", "richsp_dR", 600, 600);
639  DrawH1andFitGauss(fHM->H1("fh_dR"), true, false);
640  fHM->H1("fh_dR")->SetTitle("dR");
641  }
642 
643  {
644  TCanvas* c = fHM->CreateCanvas(
645  "richsp_ring_center_xy", "richsp_ring_center_xy", 600, 600);
646  c->SetLogz();
647  DrawH2(fHM->H2("fh_ring_center_xy"));
648  fHM->H1("fh_ring_center_xy")->SetTitle("Ring center");
649  }
650 
651  {
652  TCanvas* c = fHM->CreateCanvas("richsp_nonphoton_pmt_points_xy",
653  "richsp_nonphoton_pmt_points_xy",
654  600,
655  600);
656  c->SetLogz();
657  DrawH2(fHM->H2("fh_nonphoton_pmt_points_xy"));
658  fHM->H2("fh_nonphoton_pmt_points_xy")->SetTitle("Non photons on PMT plane");
659  }
660 }
661 
662 void CbmRichMCbmQa::DrawEvent() {
663  cout << "CbmRichMCbmQa::DrawEvent" << endl;
664 
665  stringstream ss;
666  ss << "richsp_event_display_event_" << fEventNum;
667  fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 705, 800);
668 
669  TH2D* pad =
670  new TH2D("event_display_pad", ";X [cm];Y [cm]", 1, -65., 5., 1, -40., 40.);
671 
672  DrawH2(pad);
673  pad->GetYaxis()->SetTitleOffset(0.9);
674  gPad->SetLeftMargin(0.13);
675  gPad->SetRightMargin(0.05);
676 
677  // Draw hits
678  int nofHits = fRichHits->GetEntriesFast();
679  for (int iH = 0; iH < nofHits; iH++) {
680  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(iH);
681  if (NULL == hit) continue;
682  TEllipse* hitDr = new TEllipse(hit->GetX(), hit->GetY(), 0.25);
683  hitDr->SetFillColor(kRed);
684  hitDr->SetLineColor(kRed);
685  hitDr->Draw();
686  }
687 
688  // Draw RICH MC Points
689  int nofPoints = fRichPoints->GetEntriesFast();
690  for (int iP = 0; iP < nofPoints; iP++) {
691  CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(iP);
692  if (NULL == point) continue;
693  TEllipse* pointDr = new TEllipse(point->GetX(), point->GetY(), 0.15);
694  pointDr->SetFillColor(kBlue);
695  pointDr->SetLineColor(kBlue);
696  pointDr->Draw();
697  }
698 
699  /* //Draw proton start XY
700  for( int i = 0; i < fMCTracks->GetEntriesFast(); i++) {
701  CbmMCTrack* mctrack= static_cast<CbmMCTrack*>(fMCTracks->At(i));
702  if (mctrack == NULL) continue;
703  Double_t pdg = mctrack->GetPdgCode();
704  Double_t motherId = mctrack->GetMotherId();
705  if (pdg == 2212 && motherId < 0) {
706  TEllipse* pointDr = new TEllipse(mctrack->GetStartX(), mctrack->GetStartY(), 0.35);
707  pointDr->SetFillColor(kGreen+3);
708  pointDr->SetLineColor(kGreen+3);
709  pointDr->Draw();
710 
711  }
712  }
713  */
714 
715  // Draw rings
716  int nofRings = fRichRings->GetEntriesFast();
717  for (int iR = 0; iR < nofRings; iR++) {
718  CbmRichRing* ring = (CbmRichRing*) fRichRings->At(iR);
719  if (NULL == ring) continue;
720  DrawCircle(ring);
721  }
722 }
723 
724 void CbmRichMCbmQa::DrawCircle(CbmRichRing* ring) {
725  TEllipse* circle =
726  new TEllipse(ring->GetCenterX(), ring->GetCenterY(), ring->GetRadius());
727  circle->SetFillStyle(0);
728  circle->SetLineWidth(4);
729  circle->SetLineColor(kBlack);
730  circle->Draw();
731 
732  TEllipse* center = new TEllipse(ring->GetCenterX(), ring->GetCenterY(), 0.2);
733  center->SetFillColor(kBlack);
734  center->SetLineColor(kBlack);
735  center->Draw();
736 }
737 
738 void CbmRichMCbmQa::Finish() {
739  //DrawHist();
740  //fHM->SaveCanvasToImage(fOutputDir);
741  fHM->WriteToFile();
742 }
743 
744 
745 void CbmRichMCbmQa::DrawFromFile(const string& fileName,
746  const string& outputDir) {
747  fOutputDir = outputDir;
748 
749  if (fHM != nullptr) delete fHM;
750 
751  fHM = new CbmHistManager();
752  TFile* file = new TFile(fileName.c_str());
753  fHM->ReadFromFile(file);
754 
755  DrawHist();
756 
757  fHM->SaveCanvasToImage(fOutputDir);
758 }
759 
760 
761 ClassImp(CbmRichMCbmQa)
CbmRichPoint.h
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition: CbmMCTrack.h:172
CbmRichMCbmQa::fMCTracks
TClonesArray * fMCTracks
Definition: CbmRichMCbmQa.h:97
CbmRichMCbmQa::DrawHist
void DrawHist()
Draw histograms.
Definition: CbmRichMCbmQa.cxx:556
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
CbmHistManager::ScaleByPattern
void ScaleByPattern(const std::string &pattern, Double_t scale)
Scale histograms which name matches specified pattern.
Definition: core/base/CbmHistManager.cxx:221
CbmRichMCbmQa::DrawFromFile
void DrawFromFile(const string &fileName, const string &outputDir)
Draw histogram from file.
Definition: CbmRichMCbmQa.cxx:745
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmRichMCbmQa::fRichDigis
TClonesArray * fRichDigis
Definition: CbmRichMCbmQa.h:99
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmHistManager::ReadFromFile
void ReadFromFile(TFile *file)
Read histograms from file.
Definition: core/base/CbmHistManager.cxx:110
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
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichMCbmQa::fRefPlanePoints
TClonesArray * fRefPlanePoints
Definition: CbmRichMCbmQa.h:103
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition: CbmHistManager.h:190
CbmRichMCbmQa::CbmRichMCbmQa
CbmRichMCbmQa()
Standard constructor.
Definition: CbmRichMCbmQa.cxx:44
CbmGlobalTrack.h
CbmRichMCbmQa::fEventNum
Int_t fEventNum
Definition: CbmRichMCbmQa.h:91
CbmRichRing
Definition: CbmRichRing.h:17
CbmRichRing.h
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichDraw.h
CbmRichMCbmQa::fTofHitMatches
TClonesArray * fTofHitMatches
Definition: CbmRichMCbmQa.h:108
CbmHistManager.h
Histogram manager.
CbmMatchRecoToMC.h
FairTask for matching RECO data to MC.
CbmRichMCbmQa
Definition: CbmRichMCbmQa.h:16
CbmRichMCbmQa::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichMCbmQa.cxx:738
CbmRichGeoManager.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
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmRichMCbmQa::fHM
CbmHistManager * fHM
Definition: CbmRichMCbmQa.h:88
CbmRichMCbmQa::InitHistograms
void InitHistograms()
Initialize histograms.
Definition: CbmRichMCbmQa.cxx:153
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmTrackMatchNew.h
CbmRichMCbmQa::fOutputDir
string fOutputDir
Definition: CbmRichMCbmQa.h:94
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
CbmTofHit::GetR
Double_t GetR() const
Definition: core/data/tof/CbmTofHit.h:94
CbmRichMCbmQa::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichMCbmQa.cxx:280
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmRichMCbmQa.h
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition: CbmHistManager.h:170
CbmRichMCbmQa::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichMCbmQa.cxx:62
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
CbmRichMCbmQa::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmRichMCbmQa.h:102
CbmMCTrack.h
CbmRichMCbmQa::DrawCircle
void DrawCircle(CbmRichRing *ring)
Definition: CbmRichMCbmQa.cxx:724
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmTofPoint.h
CbmRichRing::GetRadius
Float_t GetRadius() const
Definition: CbmRichRing.h:82
CbmRichRing::GetCenterY
Float_t GetCenterY() const
Definition: CbmRichRing.h:81
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmTofPoint
Geometric intersection of a MC track with a TOFb detector.
Definition: CbmTofPoint.h:40
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
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
CbmTrdTrack.h
CbmRichMCbmQa::fRichPoints
TClonesArray * fRichPoints
Definition: CbmRichMCbmQa.h:98
CbmRichMCbmQa::fTofPoints
TClonesArray * fTofPoints
Definition: CbmRichMCbmQa.h:107
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmRichMCbmQa::fRichHits
TClonesArray * fRichHits
Definition: CbmRichMCbmQa.h:100
CbmRichMCbmQa::DrawEvent
void DrawEvent()
Definition: CbmRichMCbmQa.cxx:662
CbmRichRing::GetCenterX
Float_t GetCenterX() const
Definition: CbmRichRing.h:80
CbmRichMCbmQa::fRichRings
TClonesArray * fRichRings
Definition: CbmRichMCbmQa.h:101
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
CbmRichMCbmQa::fTofHits
TClonesArray * fTofHits
Definition: CbmRichMCbmQa.h:106