CbmRoot
CbmL1RichRingQa.cxx
Go to the documentation of this file.
1 #include "CbmL1RichRingQa.h"
2 
3 #include "CbmL1Def.h"
4 
5 #include "CbmMCTrack.h"
6 #include "CbmRichHit.h"
7 #include "CbmRichPoint.h"
8 #include "CbmRichRing.h"
9 #include "CbmRichRingFitterEllipseTau.h"
10 #include "CbmRichRingLight.h"
11 #include "FairRootManager.h"
12 #include "FairTask.h"
13 
14 #include "TArc.h"
15 #include "TClonesArray.h"
16 #include "TEllipse.h"
17 #include "TFile.h"
18 #include "TH1.h"
19 #include "TH1F.h"
20 #include "TH2.h"
21 #include "TH2F.h"
22 #include "TLatex.h"
23 #include "TProfile.h"
24 #include "TROOT.h"
25 #include "TStyle.h"
26 #include "TText.h"
27 
28 
29 #include <cmath>
30 #include <iostream>
31 #include <list>
32 #include <map>
33 #include <vector>
34 
35 #include <ctype.h>
36 #include <sstream>
37 #include <stdio.h>
38 
39 using namespace std;
40 
41 using std::cout;
42 using std::endl;
43 using std::fabs;
44 using std::ios;
45 using std::list;
46 using std::map;
47 using std::pair;
48 using std::vector;
49 
50 ClassImp(CbmL1RichRingQa)
51 
52 #define DRAW
53 
54  //------------ standard constructor (with verbosity level) ---------------------------------
55  CbmL1RichRingQa::CbmL1RichRingQa(const char* name, const char*, Int_t verbose)
56  : FairTask(name)
57  , fRingArray(0)
58  , // Array of CbmRichRings
59  fMCPointArray(0)
60  , // Array of FairMCPoints
61  fMCTrackArray(0)
62  , // Array of CbmMCTracks
63  fHitArray(0)
64  , // Array of CbmRichHits
65 
66  Chi2Ghost(0)
67  , Chi2Ref(0)
68  , Chi2All(0)
69  , Chi2Clone(0)
70  , Chi2NhitsGhost(0)
71  , Chi2NhitsAll(0)
72  , RGhost(0)
73  , REl(0)
74  , RPi(0)
75  , NHitsMC(0)
76  , NSameHits(0)
77  ,
78 
79  //Chi2NhitsGhost(0),
80  Chi2NhitsPi(0)
81  , Chi2NhitsEll(0)
82  , RNhitsGhost(0)
83  , RNhitsPi(0)
84  , RNhitsEll(0)
85  , RChi2Ghost(0)
86  , RChi2Pi(0)
87  , RChi2Ell(0)
88  , NSameHitsVsP(0)
89  , NHitsVsMCP(0)
90  , RadiusVsPForClone(0)
91  , DistanceVsPClone(0)
92  , Chi2VsPClone(0)
93  , RadiusVsDistanceClone(0)
94  , NHitsRecoVsNHitsMC(0)
95 
96 
97 {
98  Chi2Ghost = new TH1F("Chi2 for Ghost", "Chi2 for Ghost", 100, 0.f, 1.f);
99  Chi2Ghost->SetXTitle("Chi2 Ghost");
100  Chi2Ghost->SetYTitle("Enteric");
101 
102  Chi2Ref = new TH1F("Ref Chi2", "Ref Chi2", 100, 0.f, 1.f);
103  Chi2Ref->SetXTitle("Chi2 Ref");
104  Chi2Ref->SetYTitle("Enteric");
105 
106  Chi2All = new TH1F("Chi2 for reconstructed rings",
107  "Chi2 for reconstructed rings",
108  100,
109  0.f,
110  1.f);
111  Chi2All->SetXTitle("Chi2 reco rings");
112  Chi2All->SetYTitle("Enteric");
113 
114  Chi2Clone = new TH1F("Clone Chi2", "Clone Chi2", 100, 0.f, 1.f);
115  Chi2Clone->SetXTitle("Chi2 Clone");
116  Chi2Clone->SetYTitle("Enteric");
117 
118  Chi2NhitsGhost = new TH2F(
119  "Chi2 Ghost vs Nhits", "Chi2 vs Nhits", 30, 4.5f, 34.5f, 100, 0.f, 1.f);
120  Chi2NhitsGhost->SetXTitle("Ghost Nhits");
121  Chi2NhitsGhost->SetYTitle("Ghost Chi2");
122 
123  Chi2NhitsAll = new TH2F(
124  "Chi2 All vs Nhits", "Chi2 vs Nhits", 30, 4.5f, 34.5f, 100, 0.f, 1.f);
125  Chi2NhitsAll->SetXTitle("Nhits All");
126  Chi2NhitsAll->SetYTitle("Chi2 All");
127 
128  REl = new TH1F("REl", "REl", 100, 0.f, 7.f);
129  RPi = new TH1F("RPi", "RPi", 100, 0.f, 7.f);
130  RGhost = new TH1F("RGhost", "RGhost", 100, 0.f, 7.f);
131  // Chi2NhitsGhost = new TH2F("Chi2 Ghost vs Nhits", "Chi2 Ghost vs Nhits", 30, 4.5f, 34.5f,100, 0.f, 1.f);
132  Chi2NhitsPi = new TH2F(
133  "Chi2 Pi vs Nhits", "Chi2 Pi vs Nhits", 30, 4.5f, 34.5f, 100, 0.f, 1.f);
134  Chi2NhitsPi->SetXTitle("Nhits Pi");
135  Chi2NhitsPi->SetYTitle("Chi2 Pi");
136 
137  Chi2NhitsEll = new TH2F(
138  "Chi2 Ell vs Nhits", "Chi2 Ell vs Nhits", 30, 4.5f, 34.5f, 100, 0.f, 1.f);
139  Chi2NhitsEll->SetXTitle("Nhits Ell");
140  Chi2NhitsEll->SetYTitle("Chi2 Ell");
141 
142  RNhitsGhost = new TH2F(
143  "R Ghost vs Nhits", "R Ghost vs Nhits", 30, 4.5f, 34.5f, 100, 2.f, 8.f);
144  RNhitsGhost->SetXTitle("Nhits Ghost");
145  RNhitsGhost->SetYTitle("R Ghost");
146 
147  RNhitsPi =
148  new TH2F("R Pi vs Nhits", "R Pi vs Nhits", 30, 4.5f, 34.5f, 100, 2.f, 8.f);
149  RNhitsPi->SetXTitle("Nhits Pi");
150  RNhitsPi->SetYTitle("R Pi");
151 
152  RNhitsEll = new TH2F(
153  "R Ell vs Nhits", "R Ell vs Nhits", 30, 4.5f, 34.5f, 100, 2.f, 8.f);
154  RNhitsEll->SetXTitle("Nhits electrons");
155  RNhitsEll->SetYTitle("R electrons");
156 
157  RChi2Ghost = new TH2F(
158  "R Ghost vs Chi2", "R Ghost vs Chi2", 100, 0.f, 1.f, 100, 2.f, 8.f);
159  RChi2Ghost->SetXTitle("Chi2 Ghost");
160  RChi2Ghost->SetYTitle("R Ghost");
161 
162  RChi2Pi =
163  new TH2F("R Pi vs Chi2", "R Pi vs Chi2", 100, 0.f, 1.f, 100, 2.f, 8.f);
164  RChi2Pi->SetXTitle("Chi2 Pi");
165  RChi2Pi->SetYTitle("R Pi");
166 
167  RChi2Ell =
168  new TH2F("R Ell vs Chi2", "R Ell vs Chi2", 100, 0.f, 1.f, 100, 2.f, 8.f);
169  RChi2Ell->SetXTitle("Chi2 Electrons");
170  RChi2Ell->SetYTitle("R Electrons");
171 
172  NHitsMC = new TH1F("MC Hits %", "MC Hits %", 100, 0.f, 1.f);
173  NHitsMC->SetXTitle("MC hits found, %");
174  NHitsMC->SetYTitle("Entries");
175 
176  NSameHits = new TH1F("Same Hits", "Same Hits", 100, 0.f, 30.f);
177  NSameHits->SetXTitle("N Same Hits");
178  NSameHits->SetYTitle("Enteric");
179 
180  NSameHitsVsP = new TH2F(
181  "MC Hits % vs P(MC)", "MC Hits % vs P(MC)", 100, 0.f, 10.f, 100, 0.f, 1.f);
182  NSameHitsVsP->SetXTitle("P MC");
183  NSameHitsVsP->SetYTitle("MC Hits, %");
184 
185  NHitsVsMCP = new TH2F(
186  "Same Hits vs P(MC)", "Same Hits vs P(MC)", 100, 0.f, 20.f, 100, 0.f, 30.f);
187  NHitsVsMCP->SetXTitle("P MC");
188  NHitsVsMCP->SetYTitle("Same Hits");
189 
190  RadiusVsPForClone =
191  new TH2F(" Radius Vs P ", " Radius Vs P ", 100, 0.f, 12.f, 100, 0.f, 3.f);
192  RadiusVsPForClone->SetXTitle("P Clone");
193  RadiusVsPForClone->SetYTitle("R1-R2");
194 
195  DistanceVsPClone = new TH2F(
196  " Distance Vs P ", " Distance Vs P ", 100, 0.f, 12.f, 100, 0.f, 5.f);
197  DistanceVsPClone->SetXTitle("P Clone");
198  DistanceVsPClone->SetYTitle("Distance");
199 
200  Chi2VsPClone =
201  new TH2F(" Chi2 Vs P ", " Chi2 Vs P ", 100, 0.f, 12.f, 100, 0.f, 1.f);
202  Chi2VsPClone->SetXTitle("P Clone");
203  Chi2VsPClone->SetYTitle("Chi2");
204 
205  RadiusVsDistanceClone = new TH2F(" Radius Vs Distance ",
206  " Radius Vs Distance ",
207  100,
208  0.f,
209  5.f,
210  100,
211  0.f,
212  3.f);
213  RadiusVsDistanceClone->SetXTitle("Distance");
214  RadiusVsDistanceClone->SetYTitle("R1-R2");
215 
216  NHitsRecoVsNHitsMC = new TH2F(" N Hits Reco Vs N Hits MC ",
217  " N Hits Reco Vs N Hits MC ",
218  100,
219  0.f,
220  35.f,
221  100,
222  0.f,
223  35.f);
224  NHitsRecoVsNHitsMC->SetXTitle("N Hits MC");
225  NHitsRecoVsNHitsMC->SetYTitle("N Hits Reco");
226 
227  fVerbose = verbose;
228 }
229 
230 // ----- Destructor ----------------------------------------------------
231 CbmL1RichRingQa::~CbmL1RichRingQa() {}
232 
233 InitStatus CbmL1RichRingQa::Init() {
234 
235  // Get and check FairRootManager
236  FairRootManager* ioman = FairRootManager::Instance();
237  if (!ioman) {
238  cout << "-E- CbmL1RichRingQa::Init: "
239  << "RootManager not instantised!" << endl;
240  return kERROR;
241  }
242 
243  // Get hit Array
244  fHitArray = L1_DYNAMIC_CAST<TClonesArray*>(ioman->GetObject("RichHit"));
245  if (!fHitArray) {
246  cout << "-W- CbmL1RichRingQa::Init: No RichHit array!" << endl;
247  }
248 
249  // Get RichRing Array
250  fRingArray = L1_DYNAMIC_CAST<TClonesArray*>(ioman->GetObject("RichRing"));
251  if (!fRingArray) {
252  cout << "-E- CbmL1RichRingQa::Init: No RichRing array!" << endl;
253  return kERROR;
254  }
255 
256  // Get MC Point array
257  fMCPointArray = L1_DYNAMIC_CAST<TClonesArray*>(ioman->GetObject("RichPoint"));
258  if (!fMCPointArray) {
259  cout << "-E- CbmL1RichRingQa::Init: No RichPoint array!" << endl;
260  return kERROR;
261  }
262 
263  // Get MC Point array
264  fMCTrackArray = L1_DYNAMIC_CAST<TClonesArray*>(ioman->GetObject("MCTrack"));
265  if (!fMCTrackArray) {
266  cout << "-E- CbmL1RichRingQa::Init: No MCTrack array!" << endl;
267  return kERROR;
268  }
269 
270  return kSUCCESS;
271 }
272 
273 void CbmL1RichRingQa::CirFit(list<pair<Double_t, Double_t>>& P,
274  Double_t* X,
275  Double_t* Y,
276  Double_t* R) {
277  Double_t S[8] = {0, 0, 0, 0, 0, 0, 0, 0};
278  Int_t n = 0;
279  for (list<pair<Double_t, Double_t>>::iterator i = P.begin(); i != P.end();
280  ++i) {
281  Double_t& x = i->first;
282  Double_t& y = i->second;
283  Double_t r2 = x * x + y * y;
284  S[0] += x * x;
285  S[1] += y * y;
286  S[2] += x * y;
287  S[3] += x * r2;
288  S[4] += y * r2;
289  S[5] += x;
290  S[6] += y;
291  S[7] += r2;
292  n++;
293  }
294  Double_t s0 = S[6] * S[0] - S[2] * S[5];
295  Double_t s1 = S[0] * S[1] - S[2] * S[2];
296  Double_t s2 = S[0] * S[4] - S[2] * S[3];
297  *X = *Y = *R = 0.;
298  if (fabs(s0) < 1.E-6 || fabs(s1) < 1.E-6) return;
299  Double_t tmp = s1 * (S[5] * S[5] - n * S[0]) + s0 * s0;
300  Double_t A = ((S[0] * S[7] - S[3] * S[5]) * s1 - s2 * s0) / tmp;
301  *Y = (s2 + s0 * A) / s1 / 2;
302  *X = (S[3] + S[5] * A - S[2] * (*Y) * 2) / S[0] / 2;
303  Double_t R2 = (*X) * (*X) + (*Y) * (*Y) - A;
304  if (R2 < 0) return;
305  *R = sqrt(R2);
306 }
307 
308 void CbmL1RichRingQa::Exec(Option_t* /*option*/) {
309  // histogramms
310  //anagai
311  map<void*, int> pHitIndex;
312  map<Int_t, MCRing> MCRingMap;
313 
314  Int_t NRecoEx = 0, NRecoG = 0, NGostmy = 0, NMC = 0, NMC37 = 0;
315 
316  static Int_t niv = 0;
317  Double_t ArrRingX[fRingArray->GetEntriesFast()];
318  Double_t ArrRingY[fRingArray->GetEntriesFast()];
319  Double_t ArrRingR[fRingArray->GetEntriesFast()];
320 
321  Int_t my_NHits[fRingArray->GetEntriesFast()];
322 
323 
324  Double_t HitArrX[fHitArray->GetEntriesFast()];
325  Double_t HitArrY[fHitArray->GetEntriesFast()];
326  Double_t NoiseHitArrX[fHitArray->GetEntriesFast()];
327  Double_t NoiseHitArrY[fHitArray->GetEntriesFast()];
328 
329  Int_t NpiMC = 0, NeMC = 0;
330  niv++;
331  std::string s;
332  std::stringstream out;
333  out << niv;
334  s = out.str();
335 
336 #ifdef DRAW
337  TCanvas* Up = new TCanvas("Up", "Up", 0, 0, 240, 180);
338  Up->SetFillColor(kWhite);
339 
340  Up->Divide(1, 2);
341  Up->cd(1);
342  Up->Range(-110, 110, 110, 180);
343  gPad->DrawFrame(-110, 110, 110, 180);
344  Up->cd(2);
345  Up->Range(-110, -180, 110, -110);
346  gPad->DrawFrame(-110, -180, 110, -110);
347 #endif //DRAW
348  Int_t NfakeHits = 0;
349  for (Int_t hit = 0; hit < fHitArray->GetEntriesFast(); hit++) {
350 
351  CbmRichHit* phit = L1_DYNAMIC_CAST<CbmRichHit*>(fHitArray->At(hit));
352  HitArrX[hit] = (phit->GetX());
353  HitArrY[hit] = -1 * (phit->GetY());
354  if (phit->GetRefId() == -1) {
355  NoiseHitArrX[NfakeHits] = (phit->GetX());
356  NoiseHitArrY[NfakeHits] = -1 * (phit->GetY());
357  NfakeHits++;
358  }
359  }
360  for (Int_t i = 0; i < fRingArray->GetEntriesFast(); i++) {
361  CbmRichRing* ring = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(i));
362  //CbmRichRingLight *my_ring = (CbmRichRingLight*) fRingArray->At( i );
363  ArrRingX[i] = ring->GetCenterX();
364  ArrRingY[i] = -1 * (ring->GetCenterY());
365  ArrRingR[i] = ring->GetRadius();
366  my_NHits[i] = 0;
367  my_NHits[i] = ring->GetNofHits();
368  TString st = "";
369  st += my_NHits[i];
370  NGostmy++;
371 #ifdef DRAW
372  //Draw All Rings
373  Up->cd(1);
374  TArc* AllRingsUp = new TArc(ArrRingX[i], ArrRingY[i], ArrRingR[i], 0, 360);
375  AllRingsUp->SetLineWidth(1);
376  AllRingsUp->SetLineColor(1);
377  AllRingsUp->SetFillStyle(0);
378  AllRingsUp->Draw("*");
379  // Draw N Hits for rings
380  /*
381  TText* textn = new TText( ArrRingX[i] - ArrRingR[i]/5 , ArrRingY[i] + ArrRingR[i]/5 , st.Data() );
382  textn->SetTextAlign(12);
383  textn->SetTextSize(0.01);
384  textn->Draw();
385  */
386  Up->cd(2);
387  TArc* AllRingsDown =
388  new TArc(ArrRingX[i], ArrRingY[i], ArrRingR[i], 0, 360);
389  AllRingsDown->SetLineWidth(1);
390  AllRingsDown->SetLineColor(1);
391  AllRingsDown->SetFillStyle(0);
392  AllRingsDown->Draw("*");
393  // Draw N Hits for rings
394  /*
395  TText* text1 = new TText( ArrRingX[i] - ArrRingR[i]/5 , ArrRingY[i] + ArrRingR[i]/5 , st.Data() );
396  text1->SetTextAlign(12);
397  text1->SetTextSize(0.01);
398  text1->Draw();
399 */
400 #endif // DRAW
401  }
402 
403  //anagai
404 
405  static TH1F *h_MC_radius, *h_MC_nhits, *h_MC_primary_nhits, *h_MC_momentum,
406  *h_MC_primary_momentum, *h_MC_resolution, *h_MC_ref_resolution,
407  *h_MC_extra_resolution, *h_ghost_nhits;
408 
409  static TH2F* h_MC_primary_res_vs_momentum;
410 
411  static TProfile *p_ref_eff_vs_nhits, *p_extra_eff_vs_nhits;
412 
413  static TList* listHisto;
414  static bool first_call_performance = 1;
415 
416  if (first_call_performance) {
417  first_call_performance = 0;
418  TDirectory* curdir = gDirectory;
419  TDirectory* histodir = gROOT->mkdir("CbmL1RichRingQaHisto");
420  histodir->cd();
421 
422  h_MC_radius = new TH1F("h_MC_radius", "MC ring radius (cm)", 100, 0.0, 10.);
423  h_MC_nhits = new TH1F("h_MC_nhits", "Hits per MC ring", 50, 0.0, 50.);
424  h_MC_primary_nhits =
425  new TH1F("h_MC_primary_nhits", "Hits per primary MC ring", 50, 0.0, 50.);
426  h_MC_momentum =
427  new TH1F("h_MC_momentum", "MC track momentum (GeV)", 100, 0.0, 15.);
428  h_MC_primary_momentum = new TH1F("h_MC_primary_momentum",
429  "MC primary track momentum (GeV)",
430  100,
431  0.0,
432  15.);
433  h_MC_resolution = new TH1F(
434  "h_MC_resolution", "Hit deviation from MC ring (cm)", 500, -5.0, 5.0);
435  h_MC_ref_resolution = new TH1F("h_MC_ref_resolution",
436  "Hit deviation from REF MC ring (cm)",
437  500,
438  -5.0,
439  5.0);
440  h_MC_extra_resolution = new TH1F("h_MC_extra_resolution",
441  "Hit deviation from EXTRA MC ring (cm)",
442  500,
443  -5.0,
444  5.0);
445 
446  h_ghost_nhits =
447  new TH1F("h_ghost_nhits", "Hits per ghost ring", 50, 0.0, 50.);
448 
449  h_MC_primary_res_vs_momentum =
450  new TH2F("h_MC_primary_res_vs_momentum",
451  "Hit deviation from ptimary MC ring (cm) vs P",
452  100,
453  0.,
454  15.,
455  500,
456  -5.0,
457  5.0);
458 
459  p_ref_eff_vs_nhits = new TProfile("p_ref_eff_vs_nhits",
460  "Refset efficiency vs N Hits",
461  100,
462  0.0,
463  50.0,
464  0.0,
465  1.0);
466 
467  p_extra_eff_vs_nhits = new TProfile("p_extra_eff_vs_nhits",
468  "Extraset efficiency vs N Hits",
469  100,
470  0.0,
471  50.0,
472  0.0,
473  1.0);
474 
475  // ----- Create list of all histogram pointers
476  listHisto = gDirectory->GetList();
477 
478  curdir->cd();
479  }
480 
481  // Create hit vector
482 
483  if (!fHitArray || !fMCTrackArray || !fMCPointArray || !fRingArray) return;
484 
485  int NHits = fHitArray->GetEntriesFast();
486  PerfHit Hits[NHits];
487 
488  // map < void*, int > pHitIndex;
489 
490  for (Int_t i = 0; i < NHits; i++) {
491 
492  PerfHit& hit = Hits[i];
493  hit.index = i;
494  hit.x = 0;
495  hit.y = 0;
496  hit.MCTrackID = -1;
497 
498  CbmRichHit* phit = L1_DYNAMIC_CAST<CbmRichHit*>(fHitArray->At(i));
499  if (!phit) continue;
500  pHitIndex.insert(pair<void*, int>(phit, i));
501  hit.x = phit->GetX();
502  hit.y = phit->GetY();
503  Int_t pointID = phit->GetRefId();
504  if (pointID < 0) continue;
505  CbmRichPoint* point =
506  L1_DYNAMIC_CAST<CbmRichPoint*>(fMCPointArray->At(pointID));
507  if (!point) continue;
508  Int_t trackID = point->GetTrackID();
509  if (trackID < 0) continue;
510  CbmMCTrack* track =
511  L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTrackArray->At(trackID));
512  if (!track || track->GetPdgCode() != 50000050)
513  continue; // select only Cherenkov photons
514  hit.MCTrackID = track->GetMotherId();
515  }
516 
517  // Create map of MC rings
518 
519  //map < Int_t, MCRing> MCRingMap;
520 
521  // match hits
522 
523  for (int ih = 0; ih < NHits; ih++) {
524  int ID = Hits[ih].MCTrackID;
525  if (ID >= 0 && MCRingMap.find(ID) == MCRingMap.end()) {
526  MCRing tmp;
527  tmp.NHits = 0;
528  MCRingMap.insert(pair<Int_t, MCRing>(ID, tmp));
529  }
530  MCRingMap[ID].NHits++;
531  MCRingMap[ID].Hits.push_back(ih);
532  MCRingMap[ID].NHitsBestReco = 0;
533  MCRingMap[ID].BestReco = 0;
534  MCRingMap[ID].NHitsBestvsNHitsMC = 0.;
535  }
536 
537  // fit MC rings & set parameters
538  Int_t NMCRings = 0;
539  for (map<Int_t, MCRing>::iterator i = MCRingMap.begin(); i != MCRingMap.end();
540  ++i) {
541  NMCRings++;
542  }
543 
544  TH2F* h_e = new TH2F("MC ring radius (cm) vs MC track momentum (GeV)",
545  "",
546  100,
547  0.f,
548  12.f,
549  100,
550  0.f,
551  6.f);
552  h_e->SetXTitle("MC track momentum (GeV)");
553  h_e->SetYTitle("MC ring R1/R2");
554 
555  for (map<Int_t, MCRing>::iterator i = MCRingMap.begin(); i != MCRingMap.end();
556  ++i) {
557  vector<Double_t> hitX;
558  vector<Double_t> hitY;
559  // CbmRichRing *Ellipse = new CbmRichRing();
560 
561  MCRing& ring = i->second;
562  ring.MCTrackID = i->first;
563  ring.primary = 0;
564  ring.P = 0;
565  ring.PDG = 0;
566  ring.Reconstructed = 0;
567  ring.kind = 0;
568  ring.k = 0;
569 
570  // find momentum
571  if (!fMCTrackArray || ring.MCTrackID < 0) continue;
572 
573  CbmMCTrack* pm =
574  L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTrackArray->At(ring.MCTrackID));
575  if (pm) {
576  ring.PDG = pm->GetPdgCode(); // get PDG code of mother
577 
578  Double_t vx = pm->GetStartX();
579  Double_t vy = pm->GetStartY();
580  Double_t vz = pm->GetStartZ();
581 
582  if (fabs(vx) < 0.1 && fabs(vy) < 0.1 && fabs(vz) < 0.1) ring.primary = 1;
583 
584  ring.P = pm->GetP();
585  }
586 
587  // fit the MC ring
588 
589  list<pair<Double_t, Double_t>> L;
590  Int_t in = 0;
591  for (vector<int>::iterator j = ring.Hits.begin(); j != ring.Hits.end();
592  ++j) {
593  L.push_back(pair<Double_t, Double_t>(Hits[(*j)].x, Hits[(*j)].y));
594  hitX.push_back(Hits[(*j)].x);
595  hitY.push_back(Hits[(*j)].y);
596  in++;
597  }
598  CirFit(L, &ring.x, &ring.y, &ring.r);
599 #ifdef DRAW
600  // Use Fitter EllipseTau to fit MC Rings
601 /* CbmRichRingFitterEllipseTau *FitterEllipse = new CbmRichRingFitterEllipseTau();
602  FitterEllipse -> DoFit1(Ellipse, hitX , hitY );
603  Double_t x0, y0, d, delta ,A=Ellipse->GetAPar()
604  ,B=Ellipse->GetBPar()
605  ,C=Ellipse->GetCPar()
606  ,D=Ellipse->GetDPar()
607  ,E=Ellipse->GetEPar()
608  ,F=Ellipse->GetFPar()
609  ,l1,l2, r1, r2, theta;
610  d = A*C-0.25*B*B;
611  x0 = -(0.5*D*C-0.25*B*E)/d;
612  y0 = -(0.5*A*E-0.25*B*D)/d;
613  delta = A*(C*F-0.25*E*E) - 0.5*B*(0.5*B*F-0.25*E*D) + 0.5*D*(0.25*B*E-0.5*D*C);
614  l1 = 0.5*(A+C-sqrt((A+C)*(A+C)-4*A*C+B*B));
615  l2 = 0.5*(A+C+sqrt((A+C)*(A+C)-4*A*C+B*B));
616  r1 = sqrt(-delta/(l2*d));
617  r2 = sqrt(-delta/(l1*d));
618  theta = -(90*atan(B/(A-C)))/3.1415;
619  Up->cd(1);
620  TEllipse *el1 = new TEllipse( x0 , -y0 , r1 , r2, 0, 360, theta );
621  el1->SetFillStyle(0);
622  el1->SetLineColor(2);
623  el1->Draw();
624  Up->cd(2);
625  TEllipse *el2 = new TEllipse( x0 , -y0 , r1 , r2, 0, 360, theta );
626  el2->SetFillStyle(0);
627  el2->SetLineColor(2);
628  el2->Draw();
629  delete FitterEllipse;*/
630 //End of Use Fitter EllipseTau to fit MC Rings
631 #endif //DRAW
632  if (ring.r > 3. && ring.r < 7. && ring.NHits >= 7 /* && ring.P > 0.2 */) {
633  if (ring.NHits >= 15 && ring.primary) {
634  ring.kind = 2;
635  } else
636  ring.kind = 1;
637  } else
638  ring.kind = 0;
639 
640 #ifdef DRAW
641  Up->cd(1);
642  TArc* MCUp = new TArc(ring.x, -ring.y, ring.r, 0, 360);
643  MCUp->SetLineWidth(2);
644  if (ring.PDG == 11 || ring.PDG == -11) {
645  MCUp->SetLineColor(2); // electron MC Ring
646  if (ring.kind != 0) NeMC++;
647  } else if (ring.PDG == 211 || ring.PDG == -211 || ring.PDG == 111) {
648  MCUp->SetLineColor(28); // pion MC Ring
649  if (ring.kind != 0) NpiMC++;
650  } else
651  MCUp->SetLineColor(5); // other MC Ring
652  MCUp->SetFillStyle(0);
653 
654  if (ring.kind == 0) {
655  MCUp->SetLineStyle(3);
656  MCUp->SetLineWidth(1);
657  } else
658  MCUp->SetLineStyle(1);
659 
660  // Draw NHits and P for MC rings
661  /*
662  TString st1 = "";
663  st1 += ring.NHits;
664  st1 += " P = ";
665  st1 += ring.P;
666  TText* textMc1 = new TText( ring.x + ring.r/5 , -ring.y - ring.r/5 , st1.Data() );
667  textMc1->SetTextAlign(12);
668  textMc1->SetTextSize(0.01);
669  textMc1->SetTextColor(1);
670  textMc1->Draw();
671  */
672  MCUp->Draw("*");
673 
674  Up->cd(2);
675  // Draw NHits and P for MC rings
676  /*
677  TText* textMc2 = new TText( ring.x + ring.r/5 , -ring.y - ring.r/5 , st1.Data() );
678  textMc2->SetTextAlign(12);
679  textMc2->SetTextSize(0.01);
680  textMc2->SetTextColor(1);
681  textMc2->Draw();
682  */
683 
684  TArc* MCDown = new TArc(ring.x, -ring.y, ring.r, 0, 360);
685  MCDown->SetLineWidth(2);
686  if (ring.PDG == 11 || ring.PDG == -11)
687  MCDown->SetLineColor(2); //electron MC Ring
688  else if (ring.PDG == 211 || ring.PDG == -211)
689  MCDown->SetLineColor(28); //pion MC Ring
690  else
691  MCDown->SetLineColor(5); // other MC Ring
692  MCDown->SetFillStyle(0);
693  if (ring.kind == 0) {
694  MCDown->SetLineStyle(3);
695  MCDown->SetLineWidth(1);
696  } else
697  MCDown->SetLineStyle(1);
698  MCDown->Draw("*");
699 
700  for (Int_t i1 = 0; i1 < fRingArray->GetEntriesFast(); i1++) {
701  // Reco Rings
702  if (((sqrt((ArrRingX[i1] - ring.x) * (ArrRingX[i1] - ring.x)
703  + (ArrRingY[i1] + ring.y) * (ArrRingY[i1] + ring.y))
704  < 0.15 * ring.r)
705  && (sqrt((ArrRingR[i1] - ring.r) * (ArrRingR[i1] - ring.r))
706  < 0.15 * ring.r)
707  && (ring.k == 0))
708  || ((sqrt((ArrRingX[i1] - ring.x) * (ArrRingX[i1] - ring.x)
709  + (ArrRingY[i1] + ring.y) * (ArrRingY[i1] + ring.y))
710  < 0.2 * ring.r)
711  && (sqrt((ArrRingR[i1] - ring.r) * (ArrRingR[i1] - ring.r))
712  < 0.05 * ring.r)
713  && (ring.k == 0))) {
714  Up->cd(1);
715  TArc* RecoRingUp =
716  new TArc(ArrRingX[i1], ArrRingY[i1], ArrRingR[i1], 0, 360);
717  NGostmy--;
718  ring.k = 1;
719  if (ring.r > 3 && ring.r < 7 && ring.NHits >= 5) NRecoG++;
720  RecoRingUp->SetLineWidth(1);
721  RecoRingUp->SetLineColor(4);
722  if (ring.r < 3 || ring.r > 7 || ring.NHits < 5)
723  RecoRingUp->SetLineColor(5);
724  RecoRingUp->SetFillStyle(0);
725  RecoRingUp->Draw("*");
726  Up->cd(2);
727  TArc* RecoRingDown =
728  new TArc(ArrRingX[i1], ArrRingY[i1], ArrRingR[i1], 0, 360);
729  RecoRingDown->SetLineWidth(1);
730  RecoRingDown->SetLineColor(4);
731  if (ring.r < 3 || ring.r > 7 || ring.NHits < 5)
732  RecoRingDown->SetLineColor(5);
733  RecoRingDown->SetFillStyle(0);
734  RecoRingDown->Draw("*");
735  }
736  // Good Reco Rings
737  if (sqrt((ArrRingX[i1] - ring.x) * (ArrRingX[i1] - ring.x)
738  + (ArrRingY[i1] + ring.y) * (ArrRingY[i1] + ring.y))
739  < 0.05 * ring.r
740  && sqrt((ArrRingR[i1] - ring.r) * (ArrRingR[i1] - ring.r))
741  < 0.05 * ring.r
742  && ring.k != 2) {
743  Up->cd(1);
744  if (ring.k == 1) NRecoG--;
745  if (ring.r > 3 && ring.r < 7 && ring.NHits >= 5) NRecoEx++;
746  TArc* GoodRecoRingUp =
747  new TArc(ArrRingX[i1], ArrRingY[i1], ArrRingR[i1], 0, 360);
748  GoodRecoRingUp->SetLineWidth(1);
749  GoodRecoRingUp->SetLineColor(3);
750  if (ring.r < 3 || ring.r > 7 || ring.NHits < 5)
751  GoodRecoRingUp->SetLineColor(5);
752  GoodRecoRingUp->SetFillStyle(0);
753  ring.k = 2;
754  GoodRecoRingUp->Draw("*");
755  Up->cd(2);
756  TArc* GoodRecoRingDown =
757  new TArc(ArrRingX[i1], ArrRingY[i1], ArrRingR[i1], 0, 360);
758  GoodRecoRingDown->SetLineWidth(1);
759  GoodRecoRingDown->SetLineColor(3);
760  if (ring.r < 3 || ring.r > 7 || ring.NHits < 5)
761  GoodRecoRingDown->SetLineColor(5);
762  GoodRecoRingDown->SetFillStyle(0);
763  GoodRecoRingDown->Draw("*");
764  }
765  }
766 #endif //DRAW
767  }
768  // match reconstructed->simulated rings, calc. ghost
769 
770  Int_t NGhost = 0;
771  Int_t NReco = fRingArray->GetEntriesFast();
772  // Double_t local_x[fRingArray->GetEntriesFast()];
773  // Double_t local_y[fRingArray->GetEntriesFast()];
774  // Double_t local_r[fRingArray->GetEntriesFast()];
775 
776  Int_t MCRecoCor[NReco];
777 
778  for (Int_t ir = 0; ir < NReco; ir++) {
779  MCRecoCor[ir] = -1;
780  map<Int_t, Int_t> hitmap;
781  CbmRichRing* r = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(ir));
782  // local_x[ir] = r->GetCenterX();
783  // local_y[ir] = r->GetCenterY();
784  // local_r[ir] = r->GetRadius();
785  Int_t nh = r->GetNofHits();
786  for (int ih = 0; ih < nh; ih++) {
787  CbmRichHit* h =
788  L1_DYNAMIC_CAST<CbmRichHit*>(fHitArray->At(r->GetHit(ih)));
789  Int_t jh = pHitIndex[h];
790  int ID = Hits[jh].MCTrackID;
791  if (hitmap.find(ID) == hitmap.end())
792  hitmap.insert(pair<Int_t, Int_t>(ID, 0));
793  hitmap[ID]++;
794  }
795  Bool_t ghost = 1;
796  Int_t CurrentMCTrack;
797  for (map<Int_t, Int_t>::iterator j = hitmap.begin(); j != hitmap.end();
798  ++j) {
799  if ((static_cast<Double_t>(j->second)) < 0.7 * nh) continue;
800  CurrentMCTrack = j->first;
801  MCRecoCor[ir] = MCRingMap[j->first].MCTrackID;
802  MCRingMap[j->first].Reconstructed++;
803  if ((static_cast<Int_t>(j->second)) > MCRingMap[j->first].NHitsBestReco) {
804  MCRingMap[j->first].NHitsBestReco = (static_cast<Int_t>(j->second));
805  MCRingMap[j->first].NHitsBestvsNHitsMC =
806  (static_cast<Double_t>(j->second)) / MCRingMap[j->first].NHits;
807  }
808  ghost = 0;
809  break;
810  }
811  if (ghost) {
812  h_ghost_nhits->Fill(nh);
813  NGhost++;
814  Chi2Ghost->Fill(r->GetChi2());
815  Chi2NhitsGhost->Fill(r->GetNofHits(), r->GetChi2());
816  RGhost->Fill(r->GetRadius());
817  RNhitsGhost->Fill(r->GetNofHits(), r->GetRadius());
818  RChi2Ghost->Fill(r->GetChi2(), r->GetRadius());
819  } else {
820  Chi2All->Fill(r->GetChi2());
821  Chi2NhitsAll->Fill(r->GetNofHits(), r->GetChi2());
822  MCRing& ring_mc = MCRingMap[CurrentMCTrack];
823  NHitsRecoVsNHitsMC->Fill(ring_mc.NHits, r->GetNofHits());
824  if (ring_mc.PDG == 11 || ring_mc.PDG == -11) {
825  REl->Fill(ring_mc.r);
826  Chi2NhitsEll->Fill(r->GetNofHits(), r->GetChi2());
827  RNhitsEll->Fill(r->GetNofHits(), r->GetRadius());
828  RChi2Ell->Fill(r->GetChi2(), r->GetRadius());
829  } else if (ring_mc.PDG == 211 || ring_mc.PDG == -211
830  || ring_mc.PDG == 111) {
831  RPi->Fill(ring_mc.r);
832  Chi2NhitsPi->Fill(r->GetNofHits(), r->GetChi2());
833  RNhitsPi->Fill(r->GetNofHits(), r->GetRadius());
834  RChi2Pi->Fill(r->GetChi2(), r->GetRadius());
835  }
836  }
837  }
838  // Int_t IfClone[NReco];
839 
840  for (Int_t i1 = 0; i1 < NReco; i1++) {
841  // IfClone[i1] = 0;
842  CbmRichRing* r1 = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(i1));
843  Int_t nh1 = r1->GetNofHits();
844  for (Int_t i2 = 0; i2 < NReco; i2++) {
845  if (MCRecoCor[i1] != MCRecoCor[i2]) continue;
846  CbmRichRing* r2 = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(i2));
847  Int_t nh2 = r2->GetNofHits();
848  Int_t NSame = 0;
849  for (Int_t j1 = 0; j1 < nh1; j1++) {
850  for (Int_t j2 = 0; j2 < nh2; j2++) {
851  if (r1->GetHit(j1) == r2->GetHit(j2)) NSame++;
852  }
853  }
854  if (NSame != 0) {
855  NSameHits->Fill(NSame);
856  NHitsVsMCP->Fill(MCRingMap[MCRecoCor[i1]].P, NSame);
857  }
858  }
859  }
860  Int_t myNClone = 0;
861  Int_t CloneFlag[NReco];
862  for (Int_t i = 0; i < NReco; i++) {
863  CloneFlag[i] = 0;
864  }
865  for (Int_t i = 0; i < NReco; i++) {
866  CbmRichRing* r1 = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(i));
867  if (MCRecoCor[i] == -1) continue;
868  for (Int_t j = i + 1; j < NReco; j++) {
869  CbmRichRing* r2 = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(j));
870  if (MCRecoCor[j] == -1) continue;
871  if (MCRecoCor[i] == MCRecoCor[j]) {
872  CloneFlag[i] = MCRecoCor[i];
873  CloneFlag[j] = MCRecoCor[j];
874  myNClone++;
875  double dist = 0;
876  dist = sqrt((r1->GetCenterX() - r2->GetCenterX())
877  * (r1->GetCenterX() - r2->GetCenterX())
878  + (r1->GetCenterY() - r2->GetCenterY())
879  * (r1->GetCenterY() - r2->GetCenterY()));
880  double dr;
881  if (r1->GetRadius() >= r2->GetRadius())
882  dr = r1->GetRadius() - r2->GetRadius();
883  else
884  dr = r2->GetRadius() - r1->GetRadius();
885  RadiusVsDistanceClone->Fill(dist, dr);
886  for (map<Int_t, MCRing>::iterator MC = MCRingMap.begin();
887  MC != MCRingMap.end();
888  ++MC) {
889  MCRing& ring = MC->second;
890  if (ring.MCTrackID == -1) continue;
891  if (ring.MCTrackID == MCRecoCor[j]) {
892  RadiusVsPForClone->Fill(ring.P, dr);
893  DistanceVsPClone->Fill(ring.P, dist);
894  Chi2VsPClone->Fill(ring.P, r1->GetChi2());
895  continue;
896  }
897  }
898  MCRecoCor[j] = -1;
899  }
900  }
901  MCRecoCor[i] = -1;
902  }
903  //Draw Clones
904  for (Int_t i = 0; i < NReco; i++) {
905  if (CloneFlag[i] == 0) continue;
906  CbmRichRing* r = L1_DYNAMIC_CAST<CbmRichRing*>(fRingArray->At(i));
907  for (map<Int_t, MCRing>::iterator MC = MCRingMap.begin();
908  MC != MCRingMap.end();
909  ++MC) {
910  MCRing& ring = MC->second;
911  if (CloneFlag[i] != 0) {
912  if (CloneFlag[i] == ring.MCTrackID) {
913  Up->cd(1);
914  TArc* MCringUp = new TArc(ring.x, -ring.y, ring.r, 0, 360);
915  MCringUp->SetLineWidth(1);
916  MCringUp->SetLineColor(2);
917  MCringUp->SetFillStyle(0);
919  Up->cd(2);
920  TArc* MCringDown = new TArc(ring.x, -ring.y, ring.r, 0, 360);
921  MCringDown->SetLineWidth(1);
922  MCringDown->SetLineColor(2);
923  MCringDown->SetFillStyle(0);
925  }
926  Up->cd(1);
927  TArc* CloneUp =
928  new TArc(r->GetCenterX(), -r->GetCenterY(), r->GetRadius(), 0, 360);
929  CloneUp->SetLineWidth(1);
930  CloneUp->SetLineColor(7);
931  CloneUp->SetFillStyle(0);
932  CloneUp->Draw("*");
933  Up->cd(2);
934  TArc* CloneDown =
935  new TArc(r->GetCenterX(), -r->GetCenterY(), r->GetRadius(), 0, 360);
936  CloneDown->SetLineWidth(1);
937  CloneDown->SetLineColor(7);
938  CloneDown->SetFillStyle(0);
939  CloneDown->Draw("*");
940  }
941  }
942  }
943  // End of Draw Clones
944  // get statistics from MC rings
945  Int_t NAll = 0, NRef = 0, NExtra = 0, NAllReco = 0, NRefReco = 0,
946  NExtraReco = 0, NClone = 0;
947  for (map<Int_t, MCRing>::iterator i = MCRingMap.begin(); i != MCRingMap.end();
948  ++i) {
949  MCRing& ring = i->second;
950  if (ring.NHitsBestvsNHitsMC != 0.0) {
951  NHitsMC->Fill(ring.NHitsBestvsNHitsMC);
952  NSameHitsVsP->Fill(ring.P, ring.NHitsBestvsNHitsMC);
953  }
954  if (ring.kind == 0) {
955  NMC++;
956  continue;
957  }
958  NAll++;
959  NMC37++;
960  // Int_t fl=0;
961  if (ring.kind == 1) NExtra++;
962  if (ring.kind == 2) NRef++;
963  if (ring.Reconstructed) {
964  // fl = NClone;
965  NAllReco++;
966  NClone += i->second.Reconstructed - 1;
967  if (ring.kind == 2)
968  NRefReco++;
969  else
970  NExtraReco++;
971  }
972  }
973 
974 #ifdef DRAW
975  // Draw Hits
976  Up->cd(1);
977  TGraph* U1 = new TGraph(fHitArray->GetEntriesFast(), HitArrX, HitArrY);
978  U1->SetMarkerColor(kBlue);
979  U1->SetMarkerStyle(20);
980  U1->SetMarkerSize(0.3);
981  U1->Draw("P");
982  TGraph* U2 = new TGraph(NfakeHits, NoiseHitArrX, NoiseHitArrY);
983  U2->SetMarkerColor(7);
984  U2->SetMarkerStyle(20);
985  U2->SetMarkerSize(0.3);
986  U2->Draw("P,Same");
987  Up->cd(2);
988  TGraph* U = new TGraph(fHitArray->GetEntriesFast(), HitArrX, HitArrY);
989  U->SetMarkerColor(kBlue);
990  U->SetMarkerStyle(20);
991  U->SetMarkerSize(0.3);
992  U->Draw("P");
993  TGraph* U3 = new TGraph(NfakeHits, NoiseHitArrX, NoiseHitArrY);
994  U3->SetMarkerColor(7);
995  U3->SetMarkerStyle(20);
996  U3->SetMarkerSize(0.3);
997  U3->Draw("P,Same");
998 #endif //DRAW
999 
1000  // accumulated statistics
1001 
1002  static Int_t S_NAll = 0, S_NRef = 0, S_NExtra = 0, S_NReco = 0,
1003  S_NAllReco = 0, S_NRefReco = 0, S_NExtraReco = 0, S_NClone = 0,
1004  S_NGhost = 0, S_NEvents = 0;
1005 
1006  S_NAll += NAll;
1007  S_NRef += NRef;
1008  S_NExtra += NExtra;
1009  S_NReco += NReco;
1010  S_NAllReco += NAllReco;
1011  S_NRefReco += NRefReco;
1012  S_NExtraReco += NExtraReco;
1013  S_NClone += NClone;
1014  S_NGhost += NGhost;
1015 
1016 
1017  S_NEvents++;
1018 
1019  // write statistics
1020 
1021  cout.setf(ios::fixed);
1022  cout.setf(ios::showpoint);
1023  cout.precision(4);
1024 
1025 
1026  Double_t p_all = (NAll > 0) ? static_cast<Double_t>(NAllReco) / NAll : 0.;
1027  Double_t p_ref = (NRef > 0) ? static_cast<Double_t>(NRefReco) / NRef : 0.;
1028  Double_t p_extra =
1029  (NExtra > 0) ? static_cast<Double_t>(NExtraReco) / NExtra : 0.;
1030  Double_t p_clone = (NReco > 0) ? static_cast<Double_t>(NClone) / NReco : 0.;
1031  Double_t p_ghost = (NReco > 0) ? static_cast<Double_t>(NGhost) / NReco : 0.;
1032  // Double_t N_all = double(NRecoG+NRecoEx)/double(NMC37);
1033 
1034  cout << endl;
1035  cout << "L1ENN PER EVENT STAT : " << endl;
1036  cout << "MC Refset : " << NRef << endl;
1037  cout << "MC Extras : " << NExtra << endl;
1038  cout << "ALL SIMULATED : " << NAll << endl;
1039  cout << endl;
1040  cout << "RC Refset : " << NRefReco << endl;
1041  cout << "RC Extras : " << NExtraReco << endl;
1042  cout << "clones : " << NClone << endl;
1043  cout << "ghosts : " << NGhost << endl;
1044  cout << "ALL RECONSTRUCTED : " << NAllReco << endl;
1045  cout << endl;
1046  cout << "Refset efficiency : " << p_ref << endl;
1047  cout << "Allset efficiency : " << p_all << endl;
1048  cout << "Extra efficiency : " << p_extra << endl;
1049  cout << "clone probability : " << p_clone << endl;
1050  cout << "ghost probability : " << p_ghost << endl;
1051  cout << endl;
1052  // cout << "Reco time (ms) : " << fRecoTime/1000. << endl;
1053 
1054 
1055  std::string MC_Refset;
1056  std::stringstream out1;
1057  out1 << NRef;
1058  MC_Refset = out.str();
1059 
1060 
1061  TString name = "2_" + s + ".pdf";
1062 #ifdef DRAW
1063  Up->SaveAs(name);
1064 #endif // DRAW
1065 
1066  cout << endl;
1067 
1068  Double_t S_p_all =
1069  (S_NAll > 0) ? static_cast<Double_t>(S_NAllReco) / S_NAll : 0.;
1070  Double_t S_p_ref =
1071  (S_NRef > 0) ? static_cast<Double_t>(S_NRefReco) / S_NRef : 0.;
1072  Double_t S_p_extra =
1073  (S_NExtra > 0) ? static_cast<Double_t>(S_NExtraReco) / S_NExtra : 0.;
1074  Double_t S_p_clone =
1075  (S_NReco > 0) ? static_cast<Double_t>(S_NClone) / S_NReco : 0.;
1076  Double_t S_p_ghost =
1077  (S_NReco > 0) ? static_cast<Double_t>(S_NGhost) / S_NReco : 0.;
1078 
1079  cout << "ACCUMULATED STAT : " << S_NEvents << " EVENTS " << endl << endl;
1080 
1081  cout << " "
1082  << " % "
1083  << " | "
1084  << "Reco"
1085  << " | "
1086  << "MC" << endl;
1087  cout << "Refset efficiency : " << S_p_ref << " | " << S_NRefReco << " | "
1088  << S_NRef << endl;
1089  cout << "Allset efficiency : " << S_p_all << " | " << S_NAllReco << " | "
1090  << S_NAll << endl;
1091  cout << "Extra efficiency : " << S_p_extra << " | " << S_NExtraReco
1092  << " | " << S_NExtra << endl;
1093  cout << "clone probability : " << S_p_clone << " | " << S_NClone << endl;
1094  cout << "ghost probability : " << S_p_ghost << " | " << S_NGhost << endl;
1095  cout << "MC rings/event found : "
1096  << Int_t(static_cast<Double_t>(S_NAllReco)
1097  / static_cast<Double_t>(S_NEvents))
1098  << endl;
1099  cout << endl;
1100  cout << "Reco time (ms) : " << 0. * 1000. / S_NEvents << endl;
1101 
1102  cout << endl;
1103 
1104 
1105  // fill histogramms
1106 
1107  for (map<Int_t, MCRing>::iterator i = MCRingMap.begin(); i != MCRingMap.end();
1108  ++i) {
1109  MCRing& ring = i->second;
1110 
1111  h_MC_radius->Fill(ring.r);
1112  h_MC_nhits->Fill(ring.NHits);
1113  h_MC_momentum->Fill(ring.P);
1114  if (ring.primary) {
1115  h_MC_primary_nhits->Fill(ring.NHits);
1116  h_MC_primary_momentum->Fill(ring.P);
1117  }
1118 
1119  for (vector<int>::iterator j = ring.Hits.begin(); j != ring.Hits.end();
1120  ++j) {
1121  Double_t dx = Hits[(*j)].x - ring.x;
1122  Double_t dy = Hits[(*j)].y - ring.y;
1123  Double_t res = sqrt(dx * dx + dy * dy) - ring.r;
1124  h_MC_resolution->Fill(res);
1125  if (ring.kind == 1) h_MC_extra_resolution->Fill(res);
1126  if (ring.kind == 2) h_MC_ref_resolution->Fill(res);
1127  if (ring.primary) h_MC_primary_res_vs_momentum->Fill(ring.P, res);
1128  }
1129  Double_t entry = (ring.Reconstructed) ? 1.0 : 0.0;
1130  if (ring.kind == 1) p_extra_eff_vs_nhits->Fill(ring.NHits, entry);
1131  if (ring.kind == 2) p_ref_eff_vs_nhits->Fill(ring.NHits, entry);
1132  }
1133 
1134  // Open output file and write histograms
1135  TDirectory* curr = gDirectory;
1136  TFile* outfile = new TFile("L1RingQaHisto.root", "RECREATE");
1137  outfile->cd();
1138  TIter hiter(listHisto);
1139  while (TObject* obj = hiter())
1140  obj->Write();
1141  outfile->Close();
1142  curr->cd();
1143 }
1144 
1145 void CbmL1RichRingQa::Finish() {
1146  TStyle plain("Plain", "Plain Style(no colors/fill areas)");
1147  plain.SetPadColor(0);
1148  plain.SetCanvasColor(0);
1149  plain.SetTitleColor(0);
1150  plain.SetStatColor(0);
1151  plain.SetOptFit(1111);
1152  plain.SetStatW(0.225);
1153  plain.SetOptStat(0);
1154  plain.SetPalette(1);
1155  plain.cd();
1156 
1157  TCanvas* Chi2Histos = new TCanvas("Chi2Histos", "Chi2Histos", 0, 0, 240, 270);
1158  Chi2Histos->SetFillColor(kWhite);
1159  Chi2Histos->Divide(2, 2);
1160  Chi2Histos->cd(1);
1161  Chi2Ghost->Draw();
1162  Chi2Histos->cd(2);
1163  //Chi2Ref->Draw();
1164  Chi2NhitsGhost->Draw("colz");
1165  Chi2Histos->cd(3);
1166  Chi2All->Draw();
1167  Chi2Histos->cd(4);
1168  //Chi2Clone->Draw();
1169  Chi2NhitsAll->Draw("colz");
1170  Chi2Histos->SaveAs("Chi2.pdf");
1171  TCanvas* All = new TCanvas("All", "All", 0, 0, 240, 270);
1172  All->SetFillColor(kWhite);
1173  All->Divide(3, 3);
1174  All->cd(1);
1175  Chi2NhitsGhost->Draw("colz");
1176  All->cd(2);
1177  Chi2NhitsPi->Draw("colz");
1178  All->cd(3);
1179  Chi2NhitsEll->Draw("colz");
1180  All->cd(4);
1181  RNhitsGhost->Draw("colz");
1182  All->cd(5);
1183  RNhitsPi->Draw("colz");
1184  All->cd(6);
1185  RNhitsEll->Draw("colz");
1186  All->cd(7);
1187  RChi2Ghost->Draw("colz");
1188  All->cd(8);
1189  RChi2Pi->Draw("colz");
1190  All->cd(9);
1191  RChi2Ell->Draw("colz");
1192  All->SaveAs("All.pdf");
1193 
1194  TCanvas* MCNHits = new TCanvas("MCNHits", "MCNHits", 0, 0, 240, 270);
1195  MCNHits->SetFillColor(kWhite);
1196  MCNHits->Divide(2, 2);
1197  MCNHits->cd(1);
1198  NHitsMC->Draw();
1199  MCNHits->cd(2);
1200  NSameHits->Draw();
1201  MCNHits->cd(3);
1202  NSameHitsVsP->Draw("colz");
1203  MCNHits->cd(4);
1204  NHitsVsMCP->Draw("colz");
1205 
1206  MCNHits->SaveAs("MCHits.pdf");
1207 
1208  TCanvas* CloneToEll = new TCanvas("CloneToEll", "CloneToEll", 0, 0, 240, 270);
1209  CloneToEll->SetFillColor(kWhite);
1210  CloneToEll->Divide(2, 2);
1211  CloneToEll->cd(1);
1212  RadiusVsPForClone->Draw("colz");
1213  CloneToEll->cd(2);
1214  DistanceVsPClone->Draw("colz");
1215  CloneToEll->cd(3);
1216  RadiusVsDistanceClone->Draw("colz");
1217  CloneToEll->cd(4);
1218  Chi2VsPClone->Draw("colz");
1219 
1220  CloneToEll->SaveAs("CloneToEll.pdf");
1221 
1222  TCanvas* RecoVsMC = new TCanvas("RecoVsMC", "RecoVsMC", 0, 0, 240, 270);
1223  RecoVsMC->SetFillColor(kWhite);
1224  NHitsRecoVsNHitsMC->Draw("colz");
1225 
1226  RecoVsMC->SaveAs("RecoVsMC.pdf");
1227 
1228  delete Chi2Histos;
1229  delete All;
1230  delete MCNHits;
1231  delete CloneToEll;
1232  delete RecoVsMC;
1233  if (Chi2Ghost) delete Chi2Ghost;
1234  if (Chi2Ref) delete Chi2Ref;
1235  if (Chi2All) delete Chi2All;
1236  if (Chi2Clone) delete Chi2Clone;
1237  if (Chi2NhitsGhost) delete Chi2NhitsGhost;
1238  if (Chi2NhitsAll) delete Chi2NhitsAll;
1239  if (REl) delete REl;
1240  if (RPi) delete RPi;
1241  if (RGhost) delete RGhost;
1242  if (Chi2NhitsEll) delete Chi2NhitsEll;
1243  if (Chi2NhitsPi) delete Chi2NhitsPi;
1244  if (RNhitsGhost) delete RNhitsGhost;
1245  if (RNhitsEll) delete RNhitsEll;
1246  if (RNhitsPi) delete RNhitsPi;
1247  if (RChi2Ghost) delete RChi2Ghost;
1248  if (RChi2Ell) delete RChi2Ell;
1249  if (RChi2Pi) delete RChi2Pi;
1250  if (NHitsMC) delete NHitsMC;
1251  if (NSameHits) delete NSameHits;
1252  if (NHitsVsMCP) delete NHitsVsMCP;
1253  if (NSameHitsVsP) delete NSameHitsVsP;
1254  if (RadiusVsPForClone) delete RadiusVsPForClone;
1255  if (DistanceVsPClone) delete DistanceVsPClone;
1256  if (Chi2VsPClone) delete Chi2VsPClone;
1257  if (RadiusVsDistanceClone) delete RadiusVsDistanceClone;
1258  if (NHitsRecoVsNHitsMC) delete NHitsRecoVsNHitsMC;
1259 }
CbmL1RichRingQa::RNhitsPi
TH2F * RNhitsPi
Definition: CbmL1RichRingQa.h:91
CbmRichPoint.h
CbmL1RichRingQa::fMCPointArray
TClonesArray * fMCPointArray
Definition: CbmL1RichRingQa.h:70
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition: CbmMCTrack.h:75
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
CbmL1RichRingQa::MCRing::kind
Int_t kind
Definition: CbmL1RichRingQa.h:52
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
CbmL1RichRingQa::MCRing::Reconstructed
Int_t Reconstructed
Definition: CbmL1RichRingQa.h:51
CbmL1RichRingQa::fRingArray
TClonesArray * fRingArray
Definition: CbmL1RichRingQa.h:69
CbmL1RichRingQa::RPi
TH1F * RPi
Definition: CbmL1RichRingQa.h:83
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmL1RichRingQa::REl
TH1F * REl
Definition: CbmL1RichRingQa.h:82
CbmL1RichRingQa::Init
InitStatus Init()
Definition: CbmL1RichRingQa.cxx:233
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmRichRingFitterEllipseTau.h
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmL1RichRingQa::MCRing::PDG
Int_t PDG
Definition: CbmL1RichRingQa.h:49
CbmL1RichRingQa::PerfHit::MCTrackID
Int_t MCTrackID
Definition: CbmL1RichRingQa.h:63
CbmL1RichRingQa.h
CbmL1RichRingQa::~CbmL1RichRingQa
~CbmL1RichRingQa()
Definition: CbmL1RichRingQa.cxx:231
CbmL1RichRingQa::MCRing::x
Double_t x
Definition: CbmL1RichRingQa.h:53
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
CbmL1RichRingQa::MCRing::NHitsBestvsNHitsMC
Double_t NHitsBestvsNHitsMC
Definition: CbmL1RichRingQa.h:58
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRing::GetChi2
Double_t GetChi2() const
Definition: CbmRichRing.h:95
CbmL1RichRingQa::MCRing::k
Int_t k
Definition: CbmL1RichRingQa.h:54
CbmL1RichRingQa::CirFit
void CirFit(std::list< std::pair< Double_t, Double_t >> &P, Double_t *X, Double_t *Y, Double_t *R)
Definition: CbmL1RichRingQa.cxx:273
CbmL1RichRingQa::Chi2Ref
TH1F * Chi2Ref
Definition: CbmL1RichRingQa.h:76
CbmL1RichRingQa::PerfHit::index
Int_t index
Definition: CbmL1RichRingQa.h:64
CbmL1RichRingQa::RGhost
TH1F * RGhost
Definition: CbmL1RichRingQa.h:81
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
CbmL1RichRingQa::NHitsVsMCP
TH2F * NHitsVsMCP
Definition: CbmL1RichRingQa.h:97
CbmRichRing.h
CbmL1RichRingQa::Finish
void Finish()
Definition: CbmL1RichRingQa.cxx:1145
CbmL1RichRingQa::RChi2Ell
TH2F * RChi2Ell
Definition: CbmL1RichRingQa.h:95
CbmL1RichRingQa::MCRing::NHits
Int_t NHits
Definition: CbmL1RichRingQa.h:50
CbmL1RichRingQa::DistanceVsPClone
TH2F * DistanceVsPClone
Definition: CbmL1RichRingQa.h:99
CbmL1RichRingQa::MCRing::Hits
std::vector< int > Hits
Definition: CbmL1RichRingQa.h:55
CbmL1RichRingQa::MCRing
Definition: CbmL1RichRingQa.h:29
CbmL1Def.h
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition: CbmMCTrack.h:77
CbmL1RichRingQa
Definition: CbmL1RichRingQa.h:25
CbmL1RichRingQa::RNhitsGhost
TH2F * RNhitsGhost
Definition: CbmL1RichRingQa.h:90
CbmRichRing::GetHit
UInt_t GetHit(Int_t i) const
Definition: CbmRichRing.h:42
CbmL1RichRingQa::NSameHitsVsP
TH2F * NSameHitsVsP
Definition: CbmL1RichRingQa.h:96
CbmL1RichRingQa::MCRing::r
Double_t r
Definition: CbmL1RichRingQa.h:53
CbmL1RichRingQa::MCRing::y
Double_t y
Definition: CbmL1RichRingQa.h:53
CbmL1RichRingQa::RadiusVsDistanceClone
TH2F * RadiusVsDistanceClone
Definition: CbmL1RichRingQa.h:101
CbmL1RichRingQa::Chi2Ghost
TH1F * Chi2Ghost
Definition: CbmL1RichRingQa.h:75
CbmL1RichRingQa::MCRing::primary
bool primary
Definition: CbmL1RichRingQa.h:47
CbmRichRingLight.h
CbmL1RichRingQa::NHitsRecoVsNHitsMC
TH2F * NHitsRecoVsNHitsMC
Definition: CbmL1RichRingQa.h:102
CbmL1RichRingQa::Chi2NhitsEll
TH2F * Chi2NhitsEll
Definition: CbmL1RichRingQa.h:89
CbmL1RichRingQa::RChi2Ghost
TH2F * RChi2Ghost
Definition: CbmL1RichRingQa.h:93
CbmL1RichRingQa::PerfHit::y
Double_t y
Definition: CbmL1RichRingQa.h:62
CbmL1RichRingQa::NHitsMC
TH1F * NHitsMC
Definition: CbmL1RichRingQa.h:84
CbmL1RichRingQa::RadiusVsPForClone
TH2F * RadiusVsPForClone
Definition: CbmL1RichRingQa.h:98
CbmL1RichRingQa::Chi2NhitsPi
TH2F * Chi2NhitsPi
Definition: CbmL1RichRingQa.h:88
CbmL1RichRingQa::Chi2VsPClone
TH2F * Chi2VsPClone
Definition: CbmL1RichRingQa.h:100
CbmL1RichRingQa::PerfHit::x
Double_t x
Definition: CbmL1RichRingQa.h:62
CbmL1RichRingQa::fMCTrackArray
TClonesArray * fMCTrackArray
Definition: CbmL1RichRingQa.h:71
CbmL1RichRingQa::MCRing::MCTrackID
Int_t MCTrackID
Definition: CbmL1RichRingQa.h:45
CbmL1RichRingQa::Chi2All
TH1F * Chi2All
Definition: CbmL1RichRingQa.h:77
CbmL1RichRingQa::RNhitsEll
TH2F * RNhitsEll
Definition: CbmL1RichRingQa.h:92
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition: CbmMCTrack.h:76
CbmMCTrack.h
CbmL1RichRingQa::Exec
void Exec(Option_t *option)
Definition: CbmL1RichRingQa.cxx:308
CbmMCTrack
Definition: CbmMCTrack.h:34
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
CbmL1RichRingQa::fHitArray
TClonesArray * fHitArray
Definition: CbmL1RichRingQa.h:72
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmL1RichRingQa::MCRing::P
Double_t P
Definition: CbmL1RichRingQa.h:48
CbmRichRing::GetRadius
Float_t GetRadius() const
Definition: CbmRichRing.h:82
CbmRichRing::GetCenterY
Float_t GetCenterY() const
Definition: CbmRichRing.h:81
ID
int ID
Definition: CbmMvdSensorDigiToHitTask.cxx:71
CbmL1RichRingQa::Chi2NhitsGhost
TH2F * Chi2NhitsGhost
Definition: CbmL1RichRingQa.h:79
CbmRichHit.h
CbmL1RichRingQa::Chi2NhitsAll
TH2F * Chi2NhitsAll
Definition: CbmL1RichRingQa.h:80
CbmL1RichRingQa::PerfHit
Definition: CbmL1RichRingQa.h:61
ClassImp
ClassImp(CbmL1RichRingQa) CbmL1RichRingQa
Definition: CbmL1RichRingQa.cxx:50
CbmRichPoint
Definition: CbmRichPoint.h:24
CbmL1RichRingQa::RChi2Pi
TH2F * RChi2Pi
Definition: CbmL1RichRingQa.h:94
CbmMCTrack::GetP
Double_t GetP() const
Definition: CbmMCTrack.h:100
CbmL1RichRingQa::NSameHits
TH1F * NSameHits
Definition: CbmL1RichRingQa.h:85
CbmRichRing::GetCenterX
Float_t GetCenterX() const
Definition: CbmRichRing.h:80
CbmRichHit
Definition: CbmRichHit.h:19
CbmL1RichRingQa::Chi2Clone
TH1F * Chi2Clone
Definition: CbmL1RichRingQa.h:78
CbmL1RichRingQa::CbmL1RichRingQa
CbmL1RichRingQa(const CbmL1RichRingQa &)