From 30efdc7ad5b6b1ef951ef5db4fa203192fa28d6a Mon Sep 17 00:00:00 2001
From: Florian Uhlig <f.uhlig@gsi.de>
Date: Mon, 6 Jul 2020 09:48:45 +0200
Subject: [PATCH] Add changes in Rich code and macros
The commit ports the commits 16307, 16385-16386, 16426-16429, 16452, 16454,
16500, 16548 and 16566 from the mcbm2020 development branch.
---
.../unpacker/CbmMcbm2018UnpackerTaskRich.cxx | 2 +-
macro/rich/mcbm/beamtime/getToTOffset.C | 22 +-
.../rich/mcbm/beamtime/run_aerogel_analysis.C | 124 ++
.../beamtime/run_reco_mcbm_real_wToF_dec19.C | 582 ++++++++
.../beamtime/run_reco_mcbm_real_wToF_mar20.C | 541 +++++++
.../mcbm/beamtime/run_reco_mcbm_richOnly.C | 130 ++
reco/detectors/rich/CMakeLists.txt | 5 +-
reco/detectors/rich/CbmRichRecoLinkDef.h | 5 +-
.../rich/mcbm/CbmRichMCbmAerogelAna.cxx | 487 +++++++
.../rich/mcbm/CbmRichMCbmAerogelAna.h | 150 ++
.../rich/mcbm/CbmRichMCbmHitProducer.cxx | 47 +-
.../rich/mcbm/CbmRichMCbmHitProducer.h | 14 +-
.../detectors/rich/mcbm/CbmRichMCbmQaReal.cxx | 1259 ++++++++++++-----
reco/detectors/rich/mcbm/CbmRichMCbmQaReal.h | 96 +-
.../rich/mcbm/CbmRichMCbmQaRichOnly.cxx | 757 ++++++++++
.../rich/mcbm/CbmRichMCbmQaRichOnly.h | 208 +++
.../rich/mcbm/CbmRichMCbmSEDisplay.cxx | 282 ++++
.../rich/mcbm/CbmRichMCbmSEDisplay.h | 165 +++
.../rich/mcbm/CbmRichMCbmToTShifter.cxx | 41 +-
.../rich/mcbm/CbmRichMCbmToTShifter.h | 13 +-
20 files changed, 4567 insertions(+), 363 deletions(-)
create mode 100644 macro/rich/mcbm/beamtime/run_aerogel_analysis.C
create mode 100644 macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_dec19.C
create mode 100644 macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_mar20.C
create mode 100644 macro/rich/mcbm/beamtime/run_reco_mcbm_richOnly.C
create mode 100644 reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.cxx
create mode 100644 reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.h
create mode 100644 reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.cxx
create mode 100644 reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.h
create mode 100644 reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.cxx
create mode 100644 reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.h
diff --git a/fles/mcbm2018/unpacker/CbmMcbm2018UnpackerTaskRich.cxx b/fles/mcbm2018/unpacker/CbmMcbm2018UnpackerTaskRich.cxx
index f877f6da..009300c1 100644
--- a/fles/mcbm2018/unpacker/CbmMcbm2018UnpackerTaskRich.cxx
+++ b/fles/mcbm2018/unpacker/CbmMcbm2018UnpackerTaskRich.cxx
@@ -83,7 +83,7 @@ Bool_t CbmMcbm2018UnpackerTaskRich::DoUnpack(const fles::Timeslice& ts, size_t c
/*
/// Sort the buffers of hits due to the time offsets applied
- => Done in the algo!!!
+ //=> Done in the algo!!!
sort(fpvDigiRich->begin(), fpvDigiRich->end(),
[](const CbmRichDigi & a, const CbmRichDigi & b) -> bool
{
diff --git a/macro/rich/mcbm/beamtime/getToTOffset.C b/macro/rich/mcbm/beamtime/getToTOffset.C
index 708b2111..7888f402 100644
--- a/macro/rich/mcbm/beamtime/getToTOffset.C
+++ b/macro/rich/mcbm/beamtime/getToTOffset.C
@@ -12,10 +12,22 @@
// In order to call later Finish, we make this global
FairRunOnline *run = nullptr;
-void getToTOffset(TString inFile = "/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn02_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn04_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn05_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn06_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn08_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn10_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn11_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn12_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn13_*.tsa;/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_pn15_*.tsa", UInt_t uRunId = 370, UInt_t nrEvents=10000, TString outDir="/lustre/nyx/cbm/users/adrian/data/ToTOffset/", TString inDir="") //1Event is 1TS
+void getToTOffset( UInt_t uRunId = 831, UInt_t nrEvents=10000, TString outDir="./data/ToTOffset/", TString inDir="") //1Event is 1TS
{
TString srcDir = gSystem->Getenv("VMCWORKDIR");
+ TString inputDir = "/lustre/cbm/users/ploizeau/mcbm2020/data";
+ TString inFile = Form("%s/%u_pn02_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn04_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn05_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn06_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn08_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn10_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn11_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn12_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn13_*.tsa;",inputDir.Data(),uRunId);
+ inFile += Form("%s/%u_pn15_*.tsa" ,inputDir.Data(),uRunId);
+
// --- Specify number of events to be produced.
// --- -1 means run until the end of the input file.
Int_t nEvents=-1;
@@ -33,7 +45,7 @@ void getToTOffset(TString inFile = "/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_
// --- Define parameter files
TList *parFileList = new TList();
- TString paramDir = srcDir + "/macro/beamtime/mcbm2019/";
+ TString paramDir = srcDir + "/macro/beamtime/mcbm2020/";
TString paramFileRich = paramDir + "mRichPar.par";
TObjString* parRichFileName = new TObjString(paramFileRich);
@@ -66,7 +78,7 @@ void getToTOffset(TString inFile = "/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_
source->SetFileName(inFile);
source->SetInputDir(inDir);
- source->AddUnpacker(unpacker_rich, 0x30, kRich );//RICH trb
+ source->AddUnpacker(unpacker_rich, 0x30, ECbmModuleId::kRich );//RICH trb
// --- Event header
FairEventHeader* event = new CbmTbEvent();
@@ -87,8 +99,8 @@ void getToTOffset(TString inFile = "/lustre/nyx/cbm/users/ploizeau/mcbm2019/370_
// Add ToT Correction Finder
CbmRichMCbmToTShifter *tot = new CbmRichMCbmToTShifter();
- tot->GeneratePDF();
- tot->ShowTdcId(false);
+ //tot->GeneratePDF();
+ tot->ShowTdcId(true);
run->AddTask(tot);
diff --git a/macro/rich/mcbm/beamtime/run_aerogel_analysis.C b/macro/rich/mcbm/beamtime/run_aerogel_analysis.C
new file mode 100644
index 00000000..ec54a0d1
--- /dev/null
+++ b/macro/rich/mcbm/beamtime/run_aerogel_analysis.C
@@ -0,0 +1,124 @@
+void run_aerogel_analysis(
+ const string& parFile = "/lustre/cbm/users/adrian/cbmgitnew/cbmsource/macro/beamtime/mcbm2020/data/unp_mcbm_params_598.root",
+ const string& digiFile = "/lustre/cbm/users/adrian/cbmgitnew/cbmsource/macro/beamtime/mcbm2020/data/unp_mcbm_598.root",
+ const string& recoFile = "reco_aerogel_mcbm.root",
+ int nEvents = 1000
+)
+{
+ const unsigned int runId = 598;
+
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ TTree::SetMaxTreeSize(90000000000);
+
+ TString myName = "run_reco_mcbm_real";
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ TString workDir = gSystem->Getenv("VMCWORKDIR");
+
+ remove(recoFile.c_str());
+
+ //TString TofFileFolder = Form("/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/%s",cFileId.Data());
+
+ std::cout << std::endl<< "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList *parFileList = new TList();
+
+ TStopwatch timer;
+ timer.Start();
+ gDebug = 0;
+
+
+ TString geoDir = workDir;
+ TString geoFile = srcDir + "/macro/mcbm/data/mcbm_beam_2020_03.geo.root";
+ TFile* fgeo = new TFile(geoFile);
+ TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+ if (NULL == geoMan){
+ cout << "<E> FAIRGeom not found in geoFile"<<endl;
+ return;
+ }
+
+ //-----------------------------------------------//
+
+ FairRunAna *run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(digiFile.c_str());
+ run->SetSource(inputSource);
+ run->SetOutputFile(recoFile.c_str());
+ run->SetGenerateRunInfo(kTRUE);
+
+
+ /// Add the Eventbuilder
+
+ CbmMcbm2018EventBuilder* eventBuilder = new CbmMcbm2018EventBuilder();
+ // eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::MaximumTimeGap);
+ // eventBuilder->SetMaximumTimeGap(100.);
+ eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::FixedTimeWindow);
+ eventBuilder->SetFixedTimeWindow(200.);
+ eventBuilder->SetTriggerMinNumberT0(1);
+ eventBuilder->SetTriggerMinNumberSts(0);
+ eventBuilder->SetTriggerMinNumberMuch(0);
+ eventBuilder->SetTriggerMinNumberTof(10);
+ eventBuilder->SetTriggerMinNumberRich(10);
+
+ run->AddTask(eventBuilder);
+
+ CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer();
+ hitProd->SetMappingFile("mRICH_Mapping_vert_20190318_elView.geo");
+ hitProd->setToTLimits(23.7,30.0);
+ hitProd->applyToTCut();
+ //hitProd->DoRestrictToAerogelAccDec2019();
+ run->AddTask(hitProd);
+
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ richReco->UseMCbmSetup();
+ run->AddTask(richReco);
+
+
+//######################################################################
+
+ CbmRichMCbmAerogelAna* qaTask = new CbmRichMCbmAerogelAna;
+ qaTask->SetOutputDir(Form("result_run%d_Aerogel",runId));
+ //qaTask->DoRestrictToAcc();//restrict to mRICH MAR2019 in histFilling
+ qaTask->XOffsetHistos(-2.5);
+ run->AddTask(qaTask);
+
+
+ std::cout << std::endl << std::endl << "-I- " << myName << ": Set runtime DB" << std::endl;
+
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ Bool_t kParameterMerged = kTRUE;
+ FairParRootFileIo* parIo1 = new FairParRootFileIo(kParameterMerged);
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo1->open(parFile.c_str(),"UPDATE");
+ parIo1->print();
+
+
+ parIo2->open(parFileList, "in");
+ parIo2->print();
+ rtdb->setFirstInput(parIo2);
+ rtdb->setSecondInput(parIo1);
+
+
+ rtdb->print();
+ rtdb->printParamContexts();
+
+ std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+
+ // rtdb->setOutput(parIo1);
+ rtdb->saveOutput();
+
+
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+
+
+ timer.Stop();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished succesfully." << std::endl;
+ std::cout << "Output file is " << recoFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl;
+ std::cout << "Test passed" << std::endl << "All ok" << std::endl;
+
+}
diff --git a/macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_dec19.C b/macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_dec19.C
new file mode 100644
index 00000000..a1e486f4
--- /dev/null
+++ b/macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_dec19.C
@@ -0,0 +1,582 @@
+void run_reco_mcbm_real_wToF_dec19(
+// const string& parFile = "/lustre/nyx/cbm/users/adrian/data19Dec12/testNew2/unp_mcbm_params_384.root",
+// const string& digiFile = "/lustre/nyx/cbm/users/adrian/data19Dec12/testNew2/unp_mcbm_384.root",
+ const string& parFile = "/lustre/nyx/cbm/users/adrian/cbmgitnew/cbmsource/macro/beamtime/mcbm2019/data/unp_mcbm_params_384.root",
+ const string& digiFile = "/lustre/nyx/cbm/users/adrian/cbmgitnew/cbmsource/macro/beamtime/mcbm2019/data/unp_mcbm_384.root",
+ const string& recoFile = "reco_mcbm_dec19.root",
+ const unsigned int runId = 384, // used for the output folder
+ int nEvents = 1000,
+ const int taskId = 00005
+)
+{
+ const Double_t eb_fixedTimeWindow {200.};
+ const Int_t eb_TriggerMinNumberT0 {1};
+ const Int_t eb_TriggerMinNumberSts {0};
+ const Int_t eb_TriggerMinNumberMuch{0};
+ const Int_t eb_TriggerMinNumberTof {10};
+ const Int_t eb_TriggerMinNumberRich{10};
+
+ //ToF Settings
+ Int_t calMode = 93;
+ Int_t calSel = 1;
+ Int_t calSm = 40;
+ Int_t RefSel = 1;
+ //TString cFileId = "159.200.10.0.0";
+ //TString cFileId = "385.100.-192.0"; // -> 900041901
+ TString cFileId = "385.100.5.0";
+ //Int_t iCalSet = 900041901;
+ Int_t iCalSet = 14500;//900041901;//10020500;//30040500;
+ Int_t iTrackingSetup=1;
+ Int_t iSel2 = 0;
+ Double_t dDeadtime = 50;
+ // 2020/03/23
+ // 385.100.5.0_set000014500_93_1tofClust.hst.root
+ // 385.100.5.0_tofFindTracks.hst.root
+
+
+ void save_hst(TString cstr, Bool_t bROOT);
+
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ TTree::SetMaxTreeSize(90000000000);
+
+ TString myName = "run_reco_mcbm_real";
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ TString workDir = gSystem->Getenv("VMCWORKDIR");
+
+ remove(recoFile.c_str());
+
+ //TString TofFileFolder = Form("/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/%s",cFileId.Data());
+ TString TofFileFolder = Form("/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2019/%s",cFileId.Data());
+// TString setupFile = srcDir + "/geometry/setup/setup_" + geoSetup + ".C";
+// TString setupFunct = "setup_" + geoSetup + "()";
+// gROOT->LoadMacro(setupFile);
+// gROOT->ProcessLine(setupFunct);
+
+ std::cout << std::endl<< "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList *parFileList = new TList();
+
+ TStopwatch timer;
+ timer.Start();
+ gDebug = 0;
+
+ //-----------------------------------------------//
+ TString FId = cFileId;
+ TString TofGeo = "v19b_mcbm"; //v18m_mCbm
+
+ TObjString *tofDigiFile = new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par"); // TOF digi file
+ std::cout << std::endl<< "-I- digi.par file : "<<workDir << "/parameters/tof/tof_" << TofGeo << ".digi.par"<< std::endl;
+ parFileList->Add(tofDigiFile);
+
+ TObjString *tofDigiBdfFile = new TObjString( workDir + "/parameters/tof/tof_" + TofGeo +".digibdf.par");
+ std::cout << std::endl<< "-I- digibdf.par file : "<<workDir << "/parameters/tof/tof_" << TofGeo << ".digibdf.par"<< std::endl;
+ parFileList->Add(tofDigiBdfFile);
+
+ TString geoDir = workDir; // gSystem->Getenv("VMCWORKDIR");
+ //TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+ //TString geoFile = "/lustre/nyx/cbm/users/adrian/data/sim/rich/mcbm/sis18_mcbm_25deg_long_geofile_full.root"; //18m
+ TString geoFile = srcDir + "/macro/mcbm/data/mcbm_beam_2019_12.geo.root";
+ TFile* fgeo = new TFile(geoFile);
+ TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+ if (NULL == geoMan){
+ cout << "<E> FAIRGeom not found in geoFile"<<endl;
+ return;
+ }
+
+ //-----------------------------------------------//
+
+ gROOT->LoadMacro("save_hst.C");
+
+ FairRunAna *run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(digiFile.c_str());
+ run->SetSource(inputSource);
+ run->SetOutputFile(recoFile.c_str());
+ run->SetGenerateRunInfo(kTRUE);
+
+ /// Add the Eventbuilder
+
+ CbmMcbm2018EventBuilder* eventBuilder = new CbmMcbm2018EventBuilder();
+ // eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::MaximumTimeGap);
+ // eventBuilder->SetMaximumTimeGap(100.);
+ eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::FixedTimeWindow);
+ eventBuilder->SetFixedTimeWindow( eb_fixedTimeWindow );
+ eventBuilder->SetTriggerMinNumberT0( eb_TriggerMinNumberT0 );
+ eventBuilder->SetTriggerMinNumberSts( eb_TriggerMinNumberSts );
+ eventBuilder->SetTriggerMinNumberMuch( eb_TriggerMinNumberMuch);
+ eventBuilder->SetTriggerMinNumberTof( eb_TriggerMinNumberTof );
+ eventBuilder->SetTriggerMinNumberRich( eb_TriggerMinNumberRich);
+
+ run->AddTask(eventBuilder);
+
+ CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer();
+ hitProd->SetMappingFile("mRICH_Mapping_vert_20190318_elView.geo");
+ hitProd->setToTLimits(23.7,30.0);
+ hitProd->applyToTCut();
+ run->AddTask(hitProd);
+
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ richReco->UseMCbmSetup();
+ run->AddTask(richReco);
+
+ // Starting with ToF
+
+ CbmTofEventClusterizer* tofClust = new CbmTofEventClusterizer("TOF Event Clusterizer",1,0);
+ tofClust->SetCalMode(calMode);
+ tofClust->SetCalSel(calSel);
+ tofClust->SetCaldXdYMax(10.); // geometrical matching window in cm
+ tofClust->SetCalCluMulMax(3.); // Max Counter Cluster Multiplicity for filling calib histos
+ tofClust->SetCalRpc(calSm); // select detector for calibration update
+ tofClust->SetTRefId(RefSel); // reference trigger for offset calculation
+ tofClust->SetTotMax(20.); // Tot upper limit for walk corection
+ tofClust->SetTotMin(0.01); //(12000.); // Tot lower limit for walk correction
+ tofClust->SetTotPreRange(5.); // effective lower Tot limit in ns from peak position
+ tofClust->SetTotMean(5.); // Tot calibration target value in ns
+ tofClust->SetMaxTimeDist(1.0); // default cluster range in ns
+ tofClust->SetDelTofMax(10.); // acceptance range for cluster distance in ns (!)
+ tofClust->SetSel2MulMax(3); // limit Multiplicity in 2nd selector
+ tofClust->SetChannelDeadtime(dDeadtime); // artificial deadtime in ns
+ tofClust->SetEnableAvWalk(kFALSE);
+ //tofClust->SetEnableMatchPosScaling(kFALSE); // turn off projection to nominal target
+ tofClust->SetYFitMin(1.E4);
+ tofClust->SetIdMode(1); // calibrate on module level
+ // tofClust->SetDeadStrips(15,23); // declare dead strip for T0M3,Rpc0,Strip 23
+
+
+ Int_t calSelRead = calSel;
+ if (calSel<0) calSelRead=0;
+ TString cCalibFname=Form("/%s_set%09d_%02d_%01dtofClust.hst.root",cFileId.Data(),iCalSet,calMode,calSelRead);
+ //cCalibFname = "/385.100.5.0_set000014500_93_1tofClust.hst.root";
+ tofClust->SetCalParFileName(TofFileFolder + cCalibFname);
+ TString cOutFname=Form("tofClust_%s_set%09d.hst.root",cFileId.Data(),iCalSet);
+ tofClust->SetOutHstFileName(cOutFname);
+
+switch (calMode) {
+ case -1: // initial check of raw data
+ tofClust->SetTotMax(256.); // range in bin number
+ tofClust->SetTotPreRange(256.);
+ //tofClust->SetTotMin(1.);
+ tofClust->SetTRefDifMax(26000.); // in ns
+ tofClust->PosYMaxScal(10000.); // in % of length
+ tofClust->SetMaxTimeDist(0.); // no cluster building
+ //tofClust->SetTimePeriod(25600.); // inspect coarse time
+ break;
+ case 0: // initial calibration
+ tofClust->SetTotMax(256.); // range in bin number
+ tofClust->SetTotPreRange(256.);
+ //tofClust->SetTotMin(1.);
+ tofClust->SetTRefDifMax(1000.); // in ns
+ tofClust->PosYMaxScal(10.); // in % of length
+ tofClust->SetMaxTimeDist(0.); // no cluster building
+ break;
+ case 1: // save offsets, update walks, for diamonds
+ tofClust->SetTotMax(256.); // range in bin number
+ tofClust->SetTotPreRange(256.);
+ tofClust->SetTRefDifMax(6.25); // in ns
+ //tofClust->SetTimePeriod(6.25); // inspect coarse time
+ tofClust->PosYMaxScal(10.); // in % of length
+ break;
+ case 11:
+ tofClust->SetTRefDifMax(5.); // in ns
+ tofClust->PosYMaxScal(3.0); // in % of length
+ break;
+ case 21:
+ tofClust->SetTRefDifMax(2.5); // in ns
+ tofClust->PosYMaxScal(2.0); // in % of length
+ break;
+ case 31:
+ tofClust->SetTRefDifMax(2.); // in ns
+ tofClust->PosYMaxScal(1.5); // in % of length
+ break;
+ case 41:
+ tofClust->SetTRefDifMax(1.5); // in ns
+ tofClust->PosYMaxScal(0.8); // in % of length
+ break;
+ case 51:
+ tofClust->SetTRefDifMax(1.2); // in ns
+ tofClust->PosYMaxScal(0.7); // in % of length
+ break;
+ case 61:
+ tofClust->SetTRefDifMax(1.0); // in ns
+ tofClust->PosYMaxScal(0.75); // in % of length
+ break;
+ case 71:
+ tofClust->SetTRefDifMax(0.8); // in ns
+ tofClust->PosYMaxScal(0.6); // in % of length
+ break;
+
+ case 2: // time difference calibration
+ tofClust->SetTRefDifMax(300.); // in ns
+ tofClust->PosYMaxScal(1000.); //in % of length
+ break;
+
+ case 3: // time offsets
+ tofClust->SetTRefDifMax(200.); // in ns
+ tofClust->PosYMaxScal(100.); //in % of length
+ tofClust->SetMaxTimeDist(0.); // no cluster building
+ break;
+ case 12:
+ case 13:
+ tofClust->SetTRefDifMax(100.); // in ns
+ tofClust->PosYMaxScal(10.); //in % of length
+ break;
+ case 22:
+ case 23:
+ tofClust->SetTRefDifMax(50.); // in ns
+ tofClust->PosYMaxScal(5.); //in % of length
+ break;
+ case 32:
+ case 33:
+ tofClust->SetTRefDifMax(25.); // in ns
+ tofClust->PosYMaxScal(4.); //in % of length
+ break;
+ case 42:
+ case 43:
+ tofClust->SetTRefDifMax(12.); // in ns
+ tofClust->PosYMaxScal(2.); //in % of length
+ break;
+ case 52:
+ case 53:
+ tofClust->SetTRefDifMax(5.); // in ns
+ tofClust->PosYMaxScal(1.5); //in % of length
+ break;
+ case 62:
+ case 63:
+ tofClust->SetTRefDifMax(3.); // in ns
+ tofClust->PosYMaxScal(1.); //in % of length
+ break;
+ case 72:
+ case 73:
+ tofClust->SetTRefDifMax(2.5); // in ns
+ tofClust->PosYMaxScal(0.9); //in % of length
+ break;
+ case 82:
+ case 83:
+ tofClust->SetTRefDifMax(2.0); // in ns
+ tofClust->PosYMaxScal(0.8); //in % of length
+ break;
+ case 92:
+ case 93:
+ tofClust->SetTRefDifMax(1.5); // in ns
+ tofClust->PosYMaxScal(0.75); //in % of length
+ break;
+
+ case 4: // velocity dependence (DelTOF)
+ case 14:
+ tofClust->SetTRefDifMax(25.); // in ns
+ tofClust->PosYMaxScal(2.0); //in % of length
+ break;
+ case 24:
+ tofClust->SetTRefDifMax(5.); // in ns
+ tofClust->PosYMaxScal(1.5); //in % of length
+ break;
+ case 34:
+ tofClust->SetTRefDifMax(5.); // in ns
+ tofClust->PosYMaxScal(1.2); //in % of length
+ break;
+ case 44:
+ tofClust->SetTRefDifMax(3.); // in ns
+ tofClust->PosYMaxScal(1.0); //in % of length
+ break;
+ case 54:
+ tofClust->SetTRefDifMax(2.5); // in ns
+ tofClust->PosYMaxScal(0.9); //in % of length
+ break;
+ case 64:
+ tofClust->SetTRefDifMax(2.5); // in ns
+ tofClust->PosYMaxScal(0.8); //in % of length
+ break;
+ case 74:
+ tofClust->SetTRefDifMax(1.5); // in ns
+ tofClust->PosYMaxScal(0.7); //in % of length
+ break;
+ default:
+ cout << "<E> Calib mode not implemented! stop execution of script"<<endl;
+ return;
+ }
+
+ run->AddTask(tofClust);
+
+ Int_t iBRef=iCalSet%1000;
+ Int_t iSet = (iCalSet - iBRef)/1000;
+ Int_t iRSel=0;
+ Int_t iRSelTyp=0;
+ Int_t iRSelSm=0;
+ Int_t iRSelRpc=0;
+ iRSel=iBRef; // use diamond
+ if(iSel2==0){
+ // iSel2=iBRef;
+ }else{
+ if(iSel2<0) iSel2=-iSel2;
+ }
+
+ Int_t iRSelin=iRSel;
+ iRSelRpc=iRSel%10;
+ iRSelTyp = (iRSel-iRSelRpc)/10;
+ iRSelSm=iRSelTyp%10;
+ iRSelTyp = (iRSelTyp-iRSelSm)/10;
+
+ tofClust->SetBeamRefId(iRSelTyp); // define Beam reference counter
+ tofClust->SetBeamRefSm(iRSelSm);
+ tofClust->SetBeamRefDet(iRSelRpc);
+ tofClust->SetBeamAddRefMul(-1);
+ tofClust->SetBeamRefMulMax(3);
+
+ Int_t iSel2in=iSel2;
+ Int_t iSel2Rpc= iSel2%10;
+ iSel2=(iSel2-iSel2Rpc)/10;
+ Int_t iSel2Sm=iSel2%10;
+ iSel2=(iSel2-iSel2Sm)/10;
+ if(iSel2 > 0) {
+ tofClust->SetSel2Id(iSel2);
+ tofClust->SetSel2Sm(iSel2Sm);
+ tofClust->SetSel2Rpc(iSel2Rpc);
+ }
+
+ Int_t iRef = iSet %1000;
+ Int_t iDut = (iSet - iRef)/1000;
+ Int_t iDutRpc = iDut%10;
+ iDut = (iDut - iDutRpc)/10;
+ Int_t iDutSm = iDut%10;
+ iDut = (iDut - iDutSm)/10;
+
+ tofClust->SetDutId(iDut);
+ tofClust->SetDutSm(iDutSm);
+ tofClust->SetDutRpc(iDutRpc);
+
+ Int_t iRefRpc = iRef%10;
+ iRef = (iRef - iRefRpc)/10;
+ Int_t iRefSm = iRef%10;
+ iRef = (iRef - iRefSm)/10;
+
+ tofClust->SetSelId(iRef);
+ tofClust->SetSelSm(iRefSm);
+ tofClust->SetSelRpc(iRefRpc);
+
+ cout << "Run mTof Clusterizer with iRSel = "<<iRSel<<", iSel2 = "<<iSel2in<<endl;
+
+//######################################################################
+
+ // mTof Tracker Initialization
+ TString cTrkFile=Form("/%s_tofFindTracks.hst.root",cFileId.Data());
+ CbmTofTrackFinder* tofTrackFinder= new CbmTofTrackFinderNN();
+ tofTrackFinder->SetMaxTofTimeDifference(0.2); // in ns/cm
+ tofTrackFinder->SetTxLIM(0.3); // max slope dx/dz
+ tofTrackFinder->SetTyLIM(0.3); // max dev from mean slope dy/dz
+ tofTrackFinder->SetTyMean(0.); // mean slope dy/dz
+
+ CbmTofTrackFitter* tofTrackFitter= new CbmTofTrackFitterKF(0,211);
+ TFitter *MyFit = new TFitter(1); // initialize Minuit
+ tofTrackFinder->SetFitter(tofTrackFitter);
+ CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+ tofFindTracks->UseFinder(tofTrackFinder);
+ tofFindTracks->UseFitter(tofTrackFitter);
+ Int_t iGenCor=1;
+ tofFindTracks->SetCorMode(iGenCor); // valid options: 0,1,2,3,4,5,6, 10 - 19
+ tofFindTracks->SetTtTarg(0.057); // target value for inverse velocity, > 0.033 ns/cm!
+ //tofFindTracks->SetTtTarg(0.135); // target value for inverse velocity, > 0.033 ns/cm!
+ tofFindTracks->SetCalParFileName(TofFileFolder + cTrkFile); // Tracker parameter value file name
+ std::cout<<"TrackCalParFile: "<<TofFileFolder << cTrkFile<<std::endl;
+ tofFindTracks->SetBeamCounter(5,0,0); // default beam counter
+ tofFindTracks->SetStationMaxHMul(30); // Max Hit Multiplicity in any used station
+
+ tofFindTracks->SetT0MAX(1.); // in ns
+ tofFindTracks->SetSIGT(0.08); // default in ns
+ tofFindTracks->SetSIGX(0.3); // default in cm
+ tofFindTracks->SetSIGY(0.6); // default in cm
+ tofFindTracks->SetSIGZ(0.05); // default in cm
+ tofFindTracks->SetUseSigCalib(kFALSE); // ignore resolutions in CalPar file
+ Double_t dChi2Lim2= 3.5;//100;//3.5;
+ tofTrackFinder->SetSIGLIM(dChi2Lim2*2.); // matching window in multiples of chi2
+ tofTrackFinder->SetChiMaxAccept(dChi2Lim2); // max tracklet chi2
+
+
+ Int_t iMinNofHits=-1;
+ Int_t iNStations=0;
+ Int_t iNReqStations=3;
+
+ switch (iTrackingSetup){
+ case 0: // bypass mode
+ iMinNofHits=-1;
+ iNStations=1;
+ tofFindTracks->SetStation(0, 5, 0, 0); // Diamond
+ break;
+
+ case 1: // for calibration mode of full setup
+ iMinNofHits=3;
+ iNStations=39;
+ iNReqStations=4;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 0, 4, 0);
+ tofFindTracks->SetStation(2, 0, 3, 0);
+ tofFindTracks->SetStation(3, 0, 4, 1);
+ tofFindTracks->SetStation(4, 0, 3, 1);
+ tofFindTracks->SetStation(5, 0, 4, 2);
+ tofFindTracks->SetStation(6, 0, 3, 2);
+ tofFindTracks->SetStation(7, 0, 4, 3);
+ tofFindTracks->SetStation(8, 0, 3, 3);
+ tofFindTracks->SetStation(9, 0, 4, 4);
+ tofFindTracks->SetStation(10, 0, 3, 4);
+ tofFindTracks->SetStation(11, 9, 0, 0);
+ tofFindTracks->SetStation(12, 9, 0, 1);
+ tofFindTracks->SetStation(13, 7, 0, 0);
+ tofFindTracks->SetStation(14, 6, 0, 0);
+ tofFindTracks->SetStation(15, 6, 0, 1);
+ tofFindTracks->SetStation(16, 8, 0, 0);
+ tofFindTracks->SetStation(17, 8, 0, 1);
+ tofFindTracks->SetStation(18, 8, 0, 2);
+ tofFindTracks->SetStation(19, 8, 0, 3);
+ tofFindTracks->SetStation(20, 8, 0, 4);
+ tofFindTracks->SetStation(21, 8, 0, 5);
+ tofFindTracks->SetStation(22, 8, 0, 6);
+ tofFindTracks->SetStation(23, 8, 0, 7);
+
+ tofFindTracks->SetStation(24, 0, 2, 2);
+ tofFindTracks->SetStation(25, 0, 1, 2);
+ tofFindTracks->SetStation(26, 0, 0, 2);
+ tofFindTracks->SetStation(27, 0, 2, 1);
+ tofFindTracks->SetStation(28, 0, 1, 1);
+ tofFindTracks->SetStation(29, 0, 0, 1);
+ tofFindTracks->SetStation(30, 0, 2, 3);
+ tofFindTracks->SetStation(31, 0, 1, 3);
+ tofFindTracks->SetStation(32, 0, 0, 3);
+ tofFindTracks->SetStation(33, 0, 2, 0);
+ tofFindTracks->SetStation(34, 0, 1, 0);
+ tofFindTracks->SetStation(35, 0, 0, 0);
+ tofFindTracks->SetStation(36, 0, 2, 4);
+ tofFindTracks->SetStation(37, 0, 1, 4);
+ tofFindTracks->SetStation(38, 0, 0, 4);
+ break;
+
+ case 2:
+ iMinNofHits=3;
+ iNStations=14;
+ iNReqStations=5;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 0, 4, 1);
+ tofFindTracks->SetStation(2, 0, 3, 1);
+ tofFindTracks->SetStation(3, 0, 4, 0);
+ tofFindTracks->SetStation(4, 0, 3, 0);
+ tofFindTracks->SetStation(5, 0, 4, 2);
+ tofFindTracks->SetStation(6, 0, 3, 2);
+ tofFindTracks->SetStation(7, 0, 4, 3);
+ tofFindTracks->SetStation(8, 0, 3, 3);
+ tofFindTracks->SetStation(9, 0, 4, 4);
+ tofFindTracks->SetStation(10, 0, 3, 4);
+ tofFindTracks->SetStation(11, 9, 0, 0);
+ tofFindTracks->SetStation(12, 9, 0, 1);
+ tofFindTracks->SetStation(13, 7, 0, 0);
+ break;
+
+ default:
+ ;
+ }
+
+ tofFindTracks->SetMinNofHits(iMinNofHits);
+ tofFindTracks->SetNStations(iNStations);
+ tofFindTracks->SetNReqStations(iNReqStations);
+ tofFindTracks->PrintSetup();
+ run->AddTask(tofFindTracks);
+
+//######################################################################
+
+ CbmRichMCbmQaReal* qaTask = new CbmRichMCbmQaReal();
+ if (taskId < 0) {
+ qaTask->SetOutputDir(Form("result_run%d",runId));
+ } else {
+ qaTask->SetOutputDir(Form("result_run%d_%05d",runId,taskId));
+ }
+// qaTask->DoRestrictToAcc();//restrict to mRICH MAR2019 in histFilling
+ qaTask->XOffsetHistos(12);
+ run->AddTask(qaTask);
+
+
+ std::cout << std::endl << std::endl << "-I- " << myName << ": Set runtime DB" << std::endl;
+
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ Bool_t kParameterMerged = kTRUE;
+ FairParRootFileIo* parIo1 = new FairParRootFileIo(kParameterMerged);
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo1->open(parFile.c_str(),"UPDATE");
+ parIo1->print();
+
+
+ parIo2->open(parFileList, "in");
+ parIo2->print();
+ rtdb->setFirstInput(parIo2);
+ rtdb->setSecondInput(parIo1);
+
+
+ rtdb->print();
+ rtdb->printParamContexts();
+
+ std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+
+ // rtdb->setOutput(parIo1);
+ rtdb->saveOutput();
+
+
+ // print all important infos in a file
+ std::ofstream outfile;
+ if (taskId < 0) {
+ outfile.open(Form("result_run%d/run_info.dat",runId));
+ } else {
+ outfile.open(Form("result_run%d_%05d/run_info.dat",runId,taskId));
+ }
+ // write inputted data into the file.
+ outfile << "Run: " << runId << std::endl;
+ outfile << "Events: " << nEvents << std::endl;
+ outfile << "parFile: " << parFile << std::endl;
+ outfile << "digiFile: " << digiFile << std::endl;
+ outfile << "recoFile: " << recoFile << std::endl;
+ outfile << "Geometry: " << geoFile << std::endl;
+ outfile << "TrackCalParFile: " << TofFileFolder << cTrkFile << std::endl;
+ outfile << "TofClusterFile :" << TofFileFolder + cCalibFname << std::endl;
+ outfile << "TofOutput :" << cOutFname << std::endl << std::endl;
+ outfile << "Trigger:"<< std::endl;
+ outfile << " fixedTimeWindow :" << eb_fixedTimeWindow << std::endl;
+ outfile << " MinNumberT0 :" << eb_TriggerMinNumberT0 << std::endl;
+ outfile << " MinNumberSts :" << eb_TriggerMinNumberSts << std::endl;
+ outfile << " MinNumberMuch :" << eb_TriggerMinNumberMuch << std::endl;
+ outfile << " MinNumberTof :" << eb_TriggerMinNumberTof << std::endl;
+ outfile << " MinNumberRich :" << eb_TriggerMinNumberRich << std::endl;
+ outfile.close();
+
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+
+
+ timer.Stop();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished succesfully." << std::endl;
+ std::cout << "Output file is " << recoFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl;
+ std::cout << "Test passed" << std::endl << "All ok" << std::endl;
+
+}
+
+void save_hst(TString cstr="status.hst.root", Bool_t bROOT = kFALSE){
+ cout << "save all histograms to file "<<cstr.Data()<<endl;
+ TList* tList(NULL);
+ if(bROOT){
+ tList = gROOT->GetList();
+ } else {
+ tList = gDirectory->GetList();
+ }
+ TIter next(tList);
+ // Write objects to the file
+ TFile *fHist = new TFile(cstr,"RECREATE");
+ {
+ TObject* obj;
+ while( (obj= (TObject*)next()) ){
+ if(obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(TEfficiency::Class())){
+ obj->Write();
+ }
+ }
+ }
+ // fHist->ls();
+ fHist->Close();
+}
diff --git a/macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_mar20.C b/macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_mar20.C
new file mode 100644
index 00000000..b85952c5
--- /dev/null
+++ b/macro/rich/mcbm/beamtime/run_reco_mcbm_real_wToF_mar20.C
@@ -0,0 +1,541 @@
+void run_reco_mcbm_real_wToF_mar20(
+ const string srcfolder = "/lustre/cbm/users/adrian/mcbmbeamtime/cbmsource/macro/beamtime/mcbm2020/data",
+ const unsigned int runId = 831, // used for the output folder
+ int nEvents = 5000,
+ const int taskId = 3
+)
+{
+ // ----- File names --------------------------------------------------
+ const string& parFile = Form("%s/unp_mcbm_params_%d.root",srcfolder.c_str(),runId);
+ const string& digiFile = Form("%s/unp_mcbm_%d.root",srcfolder.c_str(),runId);
+ const string& recoFile = Form("reco_mcbm_mar20_%d.root",runId);
+ TString setup = "mcbm_beam_2020_03";
+ // -----------------------------------------------------------------------
+
+
+ // ----- EventBuilder Settings-----------------------------------------
+ const Double_t eb_fixedTimeWindow {200.};
+ const Int_t eb_TriggerMinNumberT0 {1};
+ const Int_t eb_TriggerMinNumberSts {0};
+ const Int_t eb_TriggerMinNumberMuch{0};
+ const Int_t eb_TriggerMinNumberTof {10};
+ const Int_t eb_TriggerMinNumberRich{5};
+ // -----------------------------------------------------------------------
+
+
+ // ----- TOF defaults --------------------------------------------------
+ Int_t calMode = 93;
+ Int_t calSel = 1;
+ Int_t calSm = 0;
+ Int_t RefSel = 0;
+ Double_t dDeadtime = 50.;
+ Int_t iSel2 = 500;
+ // -----------------------------------------------------------------------
+
+
+ // ----- TOF Settings --------------------------------------------------
+ TString cCalId = "490.100.5.0";
+ if (runId >= 759) cCalId = "759.100.4.0";
+ if (runId >= 812) cCalId = "831.100.4.0";
+
+ Int_t iCalSet = 30040500; // calibration settings
+ if (runId >= 759) iCalSet = 10020500;
+ if (runId >= 812) iCalSet = 10020500;
+
+ Double_t Tint =100.; // coincidence time interval
+ Int_t iTrackMode=2; // 2 for TofTracker
+ const Int_t iTofCluMode=1;
+ // -----------------------------------------------------------------------
+
+
+ // ----- Fair logger ---------------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ TTree::SetMaxTreeSize(90000000000);
+ // -----------------------------------------------------------------------
+
+
+ TString myName = "run_reco_mcbm_real";
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ TString workDir = gSystem->Getenv("VMCWORKDIR");
+
+ remove(recoFile.c_str());
+
+ // ----- Load the geometry setup -------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl;
+ CbmSetup *pSetup=CbmSetup::Instance();
+ pSetup->LoadSetup(setup);
+ // You can modify the pre-defined setup by using
+ // CbmSetup::Instance()->RemoveModule(ESystemId) or
+ // CbmSetup::Instance()->SetModule(ESystemId, const char*, Bool_t) or
+ //CbmSetup::Instance()->SetActive(ESystemId, Bool_t)
+ pSetup->SetActive(ECbmModuleId::kMvd, kFALSE);
+ pSetup->SetActive(ECbmModuleId::kSts, kFALSE);
+ pSetup->SetActive(ECbmModuleId::kMuch, kFALSE);
+ pSetup->SetActive(ECbmModuleId::kRich, kTRUE);
+ pSetup->SetActive(ECbmModuleId::kTrd, kFALSE);
+ pSetup->SetActive(ECbmModuleId::kPsd, kFALSE);
+ // -----------------------------------------------------------------------
+
+
+
+ //TString TofFileFolder = Form("/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2018/%s",cFileId.Data());
+// TString TofFileFolder = Form("/lustre/nyx/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2019/%s",cFileId.Data());
+ TString TofFileFolder = Form("/lustre/cbm/users/nh/CBM/cbmroot/trunk/macro/beamtime/mcbm2020/%s",cCalId.Data());
+// TString setupFile = srcDir + "/geometry/setup/setup_" + geoSetup + ".C";
+// TString setupFunct = "setup_" + geoSetup + "()";
+// gROOT->LoadMacro(setupFile);
+// gROOT->ProcessLine(setupFunct);
+
+ std::cout << std::endl<< "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList *parFileList = new TList();
+
+
+
+ //-----------------------------------------------//
+ // TString FId = cCalId;
+ // TString TofGeo = "v19b_mcbm"; //v18m_mCbm
+ // TString TofGeo = "v20a_mcbm"; //v18m_mCbm
+
+
+ // ----- TOF digitisation parameters -------------------------------------
+ TString geoTag;
+ TString geoFile;
+ if ( pSetup->IsActive(ECbmModuleId::kTof) ) {
+ pSetup->GetGeoTag(ECbmModuleId::kTof, geoTag);
+ TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
+ parFileList->Add(tofBdfFile);
+ std::cout << "-I- " << myName << ": Using parameter file "
+ << tofBdfFile->GetString() << std::endl;
+
+ //TString geoFile = srcDir + "/geometry/tof/geofile_tof_" + geoTag + ".root";
+ geoFile = srcDir + "/macro/mcbm/data/mcbm_beam_2020_03.geo.root";
+ TFile* fgeo = new TFile(geoFile);
+ TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+ if (NULL == geoMan){
+ cout << "<E> FAIRGeom not found in geoFile "<< geoFile.Data() <<endl;
+ return;
+ }
+ }
+
+ // ------------------------------------------------------------------------
+
+// TObjString *tofDigiFile = new TObjString(workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par"); // TOF digi file
+// std::cout << std::endl<< "-I- digi.par file : "<<workDir << "/parameters/tof/tof_" << TofGeo << ".digi.par"<< std::endl;
+// parFileList->Add(tofDigiFile);
+//
+// TObjString *tofDigiBdfFile = new TObjString( workDir + "/parameters/tof/tof_" + TofGeo +".digibdf.par");
+// std::cout << std::endl<< "-I- digibdf.par file : "<<workDir << "/parameters/tof/tof_" << TofGeo << ".digibdf.par"<< std::endl;
+// parFileList->Add(tofDigiBdfFile);
+
+// TString geoDir = workDir; // gSystem->Getenv("VMCWORKDIR");
+// //TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
+// //TString geoFile = "/lustre/nyx/cbm/users/adrian/data/sim/rich/mcbm/sis18_mcbm_25deg_long_geofile_full.root"; //18m
+// // TString geoFile = srcDir + "/macro/mcbm/data/mcbm_beam_2019_11.geo.root";
+// TString geoFile = srcDir + "/macro/mcbm/data/mcbm_beam_2020_03.geo.root";
+// TFile* fgeo = new TFile(geoFile);
+// TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+// if (NULL == geoMan){
+// cout << "<E> FAIRGeom not found in geoFile"<<endl;
+// return;
+// }
+
+
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+
+
+ // ---- Debug option -------------------------------------------------
+ gDebug = 0;
+ // ------------------------------------------------------------------------
+
+ gROOT->LoadMacro("save_hst.C");
+
+
+ // ----- Input file ---------------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Using input file " << digiFile << std::endl;
+ // ------------------------------------------------------------------------
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairRunAna *run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(digiFile.c_str());
+ run->SetSource(inputSource);
+ run->SetOutputFile(recoFile.c_str());
+// run->SetGenerateRunInfo(kFALSE);
+ //if (hasFairMonitor) FairMonitor::GetMonitor()->EnableMonitor(kTRUE);
+
+
+ // ----- Cbm EventBuilder ---------------------------------------------------
+ CbmMcbm2018EventBuilder* eventBuilder = new CbmMcbm2018EventBuilder();
+// eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::MaximumTimeGap);
+ eventBuilder->SetMaximumTimeGap(200.);
+ eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::FixedTimeWindow);
+ eventBuilder->SetFixedTimeWindow( eb_fixedTimeWindow );
+ eventBuilder->SetTriggerMinNumberT0( eb_TriggerMinNumberT0 );
+ eventBuilder->SetTriggerMinNumberSts( eb_TriggerMinNumberSts );
+ eventBuilder->SetTriggerMinNumberMuch( eb_TriggerMinNumberMuch);
+ eventBuilder->SetTriggerMinNumberTof( eb_TriggerMinNumberTof );
+ eventBuilder->SetTriggerMinNumberRich( eb_TriggerMinNumberRich);
+ eventBuilder->SetFillHistos(kFALSE); // to prevent memory leak???
+
+ run->AddTask(eventBuilder);
+
+
+
+ // ----- Local reconstruction of RICH Hits ---------------
+ CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer();
+ hitProd->SetMappingFile("mRICH_Mapping_vert_20190318_elView.geo");
+ hitProd->setToTLimits(23.7,30.0);
+ hitProd->applyToTCut();
+ run->AddTask(hitProd);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in RICh -> Finding of Rings ---------------
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ richReco->UseMCbmSetup();
+ run->AddTask(richReco);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in TOF ----------------------------------
+ TString cFname;
+ switch(iTofCluMode) {
+ case 1:
+ {
+ CbmTofEventClusterizer* tofCluster = new CbmTofEventClusterizer("TOF Event Clusterizer",0,1);
+// cFname=Form("/%s_set%09d_%02d_%01dtofClust.hst.root",cCalId.Data(),iCalSet,calMode,calSel);
+ cFname=Form("/%s_set%09d_%02d_%01d_noWalk_tofClust.hst.root",cCalId.Data(),iCalSet,calMode,calSel);
+ tofCluster->SetCalParFileName(TofFileFolder +cFname);
+ tofCluster->SetCalMode(calMode);
+ tofCluster->SetCalSel(calSel);
+ tofCluster->SetCaldXdYMax(3.); // geometrical matching window in cm
+ tofCluster->SetCalCluMulMax(5.); // Max Counter Cluster Multiplicity for filling calib histos
+ tofCluster->SetCalRpc(calSm); // select detector for calibration update
+ tofCluster->SetTRefId(RefSel); // reference trigger for offset calculation
+ tofCluster->SetTotMax(20.); // Tot upper limit for walk corection
+ tofCluster->SetTotMin(0.01); //(12000.); // Tot lower limit for walk correction
+ tofCluster->SetTotPreRange(5.); // effective lower Tot limit in ns from peak position
+ tofCluster->SetTotMean(5.); // Tot calibration target value in ns
+ tofCluster->SetMaxTimeDist(1.0); // default cluster range in ns
+ tofCluster->SetDelTofMax(15.); // acceptance range for cluster distance in ns (!)
+ tofCluster->SetSel2MulMax(3); // limit Multiplicity in 2nd selector
+ tofCluster->SetChannelDeadtime(dDeadtime); // artificial deadtime in ns
+ tofCluster->SetEnableAvWalk(kFALSE);
+ //tofCluster->SetEnableMatchPosScaling(kFALSE); // turn off projection to nominal target
+ tofCluster->SetYFitMin(1.E4);
+ tofCluster->SetToDAv(0.04);
+ tofCluster->SetIdMode(1); // calibrate on module level
+ tofCluster->SetTRefDifMax(2.0); // in ns
+ tofCluster->PosYMaxScal(0.75); //in % of length
+ Int_t iBRef=iCalSet%1000;
+ Int_t iSet = (iCalSet - iBRef)/1000;
+ Int_t iRSel=0;
+ Int_t iRSelTyp=0;
+ Int_t iRSelSm=0;
+ Int_t iRSelRpc=0;
+ iRSel=iBRef; // use diamond
+ Int_t iRSelin=iRSel;
+ iRSelRpc=iRSel%10;
+ iRSelTyp = (iRSel-iRSelRpc)/10;
+ iRSelSm=iRSelTyp%10;
+ iRSelTyp = (iRSelTyp-iRSelSm)/10;
+ tofCluster->SetBeamRefId(iRSelTyp); // define Beam reference counter
+ tofCluster->SetBeamRefSm(iRSelSm);
+ tofCluster->SetBeamRefDet(iRSelRpc);
+ tofCluster->SetBeamAddRefMul(-1);
+ tofCluster->SetBeamRefMulMax(3);
+ Int_t iSel2in=iSel2;
+ Int_t iSel2Rpc= iSel2%10;
+ iSel2=(iSel2-iSel2Rpc)/10;
+ Int_t iSel2Sm=iSel2%10;
+ iSel2=(iSel2-iSel2Sm)/10;
+ if(iSel2 > -1) {
+ tofCluster->SetSel2Id(iSel2);
+ tofCluster->SetSel2Sm(iSel2Sm);
+ tofCluster->SetSel2Rpc(iSel2Rpc);
+ }
+ Int_t iRef = iSet %1000;
+ Int_t iDut = (iSet - iRef)/1000;
+ Int_t iDutRpc = iDut%10;
+ iDut = (iDut - iDutRpc)/10;
+ Int_t iDutSm = iDut%10;
+ iDut = (iDut - iDutSm)/10;
+ tofCluster->SetDutId(iDut);
+ tofCluster->SetDutSm(iDutSm);
+ tofCluster->SetDutRpc(iDutRpc);
+
+ Int_t iRefRpc = iRef%10;
+ iRef = (iRef - iRefRpc)/10;
+ Int_t iRefSm = iRef%10;
+ iRef = (iRef - iRefSm)/10;
+
+ tofCluster->SetSelId(iRef);
+ tofCluster->SetSelSm(iRefSm);
+ tofCluster->SetSelRpc(iRefRpc);
+
+ run->AddTask(tofCluster);
+ std::cout << "-I- " << myName << ": Added task "
+ << tofCluster->GetName() << std::endl;
+ }
+ break;
+
+ default:
+ {
+ ;
+ }
+ }
+ // -------------------------------------------------------------------------
+
+
+ // ----- Track reconstruction ------------------------------------------
+ Double_t beamWidthX = 0.1;
+ Double_t beamWidthY = 0.1;
+ switch(iTrackMode) {
+ case 2:
+ {
+ Int_t iGenCor=1;
+ Double_t dScalFac=1.;
+ Double_t dChi2Lim2=3.5;
+ TString cTrkFile=Form("/%s_tofFindTracks.hst.root",cCalId.Data());
+ Int_t iTrackingSetup=1;
+ CbmTofTrackFinder* tofTrackFinder= new CbmTofTrackFinderNN();
+ tofTrackFinder->SetMaxTofTimeDifference(0.2); // in ns/cm
+ tofTrackFinder->SetTxLIM(0.3); // max slope dx/dz
+ tofTrackFinder->SetTyLIM(0.3); // max dev from mean slope dy/dz
+ tofTrackFinder->SetTyMean(0.); // mean slope dy/dz
+ CbmTofTrackFitter* tofTrackFitter= new CbmTofTrackFitterKF(0,211);
+ TFitter *MyFit = new TFitter(1); // initialize Minuit
+ tofTrackFinder->SetFitter(tofTrackFitter);
+ CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+ tofFindTracks->UseFinder(tofTrackFinder);
+ tofFindTracks->UseFitter(tofTrackFitter);
+ tofFindTracks->SetCorMode(iGenCor); // valid options: 0,1,2,3,4,5,6, 10 - 19
+ tofFindTracks->SetTtTarg(0.041); // target value for inverse velocity, > 0.033 ns/cm!
+ //tofFindTracks->SetTtTarg(0.035); // target value for inverse velocity, > 0.033 ns/cm!
+ tofFindTracks->SetCalParFileName(TofFileFolder + cTrkFile); // Tracker parameter value file name
+ tofFindTracks->SetBeamCounter(5,0,0); // default beam counter
+ tofFindTracks->SetStationMaxHMul(30); // Max Hit Multiplicity in any used station
+
+
+ tofFindTracks->SetT0MAX(dScalFac); // in ns
+ tofFindTracks->SetSIGT(0.08); // default in ns
+ tofFindTracks->SetSIGX(0.3); // default in cm
+ tofFindTracks->SetSIGY(0.45); // default in cm
+ tofFindTracks->SetSIGZ(0.05); // default in cm
+ tofFindTracks->SetUseSigCalib(kFALSE); // ignore resolutions in CalPar file
+ tofTrackFinder->SetSIGLIM(dChi2Lim2*2.); // matching window in multiples of chi2
+ tofTrackFinder->SetChiMaxAccept(dChi2Lim2); // max tracklet chi2
+
+ Int_t iMinNofHits=-1;
+ Int_t iNStations=0;
+ Int_t iNReqStations=3;
+
+ switch (iTrackingSetup){
+ case 0: // bypass mode
+ iMinNofHits=-1;
+ iNStations=1;
+ tofFindTracks->SetStation(0, 5, 0, 0); // Diamond
+ break;
+
+ case 1: // for calibration mode of full setup
+ iMinNofHits=3;
+ iNStations=30;
+ iNReqStations=4;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 0, 2, 2);
+ tofFindTracks->SetStation(2, 0, 1, 2);
+ tofFindTracks->SetStation(3, 0, 0, 2);
+ tofFindTracks->SetStation(4, 0, 2, 1);
+ tofFindTracks->SetStation(5, 0, 1, 1);
+ tofFindTracks->SetStation(6, 0, 0, 1);
+ tofFindTracks->SetStation(7, 0, 2, 3);
+ tofFindTracks->SetStation(8, 0, 1, 3);
+ tofFindTracks->SetStation(9, 0, 0, 3);
+ tofFindTracks->SetStation(10, 0, 2, 0);
+ tofFindTracks->SetStation(11, 0, 1, 0);
+ tofFindTracks->SetStation(12, 0, 0, 0);
+ tofFindTracks->SetStation(13, 0, 2, 4);
+ tofFindTracks->SetStation(14, 0, 1, 4);
+ tofFindTracks->SetStation(15, 0, 0, 4);
+ tofFindTracks->SetStation(16, 0, 4, 0);
+ tofFindTracks->SetStation(17, 0, 3, 0);
+ tofFindTracks->SetStation(18, 0, 4, 1);
+ tofFindTracks->SetStation(19, 0, 3, 1);
+ tofFindTracks->SetStation(20, 0, 4, 2);
+ tofFindTracks->SetStation(21, 0, 3, 2);
+ tofFindTracks->SetStation(22, 0, 4, 3);
+ tofFindTracks->SetStation(23, 0, 3, 3);
+ tofFindTracks->SetStation(24, 0, 4, 4);
+ tofFindTracks->SetStation(25, 0, 3, 4);
+ tofFindTracks->SetStation(26, 9, 0, 0);
+ tofFindTracks->SetStation(27, 9, 0, 1);
+ tofFindTracks->SetStation(28, 6, 0, 0);
+ tofFindTracks->SetStation(29, 6, 0, 1);
+ break;
+ }
+ tofFindTracks->SetMinNofHits(iMinNofHits);
+ tofFindTracks->SetNStations(iNStations);
+ tofFindTracks->SetNReqStations(iNReqStations);
+ tofFindTracks->PrintSetup();
+ run->AddTask(tofFindTracks);
+ }
+ break;
+
+ case 1:
+ {
+
+ }
+
+ case 0:
+
+ default:
+ ;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // =========================================================================
+ // === RICH QA ===
+ // =========================================================================
+
+ CbmRichMCbmQaReal* qaTask = new CbmRichMCbmQaReal();
+ if (taskId < 0) {
+ qaTask->SetOutputDir(Form("result_run%d",runId));
+ } else {
+ qaTask->SetOutputDir(Form("result_run%d_%05d",runId,taskId));
+ }
+// qaTask->DoRestrictToAcc();//restrict to mRICH MAR2019 in histFilling
+// qaTask->XOffsetHistos(12);
+ qaTask->XOffsetHistos(+2.5);
+ qaTask->SetMaxNofDrawnEvents(100);
+ qaTask->SetTotRich(23.7,30.0);
+ qaTask->SetTriggerRichHits(eb_TriggerMinNumberRich);
+ qaTask->SetTriggerTofHits(eb_TriggerMinNumberTof);
+ run->AddTask(qaTask);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Parameter database --------------------------------------------
+ std::cout << std::endl << std::endl << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ Bool_t kParameterMerged = kTRUE;
+ FairParRootFileIo* parIo1 = new FairParRootFileIo(kParameterMerged);
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo1->open(parFile.c_str(),"UPDATE");
+ parIo2->open(parFileList, "in");
+ parIo1->print();
+ parIo2->print();
+ rtdb->setFirstInput(parIo1);
+ rtdb->setSecondInput(parIo2);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+ // ------------------------------------------------------------------------
+
+
+ // ----- Database update ----------------------------------------------
+ rtdb->setOutput(parIo1);
+ rtdb->saveOutput();
+ rtdb->print();
+ rtdb->printParamContexts();
+ // ------------------------------------------------------------------------
+
+
+ //--- House Keeping -------------------------------------------------------
+ // print all important infos in a file
+ std::ofstream outfile;
+ if (taskId < 0) {
+ outfile.open(Form("result_run%d/run_info.dat",runId));
+ } else {
+ outfile.open(Form("result_run%d_%05d/run_info.dat",runId,taskId));
+ }
+ // write inputted data into the file.
+ outfile << "Run: " << runId << std::endl;
+ outfile << "Events: " << nEvents << std::endl;
+ outfile << "parFile: " << parFile << std::endl;
+ outfile << "digiFile: " << digiFile << std::endl;
+ outfile << "recoFile: " << recoFile << std::endl;
+ outfile << "Geometry: " << geoFile << std::endl;
+ outfile << "TrackCalParFile: " << TofFileFolder << cCalId << std::endl;
+ outfile << "TofClusterFile :" << TofFileFolder +cFname << std::endl;
+ outfile << "TofOutput :" << cFname << std::endl << std::endl;
+ outfile << "Trigger:"<< std::endl;
+ outfile << " fixedTimeWindow :" << eb_fixedTimeWindow << std::endl;
+ outfile << " MinNumberT0 :" << eb_TriggerMinNumberT0 << std::endl;
+ outfile << " MinNumberSts :" << eb_TriggerMinNumberSts << std::endl;
+ outfile << " MinNumberMuch :" << eb_TriggerMinNumberMuch << std::endl;
+ outfile << " MinNumberTof :" << eb_TriggerMinNumberTof << std::endl;
+ outfile << " MinNumberRich :" << eb_TriggerMinNumberRich << std::endl;
+ outfile.close();
+ // ------------------------------------------------------------------------
+
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished succesfully." << std::endl;
+ std::cout << "Output file is " << recoFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl;
+ std::cout << "Test passed" << std::endl << "All ok" << std::endl;
+
+ // ----- Resource monitoring ------------------------------------------
+// if ( hasFairMonitor /*Has_Fair_Monitor()*/ ) { // FairRoot Version >= 15.11
+// // Extract the maximal used memory an add is as Dart measurement
+// // This line is filtered by CTest and the value send to CDash
+// FairSystemInfo sysInfo;
+// Float_t maxMemory=sysInfo.GetMaxMemory();
+// std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+// std::cout << maxMemory;
+// std::cout << "</DartMeasurement>" << std::endl;
+//
+// Float_t cpuUsage=ctime/rtime;
+// std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+// std::cout << cpuUsage;
+// std::cout << "</DartMeasurement>" << std::endl;
+//
+// FairMonitor* tempMon = FairMonitor::GetMonitor();
+// tempMon->Print();
+// }
+
+}
+
+void save_hst(TString cstr="status.hst.root", Bool_t bROOT = kFALSE){
+ cout << "save all histograms to file "<<cstr.Data()<<endl;
+ TList* tList(NULL);
+ if(bROOT){
+ tList = gROOT->GetList();
+ } else {
+ tList = gDirectory->GetList();
+ }
+ TIter next(tList);
+ // Write objects to the file
+ TFile *fHist = new TFile(cstr,"RECREATE");
+ {
+ TObject* obj;
+ while( (obj= (TObject*)next()) ){
+ if(obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(TEfficiency::Class())){
+ obj->Write();
+ }
+ }
+ }
+ // fHist->ls();
+ fHist->Close();
+}
diff --git a/macro/rich/mcbm/beamtime/run_reco_mcbm_richOnly.C b/macro/rich/mcbm/beamtime/run_reco_mcbm_richOnly.C
new file mode 100644
index 00000000..4f9eb28c
--- /dev/null
+++ b/macro/rich/mcbm/beamtime/run_reco_mcbm_richOnly.C
@@ -0,0 +1,130 @@
+
+void run_reco_mcbm_richOnly(
+ const string srcfolder = "/lustre/cbm/users/adrian/cbmgitnew/cbmsource/macro/beamtime/mcbm2020/data",
+ const unsigned int runId = 790, // used for the output folder
+ int nEvents = 0,
+ const int taskId = -1
+)
+{
+ const string& parFile = Form("%s/unp_mcbm_params_%d.root",srcfolder.c_str(),runId);
+ const string& digiFile = Form("%s/unp_mcbm_%d.root",srcfolder.c_str(),runId);
+ const string& recoFile = Form("reco_mcbm_mar20_%d.root",runId);
+
+ const Double_t eb_fixedTimeWindow {200.};
+ const Int_t eb_TriggerMinNumberT0 {0};
+ const Int_t eb_TriggerMinNumberSts {0};
+ const Int_t eb_TriggerMinNumberMuch{0};
+ const Int_t eb_TriggerMinNumberTof {0};
+ const Int_t eb_TriggerMinNumberRich{5};
+
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ TTree::SetMaxTreeSize(90000000000);
+
+ TString myName = "run_reco_mcbm_reichOnly";
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ TString workDir = gSystem->Getenv("VMCWORKDIR");
+
+ remove(recoFile.c_str());
+
+// TString setupFile = srcDir + "/geometry/setup/setup_" + geoSetup + ".C";
+// TString setupFunct = "setup_" + geoSetup + "()";
+// gROOT->LoadMacro(setupFile);
+// gROOT->ProcessLine(setupFunct);
+
+ std::cout << std::endl<< "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList *parFileList = new TList();
+
+ TString geoFile = srcDir + "/macro/mcbm/data/mcbm_beam_2020_03.geo.root";
+ TFile* fgeo = new TFile(geoFile);
+ TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+ if (geoMan == nullptr){
+ cout << "<E> FAIRGeom not found in geoFile"<<endl;
+ return;
+ }
+
+ TStopwatch timer;
+ timer.Start();
+ gDebug = 0;
+
+ FairRunAna *run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(digiFile.c_str());
+ run->SetSource(inputSource);
+ run->SetOutputFile(recoFile.c_str());
+ run->SetGenerateRunInfo(kTRUE);
+
+ CbmMcbm2018EventBuilder* eventBuilder = new CbmMcbm2018EventBuilder();
+ // eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::MaximumTimeGap);
+ // eventBuilder->SetMaximumTimeGap(100.);
+ eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::FixedTimeWindow);
+ eventBuilder->SetFixedTimeWindow(eb_fixedTimeWindow);
+ eventBuilder->SetTriggerMinNumberT0(eb_TriggerMinNumberT0);
+ eventBuilder->SetTriggerMinNumberSts(eb_TriggerMinNumberSts);
+ eventBuilder->SetTriggerMinNumberMuch(eb_TriggerMinNumberMuch);
+ eventBuilder->SetTriggerMinNumberTof(eb_TriggerMinNumberTof);
+ eventBuilder->SetTriggerMinNumberRich(eb_TriggerMinNumberRich);
+ eventBuilder->SetFillHistos(kFALSE); // to prevent memory leak???
+
+
+ run->AddTask(eventBuilder);
+
+
+ CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer();
+ hitProd->SetMappingFile("mRICH_Mapping_vert_20190318_elView.geo");
+ hitProd->setToTLimits(23.7,30.0);
+ hitProd->applyToTCut();
+ run->AddTask(hitProd);
+
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ richReco->UseMCbmSetup();
+ run->AddTask(richReco);;
+
+ CbmRichMCbmQaRichOnly* qaTask = new CbmRichMCbmQaRichOnly();
+ if (taskId < 0) {
+ qaTask->SetOutputDir(Form("result_richOnly_run%d",runId));
+ } else {
+ qaTask->SetOutputDir(Form("result_richOnly_run%d_%05d",runId,taskId));
+ }
+// qaTask->DoRestrictToAcc();//restrict to mRICH MAR2019 in histFilling
+ qaTask->XOffsetHistos(+2.5);
+ qaTask->SetMaxNofDrawnEvents(100);
+ //qaTask->SetTotRich(23.7,30.0);
+ qaTask->SetTriggerRichHits(eb_TriggerMinNumberRich);
+ //qaTask->DoWriteHistToFile(false);
+ run->AddTask(qaTask);
+
+ std::cout << std::endl << std::endl << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo1->open(parFile.c_str(),"UPDATE");
+ rtdb->setFirstInput(parIo1);
+ if ( ! parFileList->IsEmpty() ) {
+ parIo2->open(parFileList, "in");
+ rtdb->setSecondInput(parIo2);
+ }
+
+
+ std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+
+ rtdb->setOutput(parIo1);
+ rtdb->saveOutput();
+ rtdb->print();
+
+
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+
+
+ timer.Stop();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished succesfully." << std::endl;
+ std::cout << "Output file is " << recoFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl;
+ std::cout << "Test passed" << std::endl << "All ok" << std::endl;
+
+}
diff --git a/reco/detectors/rich/CMakeLists.txt b/reco/detectors/rich/CMakeLists.txt
index a50373d2..433ee0be 100644
--- a/reco/detectors/rich/CMakeLists.txt
+++ b/reco/detectors/rich/CMakeLists.txt
@@ -68,8 +68,11 @@ CbmRichHitProducer.cxx
mcbm/CbmRichMCbmQa.cxx
mcbm/CbmRichMCbmQaReal.cxx
+mcbm/CbmRichMCbmQaRichOnly.cxx
mcbm/CbmRichMCbmHitProducer.cxx
-#mcbm/CbmRichMCbmTotShifter.cxx
+mcbm/CbmRichMCbmAerogelAna.cxx
+mcbm/CbmRichMCbmToTShifter.cxx
+mcbm/CbmRichMCbmSEDisplay.cxx
qa/CbmRichRingFitterQa.cxx
qa/CbmRichGeoOpt.cxx
diff --git a/reco/detectors/rich/CbmRichRecoLinkDef.h b/reco/detectors/rich/CbmRichRecoLinkDef.h
index 3ab29496..8a5087c8 100644
--- a/reco/detectors/rich/CbmRichRecoLinkDef.h
+++ b/reco/detectors/rich/CbmRichRecoLinkDef.h
@@ -14,7 +14,10 @@
#pragma link C++ class CbmRichMCbmHitProducer;
#pragma link C++ class CbmRichMCbmQa;
#pragma link C++ class CbmRichMCbmQaReal;
-//#pragma link C++ class CbmRichMCbmToTShifter;
+#pragma link C++ class CbmRichMCbmQaRichOnly;
+#pragma link C++ class CbmRichMCbmAerogelAna;
+#pragma link C++ class CbmRichMCbmToTShifter;
+#pragma link C++ class CbmRichMCbmSEDisplay;
//qa
#pragma link C++ class CbmRichGeoTest+;
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.cxx
new file mode 100644
index 00000000..89c2cb2f
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.cxx
@@ -0,0 +1,487 @@
+#include "CbmRichMCbmAerogelAna.h"
+
+#include "TH1D.h"
+#include "TH1.h"
+#include "TCanvas.h"
+#include "TClonesArray.h"
+#include "TF1.h"
+#include "TStyle.h"
+#include "TEllipse.h"
+#include "TLine.h"
+#include "TMarker.h"
+#include "TGeoNode.h"
+#include "TGeoManager.h"
+#include "TGeoBBox.h"
+#include "TMath.h"
+#include <TFile.h>
+
+#include "CbmDigiManager.h"
+#include "CbmRichDigi.h"
+#include "CbmRichHit.h"
+#include "CbmRichRingFinderHoughImpl.h"
+#include "TLatex.h"
+#include "CbmDrawHist.h"
+#include "CbmTrackMatchNew.h"
+#include "CbmRichRing.h"
+#include "CbmRichHit.h"
+#include "CbmMatchRecoToMC.h"
+#include "CbmRichGeoManager.h"
+#include "CbmRichPoint.h"
+#include "CbmRichDraw.h"
+#include "CbmGlobalTrack.h"
+#include "CbmTrdTrack.h"
+#include "CbmTofHit.h"
+#include "CbmTofDigi.h"
+#include "CbmEvent.h"
+#include "CbmTofTracklet.h"
+#include "CbmRichUtil.h"
+
+#include "CbmRichConverter.h"
+
+#include "CbmUtils.h"
+#include "CbmHistManager.h"
+
+#include <iostream>
+#include <string>
+#include <boost/assign/list_of.hpp>
+#include <sstream>
+#include <cmath>
+
+using namespace std;
+using boost::assign::list_of;
+
+#define RichZPos 348.
+
+CbmRichMCbmAerogelAna::CbmRichMCbmAerogelAna() :
+ FairTask("CbmRichMCbmAerogelAna"),
+ fRichHits(nullptr),
+ fRichRings(nullptr),
+ fCbmEvent(nullptr),
+ fHM(nullptr),
+ fXOffsetHisto(12.),
+ fEventNum(0),
+ fNofDrawnRings(0),
+ fNofDrawnRichTofEv(0),
+ fNofDrawnEvents(0),
+ fOutputDir("result")
+{
+
+}
+
+InitStatus CbmRichMCbmAerogelAna::Init()
+{
+ cout << "CbmRichMCbmAerogelAna::Init" << endl;
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (nullptr == ioman) { Fatal("CbmRichMCbmQaReal::Init","RootManager not instantised!"); }
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+
+ if ( ! fDigiMan->IsPresent(ECbmModuleId::kRich) )
+ Fatal("CbmRichMCbmQaReal::Init", "No Rich Digis!");
+
+ if ( ! fDigiMan->IsPresent(ECbmModuleId::kTof) )
+ Fatal("CbmRichMCbmQaReal::Init", "No Tof Digis!");
+
+ fRichHits =(TClonesArray*) ioman->GetObject("RichHit");
+ if (nullptr == fRichHits) { Fatal("CbmRichMCbmAerogelAna::Init", "No Rich Hits!");}
+
+ fRichRings =(TClonesArray*) ioman->GetObject("RichRing");
+ if (nullptr == fRichRings) { Fatal("CbmRichMCbmAerogelAna::Init", "No Rich Rings!");}
+
+
+// fTofHits =(TClonesArray*) ioman->GetObject("TofHit");
+// if (nullptr == fTofHits) { Fatal("CbmRichMCbmQaReal::Init", "No Tof Hits!");}
+
+// fTofTracks =(TClonesArray*) ioman->GetObject("TofTracks");
+// if (nullptr == fTofTracks) { Fatal("CbmRichMCbmQaReal::Init", "No Tof Tracks!");}
+
+// fT0Digis =(TClonesArray*) ioman->GetObject("CbmT0Digi");
+// if (nullptr == fT0Digis) { Fatal("CbmRichMCbmQaReal::Init", "No T0 Digis!");}
+
+ //fT0Digis = ioman->InitObjectAs<std::vector<CbmTofDigi> const *>("T0Digi");
+
+ fCbmEvent =(TClonesArray*) ioman->GetObject("CbmEvent");
+ if (nullptr == fCbmEvent) { Fatal("fTofDigis::Init", "No Event!");}
+
+ InitHistograms();
+
+ return kSUCCESS;
+}
+
+void CbmRichMCbmAerogelAna::InitHistograms()
+{
+ fHM = new CbmHistManager();
+
+ fHM->Create1<TH1D>("fhNofEvents","fhNofEvents;Entries", 1, 0.5, 1.5);
+ fHM->Create1<TH1D>("fhNofCbmEvents","fhNofCbmEvents;Entries", 1, 0.5, 1.5);
+ fHM->Create1<TH1D>("fhNofCbmEventsRing","fhNofCbmEventsRing;Entries", 1, 0.5, 1.5);
+
+ fHM->Create1<TH1D>("fhHitsInTimeslice","fhHitsInTimeslice;Timeslice;#Hits", 200, 1, 200);
+
+ // nof objects per timeslice
+ fHM->Create1<TH1D>("fhNofRichDigisInTimeslice","fhNofRichDigisInTimeslice;# RICH digis / timeslice;Entries", 100, -0.5, 999.5);
+ fHM->Create1<TH1D>("fhNofRichHitsInTimeslice","fhNofRichHitsInTimeslice;# RICH hits / timeslice;Entries", 100, -0.5, 999.5);
+ fHM->Create1<TH1D>("fhNofRichRingsInTimeslice","fhNofRichRingsInTimeslice;# RICH rings / timeslice;Entries", 10, -0.5, 9.5);
+
+ // RICH hits
+ fHM->Create2<TH2D>("fhRichHitXY","fhRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhEventRichHitXY","fhEventRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create1<TH1D>("fhEventRichHitX","fhEventRichHitX;RICH hit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto);
+ fHM->Create1<TH1D>("fhEventRichHitY","fhEventRichHitY;RICH hit Y [cm];Entries", 84, -25.2, 25.2);
+ fHM->Create1<TH1D>("fhRichHitX","fhRichHitX;RICH hit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto);
+ fHM->Create1<TH1D>("fhRichHitY","fhRichHitY;RICH hit Y [cm];Entries", 84, -25.2, 25.2);
+
+ //ToT
+ fHM->Create1<TH1D>("fhRichDigisToT","fhRichDigisToT;ToT [ns];Entries", 601, 9.975, 40.025);
+ fHM->Create1<TH1D>("fhRichHitToT","fhRichHitToT;ToT [ns];Entries", 601, 9.975, 40.025);
+
+ // RICH rings
+ fHM->Create2<TH2D>("fhRichRingXY","fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries", 100, -20 + fXOffsetHisto, 20 + fXOffsetHisto, 100, -20, 20);
+ fHM->Create1<TH1D>("fhRichRingRadius","fhRichRingRadius;Ring radius [cm];Entries", 100, 0., 7.);
+ fHM->Create1<TH1D>("fhNofHitsInRing","fhNofHitsInRing;# hits in ring;Entries", 50, -0.5, 49.5);
+
+ // RICH rings aerogel/ Event
+ fHM->Create2<TH2D>("fhEventRichRingXY","fhEventRichRingXY;Ring center X [cm];Ring center Y [cm];Entries", 100, -20 + fXOffsetHisto, 20 + fXOffsetHisto, 100, -20, 20);
+ fHM->Create1<TH1D>("fhEventRichRingRadius","fhEventRichRingRadius;Ring radius [cm];Entries", 100, 0., 7.);
+ fHM->Create1<TH1D>("fhEventNofHitsInRing","fhEventNofHitsInRing;# hits in ring;Entries", 50, -0.5, 49.5);
+ fHM->Create1<TH1D>("fhEventNofHitsInRingTop","fhEventNofHitsInRingTop;# hits in ring;Entries", 50, -0.5, 49.5);
+ fHM->Create1<TH1D>("fhEventNofHitsInRingBottom","fhEventNofHitsInRingBottom;# hits in ring;Entries", 50, -0.5, 49.5);
+
+ fHM->Create1<TH1D>("fhEventRichHitToT","fhRichHitToT;ToT [ns];Entries", 601, 9.975, 40.025);
+
+ fHM->Create1<TH1D>("fhEventRichRingRadiusTop","fhEventRichRingRadiusTop;Ring radius [cm];Entries", 100, 0., 7.);
+ fHM->Create1<TH1D>("fhEventRichRingRadiusBottom","fhEventRichRingRadiusBottom;Ring radius [cm];Entries", 100, 0., 7.);
+
+ fHM->Create1<TH1D>("fhNofRingsTopBottom","fhNofRingsTopBottom;Entries", 2, -0.5, 1.5);
+
+}
+
+
+void CbmRichMCbmAerogelAna::Exec(Option_t* /*option*/)
+{
+ fEventNum++;
+ fHM->H1("fhNofEvents")->Fill(1);
+ cout << "CbmRichMCbmAerogelAna, event No. " << fEventNum << endl;
+
+ {
+ for (int i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
+ const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
+ fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
+ }
+ }
+
+ int nofRichHits = fRichHits->GetEntries();
+ fHM->H1("fhNofRichHitsInTimeslice")->Fill(nofRichHits);
+ fHM->H1("fhHitsInTimeslice")->Fill(fEventNum,nofRichHits);
+ for (int iH = 0; iH < nofRichHits; iH++) {
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
+ fHM->H2("fhRichHitXY")->Fill(richHit->GetX(), richHit->GetY());
+ fHM->H1("fhRichHitToT")->Fill(richHit->GetToT());
+ fHM->H1("fhRichHitX")->Fill(richHit->GetX());
+ fHM->H1("fhRichHitY")->Fill(richHit->GetY());
+ //printf("HitToT: %f \n", richHit->GetToT());
+ }
+
+
+ //CBMEVENT
+ auto fNCbmEvent = fCbmEvent->GetEntriesFast();
+
+ for (int i=0;i<fNCbmEvent;i++){
+ fHM->H1("fhNofCbmEvents")->Fill(1);
+ CbmEvent* ev = static_cast<CbmEvent*>(fCbmEvent->At(i));
+ std::vector<int> ringIndx; // Rings in CbmEvent
+ std::vector<int> evRichHitIndx;
+
+
+ for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
+ evRichHitIndx.push_back(iRichHit);
+
+ int nofRichRings= fRichRings->GetEntries();
+ for (int l = 0; l < nofRichRings; l++)
+ {
+ CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(l));
+ auto NofRingHits = ring->GetNofHits();
+ for (int m = 0; m < NofRingHits; m++)
+ {
+ auto RingHitIndx = ring->GetHit(m);
+ if (RingHitIndx == iRichHit) {
+ Bool_t used = false;
+ for (auto check : ringIndx) {
+ if (check == l){ used = true; break;}
+ }
+ if (used == false) ringIndx.push_back(l);
+ break;
+ }
+ }
+ }
+ }
+
+
+
+ //DrawEvent(ev, ringIndx, 1);
+
+
+
+ if (ringIndx.size() > 0) { // Ring in CbmEvent
+ //loop over all Hits in a CbmEvent with Ring:
+ for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
+ fHM->H1("fhEventRichHitToT")->Fill(richHit->GetToT());
+ fHM->H2("fhEventRichHitXY")->Fill(richHit->GetX(), richHit->GetY());
+ fHM->H1("fhEventRichHitX")->Fill(richHit->GetX());
+ fHM->H1("fhEventRichHitY")->Fill(richHit->GetY());
+
+ }
+
+
+ //loop over rings in CbmEvent
+ for (unsigned int rings=0; rings < ringIndx.size();rings++) {
+ CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(ringIndx[rings]));
+ if (ring == nullptr) continue;
+
+ fHM->H2("fhEventRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
+ fHM->H1("fhEventRichRingRadius")->Fill(ring->GetRadius());
+ fHM->H1("fhEventNofHitsInRing")->Fill(ring->GetNofHits());
+
+ // for now ignore overlap of hits in crossover region!
+ if (ring->GetCenterY() >= 0.) { // new Aerogel Block
+ fHM->H1("fhEventNofHitsInRingTop")->Fill(ring->GetNofHits());
+ fHM->H1("fhEventRichRingRadiusTop")->Fill(ring->GetRadius());
+ fHM->H1("fhNofRingsTopBottom")->Fill(0);
+ } else { // Aerogel from Mar2019
+ fHM->H1("fhEventNofHitsInRingBottom")->Fill(ring->GetNofHits());
+ fHM->H1("fhEventRichRingRadiusBottom")->Fill(ring->GetRadius());
+ fHM->H1("fhNofRingsTopBottom")->Fill(1);
+ }
+ }
+
+ //DrawRichTofEv(RichTofEv);
+ fHM->H1("fhNofCbmEventsRing")->Fill(1);
+ }
+ } //End CbmEvent loop
+
+
+ RichRings();
+
+}
+
+void CbmRichMCbmAerogelAna::RichRings()
+{
+ int nofRichRings = fRichRings->GetEntries();
+ //fHM->H1("fhNofRichRingsInTimeslice")->Fill(nofRichRings);
+ for (int i = 0; i < nofRichRings; i++) {
+ CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(i));
+ if (ring == nullptr) continue;
+ //DrawRing(ring);
+ fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
+ fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
+ fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
+ }
+}
+
+
+void CbmRichMCbmAerogelAna::DrawHist()
+{
+ cout.precision(4);
+
+ //SetDefaultDrawStyle();
+ double nofEvents = fHM->H1("fhNofCbmEvents")->GetEntries();
+ fHM->ScaleByPattern("fh_.*", 1./nofEvents);
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhNofCbmEvents","rich_mcbm_fhNofCbmEvents", 600 , 600);
+ DrawH1(fHM->H1("fhNofCbmEvents"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhNofCbmEventsRing","rich_mcbm_fhNofCbmEventsRing", 600 , 600);
+ DrawH1(fHM->H1("fhNofCbmEventsRing"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhNofEvents","rich_mcbm_fhNofEvents", 600 , 600);
+ DrawH1(fHM->H1("fhNofEvents"));
+ }
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_XY","rich_mcbm_XY", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhRichHitXY"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhRichRingXY"));
+ }
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_XY_inEvent","rich_mcbm_XY_inEvent", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhEventRichHitXY"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhEventRichRingXY"));
+ }
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_X_and_Y_inEvent","rich_mcbm_X_and_Y_inEvent", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichHitX"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhRichHitY"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_X_and_Y_inEventWithRing","rich_mcbm_X_and_Y_inEventWithRing", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhEventRichHitX"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhEventRichHitY"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_ToT","rich_ToT", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichDigisToT"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhRichHitToT"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_ToT_Event","rich_ToT_Event", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichHitToT"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhEventRichHitToT"));
+ }
+
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_rings","rich_mcbm_rings", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichRingRadius"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhNofHitsInRing"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_rings_inEvent","rich_mcbm_rings_inEvent", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhEventRichRingRadius"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhEventNofHitsInRing"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_Aerogel_Top_Bottom_Hits","rich_Aerogel_Top_Bottom_Hits", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhEventNofHitsInRingTop"));
+ unsigned int sumHitsTop = 0;
+ for (unsigned int i=1;i<=50; i++){
+ sumHitsTop += fHM->H1("fhEventNofHitsInRingTop")->GetBinContent(i)*(i-1);
+ }
+ std::cout<<"Sum Hits Top: "<< sumHitsTop <<std::endl;
+
+ c->cd(2);
+ DrawH1(fHM->H1("fhEventNofHitsInRingBottom"));
+ unsigned int sumHitsBottom = 0;
+ for (unsigned int i=1;i<=50; i++){
+ sumHitsBottom += fHM->H1("fhEventNofHitsInRingBottom")->GetBinContent(i)*(i-1);
+ }
+ std::cout<<"Sum Hits Bottom: "<< sumHitsBottom <<std::endl;
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_Aerogel_Top_Bottom_Radius","rich_Aerogel_Top_Bottom_Radius", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhEventRichRingRadiusTop"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhEventRichRingRadiusBottom"));
+ }
+
+
+ {
+ fHM->CreateCanvas("rich_Aerogel_#RingsTopVsBottom","rich_Aerogel_#RingsTopVsBottom", 1200 , 600);
+ fHM->H1("fhNofRingsTopBottom")->Draw("HIST TEXT");
+ }
+}
+
+
+
+
+
+void CbmRichMCbmAerogelAna::Finish()
+{
+ //std::cout<<"Tracks: "<< fTofTracks->GetEntriesFast() <<std::endl;
+ std::cout<<"Drawing Hists...";
+ DrawHist();
+ std::cout<<"DONE!"<<std::endl;
+
+ if (this->fDoDrawCanvas){
+ fHM->SaveCanvasToImage(fOutputDir,"png");
+ std::cout<<"Canvas saved to Images!"<<std::endl;
+ }
+
+ if (this->fDoWriteHistToFile){
+ TDirectory * oldir = gDirectory;
+ std::string s = fOutputDir + "/RecoHists.root";
+ TFile *outFile = new TFile(s.c_str(),"RECREATE");
+ if (outFile->IsOpen()) {
+ fHM->WriteToFile();
+ std::cout<<"Written to Root-file \""<< s << "\" ...";
+ outFile->Close();
+ std::cout<<"Done!"<<std::endl;
+ }
+ gDirectory->cd( oldir->GetPath() );
+ }
+}
+
+
+void CbmRichMCbmAerogelAna::DrawFromFile(
+ const string& fileName,
+ const string& outputDir)
+{
+ fOutputDir = outputDir;
+
+ if (fHM != nullptr) delete fHM;
+
+ fHM = new CbmHistManager();
+ TFile* file = new TFile(fileName.c_str());
+ fHM->ReadFromFile(file);
+ DrawHist();
+
+ fHM->SaveCanvasToImage(fOutputDir);
+}
+
+
+bool CbmRichMCbmAerogelAna::doToT(CbmRichHit * hit){
+ bool check = false;
+ if ((hit->GetToT() > 23.7) && (hit->GetToT() < 30.0)) check = true;
+
+ return check;
+}
+
+
+Bool_t CbmRichMCbmAerogelAna::cutRadius(CbmRichRing *ring){
+ if (ring->GetRadius() > 1. && ring->GetRadius() < 100. ) return true;
+
+ return false;
+}
+
+
+ClassImp(CbmRichMCbmAerogelAna)
+
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.h b/reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.h
new file mode 100644
index 00000000..a8adaecc
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmAerogelAna.h
@@ -0,0 +1,150 @@
+#ifndef MCBM_RICH_Aerogel
+#define MCBM_RICH_Aerogel
+
+#include "FairTask.h"
+#include "CbmRichRingFinderHoughImpl.h"
+#include "CbmEvent.h"
+class TClonesArray;
+class CbmRichRing;
+class CbmRichHit;
+class CbmTofTracklet;
+class CbmHistManager;
+class TVector3;
+class CbmDigiManager;
+
+#include <vector>
+#include <map>
+
+using namespace std;
+
+
+
+class CbmRichMCbmAerogelAna : public FairTask
+{
+
+public:
+ /**
+ * \brief Standard constructor.
+ */
+ CbmRichMCbmAerogelAna();
+
+ /**
+ * \brief Standard destructor.
+ */
+ virtual ~CbmRichMCbmAerogelAna() {};
+
+ /**
+ * \brief Inherited from FairTask.
+ */
+ virtual InitStatus Init();
+
+ /**
+ * \brief Inherited from FairTask.
+ */
+ virtual void Exec(Option_t* option);
+
+ /**
+ * \brief Inherited from FairTask.
+ */
+ virtual void Finish();
+
+ /**
+ * \brief Set output directory where you want to write results (figures and json).
+ * \param[in] dir Path to the output directory.
+ */
+ void SetOutputDir(const string& dir) {fOutputDir = dir;}
+
+
+ /**
+ * \brief Draw histogram from file
+ */
+ void DrawFromFile(
+ const string& fileName,
+ const string& outputDir);
+
+
+ /**
+ * Apply restriction to full mRICH Acceptance (for Simulations)
+ */
+ void DoDrawCanvas (bool val=true){
+ fDoDrawCanvas = val;
+ }
+
+ /**
+ * Apply restriction to full mRICH Acceptance (for Simulations)
+ */
+ void DoWriteHistToFile (bool val=true){
+ fDoWriteHistToFile = val;
+ }
+
+
+ /**
+ * Move X-Position of mRICH in Histograms (e.g. for Geometry changes)
+ */
+ void XOffsetHistos (Double_t offset = 0.){
+ fXOffsetHisto = offset;
+ }
+
+private:
+
+ CbmDigiManager* fDigiMan = nullptr;
+
+ TClonesArray* fRichHits;
+
+ TClonesArray* fRichRings;
+
+ TClonesArray* fCbmEvent;
+
+ CbmHistManager* fHM;
+
+ Double_t fXOffsetHisto;
+
+ Int_t fEventNum;
+
+ Int_t fNofDrawnRings;
+
+ Int_t fNofDrawnRichTofEv;
+
+ Int_t fNofDrawnEvents;
+
+
+ string fOutputDir; // output dir for results
+
+ bool fDoWriteHistToFile = true;
+ bool fDoDrawCanvas = true;
+
+ /**
+ * \brief Initialize histograms.
+ */
+ void InitHistograms();
+
+ /**
+ * \brief Draw histograms.
+ */
+ void DrawHist();
+
+ void RichRings();;
+
+ bool doToT(CbmRichHit* hit);
+
+ Bool_t cutRadius(CbmRichRing *ring);
+
+
+ /**
+ * \brief Copy constructor.
+ */
+ CbmRichMCbmAerogelAna(const CbmRichMCbmAerogelAna&);
+
+ /**
+ * \brief Assignment operator.
+ */
+ CbmRichMCbmAerogelAna& operator=(const CbmRichMCbmAerogelAna&);
+
+
+
+
+
+ ClassDef(CbmRichMCbmAerogelAna,1)
+};
+
+#endif
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx
index 7c06b7b5..9b4a3154 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx
@@ -1,10 +1,10 @@
#include "CbmRichMCbmHitProducer.h"
+#include "CbmRichDetectorData.h" // for CbmRichPmtData, CbmRichPixelData
#include "CbmRichPoint.h"
#include "CbmRichHit.h"
#include "CbmRichDigi.h"
-#include "CbmRichDetectorData.h" // for CbmRichPmtData, CbmRichPixelData
#include "TClonesArray.h"
#include "CbmRichGeoManager.h"
#include "CbmRichDigiMapManager.h"
@@ -126,15 +126,17 @@ void CbmRichMCbmHitProducer::ProcessData(
CbmEvent* event)
{
if (event != NULL) {
- LOG(info) << "CbmRichHitProducer CbmEvent mode. CbmEvent # " << event->GetNumber();
+ LOG(info) << "CbmRichMCbmHitProducer CbmEvent mode. CbmEvent # " << event->GetNumber();
Int_t nofDigis = event->GetNofData(ECbmDataType::kRichDigi);
- LOG(info) << "nofDigis: " << nofDigis;
+ //LOG(info) << "nofDigis: " << nofDigis;
+ fNofHits = 0;
for (Int_t iDigi = 0; iDigi < nofDigis; iDigi++) {
Int_t digiIndex = event->GetIndex(ECbmDataType::kRichDigi, iDigi);
ProcessDigi(event, digiIndex);
}
-
+ LOG(info) << "nofDigis: " <<nofDigis<< "\t\t "<<"nofHits : " << fNofHits;
+ fNofHits = 0;
}
else {
for(Int_t iDigi = 0; iDigi < fDigiMan->GetNofDigis(ECbmModuleId::kRich); iDigi++){
@@ -170,7 +172,9 @@ void CbmRichMCbmHitProducer::ProcessDigi(
//std::cout<<std::hex<<digi->GetAddress()<<std::dec<<" "<<data->fX<<" "<<data->fY<<" "<<data->fZ<<std::endl;
if (!RestrictToAcc(posPoint)) return;
- AddHit(event,posPoint, digi->GetTime(),digi->GetToT(), digiIndex);
+ if (!RestrictToFullAcc(posPoint)) return;
+ if (!RestrictToAerogelAccDec2019(posPoint)) return;
+ AddHit(event,posPoint, digi, digiIndex,data->fPmtId);
}
}
@@ -179,9 +183,10 @@ void CbmRichMCbmHitProducer::ProcessDigi(
void CbmRichMCbmHitProducer::AddHit(
CbmEvent* event,
TVector3 &posHit,
- Double_t time,
- Double_t tot,
- Int_t index)
+ const CbmRichDigi *digi,
+ Int_t index,
+ Int_t PmtId
+ )
{
Int_t nofHits = fRichHits->GetEntriesFast();
new((*fRichHits)[nofHits]) CbmRichHit();
@@ -190,11 +195,15 @@ void CbmRichMCbmHitProducer::AddHit(
hit->SetDx(fHitError);
hit->SetDy(fHitError);
hit->SetRefId(index);
- hit->SetTime(time);
- hit->SetToT(tot);
+ hit->SetTime(digi->GetTime());
+ hit->SetToT(digi->GetToT());
+ hit->SetAddress(digi->GetAddress());
+ hit->SetPmtId(PmtId);
+ // Not Implemented in Digi: hit->SetPmtId(digi->GetPmtId());
if (event != NULL) {
event->AddData(ECbmDataType::kRichHit, nofHits);
+ fNofHits++;
}
}
@@ -264,5 +273,23 @@ bool CbmRichMCbmHitProducer::RestrictToFullAcc(Double_t x, Double_t y)
return inside;
}
+bool CbmRichMCbmHitProducer::RestrictToAerogelAccDec2019(TVector3 &pos){
+ Double_t x = pos.X();
+ Double_t y = pos.Y();
+
+ return this->RestrictToAerogelAccDec2019(x, y);
+}
+
+bool CbmRichMCbmHitProducer::RestrictToAerogelAccDec2019(Double_t x, Double_t y)
+{ //check if Track is in mRICH acceptance
+ if (fRestrictToAerogelAccDec2019 == false) return true;
+ //bool inside = false;
+
+ if (y > 13.5) return false;
+ if (y > 8.0 && x < 12.5) return false;
+
+ return true;
+}
+
ClassImp(CbmRichMCbmHitProducer)
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h
index f397d6ca..4cfae85b 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h
@@ -7,9 +7,6 @@
#include "CbmDigiManager.h" // for ROOTCLING
-#include <map> // for map
-#include <string> // for string
-
class TClonesArray;
class TVector3;
class CbmEvent;
@@ -133,6 +130,7 @@ private:
bool fDoToT = false;
bool fRestrictToAcc = false;
bool fRestrictToFullAcc = false;
+ bool fRestrictToAerogelAccDec2019 = false;
double fToTLimitLow = 0.;
double fToTLimitHigh = 1000.;
@@ -140,6 +138,8 @@ private:
Int_t fEventNum; // event number
+ Int_t fNofHits = 0;
+
Double_t fHitError;
std::string fMappingFile;
@@ -154,6 +154,8 @@ private:
bool RestrictToFullAcc(TVector3 &pos);
bool RestrictToFullAcc(Double_t x, Double_t y);
+ bool RestrictToAerogelAccDec2019(TVector3 &pos);
+ bool RestrictToAerogelAccDec2019(Double_t x, Double_t y);
/**
* \brief Add hit to the output array (and) CbmEvent if it is not NULL.
@@ -162,9 +164,9 @@ private:
void AddHit(
CbmEvent* event,
TVector3 &posHit,
- Double_t time,
- Double_t tot,
- Int_t index);
+ const CbmRichDigi *digi,
+ Int_t index,
+ Int_t PmtId );
/**
* \brief Copy constructor.
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.cxx
index 6a1b3b7f..53a72735 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.cxx
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.cxx
@@ -13,10 +13,13 @@
#include "TGeoManager.h"
#include "TGeoBBox.h"
#include "TMath.h"
+#include "TSystem.h"
#include <TFile.h>
+
#include <TLegend.h>
#include <TBox.h>
+
#include "CbmRichDigi.h"
#include "CbmRichHit.h"
#include "CbmRichRingFinderHoughImpl.h"
@@ -24,7 +27,6 @@
#include "CbmDrawHist.h"
#include "CbmTrackMatchNew.h"
#include "CbmRichRing.h"
-#include "CbmRichHit.h"
#include "CbmMatchRecoToMC.h"
#include "CbmRichGeoManager.h"
#include "CbmRichPoint.h"
@@ -38,14 +40,17 @@
#include "CbmTofTracklet.h"
#include "CbmRichUtil.h"
+#include "CbmRichMCbmSEDisplay.h"
+
#include "CbmRichConverter.h"
#include "CbmUtils.h"
#include "CbmHistManager.h"
#include "CbmDigiManager.h"
+//#include <CbmSetup.h>
-#include <iomanip>
#include <iostream>
+#include <iomanip>
#include <string>
#include <boost/assign/list_of.hpp>
#include <sstream>
@@ -54,7 +59,8 @@
using namespace std;
using boost::assign::list_of;
-#define RichZPos 348.
+//double RichZPos = 348.;
+double RichZPos = 355.2;
CbmRichMCbmQaReal::CbmRichMCbmQaReal() :
FairTask("CbmRichMCbmQaReal"),
@@ -66,10 +72,14 @@ CbmRichMCbmQaReal::CbmRichMCbmQaReal() :
fCbmEvent(nullptr),
fHM(nullptr),
fXOffsetHisto(0.),
+ fTotRichMin(23.7),
+ fTotRichMax(30.0),
fEventNum(0),
fNofDrawnRings(0),
fNofDrawnRichTofEv(0),
- fNofDrawnEvents(0),
+ fMaxNofDrawnEvents(100),
+ fTriggerRichHits(0),
+ fTriggerTofHits(0),
fOutputDir("result")
{
@@ -106,12 +116,57 @@ InitStatus CbmRichMCbmQaReal::Init()
// fT0Digis =(TClonesArray*) ioman->GetObject("CbmT0Digi");
// if (nullptr == fT0Digis) { Fatal("CbmRichMCbmQaReal::Init", "No T0 Digis!");}
+// Get a pointer to the previous already existing data level
+ fT0Digis = ioman->InitObjectAs<std::vector<CbmTofDigi> const *>("T0Digi");
+ /*if ( ! fT0DigiVec ) {
+ fT0DigiArr = dynamic_cast<TClonesArray*>(ioman->GetObject("T0Digi"));
+ if ( ! fT0DigiArr ) {
+ LOG(fatal) << "No TClonesArray with T0 digis found.";
+ }
+ }*/
+
+ // ---- Init Z Positioning ---------
+/* CbmSetup *pSetup=CbmSetup::Instance();
+ TString geoTag;
+ if ( pSetup->IsActive(ECbmModuleId::kRich) ) {
+ pSetup->GetGeoTag(ECbmModuleId::kRich, geoTag);
+ }
+
+ if (geoTag == "v20d") RichZPos = 340.8;
+ std::cout<<"mRICH Geo Tag is "<< geoTag<< ". Z Position of PMT plane set to "<<RichZPos<<"."<<std::endl;
+ */
+ // ---------------------------------
fCbmEvent =(TClonesArray*) ioman->GetObject("CbmEvent");
if (nullptr == fCbmEvent) { Fatal("CbmRichMCbmQaReal::Init", "No Event!");}
InitHistograms();
-
+
+ //--- Init Single Event display -----------------------------
+ fSeDisplay = new CbmRichMCbmSEDisplay(fHM);
+ fSeDisplay->SetRichHits(fRichHits);
+ fSeDisplay->SetRichRings(fRichRings);
+ fSeDisplay->SetTofTracks(fTofTracks);
+ fSeDisplay->SetTotRich(fTotRichMin,fTotRichMax);
+ fSeDisplay->SetMaxNofDrawnEvents(fMaxNofDrawnEvents);
+ fSeDisplay->SetOutDir(fOutputDir);
+ fSeDisplay->XOffsetHistos(fXOffsetHisto);
+ //-----------------------------------------------------------
+
+
+ //--- Init Single Event display for Correlated Tracks/Rings---
+ fSeDsply_TR = new CbmRichMCbmSEDisplay(fHM);
+ fSeDsply_TR->SetRichHits(fRichHits);
+ fSeDsply_TR->SetRichRings(fRichRings);
+ fSeDsply_TR->SetTofTracks(fTofTracks);
+ fSeDsply_TR->SetTotRich(fTotRichMin,fTotRichMax);
+ fSeDsply_TR->SetMaxNofDrawnEvents(fMaxNofDrawnEvents);
+ fSeDsply_TR->SetOutDir(fOutputDir);
+ fSeDsply_TR->XOffsetHistos(fXOffsetHisto);
+ fSeDsply_TR->SetCanvasDir("SE_Corr");
+ //-----------------------------------------------------------
+
+
return kSUCCESS;
}
@@ -132,6 +187,8 @@ void CbmRichMCbmQaReal::InitHistograms()
// RICH hits
fHM->Create2<TH2D>("fhRichHitXY","fhRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhRichDigiPixelRate","fhRichDigiPixelRate;RICH Digi X [cm];RICH Digi Y [cm];Hz", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+
// RICH digis, the limits of log histograms are set in Exec method
fHM->Create1<TH1D>("fhRichDigisTimeLog", "fhNofRichDigisTimeLog;Time [ns];Entries", 400, 0., 0.);
@@ -157,15 +214,69 @@ void CbmRichMCbmQaReal::InitHistograms()
//ToT
fHM->Create1<TH1D>("fhRichDigisToT","fhRichDigisToT;ToT [ns];Entries", 601, 9.975, 40.025);
fHM->Create1<TH1D>("fhRichHitToT","fhRichHitToT;ToT [ns];Entries", 601, 9.975, 40.025);
+ fHM->Create1<TH1D>("fhRichHitToTEvent","fhRichHitToTEvent;ToT [ns];Entries", 601, 9.975, 40.025);
+ fHM->Create2<TH2D>("fhRichHitXYEvent","fhRichHitXYEvent;RICH hit X [cm];RICH hit Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+
// RICH rings
- fHM->Create2<TH2D>("fhRichRingXY","fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries", 100, -20 + fXOffsetHisto, 20 + fXOffsetHisto, 100, -20, 20);
+ fHM->Create2<TH2D>("fhRichRingXY","fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhRichRingXY_goodTrack","fhRichRingXY_goodTrack;Ring center X [cm];Ring center Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhRichRing_goodTrackXY","fhRichRing_goodTrackXY;Track center X [cm];Track center Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+
+ fHM->Create2<TH2D>("fhRichRingRadiusY","fhRichRingRadiusY;Ring Radius [cm]; Y position[cm];Entries", 70, -0.05, 6.95, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhRichHitsRingRadius","fhRichHitsRingRadius;#Rich Hits/Ring; Ring Radius [cm];Entries", 50, -0.5, 49.5, 70, -0.05, 6.95);
+
fHM->Create1<TH1D>("fhRichRingRadius","fhRichRingRadius;Ring radius [cm];Entries", 100, 0., 7.);
fHM->Create1<TH1D>("fhNofHitsInRing","fhNofHitsInRing;# hits in ring;Entries", 50, -0.5, 49.5);
+ fHM->Create1<TH1D>("fhRichRingRadius_goodRing","fhRichRingRadius_goodRing;Ring radius [cm];Entries", 100, 0., 7.);
+ fHM->Create1<TH1D>("fhNofHitsInRing_goodRing","fhNofHitsInRing_goodRing;# hits in ring;Entries", 50, -0.5, 49.5);
+
//Tof Rich correlation
fHM->Create2<TH2D>("fhTofRichX","fhTofRichX;Rich Hit X [cm];TofHit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofRichX_stack1","fhTofRichX_stack1;Rich Hit X [cm];TofHit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofRichX_stack2","fhTofRichX_stack2;Rich Hit X [cm];TofHit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofRichX_stack3","fhTofRichX_stack3;Rich Hit X [cm];TofHit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
fHM->Create2<TH2D>("fhTofRichY","fhTofRichY;Rich Hit Y [cm];TofHit Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+ //fHM->Create2<TH2D>("fhTofRichRingHitX","fhTofRichRingHitX;Rich Ring Hit X [cm];TofHit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitX","fhTofTrackRichHitX;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY","fhTofTrackRichHitY;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+
+ fHM->Create2<TH2D>("fhTofTrackHitRichHitX_oBetacuts_dtime","fhTofTrackHitRichHitX_oBetacuts_dtime;Rich Hit X [cm];TofHit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackHitRichHitY_oBetacuts_dtime","fhTofTrackHitRichHitY_oBetacuts_dtime;Rich Hit Y [cm];TofHit Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_cuts","fhTofTrackRichHitX_cuts;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_cuts","fhTofTrackRichHitY_cuts;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_oBetacuts","fhTofTrackRichHitX_oBetacuts;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_oBetacuts","fhTofTrackRichHitY_oBetacuts;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+ fHM->Create1<TH1D>("fhTofTrackRichHitTime_oBetacuts","fhTofTrackRichHitTime_oBetacuts;#Delta Time [ns];Entries", 280, -70. , 70. );
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_uBetacuts","fhTofTrackRichHitX_uBetacuts;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_uBetacuts","fhTofTrackRichHitY_uBetacuts;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_oBetacuts_dtime","fhTofTrackRichHitX_oBetacuts_dtime;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_oBetacuts_dtime","fhTofTrackRichHitY_oBetacuts_dtime;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_oBetacuts_dtime_4","fhTofTrackRichHitX_oBetacuts_dtime_4;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_oBetacuts_dtime_4","fhTofTrackRichHitY_oBetacuts_dtime_4;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_oBetacuts_dtime_6","fhTofTrackRichHitX_oBetacuts_dtime_6;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_oBetacuts_dtime_6","fhTofTrackRichHitY_oBetacuts_dtime_6;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_oBetacuts_dtime_8","fhTofTrackRichHitX_oBetacuts_dtime_8;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_oBetacuts_dtime_8","fhTofTrackRichHitY_oBetacuts_dtime_8;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichHitX_oBetacuts_dtime_10","fhTofTrackRichHitX_oBetacuts_dtime_10;RICH Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichHitY_oBetacuts_dtime_10","fhTofTrackRichHitY_oBetacuts_dtime_10;RICH Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofTrackRichRingHitX","fhTofTrackRichRingHitX;RICH Ring Hit X [cm];TofTrack X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichRingHitY","fhTofTrackRichRingHitY;RICH Ring Hit Y [cm];TofTrack Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
+
+ fHM->Create2<TH2D>("fhTofHitRichRingHitX","fhTofHitRichRingHitX;RICH Ring Hit X [cm];Tof Hit X [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofHitRichRingHitY","fhTofHitRichRingHitY;RICH Ring Hit Y [cm];Tof Hit Y [cm];Entries", 84, -25.2, 25.2, 200, -80, 80);
//Tof Rich correlation
fHM->Create2<TH2D>("fhTofRichX_zoomed","fhTofRichX_zoomed;Rich Hit X [cm];TofHit X [cm];Entries", 27, -8.1 + fXOffsetHisto, 8.1 + fXOffsetHisto, 180, -15, 75);
@@ -182,27 +293,48 @@ void CbmRichMCbmQaReal::InitHistograms()
fHM->Create2<TH2D>("fhTofClosTrackRichRingXY","fhTofClosTrackRichRingXY; X [cm]; Y [cm];Entries", 100, -20 + fXOffsetHisto, 20 + fXOffsetHisto, 120, -10, 20); //1bin == 2mm
fHM->Create3<TH3D>("fhTofXYZ","fhTofXYZ;Tof Hit X [cm];TofHit Z [cm];Tof Hit Y [cm];Entries", 100, -20, 20, 141, 230., 370., 100, -20, 20);
+ fHM->Create2<TH2D>("fhTofHitsXY","fhTofHitsXY;Tof Hit X [cm];Tof Hit Y [cm];Entries", 100, -20, 20, 200, -80, 80);
fHM->Create1<TH1D>("fhTofHitsZ","fhTofHitsZ;Z [cm];Entries", 350, -0.5, 349.5);
fHM->Create2<TH2D>("fhTofHitsXZ","fhTofHitsXZ;Z [cm];X [cm];Entries", 350, -0.5, 349.5, 400, -50, 110);
//Tof Tracks
- fHM->Create1<TH1D>("fhTofTracksPerEvent","fhTofTracksPerEvent;NofTracks/RichEvent;Entries", 20, -0.5, 19.5);
- fHM->Create2<TH2D>("fhTofTracksXY","fhTofTracksXY;X[cm];Y[cm];NofTracks/cm^2",50 , -20 + fXOffsetHisto, 30 + fXOffsetHisto, 180,-90,90); // projected in RICH Plane
+ fHM->Create1<TH1D>("fhTofTracksPerEvent","fhTofTracksPerEvent;NofTracks/Event;Entries", 20, -0.5, 19.5);
+ fHM->Create1<TH1D>("fhTofTracksPerRichEvent","fhTofTracksPerRichEvent;NofTracks/RichEvent;Entries", 20, -0.5, 19.5);
+ fHM->Create2<TH2D>("fhTofTracksXY","fhTofTracksXY;X[cm];Y[cm];NofTracks/cm^2",250 , -100, 150, 300,-150,150);//50 , -20 + fXOffsetHisto, 30 + fXOffsetHisto, 180,-90,90); // projected in RICH Plane
+
+ fHM->Create2<TH2D>("fhTofTracksXY_Target","fhTofTracksXY_Target;X[cm];Y[cm];NofTracks/cm^2",100 , -50, 50, 100,-50,50);//50 , -20 + fXOffsetHisto, 30 + fXOffsetHisto, 180,-90,90); // projected in Z=0
+ fHM->Create2<TH2D>("fhGoodRingsXY_TargetPos","fhGoodRingsXY_TargetPos;X[cm];Y[cm];NofTracks/cm^2",100 , -50, 50, 100,-50,50);//50 , -20 + fXOffsetHisto, 30 + fXOffsetHisto, 180,-90,90); // projected in Z=0
+
+ fHM->Create2<TH2D>("fhTofTrackRichRingX","fhTofTrackRichRingX;Rich Ring Center X [cm];TofTrack X [cm];Entries", 100, -20 + fXOffsetHisto, 20 + fXOffsetHisto, 400, -50, 110);
+ fHM->Create2<TH2D>("fhTofTrackRichRingY","fhTofTrackRichRingY;Ring Ring Center Y [cm];TofTrack Y [cm];Entries", 125, -25, 25, 200, -80, 80);
+
fHM->Create2<TH2D>("fhTofTracksXYRICH","fhTofTracksXYRICH;X[cm];Y[cm];NofTracks/cm^2", 50 , -20 + fXOffsetHisto, 30 + fXOffsetHisto, 180,-90,90); // projected in RICH Plane
fHM->Create1<TH1D>("fhNofTofTracks","fhNofTofTracks;X[cm];Y[cm];NofTracks", 4,0,4); // 1: All 2: left; 3: right; 4: RICH
+ fHM->Create1<TH1D>("fhRingTrackDistance","fhRingTrackDistance;RingTrackDistance [cm];Entries", 31, -0.25, 14.75);
+
fHM->Create1<TH1D>("fhTrackRingDistance","fhTrackRingDistance;TrackRingDistance [cm];Entries", 31, -0.5, 30.5);
+ fHM->Create1<TH1D>("fhTrackRingDistanceOnTarget","fhTrackRingDistanceOnTarget;TrackRingDistance [cm];Entries", 31, -0.5, 30.5);
+ fHM->Create1<TH1D>("fhTrackRingDistanceOffTarget","fhTrackRingDistanceOffTarget;TrackRingDistance [cm];Entries", 31, -0.5, 30.5);
fHM->Create2<TH2D>("fhTrackRingDistanceVSRingradius","fhTrackRingDistanceVSRingradius;TrackRingDistance [cm];RingRadius [cm];Entries", 81, -0.5, 80.5, 100, 0., 7.);
+
+ fHM->Create2<TH2D>("fhTrackRingDistanceVSRingChi2","fhTrackRingDistanceVSRingChi2;TrackRingDistance [cm];\\Chi^2;Entries", 60, 0, 30, 101, 0.,10.1);
+ fHM->Create2<TH2D>("fhTrackRingDistanceVSRingChi2_goodRing","fhTrackRingDistanceVSRingChi2_goodRing;TrackRingDistance [cm];\\Chi^2;Entries", 40, 0, 10.0, 101, 0.,10.1);
+
+ fHM->Create1<TH1D>("fhTrackRingDistance_corr","fhTrackRingDistance_corr;TrackRingDistance [cm];Entries", 31, -0.5, 30.5);
- fHM->Create1<TH1D>("fhTofBetaTracksWithHitsNoRing","fhTofBetaTracksWithHitsNoRing; \\beta;Entries", 111, -0.005, 1.105);
- fHM->Create1<TH1D>("fhTofBetaTracksWithHits","fhTofBetaTracksWithHits; \\beta;Entries", 111, -0.005, 1.105);
- fHM->Create1<TH1D>("fhTofBetaTracksNoRing","fhTofBetaTracksNoRing; \\beta;Entries", 111, -0.005, 1.105);
- fHM->Create1<TH1D>("fhTofBetaTrackswithClosestRingInRange","fhTofBetaTrackswithClosestRingInRange; \\beta;Entries", 111, -0.005, 1.105);
+ fHM->Create1<TH1D>("fhTofBetaTracksWithHitsNoRing","fhTofBetaTracksWithHitsNoRing; \\beta;Entries", 151, -0.005, 1.505);
+ fHM->Create1<TH1D>("fhTofBetaTracksWithHits","fhTofBetaTracksWithHits; \\beta;Entries", 151, -0.005, 1.505);
+ fHM->Create1<TH1D>("fhTofBetaTracksNoRing","fhTofBetaTracksNoRing; \\beta;Entries", 151, -0.005, 1.505);
+ fHM->Create1<TH1D>("fhTofBetaTrackswithClosestRingInRange","fhTofBetaTrackswithClosestRingInRange; \\beta;Entries", 151, -0.005, 1.505);
+
+ fHM->Create1<TH1D>("fhRichRingBeta","fhRichRingBeta; \\beta;Entries", 151, -0.005, 1.505);
+ fHM->Create1<TH1D>("fhRichRingBeta_GoodRing","fhRichRingBeta_GoodRing; \\beta;Entries", 151, -0.005, 1.505);
- fHM->Create1<TH1D>("fhTofBetaRing","fhTofBetaRing; \\beta;Entries", 111, -0.005, 1.105);
- fHM->Create1<TH1D>("fhTofBetaAll","fhTofBetaAll; \\beta;Entries", 111, -0.005, 1.105);
- fHM->Create2<TH2D>("fhTofBetaVsRadius","fhTofBetaVsRadius; \\beta;ring radius [cm];Entries", 102, -0.005, 1.015, 100, 0., 7.);
- fHM->Create2<TH2D>("fhTofBetaRingDist","fhTofBetaRingDist; \\beta;ring Dist [cm];Entries", 102, -0.005, 1.015, 100, 0., 20.);
- fHM->Create1<TH1D>("fhTofBetaAllFullAcc","fhTofBetaAllFullAcc; \\beta;Entries", 221, -1.105, 1.105);
+ fHM->Create1<TH1D>("fhTofBetaRing","fhTofBetaRing; \\beta;Entries", 151, -0.005, 1.505);
+ fHM->Create1<TH1D>("fhTofBetaAll","fhTofBetaAll; \\beta;Entries", 151, -0.005, 1.505);
+ fHM->Create2<TH2D>("fhTofBetaVsRadius","fhTofBetaVsRadius; \\beta;ring radius [cm];Entries", 151, -0.005, 1.505, 100, 0., 7.);
+ fHM->Create2<TH2D>("fhTofBetaRingDist","fhTofBetaRingDist; \\beta;ring Dist [cm];Entries", 151, -0.005, 1.505, 100, 0., 20.);
+ fHM->Create1<TH1D>("fhTofBetaAllFullAcc","fhTofBetaAllFullAcc; \\beta;Entries", 301, -1.505, 1.505);
fHM->Create1<TH1D>("fhRingDeltaTime","fhRingDeltaTime; \\Delta Time/ns;Entries", 101, -10.1, 10.1);
fHM->Create1<TH1D>("fhRingToTs","fhRingToTs; ToT/ns;Entries", 601, 9.975, 40.025);
@@ -222,6 +354,25 @@ void CbmRichMCbmQaReal::InitHistograms()
fHM->Create1<TH1D>("fhNofInnerHits","fhNofInnerHits;#Hits;Entries",31,-0.5,30.5);
fHM->Create1<TH1D>("fhDiRICHsInRegion","fhNofInnerHits;#Hits;DiRICH",4096,28672,32767);
+
+ fHM->Create1<TH1D>("fhBlobTrackDistX","fhBlobTrackDistX; |TofTrackX - MAPMT center X| [cm];Entries", 30, -0.5, 29.5);
+ fHM->Create1<TH1D>("fhBlobTrackDistY","fhBlobTrackDistY; |TofTrackY - MAPMT center Y| [cm];Entries", 30, -0.5, 29.5);
+ fHM->Create1<TH1D>("fhBlobTrackDist","fhBlobTrackDist; |TofTrack - MAPMT center dist| [cm];Entries", 30, -0.5, 29.5);
+
+ fHM->Create1<TH1D>("fhNofBlobEvents","fhNofBlobEvents;;#Events with min. one Blob", 1, 0.5, 1.5);
+ fHM->Create1<TH1D>("fhNofBlobsInEvent","fhNofBlobsInEvent;#Blobs in Event;Entries", 36, 0.5, 36.5);
+
+ fHM->Create1<TH1D>("fhRichDigisConsecTime","fhRichDigisConsecTime;consecutive time [ns];Entries", 500, -0.5, 499.5);
+ fHM->Create1<TH1D>("fhRichDigisConsecTimeTOT","fhRichDigisConsecTimeTOT;consecutive time [ns];Entries", 500, -0.5, 499.5);
+
+ fHM->Create1<TH1D>("fhNofHitsInGoodRing","fhNofHitsInGoodRing;# hits in ring;Entries", 50, -0.5, 49.5);
+ fHM->Create1<TH1D>("fhTracksWithRings","fhTracksWithRings;Scenarios;Entries", 5, -0.5, 4.5);
+
+ fHM->Create1<TH1D>("fhRichRingChi2","fhRichRingChi2;\\Chi^2;Entries", 101, 0.,10.1);
+ fHM->Create1<TH1D>("fhRichRingChi2_goodRing","fhRichRingChi2_goodRing;\\Chi^2;Entries", 101, 0.,10.1);
+
+ fHM->Create2<TH2D>("fhTofTracksXYRICH_Accectance","fhTofTracksXYRICH_Accectance;X[cm];Y[cm];NofTracks/cm^2", 50 , -20 + fXOffsetHisto, 30 + fXOffsetHisto, 180,-90,90); // projected in RICH Plane
+
}
@@ -231,117 +382,181 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
fHM->H1("fhNofEvents")->Fill(1);
cout << "CbmRichMCbmQaReal, event No. " << fEventNum << endl;
- if (fEventNum == 2) {
- double minTime = std::numeric_limits<double>::max();
- for (int i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
- const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
- // fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
- if (richDigi->GetTime() < minTime) minTime = richDigi->GetTime();
- }
+ if (fDigiHitsInitialized == false) {
+ auto nOfCbmEvent = fCbmEvent->GetEntriesFast();
+ if (nOfCbmEvent > 0){
+ CbmEvent* ev = static_cast<CbmEvent*>(fCbmEvent->At(0));
+ /*if ((fDigiMan->GetNofDigis(ECbmModuleId::kRich)> 10) &&
+ //(fDigiMan->GetNofDigis(ECbmModuleId::kSts) > 30) ||
+ //(fRichRings->GetEntries()> 0) ||
+ (fDigiMan->GetNofDigis(ECbmModuleId::kTof) > 10)
+ )*/
+ if (ev != nullptr)
+ {
+ double minTime = ev->GetStartTime();//std::numeric_limits<double>::max();
+ /* for (int i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
+ const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
+ // fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
+ if (richDigi->GetTime() < minTime) minTime = richDigi->GetTime();
+ }*/
+
+ double dT = 40e9;
+ double dTZoom1 = 0.8e9;
+ double dTZoom2 = 0.008e9;
+ fHM->H1("fhRichDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
+ fHM->H1("fhRichDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
+ fHM->H1("fhRichDigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
+
+ fHM->H1("fhRichRingsTimeLog")->GetXaxis()->SetLimits(minTime, minTime + dT);
+ fHM->H1("fhRichRingsTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime + dTZoom1);
+ fHM->H1("fhRichRingsTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
+
+ fHM->H1("fhTofDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
+ fHM->H1("fhTofDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
+ fHM->H1("fhTofDigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
- double dT = 40e9;
- double dTZoom1 = 0.8e9;
- double dTZoom2 = 0.008e9;
- fHM->H1("fhRichDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
- fHM->H1("fhRichDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
- fHM->H1("fhRichDigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
+ fHM->H1("fhStsDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
+ fHM->H1("fhStsDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
+ fHM->H1("fhStsDigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
- fHM->H1("fhRichRingsTimeLog")->GetXaxis()->SetLimits(minTime, minTime + dT);
- fHM->H1("fhRichRingsTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime + dTZoom1);
- fHM->H1("fhRichRingsTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
-
- fHM->H1("fhTofDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
- fHM->H1("fhTofDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
- fHM->H1("fhTofDigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
+ fHM->H1("fhT0DigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
+ fHM->H1("fhT0DigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
+ fHM->H1("fhT0DigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
- fHM->H1("fhStsDigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
- fHM->H1("fhStsDigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
- fHM->H1("fhStsDigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
+ fDigiHitsInitialized = true;
+ }
+ }
+ } // if (fDigiHitsInitialized == false)
+
+ if (fDigiHitsInitialized == true) {
- fHM->H1("fhT0DigisTimeLog")->GetXaxis()->SetLimits(minTime, minTime+dT);
- fHM->H1("fhT0DigisTimeLogZoom")->GetXaxis()->SetLimits(minTime, minTime+dTZoom1);
- fHM->H1("fhT0DigisTimeLogZoom2")->GetXaxis()->SetLimits(minTime, minTime + dTZoom2);
- }
+ int nofRichDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
+ fHM->H1("fhNofRichDigisInTimeslice")->Fill(nofRichDigis);
+ for (int i = 0; i < nofRichDigis; i++) {
+ const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(i);
+ fHM->H1("fhRichDigisTimeLog")->Fill(digi->GetTime() );
+ fHM->H1("fhRichDigisTimeLogZoom")->Fill(digi->GetTime() );
+ fHM->H1("fhRichDigisTimeLogZoom2")->Fill(digi->GetTime());
+ }
- {
- for (int i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
- const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
- fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
+ int nofRichRings= fRichRings->GetEntries();
+ for (int i = 0; i < nofRichRings; i++)
+ {
+ CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(i));
+ fHM->H1("fhRichRingsTimeLog")->Fill(ring->GetTime());
+ fHM->H1("fhRichRingsTimeLogZoom")->Fill(ring->GetTime());
+ fHM->H1("fhRichRingsTimeLogZoom2")->Fill(ring->GetTime());
}
- }
- int nofRichDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
- fHM->H1("fhNofRichDigisInTimeslice")->Fill(nofRichDigis);
- for (int i = 0; i < nofRichDigis; i++) {
- const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(i);
- fHM->H1("fhRichDigisTimeLog")->Fill(digi->GetTime() );
- fHM->H1("fhRichDigisTimeLogZoom")->Fill(digi->GetTime() );
- fHM->H1("fhRichDigisTimeLogZoom2")->Fill(digi->GetTime());
- }
+ int nofTofDigis = fDigiMan->GetNofDigis(ECbmModuleId::kTof);
+ for (int i = 0; i < nofTofDigis; i++) {
+ const CbmTofDigi* digi = fDigiMan->Get<CbmTofDigi>(i);
+ fHM->H1("fhTofDigisTimeLog")->Fill(digi->GetTime() );
+ fHM->H1("fhTofDigisTimeLogZoom")->Fill(digi->GetTime() );
+ fHM->H1("fhTofDigisTimeLogZoom2")->Fill(digi->GetTime());
+ }
- int nofRichRings= fRichRings->GetEntries();
- for (int i = 0; i < nofRichRings; i++)
- {
- CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(i));
- fHM->H1("fhRichRingsTimeLog")->Fill(ring->GetTime());
- fHM->H1("fhRichRingsTimeLogZoom")->Fill(ring->GetTime());
- fHM->H1("fhRichRingsTimeLogZoom2")->Fill(ring->GetTime());
- }
-
- int nofTofDigis = fDigiMan->GetNofDigis(ECbmModuleId::kTof);
- for (int i = 0; i < nofTofDigis; i++) {
- const CbmTofDigi* digi = fDigiMan->Get<CbmTofDigi>(i);
- fHM->H1("fhTofDigisTimeLog")->Fill(digi->GetTime() );
- fHM->H1("fhTofDigisTimeLogZoom")->Fill(digi->GetTime() );
- fHM->H1("fhTofDigisTimeLogZoom2")->Fill(digi->GetTime());
- }
+ if (fDigiMan->IsPresent(ECbmModuleId::kSts)){
+ int nofStsDigis = fDigiMan->GetNofDigis(ECbmModuleId::kSts);
+ for (int i = 0; i < nofStsDigis; i++) {
+ const CbmStsDigi* digi = fDigiMan->Get<CbmStsDigi>(i);
+ fHM->H1("fhStsDigisTimeLog")->Fill(digi->GetTime() );
+ fHM->H1("fhStsDigisTimeLogZoom")->Fill(digi->GetTime() );
+ fHM->H1("fhStsDigisTimeLogZoom2")->Fill(digi->GetTime());
+ }
+ }
- if (fDigiMan->IsPresent(ECbmModuleId::kSts)){
- int nofStsDigis = fDigiMan->GetNofDigis(ECbmModuleId::kSts);
- for (int i = 0; i < nofStsDigis; i++) {
- const CbmStsDigi* digi = fDigiMan->Get<CbmStsDigi>(i);
- fHM->H1("fhStsDigisTimeLog")->Fill(digi->GetTime() );
- fHM->H1("fhStsDigisTimeLogZoom")->Fill(digi->GetTime() );
- fHM->H1("fhStsDigisTimeLogZoom2")->Fill(digi->GetTime());
+ {
+ Int_t nrT0Digis = 0;
+ if ( fT0Digis ) nrT0Digis = fT0Digis->size();
+ //else if ( fT0DigiArr ) nrT0Digis = fT0DigiArr->GetEntriesFast();
+ LOG(debug) << "T0Digis: " << nrT0Digis;
+
+ for (Int_t iT0 = 0; iT0 < nrT0Digis; ++iT0) {
+ const CbmTofDigi* T0Digi = nullptr;
+ if ( fT0Digis ) T0Digi = &(fT0Digis->at(iT0));
+ //else if ( fT0DigiArr ) T0Digi = dynamic_cast<const CbmTofDigi*>(fT0DigiArr->At(iT0));
+ assert (T0Digi);
+ fHM->H1("fhT0DigisTimeLog")->Fill(T0Digi->GetTime() );
+ fHM->H1("fhT0DigisTimeLogZoom")->Fill(T0Digi->GetTime() );
+ fHM->H1("fhT0DigisTimeLogZoom2")->Fill(T0Digi->GetTime());
+ }
}
}
-// int nofT0Digis = fT0Digis->GetEntries();
-// for (int i = 0; i < nofT0Digis; i++) {
-// CbmDigi* digi = static_cast<CbmDigi*>(fT0Digis->At(i));
-// fHM->H1("fhT0DigisTimeLog")->Fill(digi->GetTime() );
-// fHM->H1("fhT0DigisTimeLogZoom")->Fill(digi->GetTime() );
-// fHM->H1("fhT0DigisTimeLogZoom2")->Fill(digi->GetTime());
-// }
-
-
int nofRichHits = fRichHits->GetEntries();
fHM->H1("fhNofRichHitsInTimeslice")->Fill(nofRichHits);
fHM->H1("fhHitsInTimeslice")->Fill(fEventNum,nofRichHits);
+
+ {
+ int nofRichDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
+ for (int i = 0; i < nofRichDigis; i++) {
+ const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
+ fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
+
+ if (i>0){
+ const CbmRichDigi* richDigi_prev = fDigiMan->Get<CbmRichDigi>(i-1);
+ fHM->H1("fhRichDigisConsecTime")->Fill(richDigi->GetTime()-richDigi_prev->GetTime());
+ if (doToT(richDigi) && doToT(richDigi_prev)) fHM->H1("fhRichDigisConsecTimeTOT")->Fill(richDigi->GetTime()-richDigi_prev->GetTime());
+ }
+
+ //fHM->H2("fhRichDigiPixelRate")->Fill(richDigi->GetX(),richDigi->GetY());
+ }
+ }
+
+ Double_t zPos_tmp = 0.;
for (int iH = 0; iH < nofRichHits; iH++) {
CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
fHM->H2("fhRichHitXY")->Fill(richHit->GetX(), richHit->GetY());
fHM->H1("fhRichHitToT")->Fill(richHit->GetToT());
- //printf("HitToT: %f \n", richHit->GetToT());
- }
-
+ zPos_tmp = ((zPos_tmp*iH)+richHit->GetZ())/ (iH+1);
+ //printf("HitToT: %f \n", richHit->GetToT());
+ }
+ if (nofRichHits > 0) RichZPos = zPos_tmp;
+ //std::cout<<"[INFO] Mean Z-Position of RICH Hits: "<<RichZPos<<std::endl;
//CBMEVENT
auto fNCbmEvent = fCbmEvent->GetEntriesFast();
for (int i=0;i<fNCbmEvent;i++){
fHM->H1("fhNofCbmEvents")->Fill(1);
CbmEvent* ev = static_cast<CbmEvent*>(fCbmEvent->At(i));
+
+ if (fTriggerRichHits != 0 && (ev->GetNofData(ECbmDataType::kRichHit) < fTriggerRichHits)) continue;
+ if (fTriggerTofHits != 0 && (ev->GetNofData(ECbmDataType::kTofHit) < fTriggerTofHits )) continue;
+
+ //if (ev->GetNofData(ECbmDataType::kTofHit) > (fTriggerTofHits+10) ) continue;
+
std::vector<int> ringIndx;
std::vector<int> evRichHitIndx;
+ std::array<uint32_t,36> pmtHits;
+ for (auto& a : pmtHits) a = 0;
+ fEventPnt = ev;
+ fCbmEventStartTime = fEventPnt->GetStartTime();
+
+
+ // Scan Event to find first Digi that triggered.
+// std::cout<<"Sts Digis:"<< ev->GetNofData(kStsDigi)<<std::endl;
+// std::cout<<"Much Digis:"<< ev->GetNofData(kMuchDigi)<<std::endl;
+// std::cout<<"Tof Digis:"<< ev->GetNofData(kTofDigi)<<std::endl;
+// std::cout<<"Rich Digis:"<< ev->GetNofData(kRichDigi)<<std::endl;
+// std::cout<<"Psd Digis:"<< ev->GetNofData(kPsdDigi)<<std::endl;
+// unsigned int flagRich = 0;
+ Double_t startTime = std::numeric_limits<Double_t>::max();
for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
evRichHitIndx.push_back(iRichHit);
CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
- //std::cout<<"\t\t * "<<i<<". Event, Hit "<< j <<": "<< iRichHit <<"\t " << std::fixed << std::setprecision(5) << richHit->GetTime() <<std::endl;
+ fHM->H1("fhRichHitToTEvent")->Fill(richHit->GetToT());
+ fHM->H2("fhRichHitXYEvent")->Fill(richHit->GetX(), richHit->GetY());
+ //Blob finder
+ uint32_t pmtId = (((richHit->GetAddress())>>20)&0xF) + (((richHit->GetAddress())>>24)&0xF)*9;
+ pmtHits[pmtId]++;
+ //std::cout<<"\t\t * "<<i<<". Event, Hit "<< j <<": "<< iRichHit <<"\t " << std::fixed << std::setprecision(5) << richHit->GetTime() <<std::endl;
+ if (richHit->GetTime() < startTime ) {startTime = richHit->GetTime();/*flagRich = 1;*/}
int nofRichRings2= fRichRings->GetEntries();
for (int l = 0; l < nofRichRings2; l++)
{
@@ -353,7 +568,9 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
if (RingHitIndx == iRichHit) {
Bool_t used = false;
for (auto check : ringIndx) {
- if (check == l){ used = true; break;}
+ if (check == l){ used = true;
+ break;
+ }
}
if (used == false) ringIndx.push_back(l);
break;
@@ -364,17 +581,31 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
for (int k=0;k<ev->GetNofData(ECbmDataType::kTofHit);k++){
auto iTofHit = ev->GetIndex(ECbmDataType::kTofHit, k);
CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(iTofHit));
+ if (tofHit->GetTime() < startTime ) {startTime = tofHit->GetTime();/* flagRich = 0;*/}
if (tofHit->GetZ() < 2.) continue; // Cut T0 away!
fHM->H2("fhTofRichX")->Fill(richHit->GetX(),tofHit->GetX());
+
+ //Befor tof_v20b
+ /*if (tofHit->GetZ()> 230. && tofHit->GetZ() < 250) fHM->H2("fhTofRichX_stack1")->Fill(richHit->GetX(),tofHit->GetX());
+ if (tofHit->GetZ()> 250. && tofHit->GetZ() < 265) fHM->H2("fhTofRichX_stack2")->Fill(richHit->GetX(),tofHit->GetX());
+ if (tofHit->GetZ()> 265. && tofHit->GetZ() < 285) fHM->H2("fhTofRichX_stack3")->Fill(richHit->GetX(),tofHit->GetX());
+ */
+ if (tofHit->GetZ() > 230. && tofHit->GetZ() < 255) fHM->H2("fhTofRichX_stack1")->Fill(richHit->GetX(),tofHit->GetX());
+ if (tofHit->GetZ() >= 255. && tofHit->GetZ() < 272) fHM->H2("fhTofRichX_stack2")->Fill(richHit->GetX(),tofHit->GetX());
+ if (tofHit->GetZ() >= 272. && tofHit->GetZ() < 290) fHM->H2("fhTofRichX_stack3")->Fill(richHit->GetX(),tofHit->GetX());
+
fHM->H2("fhTofRichY")->Fill(richHit->GetY(),tofHit->GetY());
fHM->H2("fhTofRichX_zoomed")->Fill(richHit->GetX(),tofHit->GetX());
- fHM->H2("fhTofRichY_zoomed")->Fill(richHit->GetY(),tofHit->GetY());
+ fHM->H2("fhTofRichY_zoomed")->Fill(richHit->GetY(),tofHit->GetY());
+
+ fHM->H2("fhTofHitsXY")->Fill(tofHit->GetX(),tofHit->GetY());
}
}
-
+ //std::cout<<"Diff in StartTime DigiToHit: "<< startTime - fCbmEventStartTime << "\t" <<flagRich<<std::endl;
+ fCbmEventStartTime = startTime;
- DrawEvent(ev, ringIndx, 1);
+ fSeDisplay->DrawEvent(ev, ringIndx, 1);
// std::cout<<DrawCbmEvent<<std::endl;
// for (int j=0;j<ev->GetNofData(ECbmDataType::kRichDigi);j++){
@@ -387,14 +618,38 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
//*** Tracks in CbmEvent -> Here Trigger is on! This is seen in Track Position
//std::cout<<"TRACKS in CbmEvent: "<<ev->GetNofData(ECbmDataType::kTofTrack)<<std::endl;
- for (int j=0;j<ev->GetNofData(ECbmDataType::kTofTrack);j++){
+ auto noftofTracks = ev->GetNofData(ECbmDataType::kTofTrack);
+ fHM->H1("fhTofTracksPerEvent")->Fill(noftofTracks);
+
+ for (int j=0;j<noftofTracks;j++) {
auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, j);
CbmTofTracklet* tTrack = static_cast<CbmTofTracklet*>(fTofTracks->At(iTofTrack));
-
+ if (tTrack == nullptr) continue;
+
+ fHM->H1("fhTracksWithRings")->Fill(0);
fHM->H2("fhTofTracksXY")->Fill(tTrack->GetFitX(RichZPos),tTrack->GetFitY(RichZPos));
+ fHM->H2("fhTofTracksXY_Target")->Fill(tTrack->GetFitX(0.),tTrack->GetFitY(0.));
fHM->H1("fhNofTofTracks")->Fill(0.5); // 1: All 2: left; 3: right; 4: RICH
fHM->H1("fhTofBetaAllFullAcc")->Fill(getBeta(tTrack));
//std::cout<<"beta Track "<< j <<": "<<getBeta(tTrack)<<std::endl;
+
+ if (tTrack->GetFitX(RichZPos) > -10. && tTrack->GetFitX(RichZPos) < +10. && tTrack->GetFitY(RichZPos) > -25. && tTrack->GetFitY(RichZPos) < +25. && isOnTarget(tTrack)){
+ //Track in RICH
+ fTracksinRich++;
+ Int_t goodHit = 0;
+ for (int k=0;k<ev->GetNofData(ECbmDataType::kRichHit);k++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, k);
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
+ if (richHit == nullptr) continue;
+ if (std::fabs(richHit->GetY()-tTrack->GetFitY(RichZPos)) < 5. && std::fabs(richHit->GetX()-tTrack->GetFitX(RichZPos)) < 9.) goodHit++;
+ }
+ if (goodHit > 0) fTracksinRichWithRichHits[0]++;
+ if (goodHit > 5) fTracksinRichWithRichHits[1]++;
+ if (goodHit > 10) fTracksinRichWithRichHits[2]++;
+ if (goodHit > 15) fTracksinRichWithRichHits[3]++;
+
+ }
+
if (isAccmRICH(tTrack)){
if (RestrictToFullAcc(tTrack)){
fHM->H2("fhTofTracksXYRICH")->Fill(tTrack->GetFitX(RichZPos),tTrack->GetFitY(RichZPos)); // projected in RICH Plane
@@ -406,7 +661,36 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
if (ringIndx.size() == 0 && evRichHitIndx.size() > 0) fHM->H1("fhTofBetaTracksWithHitsNoRing")->Fill(getBeta(tTrack)); // no Ring in CbmEvent found
if (evRichHitIndx.size() > 0)fHM->H1("fhTofBetaTracksWithHits")->Fill(getBeta(tTrack));
if (ringIndx.size() == 0)fHM->H1("fhTofBetaTracksNoRing")->Fill(getBeta(tTrack));
- if (FindClosestRing(tTrack,ringIndx).first > -1) fHM->H1("fhTofBetaTrackswithClosestRingInRange")->Fill(getBeta(tTrack));
+
+ if ((tTrack->GetNofHits() == 4 && (getBeta(tTrack) > 0.90 && getBeta(tTrack) < 1.10))){
+ //tracks after cut
+ Double_t trackXpos = tTrack->GetFitX(RichZPos);
+ Double_t trackYpos = tTrack->GetFitY(RichZPos);
+
+ if ( trackXpos > -8 && trackXpos < 13 && trackYpos > -21 && trackYpos < 24 ){
+ if (!(trackXpos > -8 && trackXpos < -3 && trackYpos > 5 && trackYpos < 7.5) &&
+ !(trackXpos > 7.8 && trackXpos < 13 && trackYpos > 5 && trackYpos < 7.5) &&
+ !(trackXpos > -8 && trackXpos < 2 && trackYpos > 21 && trackYpos < 24 ) &&
+ !(trackXpos > 7.8 && trackXpos < 13 && trackYpos > 21 && trackYpos < 24 ) &&
+ !(trackXpos > 2.2 && trackXpos < 13 && trackYpos > 21 && trackYpos < 16 )
+ ){
+ fHM->H2("fhTofTracksXYRICH_Accectance")->Fill(trackXpos,trackYpos);
+ fHM->H1("fhTracksWithRings")->Fill(2);
+ if (ringIndx.size() > 0) fHM->H1("fhTracksWithRings")->Fill(4);
+ if (FindClosestRing(tTrack,ringIndx).first > -1) fHM->H1("fhTracksWithRings")->Fill(3);
+ }
+ }
+
+ // filter on Rich Acc.
+ // filter on
+ }
+
+ std::pair<int, double> closeRing = FindClosestRing(tTrack,ringIndx);
+ if (closeRing.first > -1) {
+ fHM->H1("fhTofBetaTrackswithClosestRingInRange")->Fill(getBeta(tTrack));
+ fHM->H1("fhRingTrackDistance")->Fill(closeRing.second);
+ fHM->H1("fhTracksWithRings")->Fill(1);
+ }
//if (ringIndx.size() > 0)fHM->H1("fhTofBetaTrackswithRing")->Fill(getBeta(tTrack)); // ring is somewehere in Acc
fHM->H1("fhTofBetaAll")->Fill(getBeta(tTrack));
}
@@ -427,7 +711,8 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
if (tofHit->GetZ() < 2.) continue; // Cut T0 away!
fHM->H1("fhTofHitsZ")->Fill(tofHit->GetZ());
- fHM->H1("fhTofHitsXZ")->Fill(tofHit->GetZ(),tofHit->GetX());
+ fHM->H2("fhTofHitsXZ")->Fill(tofHit->GetZ(),tofHit->GetX());
+ //fHM->H2("fhTofHitsXY")->Fill(tofHit->GetX(),tofHit->GetY());
fHM->H3("fhTofXYZ")->Fill(tofHit->GetX(),tofHit->GetZ(),tofHit->GetY());
if (ringIndx.size() > 0) {
@@ -437,6 +722,14 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(ringIndx[rings]));
fHM->H2("fhTofRichRingX")->Fill(ring->GetCenterX(),tofHit->GetX());
fHM->H2("fhTofRichRingY")->Fill(ring->GetCenterY(),tofHit->GetY());
+ for (int k = 0; k < ring->GetNofHits(); k++){
+ Int_t hitInd = ring->GetHit(k);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+ fHM->H2("fhTofHitRichRingHitX")->Fill(hit->GetX(),tofHit->GetX());
+ fHM->H2("fhTofHitRichRingHitY")->Fill(hit->GetY(),tofHit->GetY());
+ }
+ //fhTofRichRingHitX here loop over hits in ring.
}
}
}
@@ -445,34 +738,156 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
std::vector<CbmTofTracklet*> TracksOfEvnt;
//TRACKS
auto nofTofTracks = ev->GetNofData(ECbmDataType::kTofTrack);
- fHM->H1("fhTofTracksPerEvent")->Fill(nofTofTracks);
+ fHM->H1("fhTofTracksPerRichEvent")->Fill(nofTofTracks);
for (int j=0;j<nofTofTracks;j++){
auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, j);
CbmTofTracklet *track = static_cast<CbmTofTracklet *>(fTofTracks->At(iTofTrack));
- fHM->H2("fhTofTrackRichRingXY")->Fill(track->GetFitX(RichZPos),track->GetFitY(RichZPos));
+ if (track == nullptr) continue;
+ if (track->GetNofHits() <= 3) continue;
TracksOfEvnt.emplace_back(track);
+ if (!isOnTarget(track)) continue;
+ fHM->H2("fhTofTrackRichRingXY")->Fill(track->GetFitX(RichZPos),track->GetFitY(RichZPos));
+ for (int k=0;k<ev->GetNofData(ECbmDataType::kRichHit);k++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, k);
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
+ if (richHit == nullptr) continue;
+ fHM->H2("fhTofTrackRichHitX")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+
+ //cuts:
+ // >3 hits /track
+ // if (track->GetNofHits() > 3) {
+ fHM->H2("fhTofTrackRichHitX_cuts")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_cuts")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ if (getBeta(track)> 0.90 && getBeta(track)< 1.10) {
+ fHM->H2("fhTofTrackRichHitX_oBetacuts")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_oBetacuts")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ double deltatime = richHit->GetTime() - track->GetTime();
+ fHM->H1("fhTofTrackRichHitTime_oBetacuts")->Fill(deltatime);
+ if (deltatime > -10 && deltatime < +40) {
+ if (track->GetNofHits() == 4) {
+ fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ }
+
+ if (track->GetNofHits() > 4){
+ fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_4")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_4")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ } else if (track->GetNofHits() > 6){
+ fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_6")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_6")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ } else if (track->GetNofHits() > 8){
+ fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_8")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_8")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ } else if (track->GetNofHits() > 10){
+ fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_10")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_10")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ }
+
+ if (track->GetNofHits() == 4) {
+ for (int l=0; l<track->GetNofHits();++l){
+ auto hitIndex = track->GetHitIndex(l);
+ auto iTofHit = ev->GetIndex(ECbmDataType::kTofHit, hitIndex);
+ const CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(iTofHit));
+ if (tofHit->GetZ() < 2.) continue; // Cut T0 away!
+ fHM->H2("fhTofTrackHitRichHitX_oBetacuts_dtime")->Fill(richHit->GetX(),tofHit->GetX());
+ fHM->H2("fhTofTrackHitRichHitY_oBetacuts_dtime")->Fill(richHit->GetY(),tofHit->GetY());
+ }
+ }
+ }
+ } else {
+ fHM->H2("fhTofTrackRichHitX_uBetacuts")->Fill(richHit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichHitY_uBetacuts")->Fill(richHit->GetY(),track->GetFitY(RichZPos));
+ }
+ //}
+ }
+
}
//rings
+ //std::cout<< "Found Rings in CbmEvent: "<< ringIndx.size() <<std::endl;
+ fRingsWithTrack[1] += ringIndx.size();
+ fRingsWithTrack[2] += nofTofTracks;
+ fRingsWithTrack[5] += ringIndx.size()*nofTofTracks;
+
for (unsigned int rings=0; rings < ringIndx.size();rings++) {
CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(ringIndx[rings]));
+ if (nullptr == ring) continue;
+
+ fHM->H1("fhRichRingChi2")->Fill(ring->GetChi2());
+
+ if (!cutRadius(ring)) continue;
+ fRingsWithTrack[3]++; // Rings After Cut
+ //DEFAULT: //DrawRing(ring,TracksOfEvnt,true); // here normally; now only Rings with high BETA
+
+ fHM->H1("fhRichRingBeta")->Fill(getBeta(ring));
- //DrawRing(ring,TracksOfEvnt,true); // here normally; now only Rings with high BETA
auto clTrack = FindClosestTrack(ring,TracksOfEvnt); // has no cut on distance
analyseRing(ring,ev,clTrack);
+
+ for (int j=0;j<nofTofTracks;j++){
+ auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, j);
+ CbmTofTracklet *track = static_cast<CbmTofTracklet *>(fTofTracks->At(iTofTrack));
+ if (track == nullptr) continue;
+ if (!(track->GetNofHits() == 4 && (getBeta(track) > 0.90 && getBeta(track) < 1.10))) continue;
+ if (rings == 0) {
+ fRingsWithTrack[4]++;
+ if (ring->GetChi2()<4.) fSeDsply_TR->DrawEvent(ev, ringIndx, 1); //Some will be drawn double, but for now ok
+ }
+ fHM->H2("fhTofTrackRichRingX")->Fill(ring->GetCenterX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichRingY")->Fill(ring->GetCenterY(),track->GetFitY(RichZPos));
+
+ const double xDist = (track->GetFitX(RichZPos) - ring->GetCenterX());
+ const double yDist = (track->GetFitY(RichZPos) - ring->GetCenterY());
+ const double rDist = std::sqrt(xDist*xDist + yDist*yDist);
+ fHM->H1("fhTrackRingDistance_corr")->Fill(rDist);
+
+ fRingsWithTrack[0]++;
+
+ for (int k = 0; k < ring->GetNofHits(); k++){
+ Int_t hitInd = ring->GetHit(k);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+ fHM->H2("fhTofTrackRichRingHitX")->Fill(hit->GetX(),track->GetFitX(RichZPos));
+ fHM->H2("fhTofTrackRichRingHitY")->Fill(hit->GetY(),track->GetFitY(RichZPos));
+ }
+ }
+
if (clTrack.first > -1){
// first: TrackIndex; second: Distance;
//FIXME
//if (getBeta(TracksOfEvnt[clTrack.first]) > 0.9) DrawRing(ring,TracksOfEvnt,true);
fHM->H1("fhTrackRingDistance")->Fill(clTrack.second);
+ CbmTofTracklet *track = TracksOfEvnt[clTrack.first];//static_cast<CbmTofTracklet *>(fTofTracks->At(clTrack.first));
+ if (track == nullptr) continue;
+ if (isOnTarget(track)){
+ fHM->H1("fhTrackRingDistanceOnTarget")->Fill(clTrack.second);
+
+ } else {
+ fHM->H1("fhTrackRingDistanceOffTarget")->Fill(clTrack.second);
+
+ }
fHM->H2("fhTrackRingDistanceVSRingradius")->Fill(clTrack.second,ring->GetRadius());
- fHM->H2("fhClosTrackRingX")->Fill(ring->GetCenterX(),TracksOfEvnt[clTrack.first]->GetFitX(RichZPos));
- fHM->H2("fhClosTrackRingY")->Fill(ring->GetCenterY(),TracksOfEvnt[clTrack.first]->GetFitY(RichZPos));
- fHM->H2("fhTofClosTrackRichRingXY")->Fill(TracksOfEvnt[clTrack.first]->GetFitX(RichZPos),TracksOfEvnt[clTrack.first]->GetFitY(RichZPos));
- if (cutDistance(clTrack) && cutRadius(ring)) { //Good Ring
- fHM->H1("fhTofBetaRing")->Fill(getBeta(TracksOfEvnt[clTrack.first]));
- fHM->H2("fhTofBetaRingDist")->Fill(getBeta(TracksOfEvnt[clTrack.first]),clTrack.second);
- fHM->H2("fhTofBetaVsRadius")->Fill(getBeta(TracksOfEvnt[clTrack.first]),ring->GetRadius());
+ fHM->H2("fhTrackRingDistanceVSRingChi2")->Fill(clTrack.second,ring->GetChi2());
+ //if ( (clTrack.second < 20.0 )) {
+ fHM->H2("fhRichRingXY_goodTrack")->Fill(ring->GetCenterX(),ring->GetCenterY());
+ fHM->H2("fhRichRing_goodTrackXY")->Fill(track->GetFitX(RichZPos),track->GetFitY(RichZPos));
+ //}
+ fHM->H2("fhClosTrackRingX")->Fill(ring->GetCenterX(),track->GetFitX(RichZPos));
+ fHM->H2("fhClosTrackRingY")->Fill(ring->GetCenterY(),track->GetFitY(RichZPos));
+ fHM->H2("fhTofClosTrackRichRingXY")->Fill(track->GetFitX(RichZPos),track->GetFitY(RichZPos));
+ if ((clTrack.second < (ring->GetRadius()*1.2))) { //Good Ring
+ fHM->H2("fhGoodRingsXY_TargetPos")->Fill(track->GetFitX(0.),track->GetFitY(0.));
+ fHM->H1("fhRichRingChi2_goodRing")->Fill(ring->GetChi2());
+ fHM->H2("fhTrackRingDistanceVSRingChi2_goodRing")->Fill(clTrack.second,ring->GetChi2());
+ fHM->H1("fhTofBetaRing")->Fill(getBeta(track));
+ fHM->H2("fhTofBetaRingDist")->Fill(getBeta(track),clTrack.second);
+ fHM->H2("fhTofBetaVsRadius")->Fill(getBeta(track),ring->GetRadius());
+ //Ring properties of "Good rings"
+ fHM->H1("fhRichRingRadius_goodRing")->Fill(ring->GetRadius());
+ fHM->H1("fhNofHitsInRing_goodRing")->Fill(ring->GetNofHits());
+ fHM->H1("fhRichRingBeta_GoodRing")->Fill(getBeta(ring));
+
}
}
@@ -491,7 +906,33 @@ void CbmRichMCbmQaReal::Exec(Option_t* /*option*/)
// std::cout<< "TOF Hits:\t" << ev->GetNofData(ECbmDataType::kTofHit) <<std::endl;
// std::cout<< "TOF Tracks:\t" << ev->GetNofData(ECbmDataType::kTofTrack) <<std::endl;
-
+
+ //Select and Analyse Blobs
+ uint32_t blob = 0;
+ for (unsigned int j=0; j< pmtHits.size(); ++j ){
+ if (pmtHits[j] > 30) {
+ blob++;
+ // Blob found
+ double xBlob = -7.95 + (int((j/9))*5.3) + fXOffsetHisto;
+ double yBlob = 21.2 - ((j%9)*5.3);
+ //std::cout<<"BLOB X:" << i << " "<< xBlob << std::endl;
+ for (int k=0;k<ev->GetNofData(ECbmDataType::kTofTrack);k++){
+ auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, k);
+ CbmTofTracklet *track = static_cast<CbmTofTracklet *>(fTofTracks->At(iTofTrack));
+ if (track == nullptr) continue;
+ double xBlobTrack = track->GetFitX(RichZPos) - xBlob;
+ double yBlobTrack = track->GetFitY(RichZPos) - yBlob;
+ fHM->H1("fhBlobTrackDistX")->Fill(std::fabs(xBlobTrack));
+ fHM->H1("fhBlobTrackDistY")->Fill(std::fabs(yBlobTrack));
+ fHM->H1("fhBlobTrackDist" )->Fill(std::sqrt(xBlobTrack*xBlobTrack + yBlobTrack*yBlobTrack));
+ }
+ }
+ }
+ if (blob > 0) {
+ fHM->H1("fhNofBlobEvents")->Fill(1);
+ fHM->H1("fhNofBlobsInEvent")->Fill(blob);
+ }
+
} //End CbmEvent loop
@@ -510,6 +951,8 @@ void CbmRichMCbmQaReal::RichRings()
fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
+ fHM->H2("fhRichRingRadiusY")->Fill(ring->GetRadius(),ring->GetCenterY());
+ fHM->H2("fhRichHitsRingRadius")->Fill(ring->GetNofHits(),ring->GetRadius());
}
}
@@ -523,14 +966,16 @@ std::pair<int, double> CbmRichMCbmQaReal::FindClosestTrack(const CbmRichRing* ri
for (unsigned int indx = 0; indx < track.size(); ++indx){
//Calc if Track is in Ring (+20% )
+ if (track[indx]->GetNofHits() <= 3) continue;
const double xDist = (track[indx]->GetFitX(RichZPos) - ringX);
const double yDist = (track[indx]->GetFitY(RichZPos) - ringY);
const double rDist = std::sqrt(xDist*xDist + yDist*yDist);
const double RadiusFactor = 1.2; // Factor of how big radius of acceptance should
- if (rDist < ring->GetRadius()*RadiusFactor) {
- std::cout<<"Track in defined Ring range ("<<ring->GetRadius()*RadiusFactor<<"cm) (RingRadius: "<<ring->GetRadius()<<"cm). ";
+ if (rDist < ring->GetRadius()*RadiusFactor ) {
+ //std::cout<<"Track in defined Ring range ("<<ring->GetRadius()*RadiusFactor<<"cm) (RingRadius: "<<ring->GetRadius()<<"cm). ";
}
+ //if (rDist < ring->GetRadius()*RadiusFactor && ring->GetRadius() > 2. && ring->GetRadius() < 4.2) {
if (indx== 0) {
closDist = rDist;
closTrack = indx;
@@ -540,8 +985,9 @@ std::pair<int, double> CbmRichMCbmQaReal::FindClosestTrack(const CbmRichRing* ri
closTrack = indx;
}
}
+ //}
}
-
+ //if (closTrack > -1 ) std::cout<< "closestTrack to ring "<< ring <<": "<<closTrack<<" " << static_cast<CbmTofTracklet *>(fTofTracks->At(closTrack)) <<std::endl;
std::pair<int, double> p;
p.first = closTrack;
p.second = closDist;
@@ -569,7 +1015,7 @@ std::pair<int, double> CbmRichMCbmQaReal::FindClosestRing(CbmTofTracklet* track,
const double rDist = std::sqrt(xDist*xDist + yDist*yDist);
const double RadiusFactor = 1.2; // Factor of how big radius of acceptance should
- if (rDist < ring->GetRadius()*RadiusFactor) {
+ if (rDist < ring->GetRadius()*RadiusFactor && cutRadius(ring) ) {
//std::cout<<"Track in defined Ring range ("<<ring->GetRadius()*RadiusFactor<<"cm) (RingRadius: "<<ring->GetRadius()<<"cm). ";
if (indx== 0) {
@@ -583,6 +1029,8 @@ std::pair<int, double> CbmRichMCbmQaReal::FindClosestRing(CbmTofTracklet* track,
}
}
}
+
+ //if (closTrack > -1 ) std::cout<< "closestRing to Track "<< track <<": "<<closTrack <<" "<< static_cast<CbmRichRing *>(fRichRings->At(ringIndx[closTrack])) <<std::endl;
std::pair<int, double> p;
p.first = closTrack;
p.second = closDist;
@@ -592,7 +1040,7 @@ std::pair<int, double> CbmRichMCbmQaReal::FindClosestRing(CbmTofTracklet* track,
void CbmRichMCbmQaReal::DrawHist()
{
- cout.precision(4);
+ cout.precision(4);
//SetDefaultDrawStyle();
double nofEvents = fHM->H1("fhNofCbmEvents")->GetEntries();
@@ -617,7 +1065,36 @@ void CbmRichMCbmQaReal::DrawHist()
fHM->CreateCanvas("HitsInTimeslice","HitsInTimeslice", 600 , 600);
DrawH1(fHM->H1("fhHitsInTimeslice"));
}
+
+ {
+ fHM->CreateCanvas("RichDigisConsecTime","RichDigisConsecTime", 600 , 600);
+ DrawH1(fHM->H1("fhRichDigisConsecTime"), kLinear, kLog);
+ }
+
+ {
+ fHM->CreateCanvas("RichDigisConsecTimeTOT","RichDigisConsecTimeTOT", 600 , 600);
+ DrawH1(fHM->H1("fhRichDigisConsecTimeTOT"), kLinear, kLog);
+ }
+
+
+ /*{ // NO position info for Digis :(
+ fHM->CreateCanvas("Pixelrate","Pixelrate", 600 , 600);
+ double time = nofEvents * 0.0128; //seconds
+ fHM->H1("fhRichDigiPixelRate")->Scale(1./time);
+ DrawH1(fHM->H1("fhRichDigiPixelRate"));
+ }*/
+
+ {
+ TCanvas* c = fHM->CreateCanvas("RichRingXY_goodTrack","RichRingXY_goodTrack", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhRichRingXY_goodTrack"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhRichRing_goodTrackXY"));
+ }
+
+
{
TCanvas* c = fHM->CreateCanvas("rich_mcbm_nofObjectsInTimeslice","rich_mcbm_nofObjectsInTimeslice", 1500 , 500);
c->Divide(3,1);
@@ -663,15 +1140,15 @@ void CbmRichMCbmQaReal::DrawHist()
TCanvas* c = fHM->CreateCanvas("rich_mcbm_richDigisTimeLog","rich_mcbm_richDigisTimeLog", 1200 , 1200);
c->Divide(1,2);
c->cd(1);
- DrawH1({fHM->H1("fhStsDigisTimeLog"), fHM->H1("fhTofDigisTimeLog"), fHM->H1("fhT0DigisTimeLog"), fHM->H1("fhRichDigisTimeLog")},
- {"STS", "TOF", "T0", "RICH"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
+ DrawH1({fHM->H1("fhRichDigisTimeLog"), fHM->H1("fhTofDigisTimeLog"), fHM->H1("fhT0DigisTimeLog"),fHM->H1("fhStsDigisTimeLog") },
+ {"RICH", "TOF", "T0", "STS"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.10);
fHM->H1("fhStsDigisTimeLog")->GetYaxis()->SetTitleOffset(0.7);
fHM->H1("fhStsDigisTimeLog")->SetMinimum(0.9);
c->cd(2);
- DrawH1({fHM->H1("fhStsDigisTimeLogZoom"), fHM->H1("fhTofDigisTimeLogZoom"), fHM->H1("fhT0DigisTimeLogZoom"), fHM->H1("fhRichDigisTimeLogZoom")},
- {"STS", "TOF", "T0", "RICH"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
+ DrawH1({fHM->H1("fhRichDigisTimeLogZoom"), fHM->H1("fhTofDigisTimeLogZoom"), fHM->H1("fhT0DigisTimeLogZoom"), fHM->H1("fhStsDigisTimeLogZoom")},
+ {"RICH", "TOF", "T0", "STS"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.1);
fHM->H1("fhStsDigisTimeLogZoom")->GetYaxis()->SetTitleOffset(0.7);
@@ -680,13 +1157,25 @@ void CbmRichMCbmQaReal::DrawHist()
{
fHM->CreateCanvas("rich_mcbm_richDigisTimeLog2", "rich_mcbm_richDigisTimeLog2", 1200, 600);
- DrawH1({fHM->H1("fhStsDigisTimeLogZoom2"), fHM->H1("fhTofDigisTimeLogZoom2"), fHM->H1("fhT0DigisTimeLogZoom2"), fHM->H1("fhRichDigisTimeLogZoom2")},
- {"STS", "TOF", "T0", "RICH"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
+ DrawH1({fHM->H1("fhRichDigisTimeLogZoom2"), fHM->H1("fhTofDigisTimeLogZoom2"), fHM->H1("fhT0DigisTimeLogZoom2"), fHM->H1("fhStsDigisTimeLogZoom2")},
+ {"RICH", "TOF", "T0", "STS"}, kLinear, kLog, true, 0.87, 0.75, 0.99, 0.99);
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.1);
fHM->H1("fhStsDigisTimeLogZoom2")->GetYaxis()->SetTitleOffset(0.7);
fHM->H1("fhStsDigisTimeLogZoom2")->SetMinimum(0.9);
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_richDigisRingTimeLog", "rich_mcbm_richDigisRingTimeLog", 1200, 600);
+ TH1D* copyRichDigi = (TH1D*) fHM->H1("fhRichDigisTimeLog")->Clone();
+ TH1D* copyRichRing = (TH1D*) fHM->H1("fhRichRingsTimeLog")->Clone();
+ DrawH1({copyRichDigi, fHM->H1("fhTofDigisTimeLog"), fHM->H1("fhT0DigisTimeLog"), copyRichRing},
+ {"RICH", "TOF", "T0", "RICH RING"}, kLinear, kLog, true, 0.83, 0.75, 0.99, 0.99);
+ gPad->SetLeftMargin(0.1);
+ gPad->SetRightMargin(0.1);
+ copyRichDigi->GetYaxis()->SetTitleOffset(0.7);
+ copyRichDigi->SetMinimum(0.9);
}
{
@@ -726,6 +1215,17 @@ void CbmRichMCbmQaReal::DrawHist()
c->cd(2);
DrawH1(fHM->H1("fhRichHitToT"));
}
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_BlobTrackDist","rich_BlobTrackDist", 1800 , 600);
+ c->Divide(3,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhBlobTrackDistX"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhBlobTrackDistY"));
+ c->cd(3);
+ DrawH1(fHM->H1("fhBlobTrackDist"));
+ }
{
fHM->CreateCanvas("ToF_XYZ","ToF_XYZ", 1200 , 1200);
@@ -749,18 +1249,34 @@ void CbmRichMCbmQaReal::DrawHist()
{
fHM->CreateCanvas("TofTracksPerEvent","TofTracksPerEvent", 1200 , 1200);
- fHM->H1("fhTofTracksPerEvent")->Draw();
+ DrawH1(fHM->H1("fhTofTracksPerEvent"),kLinear,kLog);
}
{
- TCanvas* c = fHM->CreateCanvas("TofRichRingX","TofRichRingX", 1200 , 1200);
+ fHM->CreateCanvas("TofTracksPerRichEvent","TofTracksPerRichEvent", 1200 , 1200);
+ DrawH1(fHM->H1("fhTofTracksPerRichEvent"),kLinear,kLog);
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofRichRingX_Y","TofRichRingX_Y", 1200 , 1200);
c->Divide(2,1);
c->cd(1);
- DrawH1(fHM->H1("fhTofRichRingX"));
+ DrawH2(fHM->H2("fhTofRichRingX"));
c->cd(2);
- DrawH1(fHM->H1("fhTofRichRingY"));
+ DrawH2(fHM->H2("fhTofRichRingY"));
}
+ {
+ TCanvas* c = fHM->CreateCanvas("TofRichX_Stacks","TofRichX_Stacks", 1800 , 600);
+ c->Divide(3,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofRichX_stack1"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofRichX_stack2"));
+ c->cd(3);
+ DrawH2(fHM->H2("fhTofRichX_stack3"));
+ }
+
{
fHM->CreateCanvas("TofRichRingXZ","TofRichRingXZ", 1200 , 1200);
@@ -773,6 +1289,12 @@ void CbmRichMCbmQaReal::DrawHist()
DrawH2(fHM->H2("fhTofHitsXZ"));
}
+ {
+ fHM->CreateCanvas("TofHitsXY","TofHitsXY", 1200 , 1200);
+
+ DrawH2(fHM->H2("fhTofHitsXY"));
+ }
+
{
fHM->CreateCanvas("TofTrackRichRingXY","TofTrackRichRingXY", 1200 , 1200);
@@ -794,6 +1316,11 @@ void CbmRichMCbmQaReal::DrawHist()
DrawH2(fHM->H2("fhTrackRingDistanceVSRingradius"));
}
+ {
+ fHM->CreateCanvas("RingTrackDist","RingTrackDist", 1200 , 800);
+ DrawH1(fHM->H1("fhRingTrackDistance"));
+ }
+
{
TCanvas* c = fHM->CreateCanvas("ClosTrackRingXY","ClosTrackRingXY", 1200 , 800);
c->Divide(2,1);
@@ -803,6 +1330,15 @@ void CbmRichMCbmQaReal::DrawHist()
DrawH2(fHM->H2("fhClosTrackRingY"));
}
+ {
+ TCanvas* c = fHM->CreateCanvas("TrackRingXY","TrackRingXY", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichRingX"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichRingY"));
+ }
+
{
TCanvas* c = fHM->CreateCanvas("TofTracksXY","TofTracksXY", 1200 , 800);
c->Divide(2,1);
@@ -812,6 +1348,16 @@ void CbmRichMCbmQaReal::DrawHist()
DrawH2(fHM->H2("fhTofTracksXYRICH"));
}
+ {
+ fHM->CreateCanvas("TofTracksXY_Target","TofTracksXY_Target", 800 , 800);
+ DrawH2(fHM->H2("fhTofTracksXY_Target"));
+ }
+
+ {
+ fHM->CreateCanvas("GoodRingsXY_TargetPos","GoodRingsXY_TargetPos", 800 , 800);
+ DrawH2(fHM->H2("fhGoodRingsXY_TargetPos"));
+ }
+
{
fHM->CreateCanvas("TofBeta","TofBeta", 800 , 800);
gStyle->SetOptStat(0);
@@ -848,9 +1394,13 @@ void CbmRichMCbmQaReal::DrawHist()
}
{
- TCanvas* c = fHM->CreateCanvas("TofBetaLog","TofBetaLog", 800 , 800);
+ TCanvas* c = fHM->CreateCanvas("TofBetaLog","TofBetaLog", 1000 , 1000);
c->SetLogy();
+
+ Double_t max = fHM->H1("fhTofBetaAll")->GetMaximum();
fHM->H1("fhTofBetaAll")->Draw("HIST");
+
+ if ( fHM->H1("fhTofBetaTracksWithHitsNoRing")->GetMaximum() > max) max = fHM->H1("fhTofBetaTracksWithHitsNoRing")->GetMaximum();
fHM->H1("fhTofBetaTracksWithHitsNoRing")->SetLineColorAlpha(kGreen, 1);
fHM->H1("fhTofBetaTracksWithHitsNoRing")->Draw("HIST SAME");
@@ -860,13 +1410,17 @@ void CbmRichMCbmQaReal::DrawHist()
//fHM->H1("fhTofBetaTracksWithHits")->SetLineColorAlpha(44, 1);
//fHM->H1("fhTofBetaTracksWithHits")->Draw("HIST SAME");
- fHM->H1("fhTofBetaTrackswithClosestRingInRange")->SetLineColorAlpha(44, 1);
+ if ( fHM->H1("fhTofBetaTrackswithClosestRingInRange")->GetMaximum() > max) max = fHM->H1("fhTofBetaTrackswithClosestRingInRange")->GetMaximum();
+ fHM->H1("fhTofBetaTrackswithClosestRingInRange")->SetLineColorAlpha(12, 1);
fHM->H1("fhTofBetaTrackswithClosestRingInRange")->Draw("HIST SAME");
+ if ( fHM->H1("fhTofBetaRing")->GetMaximum() > max) max = fHM->H1("fhTofBetaRing")->GetMaximum();
fHM->H1("fhTofBetaRing")->SetLineColorAlpha(kRed, 1);
fHM->H1("fhTofBetaRing")->Draw("HIST SAME");
- auto legend = new TLegend(0.1,0.75,0.4,0.9);
+ fHM->H1("fhTofBetaAll")->SetAxisRange(1., max*1.8,"Y");
+
+ auto legend = new TLegend(0.75,0.77,0.99,0.96);
legend->AddEntry(fHM->H1("fhTofBetaTracksWithHitsNoRing"),"Tracks with RichHits and no Ring","l");
//legend->AddEntry(fHM->H1("fhTofBetaTracksWithHits"),"Tracks with RichHits","l");
//legend->AddEntry(fHM->H1("fhTofBetaTracksNoRing"),"Tracks with no Ring","l");
@@ -940,6 +1494,203 @@ void CbmRichMCbmQaReal::DrawHist()
//c->SetLogy();
DrawH1(fHM->H1("fhDiRICHsInRegion"));
}
+
+ {
+ fHM->CreateCanvas("RichRingRadiusVsY","RichRingRadiusVsY", 800 , 800);
+ //c->SetLogy();
+ DrawH2(fHM->H2("fhRichRingRadiusY"));
+ }
+
+ {
+ fHM->CreateCanvas("RichRingHitsVsRadius","RichRingHitsVsRadius", 800 , 800);
+ //c->SetLogy();
+ DrawH2(fHM->H2("fhRichHitsRingRadius"));
+ }
+
+ {
+ fHM->CreateCanvas("RichHitToTEvent","RichHitToTEvent", 800 , 800);
+ DrawH1(fHM->H1("fhRichHitToTEvent"));
+ }
+
+ {
+ fHM->CreateCanvas("RichHitXYEvent","RichHitXYEvent", 800 , 800);
+ DrawH2(fHM->H2("fhRichHitXYEvent"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhBlobEvents","rich_mcbm_fhBlobEvents", 600 , 600);
+ DrawH1(fHM->H1("fhNofBlobEvents"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhBlobsInCbmEvent","rich_mcbm_fhBlobsInCbmEvent", 600 , 600);
+ DrawH1(fHM->H1("fhNofBlobsInEvent"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackRichHitXY","TofTrackRichHitXY", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichHitX"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichHitY"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackRichRingHitXY","TofTrackRichRingHitXY", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichRingHitX"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichRingHitY"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofHitRichRingHitXY","TofHitRichRingHitXY", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofHitRichRingHitX"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofHitRichRingHitY"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TrackRingDistance_Target","TrackRingDistance_Target", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhTrackRingDistanceOnTarget"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhTrackRingDistanceOffTarget"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackRichHitXY_cut","TofTrackRichHitXY_cut", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_cuts"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_cuts"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackRichHitXY_overBetaCut","TofTrackRichHitXY_overBetaCut", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_oBetacuts"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_oBetacuts"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackRichHitXY_underBetaCut","TofTrackRichHitXY_underBetaCut", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_uBetacuts"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_uBetacuts"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackRichHitXY_oBetacuts_dtime","TofTrackRichHitXY_oBetacuts_dtime", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackHitRichHitXY_oBetacuts_dtime","TofTrackHitRichHitXY_oBetacuts_dtime", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackHitRichHitX_oBetacuts_dtime"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackHitRichHitY_oBetacuts_dtime"));
+ }
+
+ {
+ fHM->CreateCanvas("TofTrackRichHitTime_oBetacuts","TofTrackRichHitTime_oBetacuts", 600 , 600);
+ DrawH1(fHM->H1("fhTofTrackRichHitTime_oBetacuts"));
+ }
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TofTrackHitRichHitXY_oBetacuts_dtime_evo","TofTrackHitRichHitXY_oBetacuts_dtime_evo", 1200 , 3200);
+ c->Divide(2,4);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_4"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_4"));
+
+ c->cd(3);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_6"));
+ c->cd(4);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_6"));
+
+ c->cd(5);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_8"));
+ c->cd(6);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_8"));
+
+ c->cd(7);
+ DrawH2(fHM->H2("fhTofTrackRichHitX_oBetacuts_dtime_10"));
+ c->cd(8);
+ DrawH2(fHM->H2("fhTofTrackRichHitY_oBetacuts_dtime_10"));
+ }
+
+ {
+ fHM->CreateCanvas("TrackRingDistance_AfterCorr","TrackRingDistance_AfterCorr", 600 , 600);
+ DrawH1(fHM->H1("fhTrackRingDistance_corr"));
+ }
+
+ {
+ fHM->CreateCanvas("TracksWithRings","TracksWithRings", 600 , 600);
+ DrawH1(fHM->H1("fhTracksWithRings"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("GoodRings","GoodRings", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichRingRadius_goodRing"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhNofHitsInRing_goodRing"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("richRingBeta","richRingBeta", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichRingBeta"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhRichRingBeta_GoodRing"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("TrackRingDistanceVSRingChi2","TrackRingDistanceVSRingChi2", 1200 , 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhTrackRingDistanceVSRingChi2"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhTrackRingDistanceVSRingChi2_goodRing"));
+ }
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("RichRingChi2","RichRingChi2", 1600 , 800);
+ //c->SetLogy();
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichRingChi2"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhRichRingChi2_goodRing"));
+ }
+
+ {
+ fHM->CreateCanvas("TofTracksXYRICH_Accectance","TofTracksXYRICH_Accectance", 1200 , 1200);
+ DrawH2(fHM->H2("fhTofTracksXYRICH_Accectance"));
+ }
+
}
void CbmRichMCbmQaReal::DrawRing(
@@ -956,7 +1707,7 @@ void CbmRichMCbmQaReal::DrawRing(
bool full
)
{
- std::cout<<"#!#DRAW!!!"<<std::endl;
+ //std::cout<<"#!#DRAW!!!"<<std::endl;
if (fNofDrawnRings > 20) return;
fNofDrawnRings++;
stringstream ss;
@@ -1033,6 +1784,7 @@ void CbmRichMCbmQaReal::DrawRing(
CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
if (nullptr == hit) continue;
TEllipse* hitDr = new TEllipse(hit->GetX() - xCur, hit->GetY() - yCur, .25);
+ //std::cout<<"LE of Hit: "<< hit->GetTime()- fCbmEventStartTime << "\t" << hit->GetTime() << "\t" << fCbmEventStartTime <<std::endl;
if (doToT(hit)) { // Good ToT selection
hitDr->SetFillColor(kRed);
} else {
@@ -1072,7 +1824,7 @@ void CbmRichMCbmQaReal::DrawRing(
TLatex* latex1 = new TLatex(-4., 0.5, ss3.str().c_str());
latex1->Draw();
} else {
- std::cout<<"No Tracks to Draw."<<std::endl;
+ //std::cout<<"No Tracks to Draw."<<std::endl;
}
@@ -1101,6 +1853,8 @@ void CbmRichMCbmQaReal::DrawRichTofEv(
if (richHitIndx.size() < 1 || tofTrackIndx.size() < 1 ) return;
if (fNofDrawnRichTofEv > 30) return;
fNofDrawnRichTofEv++;
+ std::string dir = "./"+fOutputDir+"/png/TREv/";
+ gSystem->mkdir(dir.c_str(), true); // create directory if it does not exist
stringstream ss;
ss << "TREv/TofRichEvent" << fNofDrawnRichTofEv;
TCanvas *c = nullptr;
@@ -1152,6 +1906,7 @@ void CbmRichMCbmQaReal::DrawRichTofEv(
for (unsigned int i = 0; i < tofTrackIndx.size(); i++){
CbmTofTracklet* track = (CbmTofTracklet*) fTofTracks->At(tofTrackIndx[i]);
+ if (track == nullptr) continue;
TEllipse* trDr = new TEllipse(track->GetFitX(hitZ), track->GetFitY(hitZ), .25);
trDr->SetFillColor(kGreen);
trDr->Draw();
@@ -1166,185 +1921,6 @@ void CbmRichMCbmQaReal::DrawRichTofEv(
}
-
-void CbmRichMCbmQaReal::DrawEvent(CbmEvent *ev, std::vector<int> &ringIndx, bool full=true){
-
-
- double pmtWidth = 5.20;
- double pmtHeight = 5.20;
- double pmtGap = 0.10;
-
- double left = 0.0;
- //double right = 0.0;
- double top = 0.0;
- //double bottom = 0.0;
- if (full == true){
- left = -16.85;
- //right = 4.25;
- top = 23.8;
- //bottom = -23.8;
- }
-
-
- if (fNofDrawnEvents > 50) return;
- if ((ev->GetNofData(ECbmDataType::kRichHit) <= 4) /*|| ringIndx.size() == 0*/) return;
- fNofDrawnEvents++;
- stringstream ss;
- ss << "Ev/CbmEvent" << fNofDrawnEvents;
- TCanvas *c = nullptr;
- if (full == true) {
- c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 600, 1000);
- } else {
- c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 800, 800);
- }
- c->SetGrid(true, true);
- TH2D* pad = nullptr;
- if (full == true){
- pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -18., 5., 1, -26., 26.);
- } else {
- pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -15., 10., 1, -5., 20);
- }
-
- pad->SetStats(false);
- pad->Draw();
-
- if (full == true){
- for(unsigned int x=0;x<4;++x){
- for(unsigned int y=0;y<9;++y){
- double pmtLeft = left+(pmtWidth+pmtGap)*x;
- double pmtTop = top-(pmtHeight+pmtGap)*y;
- TBox *box = new TBox( pmtLeft,
- pmtTop,
- pmtLeft + pmtWidth,
- pmtTop-pmtHeight
- );
- //box->SetFillColorAlpha(8,0.2);
- box->Draw();
-
- pmtLeft += 0.175;
- pmtTop -= 0.175;
- for (unsigned int pX=0; pX<8; ++pX){
- for (unsigned int pY=0; pY<8; ++pY){
- double xStart = 0.0;
- double xEnd = 0.0;
- double yStart = 0.0;
- double yEnd = 0.0;
- if (pX == 0){
- xStart = pmtLeft;
- xEnd = pmtLeft + 0.625;
- } else if (pX < 7){
- xStart = pmtLeft + 0.625 + 0.6*(pX-1);
- xEnd = pmtLeft + 0.625 + 0.6*(pX);
- } else if (pX == 7){
- xStart = pmtLeft + 0.625 + 0.6*6;
- xEnd = pmtLeft + 0.625*2 + 0.6*6;
- }
-
- if (pY == 0){
- yStart = pmtTop;
- yEnd = pmtTop - 0.625;
- } else if (pY < 7){
- yStart = pmtTop - 0.625 - 0.6*(pY-1);
- yEnd = pmtTop - 0.625 - 0.6*(pY);
- } else if (pY == 7){
- yStart = pmtTop - 0.625 - 0.6*6;
- yEnd = pmtTop - 0.625*2 - 0.6*6;
- }
-
- TBox *box1 = new TBox( xStart,
- yStart,
- xEnd,
- yEnd
- );
- box1->SetLineWidth(1.);
- //box1->SetFillColorAlpha(30.,0.1);
- //box1->SetLineColorAlpha(30.,0.5);
-
- box1->Draw("l");
- }
- }
- }
- }
- }
-
- if (full == true){
- //rough Drawing of RichDetectorAcceptance
-/* TLine *line0 = new TLine( left , top, right, top);
- line0->Draw();
- TLine *line1 = new TLine( right, top, right, bottom);
- line1->Draw();
- TLine *line2 = new TLine( right, bottom, left, bottom);
- line2->Draw();
- TLine *line3 = new TLine( left , bottom, left, top);
- line3->Draw(); */
- }
-
- // find min and max x and y positions of the hits
- // in order to shift drawing
-
-
- double hitZ = 0;
- uint nofDrawHits = 0;
-
- for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
- auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
- CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
- if (nullptr == hit) continue;
- TEllipse* hitDr = new TEllipse(hit->GetX(), hit->GetY(), .25);
- if (doToT(hit)) { // Good ToT selection
- hitDr->SetFillColor(kCyan);
- } else {
- hitDr->SetFillColor(kBlue);
- }
- hitZ += hit->GetZ();
- nofDrawHits++;
- hitDr->Draw();
- }
-
- //Draw circle and center
- for (unsigned int rings=0; rings < ringIndx.size();rings++) {
- CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(ringIndx[rings]));
-
- TEllipse* circle = new TEllipse(ring->GetCenterX(), ring->GetCenterY(), ring->GetRadius());
- circle->SetFillStyle(0);
- circle->SetLineWidth(3);
- circle->Draw();
- TEllipse* center = new TEllipse(ring->GetCenterX(), ring->GetCenterY(), .1);
- center->Draw();
-
- // Draw hits
- for (int i = 0; i < ring->GetNofHits(); i++){
- Int_t hitInd = ring->GetHit(i);
- CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
- if (nullptr == hit) continue;
- TEllipse* hitDr = new TEllipse(hit->GetX(), hit->GetY(), .125);
- if (doToT(hit)) { // Good ToT selection
- hitDr->SetFillColor(kMagenta);
- } else {
- hitDr->SetFillColor(kRed);
- }
- //hitZ += hit->GetZ();
- //nofDrawHits++;
- hitDr->Draw();
- }
- }
- hitZ /= nofDrawHits;
-
- //Draw Tracks
- auto nofTofTracks = ev->GetNofData(ECbmDataType::kTofTrack);
- for (int j=0;j<nofTofTracks;j++){
- auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, j);
- CbmTofTracklet *track = static_cast<CbmTofTracklet *>(fTofTracks->At(iTofTrack));
-
- TEllipse* hitDr = new TEllipse(track->GetFitX(hitZ), track->GetFitY(hitZ), .25);
- hitDr->SetFillColor(kGreen);
- hitDr->Draw();
- }
-
-}
-
-
-
void CbmRichMCbmQaReal::Finish()
{
//std::cout<<"Tracks: "<< fTofTracks->GetEntriesFast() <<std::endl;
@@ -1369,6 +1945,20 @@ void CbmRichMCbmQaReal::Finish()
}
gDirectory->cd( oldir->GetPath() );
}
+
+ std::cout<< "fTracksinRich: "<<fTracksinRich<<std::endl;
+ Double_t result = (static_cast<Double_t>(fTracksinRichWithRichHits[0])/static_cast<Double_t>(fTracksinRich))*100.;
+ std::cout<< "fTracksinRichWithRichHits > 0: "<<fTracksinRichWithRichHits[0]<<"\t Percent: "<< result <<"%"<<std::endl;
+ result = (static_cast<Double_t>(fTracksinRichWithRichHits[1])/static_cast<Double_t>(fTracksinRich))*100.;
+ std::cout<< "fTracksinRichWithRichHits > 5: "<<fTracksinRichWithRichHits[1]<<"\t Percent: "<< result <<"%"<<std::endl;
+ result = (static_cast<Double_t>(fTracksinRichWithRichHits[2])/static_cast<Double_t>(fTracksinRich))*100.;
+ std::cout<< "fTracksinRichWithRichHits > 10: "<<fTracksinRichWithRichHits[2]<<"\t Percent: "<< result <<"%"<<std::endl;
+ result = (static_cast<Double_t>(fTracksinRichWithRichHits[3])/static_cast<Double_t>(fTracksinRich))*100.;
+ std::cout<< "fTracksinRichWithRichHits > 15: "<<fTracksinRichWithRichHits[3]<<"\t Percent: "<< result <<"%"<<std::endl;
+
+
+ std::cout<< "RingsWithTrack: "<<fRingsWithTrack[0]<<" Rings: "<<fRingsWithTrack[1] <<" Tracks: "<<fRingsWithTrack[2]<<std::endl;
+ std::cout<< "\t \t Rings Cutted: "<<fRingsWithTrack[3]<<" Tracks Cutted: "<<fRingsWithTrack[4] <<" Possible Combinations: "<<fRingsWithTrack[5]<<std::endl;
}
@@ -1391,10 +1981,10 @@ void CbmRichMCbmQaReal::DrawFromFile(
bool CbmRichMCbmQaReal::isAccmRICH(CbmTofTracklet *track){
//check if Track is in mRICH acceptance
- double x = track->GetFitX(RichZPos);
- double y = track->GetFitY(RichZPos);
bool inside = false;
if (!fRestrictToAcc) return true;
+ double x = track->GetFitX(RichZPos);
+ double y = track->GetFitY(RichZPos);
if ( x >= -6.25 && x <= -1.05 ){
// left part of mRICH
if ( y >8 && y < 15.9) {
@@ -1410,12 +2000,6 @@ bool CbmRichMCbmQaReal::isAccmRICH(CbmTofTracklet *track){
return inside;
}
-bool CbmRichMCbmQaReal::doToT(CbmRichHit * hit){
- bool check = false;
- if ((hit->GetToT() > 23.7) && (hit->GetToT() < 30.0)) check = true;
-
- return check;
-}
Double_t CbmRichMCbmQaReal::getBeta(CbmTofTracklet *track){
Double_t inVel = 1e7/(track->GetTt()); // in m/s
@@ -1424,6 +2008,28 @@ Double_t CbmRichMCbmQaReal::getBeta(CbmTofTracklet *track){
return beta;
}
+Double_t CbmRichMCbmQaReal::getBeta(CbmRichRing *ring){
+
+ // calculate distance of ring center to penetration point in aerogel with assumption,
+ // that particle is from target.
+ // INFO: use center of Aerogel as reference -> 11.5cm to MAPMT (TODO: get this from GEO-File)
+ const Double_t Aerogel_Dist = 11.5;
+ const Double_t aeroPenetr_Y = ring->GetCenterY() * (1. - Aerogel_Dist/(RichZPos));
+ const Double_t aeroPenetr_X = ring->GetCenterX() * (1. - Aerogel_Dist/(RichZPos));
+
+ const Double_t dist = std::sqrt( (ring->GetCenterX()-aeroPenetr_X) * (ring->GetCenterX()-aeroPenetr_X)
+ + (ring->GetCenterY()-aeroPenetr_Y) * (ring->GetCenterY()-aeroPenetr_Y)
+ + Aerogel_Dist*Aerogel_Dist );
+
+ const Double_t ioR = 1.05; // Index of Reflection
+
+ Double_t beta = std::sqrt( (((ring->GetRadius()*ring->GetRadius())/dist*dist)+1) / (ioR*ioR) );
+
+ //std::cout<<" beta:"<<beta<<std::endl;
+
+ return beta;
+}
+
bool CbmRichMCbmQaReal::RestrictToFullAcc(CbmTofTracklet *track){
//check if Track is in mRICH acceptance
@@ -1455,7 +2061,7 @@ bool CbmRichMCbmQaReal::RestrictToFullAcc(Double_t x, Double_t y)
void CbmRichMCbmQaReal::analyseRing(CbmRichRing *ring, CbmEvent *ev,std::pair<int, double> &clTrack){
- std::cout<<"analyse a Ring"<<std::endl;
+ //std::cout<<"analyse a Ring"<<std::endl;
Double_t meanTime = 0.;
unsigned int hitCnt = 0;
@@ -1484,7 +2090,7 @@ void CbmRichMCbmQaReal::analyseRing(CbmRichRing *ring, CbmEvent *ev,std::pair<in
CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
if (nullptr == hit) continue;
//std::cout<<"DeltatTime: "<< meanTime - hit->GetTime()<<std::endl;
- fHM->H1("fhRingDeltaTime")->Fill(static_cast<Double_t>(meanTime - hit->GetTime()));
+ //fHM->H1("fhRingDeltaTime")->Fill(static_cast<Double_t>(meanTime - hit->GetTime()));
fHM->H1("fhRingToTs")->Fill(hit->GetToT());
fHM->H1("fhRingLE")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()));
@@ -1493,14 +2099,15 @@ void CbmRichMCbmQaReal::analyseRing(CbmRichRing *ring, CbmEvent *ev,std::pair<in
//tmpRingIndx.push_back(ring->GetIndex);
const Double_t Tdiff_ring = (hit->GetTime()-ev->GetStartTime());
if ((Tdiff_ring > 20.) && (Tdiff_ring < 30.)){
- std::cout<<ev->GetNumber()<<" Address_ring: "<<std::hex<< CbmRichUtil::GetDirichId(hit->GetAddress())<<std::dec<<" "<< CbmRichUtil::GetDirichChannel(hit->GetAddress()) <<" "<< hit->GetToT()<<" "<<ring->GetRadius()<<std::endl;
+ // std::cout<<ev->GetNumber()<<" Address_ring: "<<std::hex<< CbmRichUtil::GetDirichId(hit->GetAddress())<<std::dec<<" "<< CbmRichUtil::GetDirichChannel(hit->GetAddress()) <<" "<< hit->GetToT()<<" "<<ring->GetRadius()<<std::endl;
//fHM->H1("fhDiRICHsInRegion")->Fill(CbmRichUtil::GetDirichId(hit->GetAddress()));
}
if (clTrack.first == -1) fHM->H1("fhRingNoClTrackLE")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()));
- if ((clTrack.first >= 0) && !(clTrack.second < 5.)) fHM->H1("fhRingClTrackFarAwayLE")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()));
+ if ((clTrack.first >= 0) && !(clTrack.second < 10.)) fHM->H1("fhRingClTrackFarAwayLE")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()));
if (cutDistance(clTrack) && cutRadius(ring)){ //Good Ring
fHM->H1("fhGoodRingLE")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()));
+ fHM->H1("fhRingDeltaTime")->Fill(static_cast<Double_t>(meanTime - hit->GetTime()));
}
}
@@ -1518,8 +2125,8 @@ void CbmRichMCbmQaReal::analyseRing(CbmRichRing *ring, CbmEvent *ev,std::pair<in
const Double_t Tdiff_inner = (richHit->GetTime()-ev->GetStartTime());
if ((Tdiff_inner > 20.) && (Tdiff_inner < 30.)){
InnerHitCnt_cut++;
- if (InnerHitCnt_cut == 1) {DrawRing(ring);}
- std::cout<<ev->GetNumber()<<" Address_inner: "<<std::hex<< CbmRichUtil::GetDirichId(richHit->GetAddress())<<std::dec<<" "<< CbmRichUtil::GetDirichChannel(richHit->GetAddress()) <<" "<< richHit->GetToT()<<" "<<ring->GetRadius()<<std::endl;
+ //if (InnerHitCnt_cut == 1) {DrawRing(ring);}
+ //std::cout<<ev->GetNumber()<<" Address_inner: "<<std::hex<< CbmRichUtil::GetDirichId(richHit->GetAddress())<<std::dec<<" "<< CbmRichUtil::GetDirichChannel(richHit->GetAddress()) <<" "<< richHit->GetToT()<<" "<<ring->GetRadius()<<std::endl;
fHM->H1("fhDiRICHsInRegion")->Fill(CbmRichUtil::GetDirichId(richHit->GetAddress()));
}
@@ -1550,7 +2157,7 @@ Bool_t CbmRichMCbmQaReal::cutRadius(CbmRichRing *ring){
Bool_t CbmRichMCbmQaReal::cutDistance(std::pair<int, double> &clTrack){
- if ( (clTrack.first >= 0) && (clTrack.second < 5.)) return true;
+ if ( (clTrack.first >= 0) && (clTrack.second < 10.)) return true;
return false;
}
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.h b/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.h
index 9aa5fc90..f564f796 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.h
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmQaReal.h
@@ -10,9 +10,11 @@
class TClonesArray;
class CbmRichRing;
class CbmRichHit;
-class CbmTofTracklet;
+#include "CbmTofTracklet.h"
class TVector3;
class CbmDigiManager;
+class CbmRichMCbmSEDisplay;
+#include "CbmTofDigi.h"
#include <vector>
#include <map>
@@ -98,31 +100,92 @@ public:
/**
* Move X-Position of mRICH in Histograms (e.g. for Geometry changes)
*/
- void XOffsetHistos (Double_t offset = 0.){
- fXOffsetHisto = offset;
+ void XOffsetHistos (Double_t val = 0.){
+ fXOffsetHisto = val;
}
+
+ /**
+ * Limit of Single Event Displays that should be drawn
+ */
+ void SetMaxNofDrawnEvents(Int_t val = 100){
+ fMaxNofDrawnEvents = val;
+ }
+
+ /**
+ * Set an trigger on the tof Hits.
+ */
+ void SetTriggerTofHits (Int_t val = 0){
+ fTriggerTofHits = val;
+ }
+
+ /**
+ * Set an trigger on the RICH Hits.
+ */
+ void SetTriggerRichHits (Int_t val = 0){
+ fTriggerRichHits = val;
+ }
+
+ /**
+ * Set an ToT cut of the RICH Hits.
+ */
+ void SetTotRich (Double_t min, Double_t max){
+ fTotRichMin = min;
+ fTotRichMax = max;
+ }
+
+ bool isOnTarget(CbmTofTracklet* tTrack){
+
+ Double_t val = std::sqrt(tTrack->GetFitX(0.)*tTrack->GetFitX(0.)+tTrack->GetFitY(0.)*tTrack->GetFitY(0.));
+ if ( val < 10.) return true;
+
+ return false;
+ }
private:
CbmDigiManager* fDigiMan = nullptr;
- TClonesArray* fT0Digis;
+
+ //TClonesArray* fT0Digis;
+ const std::vector<CbmTofDigi>* fT0Digis = nullptr;
+
TClonesArray* fRichHits;
+
TClonesArray* fRichRings;
+
TClonesArray* fTofHits;
+
TClonesArray* fTofTracks;
+
TClonesArray* fCbmEvent;
+
CbmHistManager* fHM;
+
Double_t fXOffsetHisto;
+ Double_t fTotRichMin;
+
+ Double_t fTotRichMax;
+
Int_t fEventNum;
Int_t fNofDrawnRings;
Int_t fNofDrawnRichTofEv;
+
+ Int_t fMaxNofDrawnEvents;
+
+ Int_t fTriggerRichHits;
+
+ Int_t fTriggerTofHits;
+
+ Int_t fTracksinRich = 0;
+
+ Int_t fRingsWithTrack[6] = {0,0,0,0,0,0};//rwt;ring;track;ringCut;trackCut;combinations;
- Int_t fNofDrawnEvents;
+ Int_t fTracksinRichWithRichHits[4] = {0,0,0,0};
+
string fOutputDir; // output dir for results
@@ -132,11 +195,24 @@ private:
bool fDoWriteHistToFile = true;
bool fDoDrawCanvas = true;
+
+ bool fDigiHitsInitialized = false;
bool RestrictToFullAcc(CbmTofTracklet *track);
bool RestrictToFullAcc(TVector3 &pos);
bool RestrictToFullAcc(Double_t x, Double_t y);
+ Double_t fCbmEventStartTime = 0.;
+ CbmEvent *fEventPnt = nullptr;
+
+ std::array<Double_t,2304> offset_read;
+ std::array<Double_t,2304> offset;
+ std::array<uint32_t,2304> offset_cnt;
+
+ CbmRichMCbmSEDisplay* fSeDisplay = nullptr;
+
+ CbmRichMCbmSEDisplay* fSeDsply_TR = nullptr;
+
/**
* \brief Initialize histograms.
*/
@@ -165,15 +241,21 @@ private:
bool isAccmRICH(CbmTofTracklet *track);
- bool doToT(CbmRichHit* hit);
+ template<typename T=CbmRichHit>
+ bool doToT(T* hit){
+ if ((hit->GetToT() > fTotRichMin) && (hit->GetToT() < fTotRichMax)) return true;
+ return false;
+ }
Double_t getBeta(CbmTofTracklet *track);
+ Double_t getBeta(CbmRichRing *ring);
+
void analyseRing(CbmRichRing *ring, CbmEvent *ev,std::pair<int, double> &clTrack);
Bool_t cutRadius(CbmRichRing *ring);
Bool_t cutDistance(std::pair<int, double> &clTrack);
-
+
/**
* \brief Copy constructor.
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.cxx
new file mode 100644
index 00000000..6e9cfce1
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.cxx
@@ -0,0 +1,757 @@
+#include "CbmRichMCbmQaRichOnly.h"
+
+#include "TH1D.h"
+#include "TH1.h"
+#include "TCanvas.h"
+#include "TClonesArray.h"
+#include "TF1.h"
+#include "TStyle.h"
+#include "TEllipse.h"
+#include "TLine.h"
+#include "TMarker.h"
+#include "TGeoNode.h"
+#include "TGeoManager.h"
+#include "TGeoBBox.h"
+#include "TMath.h"
+#include "TSystem.h"
+
+#include <TLegend.h>
+#include <TBox.h>
+
+#include "CbmDigiManager.h"
+#include "CbmRichDigi.h"
+#include "CbmRichHit.h"
+#include "CbmRichRingFinderHoughImpl.h"
+#include "TLatex.h"
+#include "CbmDrawHist.h"
+#include "CbmTrackMatchNew.h"
+#include "CbmRichRing.h"
+#include "CbmRichHit.h"
+#include "CbmMatchRecoToMC.h"
+#include "CbmRichGeoManager.h"
+#include "CbmRichPoint.h"
+#include "CbmRichDraw.h"
+#include "CbmGlobalTrack.h"
+#include "CbmTrdTrack.h"
+#include "CbmTofHit.h"
+#include "CbmTofDigi.h"
+#include "CbmStsDigi.h"
+#include "CbmEvent.h"
+#include "CbmTofTracklet.h"
+#include "CbmRichUtil.h"
+
+#include "CbmRichMCbmSEDisplay.h"
+
+#include "CbmEvent.h"
+
+#include "CbmRichConverter.h"
+
+#include "CbmUtils.h"
+#include "CbmHistManager.h"
+
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <boost/assign/list_of.hpp>
+#include <sstream>
+#include <cmath>
+#include <fstream>
+
+
+using namespace std;
+using boost::assign::list_of;
+
+#define RichZPos 348.
+
+CbmRichMCbmQaRichOnly::CbmRichMCbmQaRichOnly() :
+ FairTask("CbmRichMCbmQaRichOnly"),
+ fRichHits(nullptr),
+ fRichRings(nullptr),
+ fCbmEvent(nullptr),
+ fHM(nullptr),
+ fEventNum(0),
+ fNofDrawnRings(0),
+ fNofDrawnRichTofEv(0),
+ fNofDrawnEvents(0),
+ fMaxNofDrawnEvents(100),
+ fTriggerRichHits(0),
+ fOutputDir("result")
+{
+
+}
+
+InitStatus CbmRichMCbmQaRichOnly::Init()
+{
+ cout << "CbmRichMCbmQaRichOnly::Init" << endl;
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (nullptr == ioman) { Fatal("CbmRichMCbmQaRichOnly::Init","RootManager not instantised!"); }
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+
+ if ( ! fDigiMan->IsPresent(ECbmModuleId::kRich) )
+ Fatal("CbmRichMCbmQaReal::Init", "No Rich Digis!");
+
+ fRichHits =(TClonesArray*) ioman->GetObject("RichHit");
+ if (nullptr == fRichHits) { Fatal("CbmRichMCbmQaRichOnly::Init", "No Rich Hits!");}
+
+ fRichRings =(TClonesArray*) ioman->GetObject("RichRing");
+ if (nullptr == fRichRings) { Fatal("CbmRichMCbmQaRichOnly::Init", "No Rich Rings!");}
+
+ fCbmEvent =(TClonesArray*) ioman->GetObject("CbmEvent");
+ if (nullptr == fCbmEvent) { Fatal("CbmRichMCbmQaRichOnly::Init", "No Event!");}
+
+ InitHistograms();
+
+
+ fSeDisplay = new CbmRichMCbmSEDisplay(fHM);
+ fSeDisplay->SetRichHits(fRichHits);
+ fSeDisplay->SetRichRings(fRichRings);
+ ///fSeDisplay->SetTofTracks(fTofTracks);
+ fSeDisplay->SetTotRich(23.7,30.);
+ fSeDisplay->SetMaxNofDrawnEvents(fMaxNofDrawnEvents);
+ fSeDisplay->SetOutDir(fOutputDir);
+
+ //Init OffsetCorrection ICD
+ for (auto& a : ICD_offset_read) a = 0.;
+ for (auto& a : ICD_offset) a = 0.;
+ for (auto& a : ICD_offset_cnt) a = 0;
+ read_ICD(ICD_offset_read,0);
+
+ return kSUCCESS;
+}
+
+void CbmRichMCbmQaRichOnly::InitHistograms()
+{
+ fHM = new CbmHistManager();
+
+ fHM->Create1<TH1D>("fhNofEvents","fhNofEvents;Entries", 1, 0.5, 1.5);
+ fHM->Create1<TH1D>("fhNofCbmEvents","fhNofCbmEvents;Entries", 1, 0.5, 1.5);
+ fHM->Create1<TH1D>("fhNofCbmEventsRing","fhNofCbmEventsRing;Entries", 1, 0.5, 1.5);
+
+ fHM->Create1<TH1D>("fhNofBlobEvents","fhNofBlobEvents;Entries", 1, 0.5, 1.5);
+ fHM->Create1<TH1D>("fhNofBlobsInEvent","fhNofBlobsInEvent;Entries", 36, 0.5, 36.5);
+
+ // RICH hits
+ fHM->Create2<TH2D>("fhRichHitXY","fhRichHitXY;RICH hit X [cm];RICH hit Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhRichHitXY_fromRing","fhRichHitXY_fromRing;RICH hit X [cm];RICH hit Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+
+ //ToT
+ fHM->Create1<TH1D>("fhRichDigisToT","fhRichDigisToT;ToT [ns];Entries", 601, 9.975, 40.025);
+ fHM->Create1<TH1D>("fhRichHitToT","fhRichHitToT;ToT [ns];Entries", 601, 9.975, 40.025);
+
+ // RICH rings
+ fHM->Create2<TH2D>("fhRichRingXY","fhRichRingXY;Ring center X [cm];Ring center Y [cm];Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 84, -25.2, 25.2);
+ fHM->Create1<TH1D>("fhRichRingRadius","fhRichRingRadius;Ring radius [cm];Entries", 100, 0., 7.);
+ fHM->Create1<TH1D>("fhNofHitsInRing","fhNofHitsInRing;# hits in ring;Entries", 50, -0.5, 49.5);
+ fHM->Create2<TH2D>("fhICD","fhICD;channel;DeltaTime", 2305, -0.5, 2304.5, 31, -15.5,15.5);
+
+ fHM->Create2<TH2D>("fhRichRingRadiusY","fhRichRingRadiusY;Ring Radius [cm]; Y position[cm];Entries", 70, -0.05, 6.95, 84, -25.2, 25.2);
+ fHM->Create2<TH2D>("fhRichHitsRingRadius","fhRichHitsRingRadius;#Rich Hits/Ring; Ring Radius [cm];Entries", 50, -0.5, 49.5, 70, -0.05, 6.95);
+
+ fHM->Create1<TH1D>("fhRingDeltaTime","fhRingDeltaTime; \\Delta Time/ns;Entries", 101, -10.1, 10.1);
+ fHM->Create1<TH1D>("fhRingChi2","fhRingChi2; \\Chi^2 ;Entries", 101, 0.0, 10.1);
+
+ fHM->Create2<TH2D>("fhRichRingCenterXChi2","fhRichRingCenterXChi2;Ring Center X [cm];\\Chi^2 ;;Entries", 67, -20.1 + fXOffsetHisto, 20.1 + fXOffsetHisto, 101, 0.0, 10.1 );
+ fHM->Create2<TH2D>("fhRichRingCenterYChi2","fhRichRingCenterYChi2;Ring Center Y [cm];\\Chi^2 ;;Entries", 84, -25.2, 25.2, 101, 0.0, 10.1);
+ fHM->Create2<TH2D>("fhRichRingRadiusChi2" ,"fhRichRingRadiusChi2; Ring Radius [cm];\\Chi^2 ;;Entries", 70, -0.05, 6.95, 101, 0.0, 10.1);
+
+
+ // Digis
+ fHM->Create2<TH2D>("fhDigisInChnl","fhDigisInChnl;channel;#Digis;",2304 , -0.5, 2303.5, 50, -0.5, 49.5);
+ fHM->Create2<TH2D>("fhDigisInDiRICH","fhDigisInDiRICH;DiRICH;#Digis;",72 , -0.5, 71.5, 300, -0.5, 299.5);
+}
+
+
+void CbmRichMCbmQaRichOnly::Exec(Option_t* /*option*/)
+{
+ fEventNum++;
+ fHM->H1("fhNofEvents")->Fill(1);
+
+ cout << "CbmRichMCbmQaRichOnly, event No. " << fEventNum << endl;
+
+ std::array<unsigned int,2304> chnlDigis;
+ for (auto& c : chnlDigis) c = 0;
+ for (int i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
+ const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(i);
+ fHM->H1("fhRichDigisToT")->Fill(richDigi->GetToT());
+
+ uint16_t addrDiRICH = (richDigi->GetAddress()>>16) & 0xFFFF;
+ uint16_t addrChnl = richDigi->GetAddress() & 0xFFFF;
+ uint16_t dirichNmbr = ((addrDiRICH >> 8) & 0xF)*18
+ + ((addrDiRICH >> 4) & 0xF)*2
+ + ((addrDiRICH >> 0) & 0xF);
+ uint32_t fullNmbr = (dirichNmbr << 5) | (addrChnl-1);
+ chnlDigis[fullNmbr]++;
+ }
+
+ {
+ unsigned int sum = 0;
+ for (uint16_t i = 0; i< 2304; ++i) {
+ if (chnlDigis[i] != 0 )fHM->H1("fhDigisInChnl")->Fill(i,chnlDigis[i]);
+ sum += chnlDigis[i];
+ if (i%32 == 31) {
+ uint16_t dirich = i/32;
+ if (sum != 0 ) fHM->H1("fhDigisInDiRICH")->Fill(dirich,sum);
+ sum = 0;
+ }
+
+ }
+ }
+
+ int nofRichHits = fRichHits->GetEntries();
+ for ( int iH = 0; iH < nofRichHits; iH++) {
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
+ fHM->H2("fhRichHitXY")->Fill(richHit->GetX(), richHit->GetY());
+ fHM->H1("fhRichHitToT")->Fill(richHit->GetToT());
+ }
+
+
+ //CBMEVENT
+ auto fNCbmEvent = fCbmEvent->GetEntriesFast();
+
+ for (int i=0;i<fNCbmEvent;i++){
+ fHM->H1("fhNofCbmEvents")->Fill(1);
+ CbmEvent* ev = static_cast<CbmEvent*>(fCbmEvent->At(i));
+
+ if (fTriggerRichHits != 0 && (ev->GetNofData(ECbmDataType::kRichHit) < fTriggerRichHits)) continue;
+
+
+ std::vector<int> ringIndx;
+ std::vector<int> evRichHitIndx;
+ std::array<uint32_t,36> pmtHits;
+ for (auto& a : pmtHits) a = 0;
+
+ // Map Rings to CbmEvent
+ for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
+ evRichHitIndx.push_back(iRichHit);
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
+
+ uint32_t pmtId = (((richHit->GetAddress())>>20)&0xF) + (((richHit->GetAddress())>>24)&0xF)*9;
+ pmtHits[pmtId]++;
+
+ int nofRichRings= fRichRings->GetEntries();
+ for (int l = 0; l < nofRichRings; l++)
+ {
+ CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(l));
+
+ auto NofRingHits = ring->GetNofHits();
+ for (int m = 0; m < NofRingHits; m++)
+ {
+ auto RingHitIndx = ring->GetHit(m);
+ if (RingHitIndx == iRichHit) {
+ Bool_t used = false;
+ for (auto check : ringIndx) {
+ if (check == l){ used = true; break;}
+ }
+ if (used == false) ringIndx.push_back(l);
+ break;
+ }
+ }
+ }
+ }
+
+ uint16_t blob = 0;
+ for (auto a : pmtHits ) {
+ if (a > 30) {blob++;}
+ }
+ if (blob > 0) {
+ fHM->H1("fhNofBlobEvents")->Fill(1);
+ fHM->H1("fhNofBlobsInEvent")->Fill(blob);
+ }
+
+ if (ringIndx.size() != 0) fHM->H1("fhNofCbmEventsRing")->Fill(1);
+
+ //Ring Loop
+ for (unsigned int k = 0 ; k < ringIndx.size(); ++k){
+ //Rigs in this CbmEvent
+ CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(ringIndx[k]));
+ if (ring == nullptr) continue;
+ fHM->H1("fhRingChi2")->Fill(ring->GetChi2());
+ fHM->H2("fhRichRingCenterXChi2")->Fill(ring->GetCenterX(), ring->GetChi2());
+ fHM->H2("fhRichRingCenterYChi2")->Fill(ring->GetCenterY(), ring->GetChi2());
+ fHM->H2("fhRichRingRadiusChi2") ->Fill(ring->GetRadius() , ring->GetChi2());
+
+ fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
+ fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
+ fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
+ fHM->H2("fhRichRingRadiusY")->Fill(ring->GetRadius(),ring->GetCenterY());
+ fHM->H2("fhRichHitsRingRadius")->Fill(ring->GetNofHits(),ring->GetRadius());
+ for (int j = 0; j < ring->GetNofHits(); ++j){
+ Int_t hitIndx = ring->GetHit(j);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitIndx);
+ if (nullptr == hit) continue;
+ fHM->H2("fhRichHitXY_fromRing")->Fill(hit->GetX(), hit->GetY());
+ }
+
+ }
+ fSeDisplay->DrawEvent(static_cast<CbmEvent*>(fCbmEvent->At(i)), ringIndx, 1);
+
+ } //End CbmEvent loop
+
+ // Loop over all Rings
+ RichRings();
+
+}
+
+void CbmRichMCbmQaRichOnly::RichRings()
+{
+ int nofRichRings = fRichRings->GetEntries();
+ for (int i = 0; i < nofRichRings; i++) {
+ CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(i));
+ if (ring == nullptr) continue;
+
+ for (int j = 1; j < ring->GetNofHits(); ++j){
+ Int_t hitIndx = ring->GetHit(j);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitIndx);
+ if (nullptr == hit) continue;
+
+ //Read Address
+ uint32_t DiRICH_Addr = hit->GetAddress();
+ unsigned int addr = (((DiRICH_Addr>>24)&0xF)*18*32)
+ +(((DiRICH_Addr>>20)&0xF)*2*32)
+ +(((DiRICH_Addr>>16)&0xF)*32)
+ +((DiRICH_Addr&0xFFFF) - 0x1);
+ ICD_offset.at(addr) += hit->GetTime() - ring->GetTime() - ICD_offset_read.at(addr);
+ fHM->H2("fhICD")->Fill(addr,hit->GetTime() - ring->GetTime() - ICD_offset_read.at(addr));
+ fHM->H1("fhRingDeltaTime")->Fill(hit->GetTime() - ring->GetTime() - ICD_offset_read.at(addr));
+ ICD_offset_cnt.at(addr)++;
+
+
+ }
+
+ //DrawRing(ring);
+ /*fHM->H2("fhRichRingXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
+ fHM->H1("fhRichRingRadius")->Fill(ring->GetRadius());
+ fHM->H1("fhNofHitsInRing")->Fill(ring->GetNofHits());
+ */
+ }
+}
+
+
+
+void CbmRichMCbmQaRichOnly::DrawHist()
+{
+ cout.precision(4);
+
+ //SetDefaultDrawStyle();
+ double nofEvents = fHM->H1("fhNofCbmEvents")->GetEntries();
+ fHM->ScaleByPattern("fh_.*", 1./nofEvents);
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhNofCbmEvents","rich_mcbm_fhNofCbmEvents", 600 , 600);
+ DrawH1(fHM->H1("fhNofCbmEvents"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhNofCbmEventsRing","rich_mcbm_fhNofCbmEventsRing", 600 , 600);
+ DrawH1(fHM->H1("fhNofCbmEventsRing"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhNofEvents","rich_mcbm_fhNofEvents", 600 , 600);
+ DrawH1(fHM->H1("fhNofEvents"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhBlobEvents","rich_mcbm_fhBlobEvents", 600 , 600);
+ DrawH1(fHM->H1("fhNofBlobEvents"));
+ }
+
+ {
+ fHM->CreateCanvas("rich_mcbm_fhBlobsInCbmEvent","rich_mcbm_fhBlobsInCbmEvent", 600 , 600);
+ DrawH1(fHM->H1("fhNofBlobsInEvent"));
+ }
+
+ {
+ fHM->CreateCanvas("inner_channel_delay","inner_channel_delay", 1200 , 600);
+ DrawH2(fHM->H2("fhICD"));
+ }
+
+
+ {
+ fHM->CreateCanvas("RingDelta","RingDelta", 600 , 600);
+ DrawH1(fHM->H1("fhRingDeltaTime"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_ToT","rich_ToT", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichDigisToT"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhRichHitToT"));
+ }
+
+ {
+ fHM->CreateCanvas("DigisInChnl","DigisInChnl", 1200 , 600);
+ DrawH2(fHM->H2("fhDigisInChnl"));
+ }
+
+ {
+ fHM->CreateCanvas("DigisInDiRICH","DigisInDiRICH", 1200 , 600);
+ DrawH2(fHM->H2("fhDigisInDiRICH"));
+ }
+
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_Hits","rich_Hits", 1200 , 1200);
+ c->Divide(2,2);
+ c->cd(1);
+ DrawH2(fHM->H2("fhRichHitXY"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhRichHitXY_fromRing"));
+ c->cd(3);
+ TH2F *hitsBg = (TH2F*) fHM->H2("fhRichHitXY")->Clone();
+ hitsBg->SetName("fhRichHitXY_bg");
+ hitsBg->SetTitle("fhRichHitXY_bg");
+ hitsBg->Add(fHM->H2("fhRichHitXY_fromRing"),-1);
+ //hitsBg->Draw("colz");
+ DrawH2(hitsBg);
+ c->cd(4);
+ DrawH2(fHM->H2("fhRichRingXY"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("rich_mcbm_rings","rich_mcbm_rings", 1200 , 600);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRichRingRadius"));
+ c->cd(2);
+ DrawH1(fHM->H1("fhNofHitsInRing"));
+ }
+
+ {
+ fHM->CreateCanvas("RichRingRadiusVsY","RichRingRadiusVsY", 800 , 800);
+ //c->SetLogy();
+ DrawH2(fHM->H2("fhRichRingRadiusY"));
+ }
+
+ {
+ fHM->CreateCanvas("RichRingHitsVsRadius","RichRingHitsVsRadius", 800 , 800);
+ //c->SetLogy();
+ DrawH2(fHM->H2("fhRichHitsRingRadius"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("RichRingChi2","RichRingChi2", 1600 , 800);
+ //c->SetLogy();
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fHM->H1("fhRingChi2"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhRichRingRadiusChi2"));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("RichRingCenterChi2","RichRingCenterChi2", 1600 , 800);
+ //c->SetLogy();
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH2(fHM->H2("fhRichRingCenterXChi2"));
+ c->cd(2);
+ DrawH2(fHM->H2("fhRichRingCenterYChi2"));
+ }
+
+}
+
+void CbmRichMCbmQaRichOnly::DrawRing(
+ CbmRichRing* ring,
+ bool full
+ )
+{
+ //std::cout<<"#!#DRAW!!!"<<std::endl;
+ if (fNofDrawnRings > 20) return;
+ fNofDrawnRings++;
+ stringstream ss;
+ ss << "Event" << fNofDrawnRings;
+ //fNofDrawnRings++;
+ TCanvas *c = nullptr;
+ if (full == true) {
+ c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 800, 800);
+ } else {
+ c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 500, 500);
+ }
+ c->SetGrid(true, true);
+ TH2D* pad = nullptr;
+ if (full == true){
+ pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -15., 10., 1, -5., 20);
+ } else {
+ pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -5., 5., 1, -5., 5);
+ }
+
+ pad->SetStats(false);
+ pad->Draw();
+
+ if (full == true){
+ //rough Drawing of RichDetectorAcceptance
+ TLine *line0 = new TLine(-6.25,8,-6.25,15.9);
+ line0->Draw();
+ TLine *line1 = new TLine(-6.25,15.9,-1.05,15.9);
+ line1->Draw();
+ TLine *line2 = new TLine(-1.05,15.9,-1.05,13.2);
+ line2->Draw();
+ TLine *line3 = new TLine(-1.05,13.2,+4.25,13.2);
+ line3->Draw();
+ TLine *line4 = new TLine(4.25,13.2,+4.25,8);
+ line4->Draw();
+ TLine *line5 = new TLine(4.25,8,-6.25,8);
+ line5->Draw();
+ }
+
+ // find min and max x and y positions of the hits
+ // in order to shift drawing
+ double xCur = 0.;
+ double yCur = 0.;
+
+ if (full == false){
+ double xmin = 99999., xmax = -99999., ymin = 99999., ymax = -99999.;
+ for (int i = 0; i < ring->GetNofHits(); i++){
+ Int_t hitInd = ring->GetHit(i);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+ if (xmin > hit->GetX()) xmin = hit->GetX();
+ if (xmax < hit->GetX()) xmax = hit->GetX();
+ if (ymin > hit->GetY()) ymin = hit->GetY();
+ if (ymax < hit->GetY()) ymax = hit->GetY();
+ }
+ xCur = (xmin + xmax) / 2.;
+ yCur = (ymin + ymax) / 2.;
+ }
+
+ //Draw circle and center
+ TEllipse* circle = new TEllipse(ring->GetCenterX() - xCur, ring->GetCenterY() - yCur, ring->GetRadius());
+ circle->SetFillStyle(0);
+ circle->SetLineWidth(3);
+ circle->Draw();
+ TEllipse* center = new TEllipse(ring->GetCenterX() - xCur, ring->GetCenterY() - yCur, .1);
+ center->Draw();
+
+
+ double hitZ = 0;
+ uint nofDrawHits = 0;
+
+ // Draw hits
+ for (int i = 0; i < ring->GetNofHits(); i++){
+ Int_t hitInd = ring->GetHit(i);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+ TEllipse* hitDr = new TEllipse(hit->GetX() - xCur, hit->GetY() - yCur, .25);
+ //std::cout<<"LE of Hit: "<< hit->GetTime()- fCbmEventStartTime << "\t" << hit->GetTime() << "\t" << fCbmEventStartTime <<std::endl;
+ if (doToT(hit)) { // Good ToT selection
+ hitDr->SetFillColor(kRed);
+ } else {
+ hitDr->SetFillColor(kBlue);
+ }
+ hitZ += hit->GetZ();
+ nofDrawHits++;
+ hitDr->Draw();
+ }
+ hitZ /= nofDrawHits;
+
+
+
+ //Draw information
+ stringstream ss2;
+ ss2 << "(x,y,r,n)=(" << setprecision(3) << ring->GetCenterX() << ", " << ring->GetCenterY()<< ", "
+ << ring->GetRadius() << ", " << ring->GetNofHits()<<")";
+ TLatex* latex = nullptr;
+ if (full == true){
+ latex = new TLatex(ring->GetCenterX()-13., ring->GetCenterY()+5., ss2.str().c_str());
+ } else {
+ latex = new TLatex(-4., 4., ss2.str().c_str());
+ }
+ latex->Draw();
+}
+
+
+void CbmRichMCbmQaRichOnly::Finish()
+{
+
+ for (unsigned int i = 0; i < ICD_offset.size();++i) {
+ if (ICD_offset_cnt.at(i) == 0 ) {
+ ICD_offset.at(i) = 0.0;
+ } else {
+ ICD_offset.at(i) /= ICD_offset_cnt.at(i);
+ }
+ ICD_offset.at(i) += ICD_offset_read.at(i);
+ }
+ save_ICD(ICD_offset,0);
+
+ //std::cout<<"Address: "<<InterChannel[0].first << std::endl;
+ //std::cout<<"Tracks: "<< fTofTracks->GetEntriesFast() <<std::endl;
+ std::cout<<"Drawing Hists...";
+ DrawHist();
+ std::cout<<"DONE!"<<std::endl;
+
+ if (this->fDoDrawCanvas){
+ fHM->SaveCanvasToImage(fOutputDir,"png");
+ std::cout<<"Canvas saved to Images!"<<std::endl;
+ }
+
+ if (this->fDoWriteHistToFile){
+ TDirectory * oldir = gDirectory;
+ std::string s = fOutputDir + "/RecoHists.root";
+ TFile *outFile = new TFile(s.c_str(),"RECREATE");
+ if (outFile->IsOpen()) {
+ fHM->WriteToFile();
+ std::cout<<"Written to Root-file \""<< s << "\" ...";
+ outFile->Close();
+ std::cout<<"Done!"<<std::endl;
+ }
+ gDirectory->cd( oldir->GetPath() );
+ }
+}
+
+
+void CbmRichMCbmQaRichOnly::DrawFromFile(
+ const string& fileName,
+ const string& outputDir)
+{
+ fOutputDir = outputDir;
+
+ if (fHM != nullptr) delete fHM;
+
+ fHM = new CbmHistManager();
+ TFile* file = new TFile(fileName.c_str());
+ fHM->ReadFromFile(file);
+ DrawHist();
+
+ fHM->SaveCanvasToImage(fOutputDir);
+}
+
+bool CbmRichMCbmQaRichOnly::doToT(CbmRichHit * hit){
+ bool check = false;
+ if ((hit->GetToT() > 23.7) && (hit->GetToT() < 30.0)) check = true;
+
+ return check;
+}
+
+void CbmRichMCbmQaRichOnly::analyseRing(CbmRichRing *ring, CbmEvent *ev){
+
+ std::cout<<"analyse a Ring"<<std::endl;
+
+ Double_t meanTime = 0.;
+ unsigned int hitCnt = 0;
+ Double_t minRHit2 = std::numeric_limits<Double_t>::max();
+ for (int i = 0; i < ring->GetNofHits(); i++){
+ Int_t hitInd = ring->GetHit(i);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+
+ meanTime += hit->GetTime();
+ hitCnt++;
+
+ const Float_t diffX = hit->GetX() - ring->GetCenterX();
+ const Float_t diffY = hit->GetY() - ring->GetCenterY();
+ const Float_t tmpHitRadius2 = (diffX*diffX+diffY*diffY);
+
+ if (tmpHitRadius2 < minRHit2) {
+ minRHit2 = tmpHitRadius2;
+ }
+ }
+ meanTime = meanTime / hitCnt;
+
+ //std::cout<<"mean: "<<meanTime<<std::endl;
+ for (int i = 0; i < ring->GetNofHits(); i++){
+ Int_t hitInd = ring->GetHit(i);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+ //std::cout<<"DeltatTime: "<< meanTime - hit->GetTime()<<std::endl;
+ fHM->H1("fhRingDeltaTime")->Fill(static_cast<Double_t>(meanTime - hit->GetTime()));
+
+ fHM->H1("fhRingToTs")->Fill(hit->GetToT());
+ fHM->H1("fhRingLE")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()));
+ fHM->H2("fhRingLEvsToT")->Fill(static_cast<Double_t>(hit->GetTime()-ev->GetStartTime()),hit->GetToT());
+ //std::vector<int> tmpRingIndx;
+ //tmpRingIndx.push_back(ring->GetIndex);
+ const Double_t Tdiff_ring = (hit->GetTime()-ev->GetStartTime());
+ if ((Tdiff_ring > 20.) && (Tdiff_ring < 30.)){
+ std::cout<<ev->GetNumber()<<" Address_ring: "<<std::hex<< CbmRichUtil::GetDirichId(hit->GetAddress())<<std::dec<<" "<< CbmRichUtil::GetDirichChannel(hit->GetAddress()) <<" "<< hit->GetToT()<<" "<<ring->GetRadius()<<std::endl;
+ //fHM->H1("fhDiRICHsInRegion")->Fill(CbmRichUtil::GetDirichId(hit->GetAddress()));
+ }
+ }
+
+ int InnerHitCnt = 0;
+ int InnerHitCnt_cut = 0;
+ for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
+ CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
+ if (nullptr == richHit) continue;
+ const Float_t diffX = richHit->GetX() - ring->GetCenterX();
+ const Float_t diffY = richHit->GetY() - ring->GetCenterY();
+ //select inner Part of Ring
+ if ( diffX*diffX+diffY*diffY < minRHit2 ) {
+ InnerHitCnt++;
+ const Double_t Tdiff_inner = (richHit->GetTime()-ev->GetStartTime());
+ if ((Tdiff_inner > 20.) && (Tdiff_inner < 30.)){
+ InnerHitCnt_cut++;
+ //if (InnerHitCnt_cut == 1) {DrawRing(ring);}
+ std::cout<<ev->GetNumber()<<" Address_inner: "<<std::hex<< CbmRichUtil::GetDirichId(richHit->GetAddress())<<std::dec<<" "<< CbmRichUtil::GetDirichChannel(richHit->GetAddress()) <<" "<< richHit->GetToT()<<" "<<ring->GetRadius()<<std::endl;
+ fHM->H1("fhDiRICHsInRegion")->Fill(CbmRichUtil::GetDirichId(richHit->GetAddress()));
+ }
+
+ fHM->H1("fhInnerRingDeltaTime")->Fill(static_cast<Double_t>(meanTime - richHit->GetTime()));
+ fHM->H1("fhInnerRingToTs")->Fill(richHit->GetToT());
+ fHM->H1("fhInnerRingLE")->Fill(static_cast<Double_t>(richHit->GetTime()-ev->GetStartTime()));
+
+ }
+ }
+ if (InnerHitCnt == 0) {
+ fHM->H1("fhInnerRingFlag")->Fill(1);
+ } else {
+ fHM->H1("fhInnerRingFlag")->Fill(0);
+ }
+ fHM->H1("fhNofInnerHits")->Fill(InnerHitCnt);
+}
+
+
+Bool_t CbmRichMCbmQaRichOnly::cutRadius(CbmRichRing *ring){
+ if (ring->GetRadius() > 2. && ring->GetRadius() < 4.2 ) return true;
+
+ return false;
+}
+
+void CbmRichMCbmQaRichOnly::save_ICD (std::array<Double_t,2304> &ICD_offsets, unsigned int iteration) {
+ std::ofstream icd_file;
+ icd_file.open (Form("icd_offset_it_%u.data",iteration));
+ if (icd_file.is_open()) {
+ for (unsigned int i = 0; i< ICD_offsets.size();++i){
+ //loop over all Offsets
+ icd_file<<i<<"\t"<<std::setprecision(25)<<ICD_offsets.at(i)<<"\n";
+ }
+ icd_file.close();
+ }
+ else std::cout << "Unable to open inter channel delay file icd_offset_it_"<< iteration <<".data\n";
+
+}
+
+
+void CbmRichMCbmQaRichOnly::read_ICD (std::array<Double_t,2304> &ICD_offsets, unsigned int iteration) {
+ std::cout<<gSystem->pwd()<<std::endl;
+ std::string line;
+ std::ifstream icd_file (Form("icd_offset_it_%u.data",iteration));
+ unsigned int lineCnt = 0;
+ if (icd_file.is_open()) {
+ while ( getline (icd_file,line) ) {
+ if (line[0] == '#') continue; // just a comment
+ std::istringstream iss(line);
+ unsigned int addr = 0;
+ Double_t value;
+ if (!(iss >> addr >> value)) {
+ std::cout << "A Problem accured in line "<< lineCnt <<"\n";
+ break;
+ } // error
+ lineCnt++;
+ ICD_offsets.at(addr) += value;
+ }
+ icd_file.close();
+ }
+ else std::cout << "Unable to open inter channel delay file icd_offset_it_"<< iteration <<".data\n";
+
+}
+
+ClassImp(CbmRichMCbmQaRichOnly)
+
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.h b/reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.h
new file mode 100644
index 00000000..89e1b073
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmQaRichOnly.h
@@ -0,0 +1,208 @@
+#ifndef MCBM_RICH_QA_RICH_ONLY
+#define MCBM_RICH_QA_RICH_ONLY
+
+#include "FairTask.h"
+#include "CbmRichRingFinderHoughImpl.h"
+#include "CbmEvent.h"
+class TClonesArray;
+class CbmRichRing;
+class CbmRichHit;
+class CbmTofTracklet;
+class CbmHistManager;
+class TVector3;
+class CbmDigiManager;
+class CbmRichMCbmSEDisplay;
+
+#include <vector>
+#include <map>
+#include <tuple>
+#include <array>
+
+using namespace std;
+
+
+
+class CbmRichMCbmQaRichOnly : public FairTask
+{
+
+public:
+ /**
+ * \brief Standard constructor.
+ */
+ CbmRichMCbmQaRichOnly();
+
+ /**
+ * \brief Standard destructor.
+ */
+ virtual ~CbmRichMCbmQaRichOnly() {};
+
+ /**
+ * \brief Inherited from FairTask.
+ */
+ virtual InitStatus Init();
+
+ /**
+ * \brief Inherited from FairTask.
+ */
+ virtual void Exec(Option_t* option);
+
+ /**
+ * \brief Inherited from FairTask.
+ */
+ virtual void Finish();
+
+ /**
+ * \brief Set output directory where you want to write results (figures and json).
+ * \param[in] dir Path to the output directory.
+ */
+ void SetOutputDir(const string& dir) {fOutputDir = dir;}
+
+
+ /**
+ * \brief Draw histogram from file
+ */
+ void DrawFromFile(
+ const string& fileName,
+ const string& outputDir);
+
+ /**
+ * Apply restriction to current mRICH Acceptance (for Simulations)
+ */
+ void DoRestrictToAcc(){
+ fRestrictToAcc = true;
+ }
+
+
+ /**
+ * Apply restriction to full mRICH Acceptance (for Simulations)
+ */
+ void DoRestrictToFullAcc(bool val=true){
+ fRestrictToFullAcc = val;
+ }
+
+
+ /**
+ * Apply restriction to full mRICH Acceptance (for Simulations)
+ */
+ void DoDrawCanvas (bool val=true){
+ fDoDrawCanvas = val;
+ }
+
+ /**
+ * Apply restriction to full mRICH Acceptance (for Simulations)
+ */
+ void DoWriteHistToFile (bool val=true){
+ fDoWriteHistToFile = val;
+ }
+
+ /**
+ * Move X-Position of mRICH in Histograms (e.g. for Geometry changes)
+ */
+ void XOffsetHistos (Double_t val = 0.){
+ fXOffsetHisto = val;
+ }
+
+ /**
+ * Limit of Single Event Displays that should be drawn
+ */
+ void SetMaxNofDrawnEvents(Int_t val = 100){
+ fMaxNofDrawnEvents = val;
+ }
+
+ /**
+ * Set an trigger on the RICH Hits.
+ */
+ void SetTriggerRichHits (Int_t val = 0){
+ fTriggerRichHits = val;
+ }
+
+private:
+
+ CbmDigiManager* fDigiMan = nullptr;
+
+ TClonesArray* fRichHits;
+
+ TClonesArray* fRichRings;
+
+ TClonesArray* fCbmEvent;
+
+ CbmHistManager* fHM;
+
+ Double_t fXOffsetHisto;
+
+ Int_t fEventNum;
+
+ Int_t fNofDrawnRings;
+
+ Int_t fNofDrawnRichTofEv;
+
+ Int_t fNofDrawnEvents;
+
+ Int_t fMaxNofDrawnEvents;
+
+ Int_t fTriggerRichHits;
+
+
+ string fOutputDir; // output dir for results
+
+ bool fRestrictToAcc = false;
+ bool fRestrictToFullAcc = false;
+
+ bool fDoWriteHistToFile = true;
+ bool fDoDrawCanvas = true;
+
+ std::array<Double_t,2304> ICD_offset_read;
+ std::array<Double_t,2304> ICD_offset;
+ std::array<uint32_t,2304> ICD_offset_cnt;
+
+ CbmRichMCbmSEDisplay* fSeDisplay = nullptr;
+
+
+ /**
+ * \brief Initialize histograms.
+ */
+ void InitHistograms();
+
+ /**
+ * \brief Draw histograms.
+ */
+ void DrawHist();
+
+ void RichRings();
+
+ void DrawEvent(CbmEvent *ev, std::vector<int> &ringIndx, bool full);
+
+ void DrawRing(CbmRichRing* ring) { DrawRing(ring, false); };
+
+ void DrawRing(CbmRichRing* ring, bool full);
+
+
+ bool doToT(CbmRichHit* hit);
+
+
+ void analyseRing(CbmRichRing *ring, CbmEvent *ev);
+
+ Bool_t cutRadius(CbmRichRing *ring);
+
+ void save_ICD (std::array<Double_t,2304> &offsets, unsigned int iteration);
+
+ void read_ICD (std::array<Double_t,2304> &offsets, unsigned int iteration);
+
+ /**
+ * \brief Copy constructor.
+ */
+ CbmRichMCbmQaRichOnly(const CbmRichMCbmQaRichOnly&);
+
+ /**
+ * \brief Assignment operator.
+ */
+ CbmRichMCbmQaRichOnly& operator=(const CbmRichMCbmQaRichOnly&);
+
+
+
+
+
+ ClassDef(CbmRichMCbmQaRichOnly,1)
+};
+
+#endif
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.cxx
new file mode 100644
index 00000000..9699f7a6
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.cxx
@@ -0,0 +1,282 @@
+#include "CbmRichMCbmSEDisplay.h"
+
+#include "TEllipse.h"
+#include "TLine.h"
+#include "TMarker.h"
+#include "TCanvas.h"
+#include <TBox.h>
+#include "TPad.h"
+#include "TH1D.h"
+#include "TClonesArray.h"
+
+#include "CbmRichHit.h"
+#include "CbmRichRing.h"
+#include "CbmTofTracklet.h"
+
+#include "TSystem.h"
+
+#include <string>
+#include <sstream>
+
+
+CbmRichMCbmSEDisplay::CbmRichMCbmSEDisplay() :
+ fRichHits(nullptr),
+ fRichRings(nullptr),
+ fTofTracks(nullptr),
+ fXOffsetHisto(0.),
+ fTotRichMin(23.7),
+ fTotRichMax(30.),
+ fNofDrawnEvents(0),
+ fMaxNofDrawnEvents(100),
+ fOutputDir("result"),
+ fHM(nullptr)
+{
+
+}
+
+CbmRichMCbmSEDisplay::CbmRichMCbmSEDisplay(CbmHistManager* manager) :
+ fRichHits(nullptr),
+ fRichRings(nullptr),
+ fTofTracks(nullptr),
+ fXOffsetHisto(0.),
+ fTotRichMin(23.7),
+ fTotRichMax(30.),
+ fNofDrawnEvents(0),
+ fMaxNofDrawnEvents(100),
+ fOutputDir("result"),
+ fHM(manager)
+{
+
+}
+
+
+void CbmRichMCbmSEDisplay::DrawEvent(CbmEvent *ev, std::vector<int> &ringIndx, bool full=true){
+
+ // ---- General Size of PMTs [cm] -------------------------------------------------
+ double pmtWidth = 5.20;
+ double pmtHeight = 5.20;
+ double pmtGap = 0.10;
+ // --------------------------------------------------------------------------------
+
+
+ // ---- Positioning ---------------------------------------------------------------
+ double left = 0.0;
+ double top = 0.0;
+ if (full == true){
+ left = -16.85 + fXOffsetHisto + 6.225;
+ top = 23.7;
+ }
+ // -------------------------------------------------------------------------------
+
+
+ // ---- Limiting Output ----------------------------------------------------------
+ if (fNofDrawnEvents > fMaxNofDrawnEvents) return;
+ // -------------------------------------------------------------------------------
+
+
+ // ---- Throw away too clean events ----------------------------------------------
+ if (ev->GetNofData(ECbmDataType::kRichHit) <= 4) return;
+ // -------------------------------------------------------------------------------
+
+
+ // ---- Init Canvas Layout for Event ---------------------------------------------
+ fNofDrawnEvents++;
+ std::string dir = fOutputDir+"/png/"+fFileName+"/";
+ gSystem->mkdir(dir.c_str(), true); // create directory if it does not exist
+ //std::cout<<"makeDir: "<<dir.c_str() << " "<<gSystem->mkdir(dir.c_str(), true)<<std::endl;
+
+ stringstream ss;
+ ss << fFileName << "/CbmEvent" << fNofDrawnEvents;
+ TCanvas *c = nullptr;
+ if (full == true) {
+ c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 600, 1250);
+ } else {
+ c = fHM->CreateCanvas(ss.str().c_str(), ss.str().c_str(), 800, 800);
+ }
+ c->SetGrid(true, true);
+ TPad* pad_event = new TPad("pad_event","event",0,0.20,1,1);
+ TH2D* pad = nullptr;
+ if (full == true){
+ pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -18. + fXOffsetHisto + 6.225, 5. + fXOffsetHisto + 6.225, 1, -26., 26.);
+ } else {
+ pad = new TH2D(ss.str().c_str(), (ss.str() + ";X [cm];Y [cm]").c_str(), 1, -15. + fXOffsetHisto + 6.225, 10. + fXOffsetHisto + 6.225, 1, -5., 20);
+ }
+
+ TPad* pad_time = new TPad("pad_time","timeDist",0,0,1,0.20);
+ TH1D* timeDistRichHit = new TH1D((ss.str() + "timeDistRichHit").c_str() , ";LE [ns];Entries",200,0.0,200. );
+ TH1D* timeDistRichHitToT = new TH1D((ss.str() + "timeDistRichHitToT").c_str(), ";LE [ns];Entries",200,0.0,200. );
+ TH1D* timeDistTofTrack = new TH1D((ss.str() + "timeDistTofTrack").c_str() , ";LE [ns];Entries",200,0.0,200. );
+ pad_event->Draw();
+ pad_time->Draw();
+ pad_event->cd();
+ //pad_event->SetBottomMargin(0);
+ pad->SetStats(false);
+ pad->Draw();
+
+ if (full == true){
+ for(unsigned int x=0;x<4;++x){
+ for(unsigned int y=0;y<9;++y){
+ double pmtLeft = left+(pmtWidth+pmtGap)*x;
+ double pmtTop = top-(pmtHeight+pmtGap)*y;
+ TBox *box = new TBox( pmtLeft,
+ pmtTop,
+ pmtLeft + pmtWidth,
+ pmtTop-pmtHeight
+ );
+ //box->SetFillColorAlpha(8,0.2);
+ box->Draw();
+
+ pmtLeft += 0.175;
+ pmtTop -= 0.175;
+ for (unsigned int pX=0; pX<8; ++pX){
+ for (unsigned int pY=0; pY<8; ++pY){
+ double xStart = 0.0;
+ double xEnd = 0.0;
+ double yStart = 0.0;
+ double yEnd = 0.0;
+ if (pX == 0){
+ xStart = pmtLeft;
+ xEnd = pmtLeft + 0.625;
+ } else if (pX < 7){
+ xStart = pmtLeft + 0.625 + 0.6*(pX-1);
+ xEnd = pmtLeft + 0.625 + 0.6*(pX);
+ } else if (pX == 7){
+ xStart = pmtLeft + 0.625 + 0.6*6;
+ xEnd = pmtLeft + 0.625*2 + 0.6*6;
+ }
+
+ if (pY == 0){
+ yStart = pmtTop;
+ yEnd = pmtTop - 0.625;
+ } else if (pY < 7){
+ yStart = pmtTop - 0.625 - 0.6*(pY-1);
+ yEnd = pmtTop - 0.625 - 0.6*(pY);
+ } else if (pY == 7){
+ yStart = pmtTop - 0.625 - 0.6*6;
+ yEnd = pmtTop - 0.625*2 - 0.6*6;
+ }
+
+ TBox *box1 = new TBox( xStart,
+ yStart,
+ xEnd,
+ yEnd
+ );
+ box1->SetLineWidth(1.);
+ //box1->SetFillColorAlpha(30.,0.1);
+ //box1->SetLineColorAlpha(30.,0.5);
+
+ box1->Draw("l");
+ }
+ }
+ }
+ }
+ }
+
+ if (full == true){
+ //rough Drawing of RichDetectorAcceptance
+ /*
+ TLine *line0 = new TLine( left , top, right, top);
+ line0->Draw();
+ TLine *line1 = new TLine( right, top, right, bottom);
+ line1->Draw();
+ TLine *line2 = new TLine( right, bottom, left, bottom);
+ line2->Draw();
+ TLine *line3 = new TLine( left , bottom, left, top);
+ line3->Draw();
+ */
+ }
+ // -------------------------------------------------------------------------------
+
+
+ // find min and max x and y positions of the hits
+ // in order to shift drawing
+
+
+ double hitZ = 0;
+ uint nofDrawHits = 0;
+
+ // ---- Draw Rich Hits -----------------------------------------------------------
+ for (int j=0;j<ev->GetNofData(ECbmDataType::kRichHit);j++){
+ auto iRichHit = ev->GetIndex(ECbmDataType::kRichHit, j);
+ CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iRichHit));
+ if (nullptr == hit) continue;
+ TEllipse* hitDr = new TEllipse(hit->GetX(), hit->GetY(), .25);
+ timeDistRichHit->Fill(hit->GetTime() - ev->GetStartTime());
+ if (doToT(hit)) { // Good ToT selection
+ hitDr->SetFillColor(kCyan);
+ timeDistRichHitToT->Fill(hit->GetTime() - ev->GetStartTime());
+ } else {
+ hitDr->SetFillColor(kBlue);
+ }
+ hitZ += hit->GetZ();
+ nofDrawHits++;
+ hitDr->Draw();
+ }
+ // -------------------------------------------------------------------------------
+
+
+ // ---- Draw Rich Rings and Ring Hits --------------------------------------------
+ for (unsigned int rings=0; rings < ringIndx.size();rings++) {
+ CbmRichRing *ring = static_cast<CbmRichRing *>(fRichRings->At(ringIndx[rings]));
+
+ //Draw Ring
+ TEllipse* circle = new TEllipse(ring->GetCenterX(), ring->GetCenterY(), ring->GetRadius());
+ circle->SetFillStyle(0);
+ circle->SetLineWidth(3);
+ circle->Draw();
+ TEllipse* center = new TEllipse(ring->GetCenterX(), ring->GetCenterY(), .1);
+ center->Draw();
+
+ // Draw Ring Hits
+ for (int i = 0; i < ring->GetNofHits(); i++){
+ Int_t hitInd = ring->GetHit(i);
+ CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
+ if (nullptr == hit) continue;
+ TEllipse* hitDr = new TEllipse(hit->GetX(), hit->GetY(), .125);
+ if (doToT(hit)) { // Good ToT selection
+ hitDr->SetFillColor(kMagenta);
+ } else {
+ hitDr->SetFillColor(kRed);
+ }
+ //hitZ += hit->GetZ();
+ //nofDrawHits++;
+ hitDr->Draw();
+ }
+ }
+ hitZ /= nofDrawHits;
+ // -------------------------------------------------------------------------------
+
+ // ---- Draw Tracks in RICH Plane ------------------------------------------------
+ auto nofTofTracks = ev->GetNofData(ECbmDataType::kTofTrack);
+ for (int j=0;j<nofTofTracks;j++){
+ auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, j);
+ CbmTofTracklet *track = static_cast<CbmTofTracklet *>(fTofTracks->At(iTofTrack));
+ if (nullptr == track) continue;
+ timeDistTofTrack->Fill(track->GetTime() - ev->GetStartTime());
+ TEllipse* hitDr = new TEllipse(track->GetFitX(hitZ), track->GetFitY(hitZ), .25);
+ hitDr->SetFillColor(kGreen);
+ hitDr->Draw();
+ }
+ // -------------------------------------------------------------------------------
+
+
+ // ---- Draw LE time distribution ------------------------------------------------
+ pad_time->cd();
+ pad_time->SetTitle("");
+ pad_time->SetTopMargin(0);
+ pad_time->SetBottomMargin(0.25);
+ timeDistRichHit->SetFillColor( kBlue);
+ timeDistRichHit->SetStats(false);
+ timeDistRichHit->GetXaxis()->SetLabelSize(0.06);
+ timeDistRichHit->GetXaxis()->SetTitleSize(0.08);
+ timeDistRichHit->GetYaxis()->SetLabelSize(0.06);
+ timeDistRichHit->GetYaxis()->SetTitleSize(0.08);
+ timeDistRichHit->GetYaxis()->SetTitleOffset(0.52);
+ timeDistRichHit->Draw("HIST");
+ timeDistRichHitToT->SetFillColor( kCyan);
+ timeDistRichHitToT->Draw("HIST SAME");
+ timeDistTofTrack->SetFillColor(kGreen);
+ timeDistTofTrack->Draw("HIST SAME");
+ // -------------------------------------------------------------------------------
+
+}
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.h b/reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.h
new file mode 100644
index 00000000..6966286b
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmSEDisplay.h
@@ -0,0 +1,165 @@
+#ifndef MCBM_RICH_SE_DISPLAY
+#define MCBM_RICH_SE_DISPLAY
+
+#include "CbmEvent.h"
+#include "CbmHistManager.h" // for ROOTCLING
+
+class TClonesArray;
+class CbmRichHit;
+class CbmRichRing;
+class CbmTofTracklet;
+
+#include <vector>
+
+using namespace std;
+
+
+
+class CbmRichMCbmSEDisplay
+{
+
+public:
+ /**
+ * \brief Standard constructor.
+ */
+ CbmRichMCbmSEDisplay();
+
+ /**
+ * \brief constructor with HistManager.
+ */
+ CbmRichMCbmSEDisplay(CbmHistManager * manager);
+
+ /**
+ * \brief Standard destructor.
+ */
+ virtual ~CbmRichMCbmSEDisplay() {};
+
+
+ /**
+ * \brief Draw histograms.
+ */
+ void DrawEvent(CbmEvent *ev, std::vector<int> &ringIndx, bool full);
+
+
+ /**
+ * Move X-Position of mRICH in Histograms (e.g. for Geometry changes)
+ */
+ void XOffsetHistos (Double_t val = 0.){
+ fXOffsetHisto = val;
+ }
+
+ /**
+ * Set an ToT cut of the RICH Hits.
+ */
+ void SetTotRich (Double_t min, Double_t max){
+ fTotRichMin = min;
+ fTotRichMax = max;
+ }
+
+ /**
+ * Set a pointer to the loaded Rich hits
+ */
+ void SetRichHits(TClonesArray* hits = nullptr){
+ fRichHits = hits;
+ }
+
+ /**
+ * Set a pointer to the loaded Rich hits
+ */
+ void SetRichRings(TClonesArray* ring = nullptr){
+ fRichRings = ring;
+ }
+
+ /**
+ * Set a pointer to the loaded Rich hits
+ */
+ void SetTofTracks(TClonesArray* track = nullptr){
+ fTofTracks = track;
+ }
+
+ /**
+ * Limit of Single Event Displays that should be drawn
+ */
+ void SetMaxNofDrawnEvents(Int_t val = 100){
+ fMaxNofDrawnEvents = val;
+ }
+
+
+ /**
+ * Limit of Single Event Displays that should be drawn
+ */
+ void SetHistmanager(CbmHistManager* manager){
+ fHM = manager;
+ }
+
+
+ /**
+ * Set the output directory of the Analysis
+ */
+ void SetOutDir(std::string dir){
+ fOutputDir = dir;
+ }
+
+ /**
+ * Set the output directory of the Canvases
+ */
+
+ void SetCanvasDir(std::string dir){
+ fFileName = dir;
+ }
+
+private:
+
+
+ TClonesArray* fRichHits;
+
+ TClonesArray* fRichRings;
+
+ TClonesArray* fTofTracks;
+
+ Double_t fXOffsetHisto;
+
+ Double_t fTotRichMin;
+
+ Double_t fTotRichMax;
+
+ Int_t fNofDrawnEvents;
+
+ Int_t fMaxNofDrawnEvents;
+
+
+
+ std::string fOutputDir; // output dir for results
+
+ std::string fFileName = "Ev";
+
+ CbmHistManager* fHM;
+
+
+ template<typename T=CbmRichHit>
+ bool doToT(T* hit){
+ bool check = false;
+ if ((hit->GetToT() > 23.7) && (hit->GetToT() < 30.0)) check = true;
+
+ return check;
+ }
+
+
+ /**
+ * \brief Copy constructor.
+ */
+ CbmRichMCbmSEDisplay(const CbmRichMCbmSEDisplay&);
+
+ /**
+ * \brief Assignment operator.
+ */
+ CbmRichMCbmSEDisplay& operator=(const CbmRichMCbmSEDisplay&);
+
+
+
+
+
+ ClassDef(CbmRichMCbmSEDisplay,1)
+};
+
+#endif
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.cxx
index 99abf87e..83d30f9a 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.cxx
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.cxx
@@ -55,7 +55,7 @@ InitStatus CbmRichMCbmToTShifter::Init()
fDigiMan = CbmDigiManager::Instance();
fDigiMan->Init();
- if (! fDigiMan->IsPresent(kRich) ) {
+ if (! fDigiMan->IsPresent(ECbmModuleId::kRich) ) {
Fatal("CbmRichMCbmToTShifter::Init", "No Rich Digis!");
}
@@ -67,7 +67,7 @@ InitStatus CbmRichMCbmToTShifter::Init()
void CbmRichMCbmToTShifter::Exec(Option_t* /*option*/)
{
fEventNum++;
- int nofDigis = fDigiMan->GetNofDigis(kRich);
+ int nofDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
for (int i=0; i < nofDigis; ++i){
const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(i);
@@ -93,12 +93,24 @@ void CbmRichMCbmToTShifter::Finish()
if (mean_cnt != 0) mean /= mean_cnt;
+ if (fShowTdcId) s << "TDC 0x"<<std::hex<<0xc000<< " " <<std::dec<<" !";
+ s << printEmpty() << " \\"<< std::endl;
+ if (fShowTdcId) s << "TDC 0x"<<std::hex<<0xc001<< " " <<std::dec<<" !";
+ s << printEmpty() << " \\"<< std::endl;
+
auto it = std::begin(fhTotMap);
+ uint32_t cnt = 0;
for (auto const &outer : fhTotMap) {
int tdc = outer.first;
TCanvas* c = new TCanvas(Form("fhToT_%x",outer.first),Form("fhToT_%x",outer.first), 2000 , 2000);
c->Divide(6,6);
- if (fShowTdcId) s << "TDC 0x"<<std::hex<<outer.first<<std::dec<<" !";
+ while (calcDirichAddr(cnt) < static_cast<uint16_t>(tdc)) { // this dirich is not in use; fill up the parameter file with printEmpty
+ if (fShowTdcId) s << "TDC 0x"<<std::hex<<calcDirichAddr(cnt)<< " " <<std::dec<<" !";
+ s << printEmpty();
+ if (std::next(it) != fhTotMap.end()) s<< " \\"<< std::endl;
+ cnt++;
+ }
+ if (fShowTdcId) s << "TDC 0x"<<std::hex<<outer.first<< " " <<std::dec<<" !";
s<<" 0.00";
for (int i=0;i<32;++i){
c->cd(1+i);
@@ -113,20 +125,32 @@ void CbmRichMCbmToTShifter::Finish()
}
}
+ cnt++;
if (it == fhTotMap.begin()) {
if (fGeneratePDFs) c->Print("ToTs/Tdc_all.pdf(","pdf");
s<< " \\"<< std::endl;
} else if (std::next(it) == fhTotMap.end()) {
+
if (fGeneratePDFs) c->Print("ToTs/Tdc_all.pdf)","pdf");
- s<<std::endl;
+ if (cnt == 71) s<<std::endl;
} else {
if (fGeneratePDFs) c->Print("ToTs/Tdc_all.pdf","pdf");
s<< " \\"<< std::endl;
}
++it;
+
+ }
+
+ //fill up till end of SetupSize
+
+ for (uint16_t i=cnt;i<72;++i){
+ s<< " \\"<< std::endl; // from last output
+
+ if (fShowTdcId) s << "TDC 0x"<<std::hex<<calcDirichAddr(i)<< " " <<std::dec<<" !";
+ s << printEmpty();
}
std::cout<<s.str()<<std::endl;
@@ -160,5 +184,14 @@ Double_t CbmRichMCbmToTShifter::GetMaxH1(TH1 * h)
return h->GetBinCenter(b);
}
+
+std::string CbmRichMCbmToTShifter::printEmpty(){
+ std::string s = "";
+ for (uint16_t i=0; i<33;++i){
+ s += " 0.00";
+ }
+ return s;
+}
+
ClassImp(CbmRichMCbmToTShifter)
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.h b/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.h
index 6ad7888a..f7a04de5 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.h
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmToTShifter.h
@@ -65,11 +65,11 @@ public:
private:
- CbmDigiManager* fDigiMan = nullptr;;
+ CbmDigiManager* fDigiMan = nullptr;
Int_t fEventNum;
- string fOutputDir; // output dir for results
+ std::string fOutputDir; // output dir for results
std::map<Int_t, std::map<Int_t, TH1*> > fhTotMap;
@@ -83,6 +83,15 @@ private:
TH1* GetTotH1(Int_t tdc, Int_t channel);
+ uint16_t calcDirichAddr( uint32_t cnt ){
+ return (0x7<<12) | ((cnt/18)<<8) | (((cnt%18)/2)<<4) | ((cnt%2)<<0);
+ };
+
+ /**
+ * \brief Fill output lines with 0's if DiRICh Address is not in use in input file
+ */
+ std::string printEmpty();
+
/**
* \brief Initialize histograms.
*/
--
GitLab