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