CbmRoot
CbmRichMCbmQaRichOnly.cxx
Go to the documentation of this file.
2 
3 #include "TCanvas.h"
4 #include "TClonesArray.h"
5 #include "TEllipse.h"
6 #include "TF1.h"
7 #include "TFile.h"
8 #include "TGeoBBox.h"
9 #include "TGeoManager.h"
10 #include "TGeoNode.h"
11 #include "TH1.h"
12 #include "TH1D.h"
13 #include "TLine.h"
14 #include "TMarker.h"
15 #include "TMath.h"
16 #include "TStyle.h"
17 #include "TSystem.h"
18 
19 #include <TBox.h>
20 #include <TLegend.h>
21 
22 #include "CbmDigiManager.h"
23 #include "CbmDrawHist.h"
24 #include "CbmEvent.h"
25 #include "CbmGlobalTrack.h"
26 #include "CbmMatchRecoToMC.h"
27 #include "CbmRichDigi.h"
28 #include "CbmRichDraw.h"
29 #include "CbmRichGeoManager.h"
30 #include "CbmRichHit.h"
31 #include "CbmRichPoint.h"
32 #include "CbmRichRing.h"
34 #include "CbmRichUtil.h"
35 #include "CbmStsDigi.h"
36 #include "CbmTofDigi.h"
37 #include "CbmTofHit.h"
38 #include "CbmTofTracklet.h"
39 #include "CbmTrackMatchNew.h"
40 #include "CbmTrdTrack.h"
41 #include "TLatex.h"
42 
43 #include "CbmRichMCbmSEDisplay.h"
44 
45 #include "CbmEvent.h"
46 
47 #include "CbmRichConverter.h"
48 
49 #include "CbmHistManager.h"
50 #include "CbmUtils.h"
51 
52 #include <boost/assign/list_of.hpp>
53 #include <cmath>
54 #include <fstream>
55 #include <iomanip>
56 #include <iostream>
57 #include <sstream>
58 #include <string>
59 
60 
61 using namespace std;
62 using boost::assign::list_of;
63 
64 #define RichZPos 348.
65 
67  : FairTask("CbmRichMCbmQaRichOnly")
68  , fRichHits(nullptr)
69  , fRichRings(nullptr)
70  , fCbmEvent(nullptr)
71  , fHM(nullptr)
72  , fEventNum(0)
73  , fNofDrawnRings(0)
74  , fNofDrawnRichTofEv(0)
75  , fNofDrawnEvents(0)
76  , fMaxNofDrawnEvents(100)
77  , fTriggerRichHits(0)
78  , fOutputDir("result") {}
79 
81  cout << "CbmRichMCbmQaRichOnly::Init" << endl;
82 
83  FairRootManager* ioman = FairRootManager::Instance();
84  if (nullptr == ioman) {
85  Fatal("CbmRichMCbmQaRichOnly::Init", "RootManager not instantised!");
86  }
87 
89  fDigiMan->Init();
90 
92  Fatal("CbmRichMCbmQaReal::Init", "No Rich Digis!");
93 
94  fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
95  if (nullptr == fRichHits) {
96  Fatal("CbmRichMCbmQaRichOnly::Init", "No Rich Hits!");
97  }
98 
99  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
100  if (nullptr == fRichRings) {
101  Fatal("CbmRichMCbmQaRichOnly::Init", "No Rich Rings!");
102  }
103 
104  fCbmEvent = (TClonesArray*) ioman->GetObject("CbmEvent");
105  if (nullptr == fCbmEvent) {
106  Fatal("CbmRichMCbmQaRichOnly::Init", "No Event!");
107  }
108 
109  InitHistograms();
110 
111 
116  fSeDisplay->SetTotRich(23.7, 30.);
119 
120  //Init OffsetCorrection ICD
121  for (auto& a : ICD_offset_read)
122  a = 0.;
123  for (auto& a : ICD_offset)
124  a = 0.;
125  for (auto& a : ICD_offset_cnt)
126  a = 0;
128 
129  return kSUCCESS;
130 }
131 
133  fHM = new CbmHistManager();
134 
135  fHM->Create1<TH1D>("fhNofEvents", "fhNofEvents;Entries", 1, 0.5, 1.5);
136  fHM->Create1<TH1D>("fhNofCbmEvents", "fhNofCbmEvents;Entries", 1, 0.5, 1.5);
137  fHM->Create1<TH1D>(
138  "fhNofCbmEventsRing", "fhNofCbmEventsRing;Entries", 1, 0.5, 1.5);
139 
140  fHM->Create1<TH1D>("fhNofBlobEvents", "fhNofBlobEvents;Entries", 1, 0.5, 1.5);
141  fHM->Create1<TH1D>(
142  "fhNofBlobsInEvent", "fhNofBlobsInEvent;Entries", 36, 0.5, 36.5);
143 
144  // RICH hits
145  fHM->Create2<TH2D>("fhRichHitXY",
146  "fhRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries",
147  67,
148  -20.1 + fXOffsetHisto,
149  20.1 + fXOffsetHisto,
150  84,
151  -25.2,
152  25.2);
153  fHM->Create2<TH2D>(
154  "fhRichHitXY_fromRing",
155  "fhRichHitXY_fromRing;RICH hit X [cm];RICH hit Y [cm];Entries",
156  67,
157  -20.1 + fXOffsetHisto,
158  20.1 + fXOffsetHisto,
159  84,
160  -25.2,
161  25.2);
162 
163  //ToT
164  fHM->Create1<TH1D>(
165  "fhRichDigisToT", "fhRichDigisToT;ToT [ns];Entries", 601, 9.975, 40.025);
166  fHM->Create1<TH1D>(
167  "fhRichHitToT", "fhRichHitToT;ToT [ns];Entries", 601, 9.975, 40.025);
168 
169  // RICH rings
170  fHM->Create2<TH2D>(
171  "fhRichRingXY",
172  "fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries",
173  67,
174  -20.1 + fXOffsetHisto,
175  20.1 + fXOffsetHisto,
176  84,
177  -25.2,
178  25.2);
179  fHM->Create1<TH1D>("fhRichRingRadius",
180  "fhRichRingRadius;Ring radius [cm];Entries",
181  100,
182  0.,
183  7.);
184  fHM->Create1<TH1D>("fhNofHitsInRing",
185  "fhNofHitsInRing;# hits in ring;Entries",
186  50,
187  -0.5,
188  49.5);
189  fHM->Create2<TH2D>(
190  "fhICD", "fhICD;channel;DeltaTime", 2305, -0.5, 2304.5, 31, -15.5, 15.5);
191 
192  fHM->Create2<TH2D>(
193  "fhRichRingRadiusY",
194  "fhRichRingRadiusY;Ring Radius [cm]; Y position[cm];Entries",
195  70,
196  -0.05,
197  6.95,
198  84,
199  -25.2,
200  25.2);
201  fHM->Create2<TH2D>(
202  "fhRichHitsRingRadius",
203  "fhRichHitsRingRadius;#Rich Hits/Ring; Ring Radius [cm];Entries",
204  50,
205  -0.5,
206  49.5,
207  70,
208  -0.05,
209  6.95);
210 
211  fHM->Create1<TH1D>("fhRingDeltaTime",
212  "fhRingDeltaTime; \\Delta Time/ns;Entries",
213  101,
214  -10.1,
215  10.1);
216  fHM->Create1<TH1D>(
217  "fhRingChi2", "fhRingChi2; \\Chi^2 ;Entries", 101, 0.0, 10.1);
218 
219  fHM->Create2<TH2D>(
220  "fhRichRingCenterXChi2",
221  "fhRichRingCenterXChi2;Ring Center X [cm];\\Chi^2 ;;Entries",
222  67,
223  -20.1 + fXOffsetHisto,
224  20.1 + fXOffsetHisto,
225  101,
226  0.0,
227  10.1);
228  fHM->Create2<TH2D>(
229  "fhRichRingCenterYChi2",
230  "fhRichRingCenterYChi2;Ring Center Y [cm];\\Chi^2 ;;Entries",
231  84,
232  -25.2,
233  25.2,
234  101,
235  0.0,
236  10.1);
237  fHM->Create2<TH2D>(
238  "fhRichRingRadiusChi2",
239  "fhRichRingRadiusChi2; Ring Radius [cm];\\Chi^2 ;;Entries",
240  70,
241  -0.05,
242  6.95,
243  101,
244  0.0,
245  10.1);
246 
247 
248  // Digis
249  fHM->Create2<TH2D>("fhDigisInChnl",
250  "fhDigisInChnl;channel;#Digis;",
251  2304,
252  -0.5,
253  2303.5,
254  50,
255  -0.5,
256  49.5);
257  fHM->Create2<TH2D>("fhDigisInDiRICH",
258  "fhDigisInDiRICH;DiRICH;#Digis;",
259  72,
260  -0.5,
261  71.5,
262  300,
263  -0.5,
264  299.5);
265 }
266 
267 
268 void CbmRichMCbmQaRichOnly::Exec(Option_t* /*option*/) {
269  fEventNum++;
270  fHM->H1("fhNofEvents")->Fill(1);
271 
272  cout << "CbmRichMCbmQaRichOnly, event No. " << fEventNum << endl;
273 
274  std::array<unsigned int, 2304> chnlDigis;
275  for (auto& c : chnlDigis)
276  c = 0;
277  for (int i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
278  const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
279  fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
280 
281  uint16_t addrDiRICH = (richDigi->GetAddress() >> 16) & 0xFFFF;
282  uint16_t addrChnl = richDigi->GetAddress() & 0xFFFF;
283  uint16_t dirichNmbr = ((addrDiRICH >> 8) & 0xF) * 18
284  + ((addrDiRICH >> 4) & 0xF) * 2
285  + ((addrDiRICH >> 0) & 0xF);
286  uint32_t fullNmbr = (dirichNmbr << 5) | (addrChnl - 1);
287  chnlDigis[fullNmbr]++;
288  }
289 
290  {
291  unsigned int sum = 0;
292  for (uint16_t i = 0; i < 2304; ++i) {
293  if (chnlDigis[i] != 0) fHM->H1("fhDigisInChnl")->Fill(i, chnlDigis[i]);
294  sum += chnlDigis[i];
295  if (i % 32 == 31) {
296  uint16_t dirich = i / 32;
297  if (sum != 0) fHM->H1("fhDigisInDiRICH")->Fill(dirich, sum);
298  sum = 0;
299  }
300  }
301  }
302 
303  int nofRichHits = fRichHits->GetEntries();
304  for (int iH = 0; iH < nofRichHits; iH++) {
305  CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
306  fHM->H2("fhRichHitXY")->Fill(richHit->GetX(), richHit->GetY());
307  fHM->H1("fhRichHitToT")->Fill(richHit->GetToT());
308  }
309 
310 
311  //CBMEVENT
312  auto fNCbmEvent = fCbmEvent->GetEntriesFast();
313 
314  for (int i = 0; i < fNCbmEvent; i++) {
315  fHM->H1("fhNofCbmEvents")->Fill(1);
316  CbmEvent* ev = static_cast<CbmEvent*>(fCbmEvent->At(i));
317 
318  if (fTriggerRichHits != 0
320  continue;
321 
322 
323  std::vector<int> ringIndx;
324  std::vector<int> evRichHitIndx;
325  std::array<uint32_t, 36> pmtHits;
326  for (auto& a : pmtHits)
327  a = 0;
328 
329  // Map Rings to CbmEvent
330  for (int j = 0; j < ev->GetNofData(ECbmDataType::kRichHit); j++) {
331  auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
332  evRichHitIndx.push_back(iRichHit);
333  CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
334 
335  uint32_t pmtId = (((richHit->GetAddress()) >> 20) & 0xF)
336  + (((richHit->GetAddress()) >> 24) & 0xF) * 9;
337  pmtHits[pmtId]++;
338 
339  int nofRichRings = fRichRings->GetEntries();
340  for (int l = 0; l < nofRichRings; l++) {
341  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(l));
342 
343  auto NofRingHits = ring->GetNofHits();
344  for (int m = 0; m < NofRingHits; m++) {
345  auto RingHitIndx = ring->GetHit(m);
346  if (RingHitIndx == iRichHit) {
347  Bool_t used = false;
348  for (auto check : ringIndx) {
349  if (check == l) {
350  used = true;
351  break;
352  }
353  }
354  if (used == false) ringIndx.push_back(l);
355  break;
356  }
357  }
358  }
359  }
360 
361  uint16_t blob = 0;
362  for (auto a : pmtHits) {
363  if (a > 30) { blob++; }
364  }
365  if (blob > 0) {
366  fHM->H1("fhNofBlobEvents")->Fill(1);
367  fHM->H1("fhNofBlobsInEvent")->Fill(blob);
368  }
369 
370  if (ringIndx.size() != 0) fHM->H1("fhNofCbmEventsRing")->Fill(1);
371 
372  //Ring Loop
373  for (unsigned int k = 0; k < ringIndx.size(); ++k) {
374  //Rigs in this CbmEvent
375  CbmRichRing* ring =
376  static_cast<CbmRichRing*>(fRichRings->At(ringIndx[k]));
377  if (ring == nullptr) continue;
378  fHM->H1("fhRingChi2")->Fill(ring->GetChi2());
379  fHM->H2("fhRichRingCenterXChi2")
380  ->Fill(ring->GetCenterX(), ring->GetChi2());
381  fHM->H2("fhRichRingCenterYChi2")
382  ->Fill(ring->GetCenterY(), ring->GetChi2());
383  fHM->H2("fhRichRingRadiusChi2")->Fill(ring->GetRadius(), ring->GetChi2());
384 
385  fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
386  fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
387  fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
388  fHM->H2("fhRichRingRadiusY")->Fill(ring->GetRadius(), ring->GetCenterY());
389  fHM->H2("fhRichHitsRingRadius")
390  ->Fill(ring->GetNofHits(), ring->GetRadius());
391  for (int j = 0; j < ring->GetNofHits(); ++j) {
392  Int_t hitIndx = ring->GetHit(j);
393  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitIndx);
394  if (nullptr == hit) continue;
395  fHM->H2("fhRichHitXY_fromRing")->Fill(hit->GetX(), hit->GetY());
396  }
397  }
399  static_cast<CbmEvent*>(fCbmEvent->At(i)), ringIndx, 1);
400 
401  } //End CbmEvent loop
402 
403  // Loop over all Rings
404  RichRings();
405 }
406 
408  int nofRichRings = fRichRings->GetEntries();
409  for (int i = 0; i < nofRichRings; i++) {
410  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(i));
411  if (ring == nullptr) continue;
412 
413  for (int j = 1; j < ring->GetNofHits(); ++j) {
414  Int_t hitIndx = ring->GetHit(j);
415  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitIndx);
416  if (nullptr == hit) continue;
417 
418  //Read Address
419  uint32_t DiRICH_Addr = hit->GetAddress();
420  unsigned int addr = (((DiRICH_Addr >> 24) & 0xF) * 18 * 32)
421  + (((DiRICH_Addr >> 20) & 0xF) * 2 * 32)
422  + (((DiRICH_Addr >> 16) & 0xF) * 32)
423  + ((DiRICH_Addr & 0xFFFF) - 0x1);
424  ICD_offset.at(addr) +=
425  hit->GetTime() - ring->GetTime() - ICD_offset_read.at(addr);
426  fHM->H2("fhICD")->Fill(
427  addr, hit->GetTime() - ring->GetTime() - ICD_offset_read.at(addr));
428  fHM->H1("fhRingDeltaTime")
429  ->Fill(hit->GetTime() - ring->GetTime() - ICD_offset_read.at(addr));
430  ICD_offset_cnt.at(addr)++;
431  }
432 
433  //DrawRing(ring);
434  /*fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
435  fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
436  fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
437  */
438  }
439 }
440 
441 
443  cout.precision(4);
444 
445  //SetDefaultDrawStyle();
446  double nofEvents = fHM->H1("fhNofCbmEvents")->GetEntries();
447  fHM->ScaleByPattern("fh_.*", 1. / nofEvents);
448 
449  {
450  fHM->CreateCanvas(
451  "rich_mcbm_fhNofCbmEvents", "rich_mcbm_fhNofCbmEvents", 600, 600);
452  DrawH1(fHM->H1("fhNofCbmEvents"));
453  }
454 
455  {
456  fHM->CreateCanvas(
457  "rich_mcbm_fhNofCbmEventsRing", "rich_mcbm_fhNofCbmEventsRing", 600, 600);
458  DrawH1(fHM->H1("fhNofCbmEventsRing"));
459  }
460 
461  {
462  fHM->CreateCanvas(
463  "rich_mcbm_fhNofEvents", "rich_mcbm_fhNofEvents", 600, 600);
464  DrawH1(fHM->H1("fhNofEvents"));
465  }
466 
467  {
468  fHM->CreateCanvas(
469  "rich_mcbm_fhBlobEvents", "rich_mcbm_fhBlobEvents", 600, 600);
470  DrawH1(fHM->H1("fhNofBlobEvents"));
471  }
472 
473  {
474  fHM->CreateCanvas(
475  "rich_mcbm_fhBlobsInCbmEvent", "rich_mcbm_fhBlobsInCbmEvent", 600, 600);
476  DrawH1(fHM->H1("fhNofBlobsInEvent"));
477  }
478 
479  {
480  fHM->CreateCanvas("inner_channel_delay", "inner_channel_delay", 1200, 600);
481  DrawH2(fHM->H2("fhICD"));
482  }
483 
484 
485  {
486  fHM->CreateCanvas("RingDelta", "RingDelta", 600, 600);
487  DrawH1(fHM->H1("fhRingDeltaTime"));
488  }
489 
490  {
491  TCanvas* c = fHM->CreateCanvas("rich_ToT", "rich_ToT", 1200, 600);
492  c->Divide(2, 1);
493  c->cd(1);
494  DrawH1(fHM->H1("fhRichDigisToT"));
495  c->cd(2);
496  DrawH1(fHM->H1("fhRichHitToT"));
497  }
498 
499  {
500  fHM->CreateCanvas("DigisInChnl", "DigisInChnl", 1200, 600);
501  DrawH2(fHM->H2("fhDigisInChnl"));
502  }
503 
504  {
505  fHM->CreateCanvas("DigisInDiRICH", "DigisInDiRICH", 1200, 600);
506  DrawH2(fHM->H2("fhDigisInDiRICH"));
507  }
508 
509 
510  {
511  TCanvas* c = fHM->CreateCanvas("rich_Hits", "rich_Hits", 1200, 1200);
512  c->Divide(2, 2);
513  c->cd(1);
514  DrawH2(fHM->H2("fhRichHitXY"));
515  c->cd(2);
516  DrawH2(fHM->H2("fhRichHitXY_fromRing"));
517  c->cd(3);
518  TH2F* hitsBg = (TH2F*) fHM->H2("fhRichHitXY")->Clone();
519  hitsBg->SetName("fhRichHitXY_bg");
520  hitsBg->SetTitle("fhRichHitXY_bg");
521  hitsBg->Add(fHM->H2("fhRichHitXY_fromRing"), -1);
522  //hitsBg->Draw("colz");
523  DrawH2(hitsBg);
524  c->cd(4);
525  DrawH2(fHM->H2("fhRichRingXY"));
526  }
527 
528  {
529  TCanvas* c =
530  fHM->CreateCanvas("rich_mcbm_rings", "rich_mcbm_rings", 1200, 600);
531  c->Divide(2, 1);
532  c->cd(1);
533  DrawH1(fHM->H1("fhRichRingRadius"));
534  c->cd(2);
535  DrawH1(fHM->H1("fhNofHitsInRing"));
536  }
537 
538  {
539  fHM->CreateCanvas("RichRingRadiusVsY", "RichRingRadiusVsY", 800, 800);
540  //c->SetLogy();
541  DrawH2(fHM->H2("fhRichRingRadiusY"));
542  }
543 
544  {
545  fHM->CreateCanvas("RichRingHitsVsRadius", "RichRingHitsVsRadius", 800, 800);
546  //c->SetLogy();
547  DrawH2(fHM->H2("fhRichHitsRingRadius"));
548  }
549 
550  {
551  TCanvas* c = fHM->CreateCanvas("RichRingChi2", "RichRingChi2", 1600, 800);
552  //c->SetLogy();
553  c->Divide(2, 1);
554  c->cd(1);
555  DrawH1(fHM->H1("fhRingChi2"));
556  c->cd(2);
557  DrawH2(fHM->H2("fhRichRingRadiusChi2"));
558  }
559 
560  {
561  TCanvas* c =
562  fHM->CreateCanvas("RichRingCenterChi2", "RichRingCenterChi2", 1600, 800);
563  //c->SetLogy();
564  c->Divide(2, 1);
565  c->cd(1);
566  DrawH2(fHM->H2("fhRichRingCenterXChi2"));
567  c->cd(2);
568  DrawH2(fHM->H2("fhRichRingCenterYChi2"));
569  }
570 }
571 
573  //std::cout<<"#!#DRAW!!!"<<std::endl;
574  if (fNofDrawnRings > 20) return;
575  fNofDrawnRings++;
576  stringstream ss;
577  ss << "Event" << fNofDrawnRings;
578  //fNofDrawnRings++;
579  TCanvas* c = nullptr;
580  if (full == true) {
581  c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 800, 800);
582  } else {
583  c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 500, 500);
584  }
585  c->SetGrid(true, true);
586  TH2D* pad = nullptr;
587  if (full == true) {
588  pad = new TH2D(ss.str().c_str(),
589  (ss.str() + ";X [cm];Y [cm]").c_str(),
590  1,
591  -15.,
592  10.,
593  1,
594  -5.,
595  20);
596  } else {
597  pad = new TH2D(ss.str().c_str(),
598  (ss.str() + ";X [cm];Y [cm]").c_str(),
599  1,
600  -5.,
601  5.,
602  1,
603  -5.,
604  5);
605  }
606 
607  pad->SetStats(false);
608  pad->Draw();
609 
610  if (full == true) {
611  //rough Drawing of RichDetectorAcceptance
612  TLine* line0 = new TLine(-6.25, 8, -6.25, 15.9);
613  line0->Draw();
614  TLine* line1 = new TLine(-6.25, 15.9, -1.05, 15.9);
615  line1->Draw();
616  TLine* line2 = new TLine(-1.05, 15.9, -1.05, 13.2);
617  line2->Draw();
618  TLine* line3 = new TLine(-1.05, 13.2, +4.25, 13.2);
619  line3->Draw();
620  TLine* line4 = new TLine(4.25, 13.2, +4.25, 8);
621  line4->Draw();
622  TLine* line5 = new TLine(4.25, 8, -6.25, 8);
623  line5->Draw();
624  }
625 
626  // find min and max x and y positions of the hits
627  // in order to shift drawing
628  double xCur = 0.;
629  double yCur = 0.;
630 
631  if (full == false) {
632  double xmin = 99999., xmax = -99999., ymin = 99999., ymax = -99999.;
633  for (int i = 0; i < ring->GetNofHits(); i++) {
634  Int_t hitInd = ring->GetHit(i);
635  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
636  if (nullptr == hit) continue;
637  if (xmin > hit->GetX()) xmin = hit->GetX();
638  if (xmax < hit->GetX()) xmax = hit->GetX();
639  if (ymin > hit->GetY()) ymin = hit->GetY();
640  if (ymax < hit->GetY()) ymax = hit->GetY();
641  }
642  xCur = (xmin + xmax) / 2.;
643  yCur = (ymin + ymax) / 2.;
644  }
645 
646  //Draw circle and center
647  TEllipse* circle = new TEllipse(
648  ring->GetCenterX() - xCur, ring->GetCenterY() - yCur, ring->GetRadius());
649  circle->SetFillStyle(0);
650  circle->SetLineWidth(3);
651  circle->Draw();
652  TEllipse* center =
653  new TEllipse(ring->GetCenterX() - xCur, ring->GetCenterY() - yCur, .1);
654  center->Draw();
655 
656 
657  double hitZ = 0;
658  uint nofDrawHits = 0;
659 
660  // Draw hits
661  for (int i = 0; i < ring->GetNofHits(); i++) {
662  Int_t hitInd = ring->GetHit(i);
663  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
664  if (nullptr == hit) continue;
665  TEllipse* hitDr = new TEllipse(hit->GetX() - xCur, hit->GetY() - yCur, .25);
666  //std::cout<<"LE of Hit: "<< hit->GetTime()- fCbmEventStartTime << "\t" << hit->GetTime() << "\t" << fCbmEventStartTime <<std::endl;
667  if (doToT(hit)) { // Good ToT selection
668  hitDr->SetFillColor(kRed);
669  } else {
670  hitDr->SetFillColor(kBlue);
671  }
672  hitZ += hit->GetZ();
673  nofDrawHits++;
674  hitDr->Draw();
675  }
676  hitZ /= nofDrawHits;
677 
678 
679  //Draw information
680  stringstream ss2;
681  ss2 << "(x,y,r,n)=(" << setprecision(3) << ring->GetCenterX() << ", "
682  << ring->GetCenterY() << ", " << ring->GetRadius() << ", "
683  << ring->GetNofHits() << ")";
684  TLatex* latex = nullptr;
685  if (full == true) {
686  latex = new TLatex(
687  ring->GetCenterX() - 13., ring->GetCenterY() + 5., ss2.str().c_str());
688  } else {
689  latex = new TLatex(-4., 4., ss2.str().c_str());
690  }
691  latex->Draw();
692 }
693 
694 
696 
697  for (unsigned int i = 0; i < ICD_offset.size(); ++i) {
698  if (ICD_offset_cnt.at(i) == 0) {
699  ICD_offset.at(i) = 0.0;
700  } else {
701  ICD_offset.at(i) /= ICD_offset_cnt.at(i);
702  }
703  ICD_offset.at(i) += ICD_offset_read.at(i);
704  }
705  save_ICD(ICD_offset, 0);
706 
707  //std::cout<<"Address: "<<InterChannel[0].first << std::endl;
708  //std::cout<<"Tracks: "<< fTofTracks->GetEntriesFast() <<std::endl;
709  std::cout << "Drawing Hists...";
710  DrawHist();
711  std::cout << "DONE!" << std::endl;
712 
713  if (this->fDoDrawCanvas) {
715  std::cout << "Canvas saved to Images!" << std::endl;
716  }
717 
718  if (this->fDoWriteHistToFile) {
719  TDirectory* oldir = gDirectory;
720  std::string s = fOutputDir + "/RecoHists.root";
721  TFile* outFile = new TFile(s.c_str(), "RECREATE");
722  if (outFile->IsOpen()) {
723  fHM->WriteToFile();
724  std::cout << "Written to Root-file \"" << s << "\" ...";
725  outFile->Close();
726  std::cout << "Done!" << std::endl;
727  }
728  gDirectory->cd(oldir->GetPath());
729  }
730 }
731 
732 
733 void CbmRichMCbmQaRichOnly::DrawFromFile(const string& fileName,
734  const string& outputDir) {
735  fOutputDir = outputDir;
736 
737  if (fHM != nullptr) delete fHM;
738 
739  fHM = new CbmHistManager();
740  TFile* file = new TFile(fileName.c_str());
741  fHM->ReadFromFile(file);
742  DrawHist();
743 
745 }
746 
748  bool check = false;
749  if ((hit->GetToT() > 23.7) && (hit->GetToT() < 30.0)) check = true;
750 
751  return check;
752 }
753 
755 
756  std::cout << "analyse a Ring" << std::endl;
757 
758  Double_t meanTime = 0.;
759  unsigned int hitCnt = 0;
760  Double_t minRHit2 = std::numeric_limits<Double_t>::max();
761  for (int i = 0; i < ring->GetNofHits(); i++) {
762  Int_t hitInd = ring->GetHit(i);
763  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
764  if (nullptr == hit) continue;
765 
766  meanTime += hit->GetTime();
767  hitCnt++;
768 
769  const Float_t diffX = hit->GetX() - ring->GetCenterX();
770  const Float_t diffY = hit->GetY() - ring->GetCenterY();
771  const Float_t tmpHitRadius2 = (diffX * diffX + diffY * diffY);
772 
773  if (tmpHitRadius2 < minRHit2) { minRHit2 = tmpHitRadius2; }
774  }
775  meanTime = meanTime / hitCnt;
776 
777  //std::cout<<"mean: "<<meanTime<<std::endl;
778  for (int i = 0; i < ring->GetNofHits(); i++) {
779  Int_t hitInd = ring->GetHit(i);
780  CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
781  if (nullptr == hit) continue;
782  //std::cout<<"DeltatTime: "<< meanTime - hit->GetTime()<<std::endl;
783  fHM->H1("fhRingDeltaTime")
784  ->Fill(static_cast<Double_t>(meanTime - hit->GetTime()));
785 
786  fHM->H1("fhRingToTs")->Fill(hit->GetToT());
787  fHM->H1("fhRingLE")
788  ->Fill(static_cast<Double_t>(hit->GetTime() - ev->GetStartTime()));
789  fHM->H2("fhRingLEvsToT")
790  ->Fill(static_cast<Double_t>(hit->GetTime() - ev->GetStartTime()),
791  hit->GetToT());
792  //std::vector<int> tmpRingIndx;
793  //tmpRingIndx.push_back(ring->GetIndex);
794  const Double_t Tdiff_ring = (hit->GetTime() - ev->GetStartTime());
795  if ((Tdiff_ring > 20.) && (Tdiff_ring < 30.)) {
796  std::cout << ev->GetNumber() << " Address_ring: " << std::hex
797  << CbmRichUtil::GetDirichId(hit->GetAddress()) << std::dec
798  << " " << CbmRichUtil::GetDirichChannel(hit->GetAddress())
799  << " " << hit->GetToT() << " " << ring->GetRadius()
800  << std::endl;
801  //fHM->H1("fhDiRICHsInRegion")->Fill(CbmRichUtil::GetDirichId(hit->GetAddress()));
802  }
803  }
804 
805  int InnerHitCnt = 0;
806  int InnerHitCnt_cut = 0;
807  for (int j = 0; j < ev->GetNofData(ECbmDataType::kRichHit); j++) {
808  auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
809  CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
810  if (nullptr == richHit) continue;
811  const Float_t diffX = richHit->GetX() - ring->GetCenterX();
812  const Float_t diffY = richHit->GetY() - ring->GetCenterY();
813  //select inner Part of Ring
814  if (diffX * diffX + diffY * diffY < minRHit2) {
815  InnerHitCnt++;
816  const Double_t Tdiff_inner = (richHit->GetTime() - ev->GetStartTime());
817  if ((Tdiff_inner > 20.) && (Tdiff_inner < 30.)) {
818  InnerHitCnt_cut++;
819  //if (InnerHitCnt_cut == 1) {DrawRing(ring);}
820  std::cout << ev->GetNumber() << " Address_inner: " << std::hex
821  << CbmRichUtil::GetDirichId(richHit->GetAddress()) << std::dec
822  << " "
824  << " " << richHit->GetToT() << " " << ring->GetRadius()
825  << std::endl;
826  fHM->H1("fhDiRICHsInRegion")
827  ->Fill(CbmRichUtil::GetDirichId(richHit->GetAddress()));
828  }
829 
830  fHM->H1("fhInnerRingDeltaTime")
831  ->Fill(static_cast<Double_t>(meanTime - richHit->GetTime()));
832  fHM->H1("fhInnerRingToTs")->Fill(richHit->GetToT());
833  fHM->H1("fhInnerRingLE")
834  ->Fill(static_cast<Double_t>(richHit->GetTime() - ev->GetStartTime()));
835  }
836  }
837  if (InnerHitCnt == 0) {
838  fHM->H1("fhInnerRingFlag")->Fill(1);
839  } else {
840  fHM->H1("fhInnerRingFlag")->Fill(0);
841  }
842  fHM->H1("fhNofInnerHits")->Fill(InnerHitCnt);
843 }
844 
845 
847  if (ring->GetRadius() > 2. && ring->GetRadius() < 4.2) return true;
848 
849  return false;
850 }
851 
852 void CbmRichMCbmQaRichOnly::save_ICD(std::array<Double_t, 2304>& ICD_offsets,
853  unsigned int iteration) {
854  std::ofstream icd_file;
855  icd_file.open(Form("icd_offset_it_%u.data", iteration));
856  if (icd_file.is_open()) {
857  for (unsigned int i = 0; i < ICD_offsets.size(); ++i) {
858  //loop over all Offsets
859  icd_file << i << "\t" << std::setprecision(25) << ICD_offsets.at(i)
860  << "\n";
861  }
862  icd_file.close();
863  } else
864  std::cout << "Unable to open inter channel delay file icd_offset_it_"
865  << iteration << ".data\n";
866 }
867 
868 
869 void CbmRichMCbmQaRichOnly::read_ICD(std::array<Double_t, 2304>& ICD_offsets,
870  unsigned int iteration) {
871  std::cout << gSystem->pwd() << std::endl;
872  std::string line;
873  std::ifstream icd_file(Form("icd_offset_it_%u.data", iteration));
874  unsigned int lineCnt = 0;
875  if (icd_file.is_open()) {
876  while (getline(icd_file, line)) {
877  if (line[0] == '#') continue; // just a comment
878  std::istringstream iss(line);
879  unsigned int addr = 0;
880  Double_t value;
881  if (!(iss >> addr >> value)) {
882  std::cout << "A Problem accured in line " << lineCnt << "\n";
883  break;
884  } // error
885  lineCnt++;
886  ICD_offsets.at(addr) += value;
887  }
888  icd_file.close();
889  } else
890  std::cout << "Unable to open inter channel delay file icd_offset_it_"
891  << iteration << ".data\n";
892 }
893 
CbmRichPoint.h
CbmRichRing::GetTime
Double_t GetTime() const
Definition: CbmRichRing.h:102
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmRichMCbmQaRichOnly
Definition: CbmRichMCbmQaRichOnly.h:24
CbmEvent::GetIndex
UInt_t GetIndex(ECbmDataType type, UInt_t iData)
Definition: CbmEvent.cxx:24
CbmRichMCbmQaRichOnly::fDoWriteHistToFile
bool fDoWriteHistToFile
Definition: CbmRichMCbmQaRichOnly.h:132
CbmRichRingFinderHoughImpl.h
Ring finder implementation based on Hough Transform method.
CbmRichMCbmQaRichOnly::ICD_offset_cnt
std::array< uint32_t, 2304 > ICD_offset_cnt
Definition: CbmRichMCbmQaRichOnly.h:137
CbmRichMCbmQaRichOnly::cutRadius
Bool_t cutRadius(CbmRichRing *ring)
Definition: CbmRichMCbmQaRichOnly.cxx:846
CbmHistManager::WriteToFile
void WriteToFile()
Write all histograms to current opened file.
Definition: core/base/CbmHistManager.cxx:103
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmRichMCbmQaRichOnly::fDigiMan
CbmDigiManager * fDigiMan
Definition: CbmRichMCbmQaRichOnly.h:102
CbmHistManager::ScaleByPattern
void ScaleByPattern(const std::string &pattern, Double_t scale)
Scale histograms which name matches specified pattern.
Definition: core/base/CbmHistManager.cxx:221
CbmRichDigi
Definition: CbmRichDigi.h:25
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmRichMCbmSEDisplay
Definition: CbmRichMCbmSEDisplay.h:17
CbmHistManager::ReadFromFile
void ReadFromFile(TFile *file)
Read histograms from file.
Definition: core/base/CbmHistManager.cxx:110
CbmRichHit::GetToT
Double_t GetToT() const
Definition: CbmRichHit.h:62
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition: CbmDigiManager.cxx:71
CbmRichMCbmQaRichOnly::fMaxNofDrawnEvents
Int_t fMaxNofDrawnEvents
Definition: CbmRichMCbmQaRichOnly.h:122
CbmHistManager::Create2
void Create2(const std::string &name, const std::string &title, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY)
Helper function for creation of 2-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:104
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRing::GetChi2
Double_t GetChi2() const
Definition: CbmRichRing.h:95
CbmRichMCbmQaRichOnly::analyseRing
void analyseRing(CbmRichRing *ring, CbmEvent *ev)
Definition: CbmRichMCbmQaRichOnly.cxx:754
CbmRichUtil::GetDirichId
static uint16_t GetDirichId(Int_t Address)
Definition: utils/CbmRichUtil.h:35
CbmRichMCbmQaRichOnly::fEventNum
Int_t fEventNum
Definition: CbmRichMCbmQaRichOnly.h:114
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition: CbmHistManager.h:190
CbmGlobalTrack.h
ECbmDataType::kRichHit
@ kRichHit
CbmRichMCbmSEDisplay::SetRichRings
void SetRichRings(TClonesArray *ring=nullptr)
Definition: CbmRichMCbmSEDisplay.h:63
CbmDigiManager::GetNofDigis
static Int_t GetNofDigis(ECbmModuleId systemId)
Definition: CbmDigiManager.cxx:62
CbmRichRing::GetNofHits
Int_t GetNofHits() const
Definition: CbmRichRing.h:40
CbmRichRing
Definition: CbmRichRing.h:17
CbmTofDigi.h
CbmRichRing.h
CbmRichMCbmQaRichOnly::save_ICD
void save_ICD(std::array< Double_t, 2304 > &offsets, unsigned int iteration)
Definition: CbmRichMCbmQaRichOnly.cxx:852
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichDraw.h
CbmHistManager.h
Histogram manager.
CbmMatchRecoToMC.h
FairTask for matching RECO data to MC.
CbmDigiManager::IsPresent
static Bool_t IsPresent(ECbmModuleId systemId)
Presence of a digi branch.
Definition: CbmDigiManager.cxx:112
CbmRichMCbmQaRichOnly::fNofDrawnRings
Int_t fNofDrawnRings
Definition: CbmRichMCbmQaRichOnly.h:116
CbmTofTracklet.h
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition: CbmDigiManager.h:93
CbmRichGeoManager.h
CbmRichDigi.h
CbmRichMCbmQaRichOnly::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichMCbmQaRichOnly.cxx:80
CbmEvent.h
CbmRichMCbmSEDisplay::SetOutDir
void SetOutDir(std::string dir)
Definition: CbmRichMCbmSEDisplay.h:85
DrawH1
void DrawH1(TH1 *hist, HistScale logx, HistScale logy, const string &drawOpt, Int_t color, Int_t lineWidth, Int_t lineStyle, Int_t markerSize, Int_t markerStyle)
Definition: CbmDrawHist.cxx:49
CbmRichMCbmQaRichOnly::fOutputDir
string fOutputDir
Definition: CbmRichMCbmQaRichOnly.h:127
CbmRichRing::GetHit
UInt_t GetHit(Int_t i) const
Definition: CbmRichRing.h:42
CbmStsDigi.h
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition: CbmDigiManager.h:52
CbmRichMCbmQaRichOnly::fDoDrawCanvas
bool fDoDrawCanvas
Definition: CbmRichMCbmQaRichOnly.h:133
CbmRichConverter.h
Convert internal data classes to cbmroot common data classes.
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
CbmTrackMatchNew.h
CbmEvent::GetNumber
Int_t GetNumber() const
Definition: CbmEvent.h:110
CbmHistManager::Create1
void Create1(const std::string &name, const std::string &title, Int_t nofBins, Double_t minBin, Double_t maxBin)
Helper function for creation of 1-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:81
CbmRichMCbmQaRichOnly::read_ICD
void read_ICD(std::array< Double_t, 2304 > &offsets, unsigned int iteration)
Definition: CbmRichMCbmQaRichOnly.cxx:869
CbmRichMCbmQaRichOnly::InitHistograms
void InitHistograms()
Initialize histograms.
Definition: CbmRichMCbmQaRichOnly.cxx:132
ECbmModuleId::kRich
@ kRich
Ring-Imaging Cherenkov Detector.
CbmRichMCbmQaRichOnly::DrawRing
void DrawRing(CbmRichRing *ring)
Definition: CbmRichMCbmQaRichOnly.h:156
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition: CbmHistManager.h:170
CbmRichMCbmQaRichOnly::fTriggerRichHits
Int_t fTriggerRichHits
Definition: CbmRichMCbmQaRichOnly.h:124
CbmRichMCbmQaRichOnly.h
CbmUtils.h
CbmHistManager::CreateCanvas
TCanvas * CreateCanvas(const std::string &name, const std::string &title, Int_t width, Int_t height)
Create and draw TCanvas and store pointer to it.
Definition: core/base/CbmHistManager.cxx:267
CbmRichMCbmSEDisplay::DrawEvent
void DrawEvent(CbmEvent *ev, std::vector< int > &ringIndx, bool full)
Draw histograms.
Definition: CbmRichMCbmSEDisplay.cxx:47
CbmEvent::GetNofData
Int_t GetNofData() const
Definition: CbmEvent.h:90
CbmRichDigi::GetAddress
Int_t GetAddress() const
Definition: CbmRichDigi.h:37
CbmRichUtil::GetDirichChannel
static uint16_t GetDirichChannel(Int_t Address)
Definition: utils/CbmRichUtil.h:39
CbmHistManager::SaveCanvasToImage
void SaveCanvasToImage(const std::string &outputDir, const std::string &options="png,eps")
Save all stored canvases to images.
Definition: core/base/CbmHistManager.cxx:276
CbmRichMCbmQaRichOnly::DrawHist
void DrawHist()
Draw histograms.
Definition: CbmRichMCbmQaRichOnly.cxx:442
CbmRichMCbmQaRichOnly::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichMCbmQaRichOnly.cxx:268
CbmRichMCbmQaRichOnly::RichRings
void RichRings()
Definition: CbmRichMCbmQaRichOnly.cxx:407
CbmRichMCbmQaRichOnly::fRichHits
TClonesArray * fRichHits
Definition: CbmRichMCbmQaRichOnly.h:104
CbmRichMCbmQaRichOnly::fRichRings
TClonesArray * fRichRings
Definition: CbmRichMCbmQaRichOnly.h:106
CbmRichMCbmQaRichOnly::ICD_offset
std::array< Double_t, 2304 > ICD_offset
Definition: CbmRichMCbmQaRichOnly.h:136
CbmDigiManager.h
CbmRichMCbmSEDisplay::SetMaxNofDrawnEvents
void SetMaxNofDrawnEvents(Int_t val=100)
Definition: CbmRichMCbmSEDisplay.h:73
CbmRichMCbmQaRichOnly::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichMCbmQaRichOnly.cxx:695
CbmRichMCbmQaRichOnly::fCbmEvent
TClonesArray * fCbmEvent
Definition: CbmRichMCbmQaRichOnly.h:108
CbmRichMCbmQaRichOnly::fXOffsetHisto
Double_t fXOffsetHisto
Definition: CbmRichMCbmQaRichOnly.h:112
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
CbmRichMCbmQaRichOnly::DrawFromFile
void DrawFromFile(const string &fileName, const string &outputDir)
Draw histogram from file.
Definition: CbmRichMCbmQaRichOnly.cxx:733
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition: CbmEvent.h:30
CbmRichRing::GetRadius
Float_t GetRadius() const
Definition: CbmRichRing.h:82
CbmRichRing::GetCenterY
Float_t GetCenterY() const
Definition: CbmRichRing.h:81
DrawH2
void DrawH2(TH2 *hist, HistScale logx, HistScale logy, HistScale logz, const string &drawOpt)
Definition: CbmDrawHist.cxx:84
CbmRichMCbmQaRichOnly::fSeDisplay
CbmRichMCbmSEDisplay * fSeDisplay
Definition: CbmRichMCbmQaRichOnly.h:139
CbmRichDigi::GetToT
Double_t GetToT() const
Definition: CbmRichDigi.h:67
CbmRichMCbmQaRichOnly::ICD_offset_read
std::array< Double_t, 2304 > ICD_offset_read
Definition: CbmRichMCbmQaRichOnly.h:135
CbmRichHit.h
CbmTrdTrack.h
CbmRichMCbmQaRichOnly::doToT
bool doToT(CbmRichHit *hit)
Definition: CbmRichMCbmQaRichOnly.cxx:747
CbmRichMCbmSEDisplay.h
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
CbmRichMCbmSEDisplay::SetRichHits
void SetRichHits(TClonesArray *hits=nullptr)
Definition: CbmRichMCbmSEDisplay.h:58
CbmRichRing::GetCenterX
Float_t GetCenterX() const
Definition: CbmRichRing.h:80
CbmRichMCbmSEDisplay::SetTotRich
void SetTotRich(Double_t min, Double_t max)
Definition: CbmRichMCbmSEDisplay.h:50
CbmRichMCbmQaRichOnly::fHM
CbmHistManager * fHM
Definition: CbmRichMCbmQaRichOnly.h:110
CbmRichHit
Definition: CbmRichHit.h:19
CbmRichMCbmQaRichOnly::CbmRichMCbmQaRichOnly
CbmRichMCbmQaRichOnly()
Standard constructor.
Definition: CbmRichMCbmQaRichOnly.cxx:66
CbmEvent::GetStartTime
Double_t GetStartTime() const
Definition: CbmEvent.h:131