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