4 #include "TClonesArray.h"
9 #include "TGeoManager.h"
34 #include "CbmRichUtil.h"
37 #include "CbmTofHit.h"
52 #include <boost/assign/list_of.hpp>
62 using boost::assign::list_of;
67 : FairTask(
"CbmRichMCbmQaRichOnly")
74 , fNofDrawnRichTofEv(0)
76 , fMaxNofDrawnEvents(100)
78 , fOutputDir(
"result") {}
81 cout <<
"CbmRichMCbmQaRichOnly::Init" << endl;
83 FairRootManager* ioman = FairRootManager::Instance();
84 if (
nullptr == ioman) {
85 Fatal(
"CbmRichMCbmQaRichOnly::Init",
"RootManager not instantised!");
92 Fatal(
"CbmRichMCbmQaReal::Init",
"No Rich Digis!");
94 fRichHits = (TClonesArray*) ioman->GetObject(
"RichHit");
96 Fatal(
"CbmRichMCbmQaRichOnly::Init",
"No Rich Hits!");
99 fRichRings = (TClonesArray*) ioman->GetObject(
"RichRing");
101 Fatal(
"CbmRichMCbmQaRichOnly::Init",
"No Rich Rings!");
104 fCbmEvent = (TClonesArray*) ioman->GetObject(
"CbmEvent");
106 Fatal(
"CbmRichMCbmQaRichOnly::Init",
"No Event!");
135 fHM->
Create1<TH1D>(
"fhNofEvents",
"fhNofEvents;Entries", 1, 0.5, 1.5);
136 fHM->
Create1<TH1D>(
"fhNofCbmEvents",
"fhNofCbmEvents;Entries", 1, 0.5, 1.5);
138 "fhNofCbmEventsRing",
"fhNofCbmEventsRing;Entries", 1, 0.5, 1.5);
140 fHM->
Create1<TH1D>(
"fhNofBlobEvents",
"fhNofBlobEvents;Entries", 1, 0.5, 1.5);
142 "fhNofBlobsInEvent",
"fhNofBlobsInEvent;Entries", 36, 0.5, 36.5);
146 "fhRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries",
154 "fhRichHitXY_fromRing",
155 "fhRichHitXY_fromRing;RICH hit X [cm];RICH hit Y [cm];Entries",
165 "fhRichDigisToT",
"fhRichDigisToT;ToT [ns];Entries", 601, 9.975, 40.025);
167 "fhRichHitToT",
"fhRichHitToT;ToT [ns];Entries", 601, 9.975, 40.025);
172 "fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries",
180 "fhRichRingRadius;Ring radius [cm];Entries",
185 "fhNofHitsInRing;# hits in ring;Entries",
190 "fhICD",
"fhICD;channel;DeltaTime", 2305, -0.5, 2304.5, 31, -15.5, 15.5);
194 "fhRichRingRadiusY;Ring Radius [cm]; Y position[cm];Entries",
202 "fhRichHitsRingRadius",
203 "fhRichHitsRingRadius;#Rich Hits/Ring; Ring Radius [cm];Entries",
212 "fhRingDeltaTime; \\Delta Time/ns;Entries",
217 "fhRingChi2",
"fhRingChi2; \\Chi^2 ;Entries", 101, 0.0, 10.1);
220 "fhRichRingCenterXChi2",
221 "fhRichRingCenterXChi2;Ring Center X [cm];\\Chi^2 ;;Entries",
229 "fhRichRingCenterYChi2",
230 "fhRichRingCenterYChi2;Ring Center Y [cm];\\Chi^2 ;;Entries",
238 "fhRichRingRadiusChi2",
239 "fhRichRingRadiusChi2; Ring Radius [cm];\\Chi^2 ;;Entries",
250 "fhDigisInChnl;channel;#Digis;",
258 "fhDigisInDiRICH;DiRICH;#Digis;",
270 fHM->
H1(
"fhNofEvents")->Fill(1);
272 cout <<
"CbmRichMCbmQaRichOnly, event No. " <<
fEventNum << endl;
274 std::array<unsigned int, 2304> chnlDigis;
275 for (
auto& c : chnlDigis)
279 fHM->
H1(
"fhRichDigisToT")->Fill(richDigi->
GetToT());
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]++;
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]);
296 uint16_t dirich =
i / 32;
297 if (sum != 0)
fHM->
H1(
"fhDigisInDiRICH")->Fill(dirich, sum);
303 int nofRichHits =
fRichHits->GetEntries();
304 for (
int iH = 0; iH < nofRichHits; iH++) {
306 fHM->
H2(
"fhRichHitXY")->Fill(richHit->
GetX(), richHit->
GetY());
312 auto fNCbmEvent =
fCbmEvent->GetEntriesFast();
314 for (
int i = 0;
i < fNCbmEvent;
i++) {
315 fHM->
H1(
"fhNofCbmEvents")->Fill(1);
323 std::vector<int> ringIndx;
324 std::vector<int> evRichHitIndx;
325 std::array<uint32_t, 36> pmtHits;
326 for (
auto& a : pmtHits)
332 evRichHitIndx.push_back(iRichHit);
335 uint32_t pmtId = (((richHit->
GetAddress()) >> 20) & 0xF)
336 + (((richHit->
GetAddress()) >> 24) & 0xF) * 9;
340 for (
int l = 0; l < nofRichRings; l++) {
344 for (
int m = 0;
m < NofRingHits;
m++) {
345 auto RingHitIndx = ring->
GetHit(
m);
346 if (RingHitIndx == iRichHit) {
348 for (
auto check : ringIndx) {
354 if (used ==
false) ringIndx.push_back(l);
362 for (
auto a : pmtHits) {
363 if (a > 30) { blob++; }
366 fHM->
H1(
"fhNofBlobEvents")->Fill(1);
367 fHM->
H1(
"fhNofBlobsInEvent")->Fill(blob);
370 if (ringIndx.size() != 0)
fHM->
H1(
"fhNofCbmEventsRing")->Fill(1);
373 for (
unsigned int k = 0; k < ringIndx.size(); ++k) {
377 if (ring ==
nullptr)
continue;
379 fHM->
H2(
"fhRichRingCenterXChi2")
381 fHM->
H2(
"fhRichRingCenterYChi2")
389 fHM->
H2(
"fhRichHitsRingRadius")
391 for (
int j = 0; j < ring->
GetNofHits(); ++j) {
392 Int_t hitIndx = ring->
GetHit(j);
394 if (
nullptr == hit)
continue;
395 fHM->
H2(
"fhRichHitXY_fromRing")->Fill(hit->
GetX(), hit->
GetY());
409 for (
int i = 0;
i < nofRichRings;
i++) {
411 if (ring ==
nullptr)
continue;
413 for (
int j = 1; j < ring->
GetNofHits(); ++j) {
414 Int_t hitIndx = ring->
GetHit(j);
416 if (
nullptr == hit)
continue;
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);
426 fHM->
H2(
"fhICD")->Fill(
428 fHM->
H1(
"fhRingDeltaTime")
446 double nofEvents =
fHM->
H1(
"fhNofCbmEvents")->GetEntries();
451 "rich_mcbm_fhNofCbmEvents",
"rich_mcbm_fhNofCbmEvents", 600, 600);
457 "rich_mcbm_fhNofCbmEventsRing",
"rich_mcbm_fhNofCbmEventsRing", 600, 600);
463 "rich_mcbm_fhNofEvents",
"rich_mcbm_fhNofEvents", 600, 600);
469 "rich_mcbm_fhBlobEvents",
"rich_mcbm_fhBlobEvents", 600, 600);
475 "rich_mcbm_fhBlobsInCbmEvent",
"rich_mcbm_fhBlobsInCbmEvent", 600, 600);
480 fHM->
CreateCanvas(
"inner_channel_delay",
"inner_channel_delay", 1200, 600);
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);
545 fHM->
CreateCanvas(
"RichRingHitsVsRadius",
"RichRingHitsVsRadius", 800, 800);
551 TCanvas* c =
fHM->
CreateCanvas(
"RichRingChi2",
"RichRingChi2", 1600, 800);
562 fHM->
CreateCanvas(
"RichRingCenterChi2",
"RichRingCenterChi2", 1600, 800);
579 TCanvas* c =
nullptr;
585 c->SetGrid(
true,
true);
588 pad =
new TH2D(ss.str().c_str(),
589 (ss.str() +
";X [cm];Y [cm]").c_str(),
597 pad =
new TH2D(ss.str().c_str(),
598 (ss.str() +
";X [cm];Y [cm]").c_str(),
607 pad->SetStats(
false);
612 TLine* line0 =
new TLine(-6.25, 8, -6.25, 15.9);
614 TLine* line1 =
new TLine(-6.25, 15.9, -1.05, 15.9);
616 TLine* line2 =
new TLine(-1.05, 15.9, -1.05, 13.2);
618 TLine* line3 =
new TLine(-1.05, 13.2, +4.25, 13.2);
620 TLine* line4 =
new TLine(4.25, 13.2, +4.25, 8);
622 TLine* line5 =
new TLine(4.25, 8, -6.25, 8);
632 double xmin = 99999., xmax = -99999., ymin = 99999., ymax = -99999.;
634 Int_t hitInd = ring->
GetHit(
i);
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();
642 xCur = (xmin + xmax) / 2.;
643 yCur = (ymin + ymax) / 2.;
647 TEllipse* circle =
new TEllipse(
649 circle->SetFillStyle(0);
650 circle->SetLineWidth(3);
658 uint nofDrawHits = 0;
662 Int_t hitInd = ring->
GetHit(
i);
664 if (
nullptr == hit)
continue;
665 TEllipse* hitDr =
new TEllipse(hit->
GetX() - xCur, hit->
GetY() - yCur, .25);
668 hitDr->SetFillColor(kRed);
670 hitDr->SetFillColor(kBlue);
681 ss2 <<
"(x,y,r,n)=(" << setprecision(3) << ring->
GetCenterX() <<
", "
684 TLatex* latex =
nullptr;
689 latex =
new TLatex(-4., 4., ss2.str().c_str());
709 std::cout <<
"Drawing Hists...";
711 std::cout <<
"DONE!" << std::endl;
715 std::cout <<
"Canvas saved to Images!" << std::endl;
719 TDirectory* oldir = gDirectory;
720 std::string s =
fOutputDir +
"/RecoHists.root";
721 TFile* outFile =
new TFile(s.c_str(),
"RECREATE");
722 if (outFile->IsOpen()) {
724 std::cout <<
"Written to Root-file \"" << s <<
"\" ...";
726 std::cout <<
"Done!" << std::endl;
728 gDirectory->cd(oldir->GetPath());
734 const string& outputDir) {
737 if (
fHM !=
nullptr)
delete fHM;
740 TFile* file =
new TFile(fileName.c_str());
749 if ((hit->
GetToT() > 23.7) && (hit->
GetToT() < 30.0)) check =
true;
756 std::cout <<
"analyse a Ring" << std::endl;
758 Double_t meanTime = 0.;
759 unsigned int hitCnt = 0;
762 Int_t hitInd = ring->
GetHit(
i);
764 if (
nullptr == hit)
continue;
771 const Float_t tmpHitRadius2 = (diffX * diffX + diffY * diffY);
773 if (tmpHitRadius2 < minRHit2) { minRHit2 = tmpHitRadius2; }
775 meanTime = meanTime / hitCnt;
779 Int_t hitInd = ring->
GetHit(
i);
781 if (
nullptr == hit)
continue;
783 fHM->
H1(
"fhRingDeltaTime")
784 ->Fill(
static_cast<Double_t
>(meanTime - hit->
GetTime()));
789 fHM->
H2(
"fhRingLEvsToT")
795 if ((Tdiff_ring > 20.) && (Tdiff_ring < 30.)) {
796 std::cout << ev->
GetNumber() <<
" Address_ring: " << std::hex
806 int InnerHitCnt_cut = 0;
810 if (
nullptr == richHit)
continue;
814 if (diffX * diffX + diffY * diffY < minRHit2) {
817 if ((Tdiff_inner > 20.) && (Tdiff_inner < 30.)) {
820 std::cout << ev->
GetNumber() <<
" Address_inner: " << std::hex
826 fHM->
H1(
"fhDiRICHsInRegion")
830 fHM->
H1(
"fhInnerRingDeltaTime")
831 ->Fill(
static_cast<Double_t
>(meanTime - richHit->
GetTime()));
832 fHM->
H1(
"fhInnerRingToTs")->Fill(richHit->
GetToT());
833 fHM->
H1(
"fhInnerRingLE")
837 if (InnerHitCnt == 0) {
838 fHM->
H1(
"fhInnerRingFlag")->Fill(1);
840 fHM->
H1(
"fhInnerRingFlag")->Fill(0);
842 fHM->
H1(
"fhNofInnerHits")->Fill(InnerHitCnt);
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) {
859 icd_file <<
i <<
"\t" << std::setprecision(25) << ICD_offsets.at(
i)
864 std::cout <<
"Unable to open inter channel delay file icd_offset_it_"
865 << iteration <<
".data\n";
870 unsigned int iteration) {
871 std::cout << gSystem->pwd() << std::endl;
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;
878 std::istringstream iss(line);
879 unsigned int addr = 0;
881 if (!(iss >> addr >> value)) {
882 std::cout <<
"A Problem accured in line " << lineCnt <<
"\n";
886 ICD_offsets.at(addr) += value;
890 std::cout <<
"Unable to open inter channel delay file icd_offset_it_"
891 << iteration <<
".data\n";