EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CbmRichUrqmdTest.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CbmRichUrqmdTest.cxx
1 
8 #include "CbmRichUrqmdTest.h"
9 
10 #include "TH1D.h"
11 #include "TH1.h"
12 #include "TCanvas.h"
13 #include "TClonesArray.h"
14 
15 #include "CbmMCTrack.h"
16 #include "FairTrackParam.h"
17 #include "CbmRichHit.h"
18 #include "FairMCPoint.h"
19 #include "CbmDrawHist.h"
20 #include "CbmTrackMatch.h"
21 #include "CbmRichRing.h"
22 #include "CbmRichHit.h"
23 
24 #include "std/utils/CbmLitUtils.h"
25 
26 #include <iostream>
27 #include <string>
28 #include <boost/assign/list_of.hpp>
29 
30 using namespace std;
31 using boost::assign::list_of;
32 
34  : FairTask("CbmRichUrqmdTest"),
35  fOutputDir(""),
36  fRichHits(NULL),
37  fRichRings(NULL),
38  fRichPoints(NULL),
39  fMcTracks(NULL),
40  fRichRingMatches(NULL),
41  fRichProjections(NULL),
42  fCanvas(),
43  fEventNum(0),
44  fMinNofHits(7),
45  fNofHitsInRingMap(),
46  fh_vertex_z(NULL),
47  fh_nof_rings_1hit(NULL),
48  fh_nof_rings_7hits(NULL),
49  fh_nof_rings_prim_1hit(NULL),
50  fh_nof_rings_prim_7hits(NULL),
51  fh_nof_rings_target_1hit(NULL),
52  fh_nof_rings_target_7hits(NULL),
53  fh_secel_mom(NULL),
54  fh_gamma_target_mom(NULL),
55  fh_gamma_nontarget_mom(NULL),
56  fh_pi_mom(NULL),
57  fh_kaon_mom(NULL),
58  fh_mu_mom(NULL),
59  fh_nof_hits_per_event(NULL),
60  fh_hits_xy_u(NULL),
61  fh_hits_xy_d(NULL),
62  fh_hitrate_xy_u(NULL),
63  fh_hitrate_xy_d(NULL),
64  fh_nof_proj_per_event(NULL),
65  fh_proj_xy_u(NULL),
66  fh_proj_xy_d(NULL),
67  fHists()
68 {
69 }
70 
72 {
73 
74 }
75 
77 {
78  cout << "CbmRichUrqmdTest::Init"<<endl;
80  if (NULL == ioman) { Fatal("CbmRichUrqmdTest::Init","RootManager not instantised!"); }
81 
82  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
83  if ( NULL == fRichHits) { Fatal("CbmRichUrqmdTest::Init","No RichHit array!"); }
84 
85  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
86  if ( NULL == fRichRings) { Fatal("CbmRichUrqmdTest::Init","No RichRing array!"); }
87 
88  fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
89  if ( NULL == fRichPoints) { Fatal("CbmRichUrqmdTest::Init","No RichPoint array!"); }
90 
91  fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
92  if ( NULL == fMcTracks) { Fatal("CbmRichUrqmdTest::Init","No MCTrack array!"); }
93 
94  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
95  if ( NULL == fRichRingMatches) { Fatal("CbmRichUrqmdTest::Init","No RichRingMatch array!"); }
96 
97  fRichProjections = (TClonesArray*) ioman->GetObject("RichProjection");
98  if ( NULL == fRichProjections) { Fatal("CbmRichUrqmdTest::Init","No fRichProjections array!"); }
99 
100  InitHistograms();
101 
102  return kSUCCESS;
103 }
104 
106  Option_t* option)
107 {
108  fEventNum++;
109 
110  cout << "CbmRichUrqmdTest, event No. " << fEventNum << endl;
111 
113  NofRings();
114  NofHits();
115  NofProjections();
116  Vertex();
117 
118 
119 }
120 
122 {
123  fh_vertex_z = new TH1D("fh_vertex_z", "fh_vertex_z;z [cm];Number of vertices per event", 200, -1., 200);
124  fh_nof_rings_1hit = new TH1D("fh_nof_rings_1hit", "fh_nof_rings_1hit;Number of detected particles/event;Yield", 100, -.5, 99.5);
125  fh_nof_rings_7hits = new TH1D("fh_nof_rings_7hits", "fh_nof_rings_7hits;Number of detected particles/event;Yield", 100, -.5, 99.5 );
126  fh_nof_rings_prim_1hit = new TH1D("fh_nof_rings_prim_1hit", "fh_nof_rings_prim_1hit;Number of detected particles/event;Yield", 50, -.5, 49.5);
127  fh_nof_rings_prim_7hits = new TH1D("fh_nof_rings_prim_7hits", "fh_nof_rings_prim_7hits;Number of detected particles/event;Yield", 50, -.5, 49.5 );
128  fh_nof_rings_target_1hit = new TH1D("fh_nof_rings_target_1hit", "fh_nof_rings_target_1hit;Number of detected particles/event;Yield", 60, -.5, 59.5);
129  fh_nof_rings_target_7hits = new TH1D("fh_nof_rings_target_7hits", "fh_nof_rings_target_7hits;Number of detected particles/event;Yield", 60, -.5, 59.5 );
130 
131  fh_secel_mom = new TH1D("fh_secel_mom", "fh_secel_mom;p [GeV/c];Number per event", 100, 0., 20);
132  fh_gamma_target_mom = new TH1D("fh_gamma_target_mom", "fh_gamma_target_mom;p [GeV/c];Number per event", 100, 0., 20);
133  fh_gamma_nontarget_mom = new TH1D("fh_gamma_nontarget_mom", "fh_gamma_nontarget_mom;p [GeV/c];Number per event", 100, 0., 20);
134  fh_pi_mom = new TH1D("fh_pi_mom", "fh_pi_mom;p [GeV/c];Number per event", 100, 0., 20);
135  fh_kaon_mom = new TH1D("fh_pi_mom", "fh_pi_mom;p [GeV/c];Number per event", 100, 0., 20);
136  fh_mu_mom = new TH1D("fh_mu_mom", "fh_mu_mom;p [GeV/c];Number per event", 100, 0., 20);
137 
138  fh_nof_hits_per_event = new TH1D("fh_nof_hits_per_event", "fh_nof_hits_per_event;Number of hits per event;Yield", 50, 0, 1500);
139  fh_hits_xy_u = new TH2D("fh_hits_xy_u", "fh_hits_xy_u;x [cm];y [cm];Number of hits/cm^{2}/event", 110, -110, 110, 45, 90, 180);
140  fh_hits_xy_d = new TH2D("fh_hits_xy_d", "fh_hits_xy_d;x [cm];y [cm];Number of hits/cm^{2}/event", 110, -110, 110, 45, -180, -90);
141 
142  // bin size is set to 1.2 cm in order to cover 4 pixels, before drawing must be normalized by 1/4
143  fh_hitrate_xy_u = new TH2D("fh_hitrate_xy_u", "fh_hitrate_xy_u;x [cm];y [cm];Number of hits/pixel/s", 184, -110, 110, 75, 90, 180);
144  fh_hitrate_xy_d = new TH2D("fh_hitrate_xy_d", "fh_hitrate_xy_d;x [cm];y [cm];Number of hits/pixel/s", 184, -110, 110, 75, -180, -90);
145 
146  fh_nof_proj_per_event = new TH1D("fh_nof_proj_per_event", "fh_nof_proj_per_event;Number of tracks per event;Yield", 50, 200, 600);
147  fh_proj_xy_u = new TH2D("fh_proj_xy_u", "fh_proj_xy_u;x [cm];y [cm];Number of tracks/cm^{2}/event", 220, -110, 110, 90, 90, 180);
148  fh_proj_xy_d = new TH2D("fh_proj_xy_d", "fh_proj_xy_d;x [cm];y [cm];Number of tracks/cm^{2}/event", 220, -110, 110, 90, -180, -90);
149 }
150 
152 {
153  fNofHitsInRingMap.clear();
154  Int_t nofRichHits = fRichHits->GetEntriesFast();
155  for (Int_t iHit=0; iHit < nofRichHits; iHit++) {
156  CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
157  if (NULL == hit) continue;
158 
159  Int_t iPoint = hit->GetRefId();
160  if (iPoint < 0) continue;
161 
162  FairMCPoint* point = static_cast<FairMCPoint*>(fRichPoints->At(iPoint));
163  if (NULL == point) continue;
164 
165  Int_t iMCTrack = point->GetTrackID();
166  CbmMCTrack* track = static_cast<CbmMCTrack*>(fMcTracks->At(iMCTrack));
167  if (NULL == track) continue;
168 
169  Int_t iMother = track->GetMotherId();
170  if (iMother == -1) continue;
171 
172  fNofHitsInRingMap[iMother]++;
173  }
174 }
175 
177 {
178  Int_t nofRings = fRichRings->GetEntriesFast();
179  int nRings1hit = 0, nRings7hits = 0;
180  int nRingsPrim1hit = 0, nRingsPrim7hits = 0;
181  int nRingsTarget1hit = 0, nRingsTarget7hits = 0;
182  for (Int_t iR = 0; iR < nofRings; iR++){
183  CbmRichRing *ring = (CbmRichRing*) fRichRings->At(iR);
184  if (NULL == ring) continue;
185  CbmTrackMatch* ringMatch = (CbmTrackMatch*) fRichRingMatches->At(iR);
186  if (NULL == ringMatch) continue;
187  Int_t mcTrackId = ringMatch->GetMCTrackId();
188  if (mcTrackId < 0) continue;
189  CbmMCTrack* mcTrack = (CbmMCTrack*)fMcTracks->At(mcTrackId);
190  if (NULL == mcTrack) continue;
191  Int_t motherId = mcTrack->GetMotherId();
192  Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
193  double mom = mcTrack->GetP();
194  TVector3 vert;
195  mcTrack->GetStartVertex(vert);
196  double dZ = vert.Z();
197 
198  if (motherId == -1 && pdg == 11) continue; // do not calculate embedded electrons
199 
200  int nofHits = ring->GetNofHits();
201  if (nofHits >= 1) nRings1hit++;
202  if (nofHits >= fMinNofHits) nRings7hits++;
203 
204  if (motherId == -1 && nofHits >= 1) nRingsPrim1hit++;
205  if (motherId == -1 && nofHits >= fMinNofHits) nRingsPrim7hits++;
206 
207  if (dZ < 0.1 && nofHits >= 1) nRingsTarget1hit++;
208  if (dZ < 0.1 && nofHits >= fMinNofHits) nRingsTarget7hits++;
209 
210  if (nofHits >= 1) {
211  if (motherId != -1) {
212  int motherPdg;
213  CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
214  if (NULL != mother) motherPdg = mother->GetPdgCode();
215  if (motherId != -1 && pdg == 11 && motherPdg != 22) fh_secel_mom->Fill(mom);
216 
217  if (motherId != -1 && pdg == 11 && motherPdg == 22){
218  if (dZ < 0.1){
219  fh_gamma_target_mom->Fill(mom);
220  } else {
221  fh_gamma_nontarget_mom->Fill(mom);
222  }
223  }
224 
225  }
226  if (pdg == 211) fh_pi_mom->Fill(mom);
227  if (pdg == 321) fh_kaon_mom->Fill(mom);
228  if (pdg == 13) fh_mu_mom->Fill(mom);
229  }
230  }
231  fh_nof_rings_1hit->Fill(nRings1hit);
232  fh_nof_rings_7hits->Fill(nRings7hits);
233 
234  fh_nof_rings_prim_1hit->Fill(nRingsPrim1hit);
235  fh_nof_rings_prim_7hits->Fill(nRingsPrim7hits);
236 
237  fh_nof_rings_target_1hit->Fill(nRingsTarget1hit);
238  fh_nof_rings_target_7hits->Fill(nRingsTarget7hits);
239 }
240 
242 {
243  int nofHits = fRichHits->GetEntriesFast();
244  int nofNoiseHits = 0;
245  int nofHitsUrqmd = 0;
246  int nofHitsEl = 0;
247  for (int i = 0; i < nofHits; i++) {
248  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(i);
249  if (NULL == hit) continue;
250 
251  Int_t iPoint = hit->GetRefId();
252  if (iPoint < 0){
253  nofNoiseHits++;
254  continue;
255  }
256 
257  FairMCPoint* point = static_cast<FairMCPoint*>(fRichPoints->At(iPoint));
258  if (NULL == point) continue;
259 
260  Int_t iMCTrack = point->GetTrackID();
261  CbmMCTrack* track = static_cast<CbmMCTrack*>(fMcTracks->At(iMCTrack));
262  if (NULL == track) continue;
263 
264  Int_t iMother = track->GetMotherId();
265  if (iMother == -1) continue;
266 
267  CbmMCTrack* track2 = static_cast<CbmMCTrack*>(fMcTracks->At(iMother));
268  if (NULL == track2) continue;
269  int pdg = TMath::Abs(track2->GetPdgCode());
270  int motherId = track2->GetMotherId();
271  if (motherId == -1 && pdg == 11){
272  nofHitsEl++;
273  continue;
274  }
275  nofHitsUrqmd++;
276  double x = hit->GetX();
277  double y = hit->GetY();
278  if (y > 0) {
279  // cout << x << " " << y << endl;
280  fh_hits_xy_u->Fill(x, y);
281  fh_hitrate_xy_u->Fill(x, y);
282  } else {
283  fh_hits_xy_d->Fill(x, y);
284  fh_hitrate_xy_d->Fill(x, y);
285  }
286  }
287 
288  fh_nof_hits_per_event->Fill(nofHitsUrqmd);
289 
290  cout << "nofHits = " << nofHits << endl;
291  cout << "nofNoiseHits = " << nofNoiseHits << endl;
292  cout << "nofHitsUrqmd = " << nofHitsUrqmd << endl;
293  cout << "nofHitsEl = " << nofHitsEl << endl;
294 }
295 
297 {
298  int nofProj = fRichProjections->GetEntriesFast();
299  int nofGoodProj = 0;
300  for (int i = 0; i < nofProj; i++){
302  if (NULL == proj) continue;
303  double x = proj->GetX();
304  double y = proj->GetY();
305  if (y > 0) {
306  fh_proj_xy_u->Fill(x, y);
307  } else {
308  fh_proj_xy_d->Fill(x, y);
309  }
310  if (proj->GetX() != 0 && proj->GetY() != 0) nofGoodProj++;
311  }
312 
313  fh_nof_proj_per_event->Fill(nofGoodProj);
314 
315 }
316 
318 {
319  int nMcTracks = fMcTracks->GetEntries();
320  for (int i = 0; i < nMcTracks; i++){
321  //At least one hit in RICH
322  if (fNofHitsInRingMap[i] < 1) continue;
323  CbmMCTrack* mctrack = (CbmMCTrack*) fMcTracks->At(i);
324  TVector3 v;
325  mctrack->GetStartVertex(v);
326  fh_vertex_z->Fill(v.Z());
327  } // nMcTracks
328 }
329 
331 {
332  cout.precision(4);
333 
335  fh_vertex_z->Scale(1./fEventNum);
336  TCanvas* cVertex = CreateCanvas("rich_urqmd_vertex_z", "rich_urqmd_vertex_z", 800, 800);
338  gPad->SetLogy(true);
339 
340  TCanvas* c2 = CreateCanvas("rich_urqmd_nof_rings", "rich_urqmd_nof_rings", 800, 800);
341  fh_nof_rings_1hit->Scale(1./fh_nof_rings_1hit->Integral());
342  fh_nof_rings_7hits->Scale(1./fh_nof_rings_7hits->Integral());
343  DrawH1(list_of(fh_nof_rings_1hit)(fh_nof_rings_7hits), list_of(string("At least 1 hit detected"))(string("At least 7 hits detected")),
344  kLinear, kLinear, true, 0.4, 0.78, 0.99, 0.99);
345  cout << "Mean nof rings per event (1 hit) = " << fh_nof_rings_1hit->GetMean() << endl;
346  cout << "Mean nof rings per event (7 hits) = " << fh_nof_rings_7hits->GetMean() << endl;
347 
348  TCanvas* c3 = CreateCanvas("rich_urqmd_nof_rings_prim", "rich_urqmd_nof_rings_prim", 800, 800);
349  fh_nof_rings_prim_1hit->Scale(1./fh_nof_rings_prim_1hit->Integral());
350  fh_nof_rings_prim_7hits->Scale(1./fh_nof_rings_prim_7hits->Integral());
351  DrawH1(list_of(fh_nof_rings_prim_1hit)(fh_nof_rings_prim_7hits), list_of("At least 1 hit detected")("At least 7 hits detected"),
352  kLinear, kLinear, true, 0.4, 0.78, 0.99, 0.99);
353  cout << "Mean nof primary rings per event (1 hit) = " << fh_nof_rings_prim_1hit->GetMean() << endl;
354  cout << "Mean nof primary rings per event (7 hits) = " << fh_nof_rings_prim_7hits->GetMean() << endl;
355 
356  TCanvas* cTarget = CreateCanvas("rich_urqmd_nof_rings_target", "rich_urqmd_nof_rings_target", 800, 800);
357  fh_nof_rings_target_1hit->Scale(1./fh_nof_rings_target_1hit->Integral());
359  DrawH1(list_of(fh_nof_rings_target_1hit)(fh_nof_rings_target_7hits), list_of("At least 1 hit detected")("At least 7 hits detected"),
360  kLinear, kLinear, true, 0.4, 0.78, 0.99, 0.99);
361  cout << "Mean nof target rings per event (1 hit) = " << fh_nof_rings_target_1hit->GetMean() << endl;
362  cout << "Mean nof target rings per event (7 hits) = " << fh_nof_rings_target_7hits->GetMean() << endl;
363 
364 
365  TCanvas* c4 = CreateCanvas("rich_urqmd_sources_mom", "rich_urqmd_sources_mom", 800, 800);
366  fh_gamma_target_mom->Scale(1./fEventNum);
368  fh_secel_mom->Scale(1./fEventNum);
369  fh_pi_mom->Scale(1./fEventNum);
370  fh_kaon_mom->Scale(1./fEventNum);
371  fh_mu_mom->Scale(1./fEventNum);
373  list_of("e^{#pm}_{target} from #gamma")("e^{#pm}_{not target} from #gamma")("e^{#pm}_{sec} other")("#pi^{#pm}")("K^{#pm}")("#mu^{#pm}"),
374  kLinear, kLog, true, 0.6, 0.7, 0.99, 0.99);
375 
376  cout << "Nof sec electrons from gamma target per event = " << fh_gamma_target_mom->GetEntries() / fEventNum << endl;
377  cout << "Nof sec electrons from gamma NOT target per event = " << fh_gamma_nontarget_mom->GetEntries() / fEventNum << endl;
378  cout << "Nof sec electrons other per event = " << fh_secel_mom->GetEntries() / fEventNum << endl;
379  cout << "Nof pions per event = " << fh_pi_mom->GetEntries() / fEventNum << endl;
380  cout << "Nof kaons per event = " << fh_kaon_mom->GetEntries() / fEventNum << endl;
381  cout << "Nof muons per event = " << fh_mu_mom->GetEntries() / fEventNum << endl;
382 
383 
384  TCanvas *c5 = CreateCanvas("rich_urqmd_hits_xy", "rich_urqmd_hits_xy", 800, 800);
385  double binArea = fh_hits_xy_u->GetXaxis()->GetBinWidth(1) * fh_hits_xy_u->GetYaxis()->GetBinWidth(1);
386  cout << "binArea = " << binArea << endl;
387  fh_hits_xy_u->Scale(1./(fEventNum * binArea));
388  fh_hits_xy_d->Scale(1./(fEventNum * binArea));
389 
390  c5->Divide(1, 2);
391  c5->cd(1);
393  fh_hits_xy_u->GetYaxis()->SetTitleOffset(0.75);
394  fh_hits_xy_u->GetZaxis()->SetTitleOffset(0.75);
395  gPad->SetLeftMargin(0.1);
396  gPad->SetRightMargin(0.15);
397  c5->cd(2);
399  fh_hits_xy_d->GetYaxis()->SetTitleOffset(0.75);
400  fh_hits_xy_d->GetZaxis()->SetTitleOffset(0.75);
401  gPad->SetLeftMargin(0.1);
402  gPad->SetRightMargin(0.15);
403 
404  TCanvas* c6 = CreateCanvas("rich_urqmd_nof_hits_per_event", "rich_urqmd_nof_hits_per_event", 800, 800);
405  fh_nof_hits_per_event->Scale(1./fh_nof_hits_per_event->Integral());
407  cout << "Mean number of hits per event = " << fh_nof_hits_per_event->GetMean() << endl;
408 
409  TCanvas *cHitRate = CreateCanvas("rich_urqmd_hitrate_xy", "rich_urqmd_hitrate_xy", 800, 800);
410  fh_hitrate_xy_u->Scale(1e7/(fEventNum * 4.));
411  fh_hitrate_xy_d->Scale(1e7/(fEventNum * 4.));
412  cHitRate->Divide(1, 2);
413  cHitRate->cd(1);
415  fh_hitrate_xy_u->GetYaxis()->SetTitleOffset(0.75);
416  fh_hitrate_xy_u->GetZaxis()->SetTitleOffset(0.87);
417  gPad->SetLeftMargin(0.1);
418  gPad->SetRightMargin(0.15);
419  cHitRate->cd(2);
421  fh_hitrate_xy_d->GetYaxis()->SetTitleOffset(0.75);
422  fh_hitrate_xy_d->GetZaxis()->SetTitleOffset(0.87);
423  gPad->SetLeftMargin(0.1);
424  gPad->SetRightMargin(0.15);
425 
426  TCanvas *c7 = CreateCanvas("rich_urqmd_proj_xy", "rich_urqmd_proj_xy", 800, 800);
427  fh_proj_xy_u->Scale(1./fEventNum);
428  fh_proj_xy_d->Scale(1./fEventNum);
429  c7->Divide(1, 2);
430  c7->cd(1);
432  fh_proj_xy_u->GetYaxis()->SetTitleOffset(0.75);
433  fh_proj_xy_u->GetZaxis()->SetTitleOffset(0.75);
434  gPad->SetLeftMargin(0.1);
435  gPad->SetRightMargin(0.15);
436  c7->cd(2);
438  fh_proj_xy_d->GetYaxis()->SetTitleOffset(0.75);
439  fh_proj_xy_d->GetZaxis()->SetTitleOffset(0.75);
440  gPad->SetLeftMargin(0.1);
441  gPad->SetRightMargin(0.15);
442 
443  TCanvas* c8 = CreateCanvas("rich_urqmd_nof_proj_per_event", "rich_urqmd_nof_proj_per_event", 800, 800);
444  fh_nof_proj_per_event->Scale(1./fEventNum);
446  cout << "Number of track projections per event = " << fh_nof_proj_per_event->GetMean() << endl;
447 }
448 
450 {
451  DrawHist();
452  for (Int_t i = 0; i < fHists.size(); i++){
453  fHists[i]->Write();
454  }
456 }
457 
458 
460  const string& name,
461  const string& title,
462  int width,
463  int height)
464 {
465  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), width, height);
466  fCanvas.push_back(c);
467  return c;
468 }
469 
471 {
472  for (int i = 0; i < fCanvas.size(); i++)
473  {
475  }
476 }
477 
479