diff --git a/core/data/CbmDefs.h b/core/data/CbmDefs.h
index 0c3b1b6f141eac03f0df4d7e8d27dd64ecaf7958..6010035a717b43a99fd4b7b2807709fb76265967 100644
--- a/core/data/CbmDefs.h
+++ b/core/data/CbmDefs.h
@@ -9,9 +9,9 @@
#define CBMDEFS_H 1
#include <RtypesCore.h> // for Double_t
+#include <type_traits> // for underlying_type
-#include <iostream> // for ostream
-#include <type_traits> // for underlying_type
+#include <iostream> // for ostream
// Convert an element of enum class to its underlying intergral type
// since with C++11 the return type can't be deduced automatically it has
@@ -21,8 +21,8 @@
// E.g. ToIntegralType(ECbmModuleId::KSts) should be evaluated at compile
// time and should not affect the run time performance at all
template<typename T>
-constexpr auto ToIntegralType(T enumerator) ->
- typename std::underlying_type<T>::type {
+constexpr auto ToIntegralType(T enumerator) -> typename std::underlying_type<T>::type
+{
return static_cast<typename std::underlying_type<T>::type>(enumerator);
}
@@ -30,7 +30,8 @@ constexpr auto ToIntegralType(T enumerator) ->
** or passive (magnet, beam pipe, target etc.)
** In order to loop over all detectors, loop until kNofSystems.
**/
-enum class ECbmModuleId {
+enum class ECbmModuleId
+{
kRef = 0, ///< Reference plane
kMvd = 1, ///< Micro-Vertex Detector
kSts = 2, ///< Silicon Tracking System
@@ -73,7 +74,8 @@ std::ostream& operator<<(std::ostream&, const ECbmModuleId&);
/** Enumerator for CBM data types **/
-enum class ECbmDataType {
+enum class ECbmDataType
+{
kUnknown = -1,
kMCTrack = 0,
kMvdPoint = ToIntegralType(ECbmModuleId::kMvd) * 100,
@@ -88,7 +90,9 @@ enum class ECbmDataType {
kRichPoint = ToIntegralType(ECbmModuleId::kRich) * 100,
kRichDigi,
kRichHit,
- kRichRing, // RICH
+ kRichRing,
+ kRichTrackParamZ,
+ kRichTrackProjection, // RICH
kMuchPoint = ToIntegralType(ECbmModuleId::kMuch) * 100,
kMuchDigi,
kMuchCluster,
@@ -127,7 +131,12 @@ std::ostream& operator<<(std::ostream&, const ECbmDataType&);
** to first entry
** kRandom: Random choice of entries between first and last one.
**/
-enum class ECbmTreeAccess { kRegular, kRepeat, kRandom };
+enum class ECbmTreeAccess
+{
+ kRegular,
+ kRepeat,
+ kRandom
+};
/** Global functions for particle masses **/
diff --git a/macro/rich/run/run_digi.C b/macro/rich/run/run_digi.C
index 78721071d63b37607c4ba7828f1b0250c9a98d5c..d2be733b3582ccf4ff8c0a28b5802189bc4ddbc1 100644
--- a/macro/rich/run/run_digi.C
+++ b/macro/rich/run/run_digi.C
@@ -1,11 +1,8 @@
-void run_digi(
- const string& mcFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/mc.00000.root",
- const string& parFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/param.00000.root",
- const string& digiFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/digi.00000.root",
- int nEvents = 10) {
+void run_digi(const string& traFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/tra.00000.root",
+ const string& parFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/par.00000.root",
+ const string& digiFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/raw.00000.root",
+ int nEvents = 3)
+{
TTree::SetMaxTreeSize(90000000000);
FairLogger::GetLogger()->SetLogScreenLevel("INFO");
FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
@@ -21,11 +18,12 @@ void run_digi(
timer.Start();
CbmDigitization run;
- run.AddInput(mcFile.c_str(), eventRate);
+ run.AddInput(traFile.c_str(), eventRate);
run.SetOutputFile(digiFile.c_str(), overwrite);
run.SetParameterRootFile(parFile.c_str());
run.SetTimeSliceLength(timeSliceLength);
run.SetEventMode(eventMode);
+ run.SetProduceNoise(false);
run.SetMonitorFile("");
CbmRichDigitizer* richDigitizer = new CbmRichDigitizer();
@@ -39,7 +37,6 @@ void run_digi(
std::cout << "Macro finished successfully." << std::endl;
std::cout << "Digi file is " << digiFile << 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 << "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/run/run_reco.C b/macro/rich/run/run_reco.C
index 7aff89be46f1ae0b599a40768c972cd7885e2f7d..582a5c07a805261932aac5909539a3a6479891eb 100644
--- a/macro/rich/run/run_reco.C
+++ b/macro/rich/run/run_reco.C
@@ -1,26 +1,33 @@
-void run_reco(
- const string& mcFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/mc.00000.root",
- const string& parFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/param.00000.root",
- const string& digiFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/digi.00000.root",
- const string& recoFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/reco.00000.root",
- const string& geoSetup = "sis100_electron",
- const string& resultDir = "results_recoqa_newqa/",
- int nEvents = 100) {
+void run_reco(const string& traFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/tra.00000.root",
+ const string& parFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/par.00000.root",
+ const string& digiFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/raw.00000.root",
+ const string& recoFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/reco.00000.root",
+ const string& geoSetup = "sis100_electron", const string& resultDir = "results_recoqa_newqa/",
+ int nofTimeSlices = 3)
+{
+
TTree::SetMaxTreeSize(90000000000);
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ remove(recoFile.c_str());
- TString myName = "run_reco";
TString srcDir = gSystem->Getenv("VMCWORKDIR");
- remove(recoFile.c_str());
+ CbmSetup* geo = CbmSetup::Instance();
+ geo->LoadSetup(geoSetup.c_str());
- CbmSetup::Instance()->LoadSetup(geoSetup.c_str());
+ TString sEvBuildRaw = "";
+ Bool_t useMC = true;
+
+ Bool_t eventBased = (sEvBuildRaw != "");
+ Bool_t useMvd = geo->IsActive(ECbmModuleId::kMvd);
+ Bool_t useSts = geo->IsActive(ECbmModuleId::kSts);
+ Bool_t useRich = geo->IsActive(ECbmModuleId::kRich);
+ Bool_t useMuch = geo->IsActive(ECbmModuleId::kMuch);
+ Bool_t useTrd = geo->IsActive(ECbmModuleId::kTrd);
+ Bool_t useTof = geo->IsActive(ECbmModuleId::kTof);
+ Bool_t usePsd = geo->IsActive(ECbmModuleId::kPsd);
- std::cout << std::endl
- << "-I- " << myName << ": Defining parameter files " << std::endl;
TList* parFileList = new TList();
TString geoTag;
@@ -29,68 +36,260 @@ void run_reco(
const Char_t* npar[4] = {"asic", "digi", "gas", "gain"};
TObjString* trdParFile(NULL);
for (Int_t i(0); i < 4; i++) {
- trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "."
- + npar[i] + ".par");
+ trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "." + npar[i] + ".par");
parFileList->Add(trdParFile);
- std::cout << "-I- " << myName << ": Using parameter file "
- << trdParFile->GetString() << std::endl;
+ std::cout << "-I- "
+ << ": Using parameter file " << trdParFile->GetString() << std::endl;
}
}
// - TOF digitisation parameters
if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
- TObjString* tofFile =
- new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par");
- parFileList->Add(tofFile);
- std::cout << "-I- " << myName << ": Using parameter file "
- << tofFile->GetString() << std::endl;
- TObjString* tofBdfFile =
- new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
+ 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;
+ std::cout << "-I- "
+ << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
}
TStopwatch timer;
timer.Start();
- gDebug = 0;
-
FairRunAna* run = new FairRunAna();
FairFileSource* inputSource = new FairFileSource(digiFile.c_str());
- inputSource->AddFriend(mcFile.c_str());
+ if (useMC) { inputSource->AddFriend(traFile.c_str()); }
run->SetSource(inputSource);
run->SetOutputFile(recoFile.c_str());
- run->SetGenerateRunInfo(kTRUE);
+ run->SetGenerateRunInfo(true);
+ if (useMC) {
+ CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
+ mcManager->AddFile(traFile.c_str());
+ run->AddTask(mcManager);
+ }
+
+ if (eventBased) {
+ if (sEvBuildRaw.EqualTo("Ideal", TString::ECaseCompare::kIgnoreCase)) {
+ FairTask* evBuildRaw = new CbmBuildEventsIdeal();
+ run->AddTask(evBuildRaw);
+ std::cout << "-I- "
+ << ": Added task " << evBuildRaw->GetName() << std::endl;
+ eventBased = kTRUE;
+ }
+ else if (sEvBuildRaw.EqualTo("Real", TString::ECaseCompare::kIgnoreCase)) {
+ CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents();
+
+ //Choose between NoOverlap, MergeOverlap, AllowOverlap
+ evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap);
+
+ // Remove detectors where digis not found
+ if (!useRich) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
+ if (!useMuch) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch);
+ if (!usePsd) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd);
+ if (!useTof) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof);
+ if (!useTrd) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd);
+ if (!useSts) {
+ std::cerr << "-E- "
+ << ": Sts must be present for raw event "
+ << "building using ``Real2019'' option. Terminating macro." << std::endl;
+ return;
+ }
+ // Set STS as reference detector
+ evBuildRaw->SetReferenceDetector(kRawEventBuilderDetSts);
+ evBuildRaw->SetTsParameters(0.0, 1.e7, 0.0);
+
+ // Use CbmMuchDigi instead of CbmMuchBeamtimeDigi
+ evBuildRaw->ChangeMuchBeamtimeDigiFlag(kFALSE);
+
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1000);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -500, 500);
+
+ run->AddTask(evBuildRaw);
+ std::cout << "-I- "
+ << ": Added task " << evBuildRaw->GetName() << std::endl;
+ eventBased = kTRUE;
+ } //? Real raw event building
+ else {
+ std::cerr << "-E- "
+ << ": Unknown option " << sEvBuildRaw << " for raw event building! Terminating macro execution."
+ << std::endl;
+ return;
+ }
+ }
- FairLogger::GetLogger()->SetLogScreenLevel("INFO");
- FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ if (eventBased && useMC) {
+ CbmBuildEventsQA* evBuildQA = new CbmBuildEventsQA();
+ run->AddTask(evBuildQA);
+ }
- CbmMCDataManager* mcManager = new CbmMCDataManager("MCManager", 1);
- mcManager->AddFile(mcFile.c_str());
- run->AddTask(mcManager);
+ if (useMvd) {
- // Reconstruction tasks
- TString macroName = srcDir + "/macro/run/modules/reconstruct.C";
- std::cout << "Loading macro " << macroName << std::endl;
- gROOT->LoadMacro(macroName);
- Bool_t recoSuccess = gROOT->ProcessLine("reconstruct(true,true)");
- if (!recoSuccess) {
- std::cerr << "-E-" << myName << ": error in executing " << macroName
- << std::endl;
- return;
+ CbmMvdClusterfinder* mvdCluster = new CbmMvdClusterfinder("MVD Cluster Finder", 0, 0);
+ run->AddTask(mvdCluster);
+ std::cout << "-I- : Added task " << mvdCluster->GetName() << std::endl;
+
+ CbmMvdHitfinder* mvdHit = new CbmMvdHitfinder("MVD Hit Finder", 0, 0);
+ mvdHit->UseClusterfinder(kTRUE);
+ run->AddTask(mvdHit);
+ std::cout << "-I- "
+ << ": Added task " << mvdHit->GetName() << std::endl;
}
- std::cout << "-I-" << myName << ": " << macroName << " excuted successfully"
- << std::endl;
- CbmMatchRecoToMC* matchRecoToMc = new CbmMatchRecoToMC();
- run->AddTask(matchRecoToMc);
- std::cout << std::endl
- << std::endl
- << "-I- " << myName << ": Set runtime DB" << std::endl;
+ if (useSts) {
+ CbmRecoSts* stsReco = new CbmRecoSts(kCbmRecoTimeslice);
+ if (eventBased) stsReco->SetMode(kCbmRecoEvent);
+ run->AddTask(stsReco);
+ std::cout << "-I- "
+ << ": Added task " << stsReco->GetName() << std::endl;
+ }
+
+
+ if (useRich) {
+ CbmRichHitProducer* richHitProd = new CbmRichHitProducer();
+ run->AddTask(richHitProd);
+ std::cout << "-I- "
+ << ": Added task " << richHitProd->GetName() << std::endl;
+ }
+
+
+ if (useTrd) {
+
+ Double_t triggerThreshold = 0.5e-6; // SIS100
+ CbmTrdClusterFinder* trdCluster = new CbmTrdClusterFinder();
+ if (eventBased) trdCluster->SetTimeBased(kFALSE);
+ else
+ trdCluster->SetTimeBased(kTRUE);
+ trdCluster->SetNeighbourEnable(true, false);
+ trdCluster->SetMinimumChargeTH(triggerThreshold);
+ trdCluster->SetRowMerger(true);
+
+ // Uncomment if you want to use all available digis.
+ // In that case clusters hits will not be added to the CbmEvent
+ // trdCluster->SetUseOnlyEventDigis(kFALSE);
+
+ run->AddTask(trdCluster);
+ std::cout << "-I- "
+ << ": Added task " << trdCluster->GetName() << std::endl;
+
+ CbmTrdHitProducer* trdHit = new CbmTrdHitProducer();
+ run->AddTask(trdHit);
+ std::cout << "-I- "
+ << ": Added task " << trdHit->GetName() << std::endl;
+ }
+
+ if (useTof) {
+ CbmTofSimpClusterizer* tofCluster = new CbmTofSimpClusterizer("TofSimpClusterizer", 0);
+ tofCluster->SetOutputBranchPersistent("TofHit", kTRUE);
+ tofCluster->SetOutputBranchPersistent("TofDigiMatch", kTRUE);
+ run->AddTask(tofCluster);
+ std::cout << "-I- "
+ << ": Added task " << tofCluster->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+ if (useMC) {
+ CbmMatchRecoToMC* match1 = new CbmMatchRecoToMC();
+ run->AddTask(match1);
+ }
+
+
+ if (useMvd || useSts) {
+ CbmKF* kalman = new CbmKF();
+ run->AddTask(kalman);
+ CbmL1* l1 = 0;
+ if (useMC) { l1 = new CbmL1("L1", 2, 3); }
+ else {
+ l1 = new CbmL1("L1", 0);
+ }
+ l1->SetDataMode(!eventBased);
+
+ // --- Material budget file names
+ TString mvdGeoTag;
+ if (geo->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile += "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetMvdMaterialBudgetFileName(parFile.Data());
+ }
+ TString stsGeoTag;
+ if (geo->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile += "/parameters/sts/sts_matbudget_" + stsGeoTag + ".root";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetStsMaterialBudgetFileName(parFile.Data());
+ }
+ run->AddTask(l1);
+ std::cout << "-I- "
+ << ": Added task " << l1->GetName() << std::endl;
+
+ CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder();
+ if (eventBased) {
+ FairTask* stsFindTracks = new CbmStsFindTracksEvents(stsTrackFinder, useMvd);
+ run->AddTask(stsFindTracks);
+ std::cout << "-I- "
+ << ": Added task " << stsFindTracks->GetName() << std::endl;
+ }
+ else {
+ FairTask* stsFindTracks = new CbmStsFindTracks(0, stsTrackFinder, useMvd);
+ run->AddTask(stsFindTracks);
+ std::cout << "-I- "
+ << ": Added task " << stsFindTracks->GetName() << std::endl;
+ }
+ }
+ // ------------------------------------------------------------------------
+
+
+ // --- Global track finding -----------------------------------------
+ CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
+ finder->SetTrackingType("branch");
+ finder->SetMergerType("nearest_hit");
+ run->AddTask(finder);
+ std::cout << "-I- : Added task " << finder->GetName() << std::endl;
+ // ----------------------------------------------------------------------
+
+
+ // ----- RICH reconstruction ----------------------------------------
+ if (useRich) {
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ run->AddTask(richReco);
+ std::cout << "-I- : Added task " << richReco->GetName() << std::endl;
+ }
+ // ----------------------------------------------------------------------
+
+ // ==== From here on, the time-based and the event-based reconstruction
+ // ==== chains differ, since time-based version of primary vertex finding
+ // ==== and global tracking are not yet available. For time-based
+ // ==== reconstruction, a track-based event finder is used; no global
+ // ==== tracks are produced.
+
+ if (eventBased) {
+
+ //----- Primary vertex finding -------------------------------------
+ CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
+ CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
+ run->AddTask(findVertex);
+ std::cout << "-I- "
+ << ": Added task " << findVertex->GetName() << std::endl;
+ // ----------------------------------------------------------------------
+
+
+ } //? event-based reco
+
+ else {
+
+ // -----Â Event building from STS tracks -----------------------------
+ run->AddTask(new CbmBuildEventsFromTracksReal());
+ // ----------------------------------------------------------------------
+
+ } //? time-based reco
+
+
+ // ----- Parameter database --------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- "
+ << ": Set runtime DB" << std::endl;
FairRuntimeDb* rtdb = run->GetRuntimeDb();
FairParRootFileIo* parIo1 = new FairParRootFileIo();
FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
@@ -100,23 +299,41 @@ void run_reco(
parIo2->open(parFileList, "in");
rtdb->setSecondInput(parIo2);
}
+ // ------------------------------------------------------------------------
- std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl;
- run->Init();
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- "
+ << ": Initialise run" << std::endl;
+ run->Init();
rtdb->setOutput(parIo1);
rtdb->saveOutput();
rtdb->print();
+ // ------------------------------------------------------------------------
+
- std::cout << "-I- " << myName << ": Starting run" << std::endl;
- run->Run(0, nEvents);
+ // ----- Register light ions (d, t, He3, He4) -------------------------
+ std::cout << std::endl;
+ TString registerLightIonsMacro = gSystem->Getenv("VMCWORKDIR");
+ registerLightIonsMacro += "/macro/KF/registerLightIons.C";
+ std::cout << "Loading macro " << registerLightIonsMacro << std::endl;
+ gROOT->LoadMacro(registerLightIonsMacro);
+ gROOT->ProcessLine("registerLightIons()");
+ // ------------------------------------------------------------------------
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- "
+ << ": Starting run" << std::endl;
+ run->Run(0, nofTimeSlices);
+ // ------------------------------------------------------------------------
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 << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl;
std::cout << std::endl << "Test passed" << std::endl << "All ok" << std::endl;
}
diff --git a/macro/rich/run/run_sim.C b/macro/rich/run/run_transport.C
similarity index 57%
rename from macro/rich/run/run_sim.C
rename to macro/rich/run/run_transport.C
index e71783ce444a024e6262921da388b85dbea859d7..a8102cf0ebada9bf7704c019d56171a474a9f0e7 100644
--- a/macro/rich/run/run_sim.C
+++ b/macro/rich/run/run_transport.C
@@ -1,28 +1,20 @@
-void run_sim(
- const string& urqmdFile =
- "/Users/slebedev/Development/cbm/data/urqmd/auau/8gev/centr/"
- "urqmd.auau.8gev.centr.00001.root", // if "", no urqmd
- const string& mcFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/mc.00000.root",
- const string& parFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/param.00000.root",
- const string& geoFile =
- "/Users/slebedev/Development/cbm/data/sim/rich/reco/geosim.00000.root",
- int nofElectrons =
- 5, // number of e- to be generated, if <=0 no electrons are embedded into event
- int nofPositrons =
- 5, // number of e- to be generated, if <=0 no positrons are embedded into event
- const string& plutoFile =
- "", // if "", no pluto particles are embedded into event
- const string& geoSetup = "sis100_electron_rich_pal_bcarb",
- int nEvents = 10) {
+void run_transport(const string& urqmdFile = "/Users/slebedev/Development/cbm/data/urqmd/auau/8gev/centr/"
+ "urqmd.auau.8gev.centr.00001.root", // if "", no urqmd
+ const string& traFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/tra.00000.root",
+ const string& parFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/par.00000.root",
+ const string& geoFile = "/Users/slebedev/Development/cbm/data/sim/rich/reco/geo.00000.root",
+ int nofElectrons = 5, // number of e- to be generated, if <=0 no e- are embedded into event
+ int nofPositrons = 5, // number of e+ to be generated, if <=0 no e+ are embedded into event
+ const string& plutoFile = "", // if "", no pluto particles are embedded into event
+ const string& geoSetup = "sis100_electron", int nEvents = 3)
+{
TTree::SetMaxTreeSize(90000000000);
FairLogger::GetLogger()->SetLogScreenLevel("INFO");
FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
remove(parFile.c_str());
- remove(mcFile.c_str());
+ remove(traFile.c_str());
remove(geoFile.c_str());
TStopwatch timer;
@@ -53,22 +45,22 @@ void run_sim(
if (plutoFile.length() > 0) { run.AddInput(plutoFile.c_str(), kPluto); }
- run.SetOutFileName(mcFile.c_str());
+ run.SetOutFileName(traFile.c_str());
run.SetParFileName(parFile.c_str());
run.SetGeoFileName(geoFile.c_str());
run.LoadSetup(geoSetup.c_str());
run.SetTarget("Gold", 0.025, 2.5);
run.SetBeamPosition(0., 0., 0.1, 0.1);
+ run.SetRandomEventPlane();
run.Run(nEvents);
timer.Stop();
std::cout << std::endl << std::endl;
std::cout << "Macro finished successfully." << std::endl;
- std::cout << "Output file is " << mcFile << std::endl;
+ std::cout << "Output file is " << traFile << std::endl;
std::cout << "Parameter file is " << parFile << std::endl;
std::cout << "Geometry file is " << geoFile << std::endl;
- std::cout << "Real time " << timer.RealTime() << " s, CPU time "
- << timer.CpuTime() << "s" << 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/L1/OffLineInterface/CbmL1RichENNRingFinder.cxx b/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.cxx
index 56e934662c07bc32a97f58ca3de24f5b6b549d1d..d9a581c230e4a1f2b4df88d0b8444c0dbe566f6e 100644
--- a/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.cxx
+++ b/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.cxx
@@ -17,20 +17,23 @@
// CBM includes
#include "CbmL1RichENNRingFinder.h"
-#include "CbmL1RichENNRingFinderParallel.h"
+#include "CbmEvent.h"
+#include "CbmL1RichENNRingFinderParallel.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "TClonesArray.h"
#include "TStopwatch.h"
-#include "assert.h"
#include <algorithm>
-#include <cmath>
#include <iostream>
#include <vector>
+#include <cmath>
+
+#include "assert.h"
+
using std::cout;
using std::endl;
using std::fabs;
@@ -44,15 +47,17 @@ ClassImp(CbmL1RichENNRingFinder)
CbmL1RichENNRingFinder::CbmL1RichENNRingFinder(Int_t verbose)
: finder(new CbmL1RichENNRingFinderParallel(verbose))
, fRecoTime(0)
- , fNEvents(0) {}
+ , fNEvents(0)
+{
+}
CbmL1RichENNRingFinder::~CbmL1RichENNRingFinder() {}
void CbmL1RichENNRingFinder::Init() {}
-Int_t CbmL1RichENNRingFinder::DoFind(TClonesArray* HitArray,
- TClonesArray* ProjArray,
- TClonesArray* RingArray) {
- return finder->DoFind(HitArray, ProjArray, RingArray);
+Int_t CbmL1RichENNRingFinder::DoFind(CbmEvent* event, TClonesArray* HitArray, TClonesArray* ProjArray,
+ TClonesArray* RingArray)
+{
+ return finder->DoFind(event, HitArray, ProjArray, RingArray);
}
diff --git a/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.h b/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.h
index a9eadb3cedd56d830eb479d070b1f720f3079af5..a540b60cbafd91594cca2a360115a67cbe4ff4e4 100644
--- a/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.h
+++ b/reco/L1/OffLineInterface/CbmL1RichENNRingFinder.h
@@ -51,9 +51,7 @@ public:
**
*@value Number of tracks created
**/
- Int_t DoFind(TClonesArray* hitArray,
- TClonesArray* projArray,
- TClonesArray* ringArray);
+ Int_t DoFind(CbmEvent* event, TClonesArray* hitArray, TClonesArray* projArray, TClonesArray* ringArray);
private:
CbmL1RichENNRingFinderParallel* finder;
diff --git a/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx b/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx
index a0bf8cee148cd0017ef0135608ce75c4cfa096c4..bb42de5bbda739f5a1a7d54dc968d8cca9d9edad 100644
--- a/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx
+++ b/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.cxx
@@ -19,18 +19,21 @@
// CBM includes
#include "CbmL1RichENNRingFinderParallel.h"
+#include "CbmEvent.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "TClonesArray.h"
#include "TStopwatch.h"
-#include "assert.h"
#include <algorithm>
-#include <cmath>
#include <iostream>
#include <vector>
+#include <cmath>
+
+#include "assert.h"
+
using std::cout;
using std::endl;
using std::fabs;
@@ -38,9 +41,8 @@ using std::ios;
using std::sqrt;
using std::vector;
-CbmL1RichENNRingFinderParallel::CbmL1RichENNRingFinderParallel(Int_t verbose)
- : fRecoTime(0), fNEvents(0) {
- fVerbose = verbose;
+CbmL1RichENNRingFinderParallel::CbmL1RichENNRingFinderParallel(Int_t verbose) : fRecoTime(0), fNEvents(0)
+{
#ifdef PRINT_TIMING
TString name_tmp[NTimers] = {
@@ -69,9 +71,9 @@ CbmL1RichENNRingFinderParallel::~CbmL1RichENNRingFinderParallel() {}
void CbmL1RichENNRingFinderParallel::Init() {}
-Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
- TClonesArray* /*ProjArray*/,
- TClonesArray* RingArray) {
+Int_t CbmL1RichENNRingFinderParallel::DoFind(CbmEvent* event, TClonesArray* HitArray, TClonesArray* /*ProjArray*/,
+ TClonesArray* RingArray)
+{
if (!RingArray || !HitArray) return 0;
RingArray->Clear();
@@ -87,9 +89,9 @@ Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
vector<ENNRingHit> Up;
vector<ENNRingHit> Down;
- const Int_t nhits = HitArray->GetEntriesFast();
-
- for (Int_t i = 0; i < nhits; ++i) {
+ const Int_t nhits = event ? event->GetNofData(ECbmDataType::kRichHit) : HitArray->GetEntriesFast();
+ for (Int_t i0 = 0; i0 < nhits; ++i0) {
+ Int_t i = event ? event->GetIndex(ECbmDataType::kRichHit, i0) : i0;
CbmRichHit* hit = L1_DYNAMIC_CAST<CbmRichHit*>(HitArray->At(i));
if (!hit) continue;
ENNRingHit tmp;
@@ -97,9 +99,8 @@ Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
tmp.y = hit->GetY();
tmp.outIndex = i;
tmp.quality = 0;
- if (tmp.y > 0.) {
- Up.push_back(tmp);
- } else {
+ if (tmp.y > 0.) { Up.push_back(tmp); }
+ else {
Down.push_back(tmp);
}
}
@@ -178,6 +179,7 @@ Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
for (vector<ENNRing>::iterator i = R.begin(); i != R.end(); ++i) {
if (!i->on) continue;
new ((*RingArray)[NRings]) CbmRichRing();
+ if (event != nullptr) event->AddData(ECbmDataType::kRichRing, NRings);
CbmRichRing* ring = L1_DYNAMIC_CAST<CbmRichRing*>(RingArray->At(NRings));
NRings++;
ring->SetCenterX(i->x);
@@ -186,8 +188,7 @@ Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
ring->SetChi2(i->chi2);
const THitIndex NHits = i->localIHits.size();
for (THitIndex j = 0; j < NHits; j++) {
- if (i->y > 0.)
- ring->AddHit(outIndicesUp[i->localIHits[j]]);
+ if (i->y > 0.) ring->AddHit(outIndicesUp[i->localIHits[j]]);
else
ring->AddHit(outIndicesDown[i->localIHits[j]]);
}
@@ -201,8 +202,8 @@ Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
cout.setf(ios::fixed);
cout.setf(ios::showpoint);
cout.precision(4);
- cout << "L1ENN Reco time stat/this [ms] : " << fRecoTime * 1000. / fNEvents
- << "/" << timer.RealTime() * 1000. << endl;
+ cout << "L1ENN Reco time stat/this [ms] : " << fRecoTime * 1000. / fNEvents << "/" << timer.RealTime() * 1000.
+ << endl;
#ifdef PRINT_TIMING
static float timeStat[NTimers];
@@ -221,23 +222,18 @@ Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
for (int i = 0; i < NTimers; i++) {
timeStat[i] += fTimers[i].RealTime();
cout.width(30);
- cout << fTimersNames[i] << " timer = " << timeStat[i] * 1000. / fNEvents
- << " / " << fTimers[i].RealTime() * 1000. << " | "
- << timeStat[i] / timeStat[1] * 100 << endl;
+ cout << fTimersNames[i] << " timer = " << timeStat[i] * 1000. / fNEvents << " / " << fTimers[i].RealTime() * 1000.
+ << " | " << timeStat[i] / timeStat[1] * 100 << endl;
}
#endif // PRINT_TIMING
return NRings;
}
-void CbmL1RichENNRingFinderParallel::ENNRingFinder(
- const int NHits,
- nsL1vector<ENNHitV>::TSimd& HitsV,
- vector<ENNRing>& Rings,
- float HitSize,
- THitIndex MinRingHits,
- fvec /*RMin*/,
- fvec RMax) {
+void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, nsL1vector<ENNHitV>::TSimd& HitsV,
+ vector<ENNRing>& Rings, float HitSize, THitIndex MinRingHits,
+ fvec /*RMin*/, fvec RMax)
+{
#ifdef PRINT_TIMING
GetTimer("All").Start(0);
#endif // PRINT_TIMING
@@ -289,8 +285,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
THitIndex i_mains[fvecLen] = {0};
- THitIndex i_main_array
- [NHits]; // need for proceed in paralled almost independent areas
+ THitIndex i_main_array[NHits]; // need for proceed in paralled almost independent areas
for (THitIndex i = 0; i < NHits; i++) {
i_main_array[i] = i; // TODO
}
@@ -315,8 +310,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
const THitIndex i_main = i_main_array[ii_main];
const int i_main_4 = i_main % fvecLen, i_main_V = i_main / fvecLen;
ENNHitV* i = &HitsV[i_main_V];
- if (i->quality[i_main_4] >= StartHitMaxQuality)
- continue; // already found hit
+ if (i->quality[i_main_4] >= StartHitMaxQuality) continue; // already found hit
i_mains[i_4] = i_main;
@@ -328,8 +322,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
// for( int j = ileft[i_4]; j < iright[i_4]; ++j ){
while ((HitsV[ileft / fvecLen].x[ileft % fvecLen] < left))
++ileft; // TODO SIMDize
- while ((iright < NHits)
- && (HitsV[iright / fvecLen].x[ileft % fvecLen] < right))
+ while ((iright < NHits) && (HitsV[iright / fvecLen].x[ileft % fvecLen] < right))
++iright;
for (int j = ileft; j < iright; ++j) {
@@ -347,10 +340,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
if (sHit.lr2[i_4] > AreaSize2) continue;
if (ISUNLIKELY(j == i_main)) continue;
- if (sHit.quality[i_4]
- >= SearchHitMaxQuality) { // CHECKME do we really need PickUpArea
- PickUpArea[static_cast<int>(PickUpAreaSize[i_4]++)].CopyHit(
- HitsV[j_V], j_4, i_4);
+ if (sHit.quality[i_4] >= SearchHitMaxQuality) { // CHECKME do we really need PickUpArea
+ PickUpArea[static_cast<int>(PickUpAreaSize[i_4]++)].CopyHit(HitsV[j_V], j_4, i_4);
continue;
}
@@ -367,12 +358,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
int MaxPickUpAreaSize = 0;
for (int i_4 = 0; i_4 < fvecLen; i_4++) {
iHit.CopyHit(HitsV[i_mains[i_4] / fvecLen], i_mains[i_4] % fvecLen, i_4);
- MaxSearchAreaSize = (MaxSearchAreaSize < SearchAreaSize[i_4])
- ? int(SearchAreaSize[i_4])
- : MaxSearchAreaSize;
- MaxPickUpAreaSize = (MaxPickUpAreaSize < PickUpAreaSize[i_4])
- ? int(PickUpAreaSize[i_4])
- : MaxPickUpAreaSize;
+ MaxSearchAreaSize = (MaxSearchAreaSize < SearchAreaSize[i_4]) ? int(SearchAreaSize[i_4]) : MaxSearchAreaSize;
+ MaxPickUpAreaSize = (MaxPickUpAreaSize < PickUpAreaSize[i_4]) ? int(PickUpAreaSize[i_4]) : MaxPickUpAreaSize;
}
#ifdef PRINT_TIMING
GetTimer("Ring finding: Prepare data").Stop();
@@ -381,9 +368,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
#ifdef PRINT_TIMING
GetTimer("Ring finding: Init of param").Start(0);
-#endif // PRINT_TIMING
- for (int isa = 0; isa < MaxSearchAreaSize;
- isa++) { // TODO don't work w\o this because of nan in wights
+#endif
+ for (int isa = 0; isa < MaxSearchAreaSize; isa++) { // TODO don't work w\o this because of nan in wights
ENNSearchHitV& sHit = SearchArea[isa];
const fvec validHit = (fvec(isa) < SearchAreaSize) & validRing;
sHit.lx = if3(validHit, sHit.lx, 0);
@@ -453,13 +439,10 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
const fvec d = fabs(sqrt(dx * dx + dy * dy) - R);
const fvec dSq = d * d;
sHit.on_ring = (d <= HitSize) & validHit;
- const fvec dp =
- if3(sHit.on_ring, -1, fabs(sHit.C + sHit.Cx * X + sHit.Cy * Y));
- Dmax = if3(((dp <= Dcut) & (dp > Dmax)), dp, Dmax);
+ const fvec dp = if3(sHit.on_ring, -1, fabs(sHit.C + sHit.Cx * X + sHit.Cy * Y));
+ Dmax = if3(((dp <= Dcut) & (dp > Dmax)), dp, Dmax);
- fvec w = if3((sHit.on_ring),
- 1. / (HitSizeSq_v + fabs(dSq)),
- 1. / (1.e-5 + fabs(dSq)));
+ fvec w = if3((sHit.on_ring), 1. / (HitSizeSq_v + fabs(dSq)), 1. / (1.e-5 + fabs(dSq)));
w = if3((dp <= Dcut) & validHit, w, 0);
S0 += w * sHit.S0;
S1 += w * sHit.S1;
@@ -507,8 +490,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
R = sqrt(R2);
} // end of the final fit
- validRing = !((NRingHits < MinRingHits_v) | (R2 > R2MaxCut)
- | (R2 < R2MinCut)); // ghost suppresion // TODO constants
+ validRing =
+ !((NRingHits < MinRingHits_v) | (R2 > R2MaxCut) | (R2 < R2MinCut)); // ghost suppresion // TODO constants
// cout << validRing << endl;
#ifdef PRINT_TIMING
GetTimer("Ring finding: Final fit").Stop();
@@ -612,8 +595,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
ring.localIHits.push_back(static_cast<THitIndex>(iHit.localIndex[i_4]));
ring.NHits = 1;
ring.chi2 = 0;
- vector<THitIndex>
- Shadow; // save hit indices of hits, which's quality will be changed
+ vector<THitIndex> Shadow; // save hit indices of hits, which's quality will be changed
Shadow.reserve(15);
Shadow.push_back(static_cast<THitIndex>(iHit.localIndex[i_4]));
for (THitIndex ih = 0; ih < SearchAreaSize[i_4]; ih++) {
@@ -625,8 +607,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
ring.chi2 += d * d;
ring.localIHits.push_back(int(sHit.localIndex[i_4]));
ring.NHits++;
- if (d <= ShadowSize)
- Shadow.push_back(static_cast<THitIndex>(sHit.localIndex[i_4]));
+ if (d <= ShadowSize) Shadow.push_back(static_cast<THitIndex>(sHit.localIndex[i_4]));
}
}
for (int ipu = 0; ipu < PickUpAreaSize[i_4]; ipu++) {
@@ -636,11 +617,9 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
const float d = fabs(sqrt(dx * dx + dy * dy) - ring.r);
if (ISLIKELY(d <= HitSize)) {
ring.chi2 += d * d;
- ring.localIHits.push_back(
- static_cast<THitIndex>(puHit.localIndex[i_4]));
+ ring.localIHits.push_back(static_cast<THitIndex>(puHit.localIndex[i_4]));
ring.NHits++;
- if (d <= ShadowSize)
- Shadow.push_back(static_cast<THitIndex>(puHit.localIndex[i_4]));
+ if (d <= ShadowSize) Shadow.push_back(static_cast<THitIndex>(puHit.localIndex[i_4]));
}
}
if (ISUNLIKELY(ring.NHits < MinRingHits)) {
@@ -661,8 +640,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
const THitIndex ih_4 = ih % fvecLen;
ENNHitV& hitV = HitsV[ih / fvecLen];
- hitV.quality[ih_4] =
- (hitV.quality[ih_4] < quality) ? quality : hitV.quality[ih_4];
+ hitV.quality[ih_4] = (hitV.quality[ih_4] < quality) ? quality : hitV.quality[ih_4];
// shHit->quality = if3( shHit->quality < quality, quality, shHit->quality );
}
} // i_4
@@ -701,9 +679,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
for (iR j = i + 1; j != Rend; ++j) {
if (j->skip) continue;
- float dist = j->r + i->r + HitSize2[0];
- float distCentr =
- sqrt((j->x - i->x) * (j->x - i->x) + (j->y - i->y) * (j->y - i->y));
+ float dist = j->r + i->r + HitSize2[0];
+ float distCentr = sqrt((j->x - i->x) * (j->x - i->x) + (j->y - i->y) * (j->y - i->y));
if (distCentr > dist) continue;
Int_t NOverlaped = 0;
@@ -772,8 +749,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
shit.S4[0] = shit.ly[0] * lr2;
shit.C[0] = -lr * 0.5;
float w;
- if (lr > 1.E-4)
- w = 1. / lr;
+ if (lr > 1.E-4) w = 1. / lr;
else
w = 1.;
shit.Cx[0] = w * shit.lx[0];
@@ -801,7 +777,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
if (shit.on_ring[0]) {
n++;
w = 1;
- } else {
+ }
+ else {
float dp = fabs(shit.C[0] + shit.Cx[0] * X + shit.Cy[0] * Y);
if (dp > Dcut) continue;
if (dp > Dmax) Dmax = dp;
@@ -874,8 +851,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
i->skip = 1;
continue;
}
- i->chi2 =
- i->chi2 / ((i->localIHits.size() - 3) * HitSize * HitSize); //.3/.3;
+ i->chi2 = i->chi2 / ((i->localIHits.size() - 3) * HitSize * HitSize); //.3/.3;
}
}
@@ -896,8 +872,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
for (iR j = i + 1; j != Rend; ++j) {
if (i->skip) continue;
float dist = j->r + best->r + HitSize2[0];
- if (fabs(j->x - best->x) > dist || fabs(j->y - best->y) > dist)
- continue;
+ if (fabs(j->x - best->x) > dist || fabs(j->y - best->y) > dist) continue;
j->NOwn = 0;
const THitIndex maxJ = j->localIHits.size();
for (THitIndex m = 0; m < maxJ; m++) {
@@ -908,7 +883,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
}
i->on = 1;
i->skip = 1;
- } else
+ }
+ else
i->skip = 1;
}
#else // NEW_SELECTION
@@ -931,17 +907,13 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
// ( ir->NHits < 10 && ir->NOwn < 1.00*ir->NHits ) ||
// ( ir->NOwn < .60*ir->NHits ) );
if (ir->skip) continue; // faster with if
- const bool skip = ((ir->NOwn < 1.0 * MinRingHits)
- || (ir->NHits < 10 && ir->NOwn < 1.00 * ir->NHits)
+ const bool skip = ((ir->NOwn < 1.0 * MinRingHits) || (ir->NHits < 10 && ir->NOwn < 1.00 * ir->NHits)
|| (ir->NOwn < .60 * ir->NHits));
- ir->skip = skip;
- const bool isBetter =
- !skip
- & ((ir->NOwn > 1.2 * bestOwn)
- || (ir->NOwn >= 1. * bestOwn && ir->chi2 < bestChi2));
- bestOwn = (isBetter) ? ir->NOwn : bestOwn; //Hits;
- bestChi2 = (isBetter) ? ir->chi2 : bestChi2;
- best = (isBetter) ? ir : best;
+ ir->skip = skip;
+ const bool isBetter = !skip & ((ir->NOwn > 1.2 * bestOwn) || (ir->NOwn >= 1. * bestOwn && ir->chi2 < bestChi2));
+ bestOwn = (isBetter) ? ir->NOwn : bestOwn; //Hits;
+ bestChi2 = (isBetter) ? ir->chi2 : bestChi2;
+ best = (isBetter) ? ir : best;
}
if (best == Rend) break;
best->skip = 1;
@@ -954,8 +926,7 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
for (iR ir = Rbeg; ir != Rend; ++ir) {
if (ir->skip) continue;
float dist = ir->r + best->r + HitSize2[0];
- if (fabs(ir->x - best->x) > dist || fabs(ir->y - best->y) > dist)
- continue;
+ if (fabs(ir->x - best->x) > dist || fabs(ir->y - best->y) > dist) continue;
ir->NOwn = 0;
const THitIndex NHitsCur = ir->localIHits.size();
for (THitIndex iih = 0; iih < NHitsCur; iih++) {
@@ -972,7 +943,8 @@ void CbmL1RichENNRingFinderParallel::ENNRingFinder(
#endif // PRINT_TIMING
}
-TStopwatch& CbmL1RichENNRingFinderParallel::GetTimer(TString t) {
+TStopwatch& CbmL1RichENNRingFinderParallel::GetTimer(TString t)
+{
int i = 0;
for (; (i < NTimers) && (fTimersNames[i] != t); i++)
;
diff --git a/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.h b/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.h
index 354518961b6dafd5b192d4c881d02bc4ac2fd771..7935155a1a95e8acbb695b32351d78cdab4d7c1c 100644
--- a/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.h
+++ b/reco/L1/OffLineInterface/CbmL1RichENNRingFinderParallel.h
@@ -43,9 +43,7 @@ class CbmL1RichENNRingFinderParallel : public CbmRichRingFinder {
int quality; // quality of the best ring with this hit
THitIndex localIndex; // index in local copy of Clone array
- static bool Compare(const ENNHit& h1, const ENNHit& h2) {
- return (h1.x < h2.x);
- };
+ static bool Compare(const ENNHit& h1, const ENNHit& h2) { return (h1.x < h2.x); };
};
@@ -62,13 +60,15 @@ class CbmL1RichENNRingFinderParallel : public CbmRichRingFinder {
fvec x, y; // coordinates
fvec quality; // quality of the best ring with this hit
fvec localIndex; // index in local copy of Clone array
- void CopyHit(ENNHit& a, int i) {
+ void CopyHit(ENNHit& a, int i)
+ {
localIndex[i] = a.localIndex;
x[i] = a.x;
y[i] = a.y;
quality[i] = a.quality;
}
- void CopyHit(ENNHitV& a, int j, int i) {
+ void CopyHit(ENNHitV& a, int j, int i)
+ {
localIndex[i] = a.localIndex[j];
x[i] = a.x[j];
y[i] = a.y[j];
@@ -81,11 +81,13 @@ class CbmL1RichENNRingFinderParallel : public CbmRichRingFinder {
fvec outIndex; // index in local copy of Clone array
- void CopyHit(ENNRingHit& a, int i) {
+ void CopyHit(ENNRingHit& a, int i)
+ {
outIndex[i] = a.outIndex;
ENNHitV::CopyHit(a, i);
}
- void CopyHit(ENNRingHitV& a, int j, int i) {
+ void CopyHit(ENNRingHitV& a, int j, int i)
+ {
outIndex[i] = a.outIndex[j];
ENNHitV::CopyHit(a, j, i);
};
@@ -144,9 +146,9 @@ class CbmL1RichENNRingFinderParallel : public CbmRichRingFinder {
bool skip; // skip the ring during selection
std::vector<THitIndex> localIHits; // indexes of hits in local array
- static bool CompareENNHRings(const ENNRing& r1, const ENNRing& r2) {
- if (r1.NHits != r2.NHits)
- return (r1.NHits > r2.NHits);
+ static bool CompareENNHRings(const ENNRing& r1, const ENNRing& r2)
+ {
+ if (r1.NHits != r2.NHits) return (r1.NHits > r2.NHits);
else
return (r1.chi2 < r2.chi2);
};
@@ -172,13 +174,8 @@ class CbmL1RichENNRingFinderParallel : public CbmRichRingFinder {
};
- void ENNRingFinder(const int NHits,
- nsL1vector<ENNHitV>::TSimd& HitsV,
- std::vector<ENNRing>& Rings,
- float HitSize = 1.,
- THitIndex MinRingHits = 5,
- fvec RMin = 2.,
- fvec RMax = 6.);
+ void ENNRingFinder(const int NHits, nsL1vector<ENNHitV>::TSimd& HitsV, std::vector<ENNRing>& Rings,
+ float HitSize = 1., THitIndex MinRingHits = 5, fvec RMin = 2., fvec RMax = 6.);
public:
/** Standard constructor **/
@@ -198,15 +195,16 @@ public:
**
*@value Number of tracks created
**/
- Int_t DoFind(TClonesArray* hitArray,
- TClonesArray* projArray,
- TClonesArray* ringArray);
+ Int_t DoFind(CbmEvent* event, TClonesArray* hitArray, TClonesArray* projArray, TClonesArray* ringArray);
private:
Float_t fRecoTime;
Int_t fNEvents;
- enum { NTimers = 11 };
+ enum
+ {
+ NTimers = 11
+ };
TStopwatch fTimers[NTimers]; // timers for different parts of algorithm
TString fTimersNames[NTimers]; // names which are correspond to the timers.
TStopwatch& GetTimer(TString t);
diff --git a/reco/base/CbmRichRingFinder.h b/reco/base/CbmRichRingFinder.h
index 34669b49c06989982258cf976450d46d0e81f335..4e6894fe86cf8338b735517385b17821a958c6b9 100644
--- a/reco/base/CbmRichRingFinder.h
+++ b/reco/base/CbmRichRingFinder.h
@@ -27,13 +27,14 @@
#include "TObject.h"
class TClonesArray;
+class CbmEvent;
class CbmRichRingFinder : public TObject {
public:
/** Default constructor **/
- CbmRichRingFinder() : TObject(), fVerbose(0) {};
+ CbmRichRingFinder() {};
/** Destructor **/
@@ -50,25 +51,15 @@ public:
** Task: Read the hit array and fill the ring array,
** pointers to which are given as arguments
**
+ *@param event CbmEvent for time-based mode, if event==nullptr then event-by-event mode.
*@param rHitArray Array of RICH hits
*@param rProjArray Array of projected tracks (for track based finders)
*@param rRingArray Array of CbmRichRing
*@value Number of rings created
**/
- virtual Int_t DoFind(TClonesArray* rHitArray,
- TClonesArray* rProjArray,
+ virtual Int_t DoFind(CbmEvent* event, TClonesArray* rHitArray, TClonesArray* rProjArray,
TClonesArray* rRingArray) = 0;
-
- /** Set verbosity
- *@param verbose Verbosity level
- **/
- void SetVerbose(Int_t verbose) { fVerbose = verbose; };
-
-
-protected:
- Int_t fVerbose; // Verbosity level
-
private:
CbmRichRingFinder(const CbmRichRingFinder&);
CbmRichRingFinder& operator=(const CbmRichRingFinder&);
diff --git a/reco/detectors/rich/CbmRichReconstruction.cxx b/reco/detectors/rich/CbmRichReconstruction.cxx
index cc3b65e0538c0e404e78cf226b513d5ac5886808..f119bc294047a0cbc9cbed3513e95d14c9e2352a 100644
--- a/reco/detectors/rich/CbmRichReconstruction.cxx
+++ b/reco/detectors/rich/CbmRichReconstruction.cxx
@@ -6,24 +6,24 @@
**/
#include "CbmRichReconstruction.h"
-#include "CbmRichRing.h"
#include "CbmRichProjectionProducerAnalytical.h"
#include "CbmRichProjectionProducerTGeo.h"
+#include "CbmRichRing.h"
//#include "prototype/CbmRichProtProjectionProducer.h"
+#include "CbmL1RichENNRingFinder.h"
+#include "CbmRichRingFinderHough.h"
+#include "CbmRichRingFinderIdeal.h"
#include "CbmRichTrackExtrapolationBase.h"
#include "CbmRichTrackExtrapolationIdeal.h"
#include "CbmRichTrackExtrapolationKF.h"
#include "CbmRichTrackExtrapolationLittrack.h"
#include "CbmRichTrackExtrapolationMirrorIdeal.h"
-
-#include "CbmL1RichENNRingFinder.h"
-#include "CbmRichRingFinderHough.h"
-#include "CbmRichRingFinderIdeal.h"
//#include "CbmL1RichENNRingFinderParallel.h"
//#include "prototype/CbmRichProtRingFinderHough.h"
+#include "CbmEvent.h"
#include "CbmGlobalTrack.h"
#include "CbmRichConverter.h"
#include "CbmRichRingFitterCOP.h"
@@ -45,81 +45,47 @@
using std::cout;
using std::endl;
-CbmRichReconstruction::CbmRichReconstruction()
- : FairTask("CbmRichReconstruction")
- , fRichHits(NULL)
- , fRichRings(NULL)
- , fRichProjections(NULL)
- , fRichTrackParamZ(NULL)
- , fGlobalTracks(NULL)
- ,
-
- fRingFinder(NULL)
- , fRingFitter(NULL)
- , fTrackExtrapolation(NULL)
- , fProjectionProducer(NULL)
- , fRingTrackAssign(NULL)
- ,
-
- fRunExtrapolation(true)
- , fRunProjection(true)
- , fRunFinder(true)
- , fRunFitter(true)
- , fRunTrackAssign(true)
- ,
-
- fUseHTAnnSelect(true)
- ,
-
- fExtrapolationName("littrack")
- , fProjectionName("analytical")
- , fFinderName("hough")
- , fFitterName("ellipse_tau")
- , fTrackAssignName("closest_distance")
- ,
-
- fZTrackExtrapolation(260.) {}
-
-CbmRichReconstruction::~CbmRichReconstruction() {
- if (NULL != fRingFinder) delete fRingFinder;
- if (NULL != fRingFitter) delete fRingFitter;
- if (NULL != fTrackExtrapolation) delete fTrackExtrapolation;
- if (NULL != fProjectionProducer) delete fProjectionProducer;
+CbmRichReconstruction::CbmRichReconstruction() : FairTask("CbmRichReconstruction") {}
+
+CbmRichReconstruction::~CbmRichReconstruction()
+{
+ if (nullptr != fRingFinder) delete fRingFinder;
+ if (nullptr != fRingFitter) delete fRingFitter;
+ if (nullptr != fTrackExtrapolation) delete fTrackExtrapolation;
+ if (nullptr != fProjectionProducer) delete fProjectionProducer;
+ if (nullptr != fRingTrackAssign) delete fRingTrackAssign;
}
-InitStatus CbmRichReconstruction::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichReconstruction::Init", "RootManager not instantised!");
+InitStatus CbmRichReconstruction::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (nullptr == manager) LOG(fatal) << "CbmRichReconstruction::Init(): FairRootManager is nullptr.";
+
+ fCbmEvents = dynamic_cast<TClonesArray*>(manager->GetObject("CbmEvent"));
+ if (fCbmEvents == nullptr) { LOG(info) << ": CbmEvent NOT found \n \n \n"; }
+ else {
+ LOG(info) << ": CbmEvent found \n \n \n";
}
- if (fRunExtrapolation && fRunProjection) {
+ if (fRunExtrapolation) {
fRichTrackParamZ = new TClonesArray("FairTrackParam", 100);
- ioman->Register("RichTrackParamZ",
- "RICH",
- fRichTrackParamZ,
- IsOutputBranchPersistent("RichTrackParamZ"));
-
- fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
- if (NULL == fGlobalTracks) {
- Fatal("CbmRichReconstruction::Init", "No GlobalTrack array!");
- }
+ manager->Register("RichTrackParamZ", "RICH", fRichTrackParamZ, IsOutputBranchPersistent("RichTrackParamZ"));
- fRichProjections = new TClonesArray("FairTrackParam");
- ioman->Register("RichProjection",
- "RICH",
- fRichProjections,
- IsOutputBranchPersistent("RichProjection"));
+ fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
+ if (fGlobalTracks == nullptr) LOG(fatal) << "CbmRichReconstruction::Init(): No GlobalTrack array.";
}
- fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
- if (NULL == fRichHits) {
- Fatal("CbmRichReconstruction::Init", "No RichHit array!");
+ if (fRunProjection) {
+ if (!fRunExtrapolation) LOG(fatal) << "CbmRichReconstruction::Init(): fRunExtrapolation must be true.";
+ fRichProjections = new TClonesArray("FairTrackParam");
+ manager->Register("RichProjection", "RICH", fRichProjections, IsOutputBranchPersistent("RichProjection"));
}
+ fRichHits = (TClonesArray*) manager->GetObject("RichHit");
+ if (fRichHits == nullptr) LOG(fatal) << "CbmRichReconstruction::Init(): No RichHit array.";
+
fRichRings = new TClonesArray("CbmRichRing", 100);
- ioman->Register(
- "RichRing", "RICH", fRichRings, IsOutputBranchPersistent("RichRing"));
+ manager->Register("RichRing", "RICH", fRichRings, IsOutputBranchPersistent("RichRing"));
if (fRunExtrapolation) InitExtrapolation();
if (fRunProjection) InitProjection();
@@ -130,115 +96,146 @@ InitStatus CbmRichReconstruction::Init() {
return kSUCCESS;
}
-void CbmRichReconstruction::Exec(Option_t* /*opt*/) {
- LOG(info) << "CbmRichReconstruction Exec";
- if (fRunExtrapolation) RunExtrapolation();
- if (fRunProjection) RunProjection();
- if (fRunFinder) RunFinder();
- if (fRunFitter) RunFitter();
- if (fRunTrackAssign) RunTrackAssign();
+void CbmRichReconstruction::Exec(Option_t* /*opt*/)
+{
+ if (fRichTrackParamZ != nullptr) fRichTrackParamZ->Delete();
+ if (fRichProjections != nullptr) fRichProjections->Delete();
+ if (fRichRings != nullptr) fRichRings->Delete();
+
+ if (fCbmEvents == nullptr) {
+ LOG(info) << "CbmRichReconstruction: Event " << ++fEventNum;
+ ProcessData(nullptr);
+ }
+ else {
+ int nEvents = fCbmEvents->GetEntriesFast();
+ LOG(info) << "CbmRichReconstruction: Timeslice " << ++fEventNum << " with " << nEvents << " events";
+ for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
+ CbmEvent* event = static_cast<CbmEvent*>(fCbmEvents->At(iEvent));
+ ProcessData(event);
+ }
+ }
+}
+
+void CbmRichReconstruction::ProcessData(CbmEvent* event)
+{
+ if (fRunExtrapolation) RunExtrapolation(event);
+ if (fRunProjection) RunProjection(event);
+ if (fRunFinder) RunFinder(event);
+ if (fRunFitter) RunFitter(event);
+ if (fRunTrackAssign) RunTrackAssign(event);
}
-void CbmRichReconstruction::InitExtrapolation() {
- if (fExtrapolationName == "ideal") {
- fTrackExtrapolation = new CbmRichTrackExtrapolationIdeal();
- } else if (fExtrapolationName == "mirror_ideal") {
+void CbmRichReconstruction::InitExtrapolation()
+{
+ if (fExtrapolationName == "ideal") { fTrackExtrapolation = new CbmRichTrackExtrapolationIdeal(); }
+ else if (fExtrapolationName == "mirror_ideal") {
fTrackExtrapolation = new CbmRichTrackExtrapolationMirrorIdeal();
- } else if (fExtrapolationName == "kf" || fExtrapolationName == "KF") {
+ }
+ else if (fExtrapolationName == "kf" || fExtrapolationName == "KF") {
fTrackExtrapolation = new CbmRichTrackExtrapolationKF();
- } else if (fExtrapolationName == "lit" || fExtrapolationName == "littrack") {
+ }
+ else if (fExtrapolationName == "lit" || fExtrapolationName == "littrack") {
fTrackExtrapolation = new CbmRichTrackExtrapolationLittrack();
- } else {
- LOG(fatal) << fExtrapolationName
- << " is not correct name for extrapolation algorithm.";
+ }
+ else {
+ LOG(fatal) << fExtrapolationName << " is not correct name for extrapolation algorithm.";
}
fTrackExtrapolation->Init();
}
-void CbmRichReconstruction::InitProjection() {
- if (fProjectionName == "analytical") {
- fProjectionProducer = new CbmRichProjectionProducerAnalytical();
- } else if (fProjectionName == "TGeo" || fProjectionName == "tgeo") {
+void CbmRichReconstruction::InitProjection()
+{
+ if (fProjectionName == "analytical") { fProjectionProducer = new CbmRichProjectionProducerAnalytical(); }
+ else if (fProjectionName == "TGeo" || fProjectionName == "tgeo") {
fProjectionProducer = new CbmRichProjectionProducerTGeo();
- } else {
- LOG(fatal) << fFinderName
- << " is not correct name for projection producer algorithm.";
+ }
+ else {
+ LOG(fatal) << fFinderName << " is not correct name for projection producer algorithm.";
}
fProjectionProducer->Init();
}
-void CbmRichReconstruction::InitFinder() {
+void CbmRichReconstruction::InitFinder()
+{
if (fFinderName == "hough") {
fRingFinder = new CbmRichRingFinderHough();
- static_cast<CbmRichRingFinderHough*>(fRingFinder)
- ->SetUseAnnSelect(fUseHTAnnSelect);
- } else if (fFinderName == "ideal") {
+ static_cast<CbmRichRingFinderHough*>(fRingFinder)->SetUseAnnSelect(fUseHTAnnSelect);
+ }
+ else if (fFinderName == "ideal") {
fRingFinder = new CbmRichRingFinderIdeal();
- } else if (fFinderName == "enn") {
+ }
+ else if (fFinderName == "enn") {
fRingFinder = new CbmL1RichENNRingFinder(0);
- } else if ((fFinderName == "enn_parallel")) {
+ }
+ else if ((fFinderName == "enn_parallel")) {
// fRingFinder = new CbmL1RichENNRingFinderParallel(0);
// } else if (fFinderName == "hough_prototype") {
// fRingFinder = new CbmRichProtRingFinderHough();
- } else {
- LOG(fatal) << fFinderName
- << " is not correct name for ring finder algorithm.";
+ }
+ else {
+ LOG(fatal) << fFinderName << " is not correct name for ring finder algorithm.";
}
fRingFinder->Init();
}
-void CbmRichReconstruction::InitFitter() {
- if (fFitterName == "circle_cop") {
- fRingFitter = new CbmRichRingFitterCOP();
- } else if (fFitterName == "circle_simple") {
+void CbmRichReconstruction::InitFitter()
+{
+ if (fFitterName == "circle_cop") { fRingFitter = new CbmRichRingFitterCOP(); }
+ else if (fFitterName == "circle_simple") {
fRingFitter = new CbmRichRingFitterCircle();
- } else if (fFitterName == "circle_tau") {
+ }
+ else if (fFitterName == "circle_tau") {
fRingFitter = new CbmRichRingFitterTAU();
- } else if (fFitterName == "circle_robust_cop") {
+ }
+ else if (fFitterName == "circle_robust_cop") {
fRingFitter = new CbmRichRingFitterRobustCOP();
- } else if (fFitterName == "ellipse_tau") {
+ }
+ else if (fFitterName == "ellipse_tau") {
fRingFitter = new CbmRichRingFitterEllipseTau();
- } else if (fFitterName == "ellipse_minuit") {
+ }
+ else if (fFitterName == "ellipse_minuit") {
fRingFitter = new CbmRichRingFitterEllipseMinuit();
- } else {
- LOG(fatal) << fFitterName
- << " is not correct name for ring fitter algorithm.";
+ }
+ else {
+ LOG(fatal) << fFitterName << " is not correct name for ring fitter algorithm.";
}
CbmRichConverter::Init();
}
-void CbmRichReconstruction::InitTrackAssign() {
- if (fTrackAssignName == "closest_distance") {
- fRingTrackAssign = new CbmRichRingTrackAssignClosestD();
- } else {
- LOG(fatal) << fTrackAssignName
- << " is not correct name for ring-track assignment algorithm.";
+void CbmRichReconstruction::InitTrackAssign()
+{
+ if (fTrackAssignName == "closest_distance") { fRingTrackAssign = new CbmRichRingTrackAssignClosestD(); }
+ else {
+ LOG(fatal) << fTrackAssignName << " is not correct name for ring-track assignment algorithm.";
}
fRingTrackAssign->Init();
}
-void CbmRichReconstruction::RunExtrapolation() {
- if (fRichTrackParamZ == NULL) LOG(info) << "fRichTrackParamZ == NULL";
- fRichTrackParamZ->Delete();
- fTrackExtrapolation->DoExtrapolation(
- fGlobalTracks, fRichTrackParamZ, fZTrackExtrapolation);
+void CbmRichReconstruction::RunExtrapolation(CbmEvent* event)
+{
+ if (fRichTrackParamZ == nullptr) LOG(info) << "fRichTrackParamZ == nullptr";
+ fTrackExtrapolation->DoExtrapolation(event, fGlobalTracks, fRichTrackParamZ, fZTrackExtrapolation);
}
-void CbmRichReconstruction::RunProjection() {
- fProjectionProducer->DoProjection(fRichProjections);
+void CbmRichReconstruction::RunProjection(CbmEvent* event)
+{
+ fProjectionProducer->DoProjection(event, fRichProjections);
}
-void CbmRichReconstruction::RunFinder() {
- fRichRings->Delete();
- fRingFinder->DoFind(fRichHits, fRichProjections, fRichRings);
+void CbmRichReconstruction::RunFinder(CbmEvent* event)
+{
+ fRingFinder->DoFind(event, fRichHits, fRichProjections, fRichRings);
}
-void CbmRichReconstruction::RunFitter() {
- int nofRings = fRichRings->GetEntriesFast();
- for (int iRing = 0; iRing < nofRings; iRing++) {
- CbmRichRing* ring = (CbmRichRing*) fRichRings->At(iRing);
- if (NULL == ring) continue;
+void CbmRichReconstruction::RunFitter(CbmEvent* event)
+{
+ const Int_t nofRings = event ? event->GetNofData(ECbmDataType::kRichRing) : fRichRings->GetEntriesFast();
+ if (nofRings <= 0) return;
+ for (Int_t iR0 = 0; iR0 < nofRings; iR0++) {
+ Int_t iR = event ? event->GetIndex(ECbmDataType::kRichRing, iR0) : iR0;
+ CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
+ if (nullptr == ring) continue;
CbmRichRingLight ringL;
CbmRichConverter::CopyHitsToRingLight(ring, &ringL);
@@ -247,8 +244,9 @@ void CbmRichReconstruction::RunFitter() {
}
}
-void CbmRichReconstruction::RunTrackAssign() {
- fRingTrackAssign->DoAssign(fRichRings, fRichProjections);
+void CbmRichReconstruction::RunTrackAssign(CbmEvent* event)
+{
+ fRingTrackAssign->DoAssign(event, fRichRings, fRichProjections);
}
void CbmRichReconstruction::Finish() {}
diff --git a/reco/detectors/rich/CbmRichReconstruction.h b/reco/detectors/rich/CbmRichReconstruction.h
index d71e3592a2dc3415304c9637a743be8df2bbcf8b..2e5d4b17631561b074392ea461623c0c72611852 100644
--- a/reco/detectors/rich/CbmRichReconstruction.h
+++ b/reco/detectors/rich/CbmRichReconstruction.h
@@ -20,6 +20,7 @@ class CbmRichRingFitterBase;
class CbmRichTrackExtrapolationBase;
class CbmRichProjectionProducerBase;
class CbmRichRingTrackAssignBase;
+class CbmEvent;
using std::string;
@@ -79,7 +80,8 @@ public:
*/
void SetZTrackExtrapolation(Double_t z) { fZTrackExtrapolation = z; }
- void UseMCbmSetup() {
+ void UseMCbmSetup()
+ {
this->SetRunExtrapolation(false);
this->SetRunProjection(false);
this->SetRunTrackAssign(false);
@@ -87,39 +89,41 @@ public:
}
private:
- TClonesArray* fRichHits;
- TClonesArray* fRichRings;
- TClonesArray* fRichProjections;
- TClonesArray* fRichTrackParamZ;
- TClonesArray* fGlobalTracks;
-
- CbmRichRingFinder* fRingFinder; // pointer to ring finder algorithm
- CbmRichRingFitterBase* fRingFitter; // pointer to ring fitting algorithm
- CbmRichTrackExtrapolationBase*
- fTrackExtrapolation; // pointer to track extrapolation algorithm
- CbmRichProjectionProducerBase*
- fProjectionProducer; // pointer to projection producer
- CbmRichRingTrackAssignBase*
- fRingTrackAssign; // pointer to track assignment algorithm
-
- // What do you wan to run.
- bool fRunExtrapolation;
- bool fRunProjection;
- bool fRunFinder;
- bool fRunFitter;
- bool fRunTrackAssign;
-
- bool fUseHTAnnSelect;
+ TClonesArray* fRichHits = nullptr;
+ TClonesArray* fRichRings = nullptr;
+ TClonesArray* fRichProjections = nullptr;
+ TClonesArray* fRichTrackParamZ = nullptr;
+ TClonesArray* fGlobalTracks = nullptr;
+ TClonesArray* fCbmEvents = nullptr;
+
+ Int_t fEventNum = 0; // event or timeslice counter
+
+ // pointers to the algorithms
+ CbmRichRingFinder* fRingFinder = nullptr;
+ CbmRichRingFitterBase* fRingFitter = nullptr;
+ CbmRichTrackExtrapolationBase* fTrackExtrapolation = nullptr;
+ CbmRichProjectionProducerBase* fProjectionProducer = nullptr;
+ CbmRichRingTrackAssignBase* fRingTrackAssign = nullptr;
+
+ // What do you want to run.
+ bool fRunExtrapolation = true;
+ bool fRunProjection = true;
+ bool fRunFinder = true;
+ bool fRunFitter = true;
+ bool fRunTrackAssign = true;
+
+ // Run ring-candidate selection algorithm based on ANN
+ bool fUseHTAnnSelect = true;
// Algorithm names for each step of reconstruction.
- string fExtrapolationName; // name of extrapolation algorithm
- string fProjectionName; // name of track projection algorithm
- string fFinderName; // name of ring finder algorithm
- string fFitterName; // name of ring fitter algorithm
- string fTrackAssignName; // name of track-ring matching algorithm
+ string fExtrapolationName = "littrack";
+ string fProjectionName = "analytical";
+ string fFinderName = "hough";
+ string fFitterName = "ellipse_tau";
+ string fTrackAssignName = "closest_distance";
- Double_t
- fZTrackExtrapolation; // Z coordinate to which one wants to extrapolate STS tracks
+ // Z coordinate where STS tracks will be extrapolated.
+ Double_t fZTrackExtrapolation = 260.;
/**
* \brief
@@ -149,27 +153,32 @@ private:
/**
* \brief
*/
- void RunExtrapolation();
+ void ProcessData(CbmEvent* event);
/**
* \brief
*/
- void RunProjection();
+ void RunExtrapolation(CbmEvent* event);
/**
* \brief
*/
- void RunFinder();
+ void RunProjection(CbmEvent* event);
/**
* \brief
*/
- void RunFitter();
+ void RunFinder(CbmEvent* event);
/**
* \brief
*/
- void RunTrackAssign();
+ void RunFitter(CbmEvent* event);
+
+ /**
+ * \brief
+ */
+ void RunTrackAssign(CbmEvent* event);
/**
* \brief Copy constructor.
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderData.h b/reco/detectors/rich/finder/CbmRichRingFinderData.h
index c3665bb0a9a5c83f754c0bdd7d1eee25bdf42700..dfeb5c6fbf22bffe59074a1ca007389b8e5ebe04 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderData.h
+++ b/reco/detectors/rich/finder/CbmRichRingFinderData.h
@@ -41,15 +41,11 @@ public:
* \author Semen Lebedev
* \date 2008
**/
-class CbmRichHoughHitCmpUp :
- public std::
- binary_function<const CbmRichHoughHit, const CbmRichHoughHit, bool> {
+class CbmRichHoughHitCmpUp : public std::binary_function<const CbmRichHoughHit, const CbmRichHoughHit, bool> {
public:
virtual ~CbmRichHoughHitCmpUp() {}
- bool operator()(const CbmRichHoughHit& m1, const CbmRichHoughHit& m2) const {
- return m1.fHit.fX < m2.fHit.fX;
- }
+ bool operator()(const CbmRichHoughHit& m1, const CbmRichHoughHit& m2) const { return m1.fHit.fX < m2.fHit.fX; }
};
@@ -61,14 +57,12 @@ public:
* \author Semen Lebedev
* \date 2008
**/
-class CbmRichRingComparatorMore :
- public std::
- binary_function<const CbmRichRingLight*, const CbmRichRingLight*, bool> {
+class CbmRichRingComparatorMore : public std::binary_function<const CbmRichRingLight*, const CbmRichRingLight*, bool> {
public:
virtual ~CbmRichRingComparatorMore() {}
- bool operator()(const CbmRichRingLight* ring1,
- const CbmRichRingLight* ring2) const {
+ bool operator()(const CbmRichRingLight* ring1, const CbmRichRingLight* ring2) const
+ {
return ring1->GetSelectionNN() > ring2->GetSelectionNN();
}
};
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx b/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx
index d1b103e8760b3be0b1dacfa499338e0b0d68fd3c..cdee4ad6bfd5dfe60c3724bb8b1beef7ace4492c 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx
@@ -6,9 +6,11 @@
**/
#include "CbmRichRingFinderHough.h"
+
#include "CbmRichRingFinderHoughImpl.h"
//#include "CbmRichRingFinderHoughSimd.h"
//#include "../../littrack/utils/CbmLitMemoryManagment.h"
+#include "CbmEvent.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
@@ -24,7 +26,7 @@ using std::endl;
using std::vector;
CbmRichRingFinderHough::CbmRichRingFinderHough()
- : fNEvent(0), fRingCount(0), fHTImpl(NULL), fUseAnnSelect(true) {
+{
#ifdef HOUGH_SERIAL
fHTImpl = new CbmRichRingFinderHoughImpl();
#endif
@@ -34,54 +36,52 @@ CbmRichRingFinderHough::CbmRichRingFinderHough()
#endif
}
-void CbmRichRingFinderHough::Init() {
+void CbmRichRingFinderHough::Init()
+{
fHTImpl->SetUseAnnSelect(fUseAnnSelect);
fHTImpl->Init();
}
CbmRichRingFinderHough::~CbmRichRingFinderHough() { delete fHTImpl; }
-Int_t CbmRichRingFinderHough::DoFind(TClonesArray* rHitArray,
- TClonesArray* /*rProjArray*/,
- TClonesArray* rRingArray) {
+Int_t CbmRichRingFinderHough::DoFind(CbmEvent* event, TClonesArray* rHitArray, TClonesArray* /*rProjArray*/,
+ TClonesArray* rRingArray)
+{
TStopwatch timer;
timer.Start();
- fNEvent++;
- LOG(info) << "-I- CbmRichRingFinderHough Event/Timeslice no. " << fNEvent;
+ fEventNum++;
vector<CbmRichHoughHit> UpH;
vector<CbmRichHoughHit> DownH;
- fRingCount = 0;
- if (NULL == rHitArray) {
- LOG(error) << "-E- CbmRichRingFinderHough::DoFind: Hit array missing!"
- << rHitArray;
+ if (rHitArray == nullptr) {
+ LOG(error) << "CbmRichRingFinderHough::DoFind(): Hit array is nullptr.";
return -1;
}
- const Int_t nhits = rHitArray->GetEntriesFast();
- if (nhits <= 0) {
- LOG(error) << "-E- CbmRichRingFinderHough::DoFind: No hits in this event!";
+
+ const Int_t nofRichHits = event ? event->GetNofData(ECbmDataType::kRichHit) : rHitArray->GetEntriesFast();
+ if (nofRichHits <= 0) {
+ LOG(debug) << "CbmRichRingFinderHough::DoFind(): No RICH hits in this event.";
return -1;
}
- UpH.reserve(nhits / 2);
- DownH.reserve(nhits / 2);
+ UpH.reserve(nofRichHits / 2);
+ DownH.reserve(nofRichHits / 2);
// convert CbmRichHit to CbmRichHoughHit and
// sort hits according to the photodetector (up or down)
- for (Int_t iHit = 0; iHit < nhits; iHit++) {
- CbmRichHit* hit = static_cast<CbmRichHit*>(rHitArray->At(iHit));
- if (hit) {
+ for (Int_t iH0 = 0; iH0 < nofRichHits; iH0++) {
+ Int_t iH = event ? event->GetIndex(ECbmDataType::kRichHit, iH0) : iH0;
+ CbmRichHit* hit = static_cast<CbmRichHit*>(rHitArray->At(iH));
+ if (hit != nullptr) {
CbmRichHoughHit tempPoint;
- tempPoint.fHit.fX = hit->GetX();
- tempPoint.fHit.fY = hit->GetY();
- tempPoint.fHit.fId = iHit;
- tempPoint.fX2plusY2 =
- hit->GetX() * hit->GetX() + hit->GetY() * hit->GetY();
- tempPoint.fTime = hit->GetTime();
- tempPoint.fIsUsed = false;
- if (hit->GetY() >= 0)
- UpH.push_back(tempPoint);
+ tempPoint.fHit.fX = hit->GetX();
+ tempPoint.fHit.fY = hit->GetY();
+ tempPoint.fHit.fId = iH;
+ tempPoint.fX2plusY2 = hit->GetX() * hit->GetX() + hit->GetY() * hit->GetY();
+ tempPoint.fTime = hit->GetTime();
+ tempPoint.fIsUsed = false;
+ if (hit->GetY() >= 0) UpH.push_back(tempPoint);
else
DownH.push_back(tempPoint);
}
@@ -91,12 +91,10 @@ Int_t CbmRichRingFinderHough::DoFind(TClonesArray* rHitArray,
Double_t dt1 = timer.RealTime();
timer.Start();
- LOG(debug) << "-I- CbmRichRingFinderHough nofHits Up:" << UpH.size();
fHTImpl->SetData(UpH);
fHTImpl->DoFind();
- if (rRingArray != NULL)
- AddRingsToOutputArray(rRingArray, rHitArray, fHTImpl->GetFoundRings());
+ if (rRingArray != nullptr) AddRingsToOutputArray(event, rRingArray, rHitArray, fHTImpl->GetFoundRings());
//for_each(UpH.begin(), UpH.end(), DeleteObject());
UpH.clear();
@@ -104,29 +102,26 @@ Int_t CbmRichRingFinderHough::DoFind(TClonesArray* rHitArray,
Double_t dt2 = timer.RealTime();
timer.Start();
- LOG(debug) << "-I- CbmRichRingFinderHough nofHits Down:" << DownH.size();
fHTImpl->SetData(DownH);
fHTImpl->DoFind();
- if (rRingArray != NULL)
- AddRingsToOutputArray(rRingArray, rHitArray, fHTImpl->GetFoundRings());
+ if (rRingArray != nullptr) AddRingsToOutputArray(event, rRingArray, rHitArray, fHTImpl->GetFoundRings());
//for_each(DownH.begin(), DownH.end(), DeleteObject());
DownH.clear();
timer.Stop();
Double_t dt3 = timer.RealTime();
- LOG(info) << "CbmRichRingFinderHough. Number of found rings "
- << rRingArray->GetEntriesFast();
-
- LOG(info) << "time1:" << dt1 << " time2:" << dt2 << " time3:" << dt3
- << " total:" << dt1 + dt2 + dt3;
+ int nofFoundRings = event ? event->GetNofData(ECbmDataType::kRichRing) : rRingArray->GetEntriesFast();
+ LOG(info) << "CbmRichRingFinderHough::DoFind(): Event:" << fEventNum << " hits:" << nofRichHits
+ << " rings:" << nofFoundRings << " ringsInTS:" << rRingArray->GetEntriesFast()
+ << " Time:" << dt1 + dt2 + dt3;
return 1;
}
-void CbmRichRingFinderHough::AddRingsToOutputArray(
- TClonesArray* rRingArray,
- TClonesArray* rHitArray,
- const vector<CbmRichRingLight*>& rings) {
+void CbmRichRingFinderHough::AddRingsToOutputArray(CbmEvent* event, TClonesArray* rRingArray, TClonesArray* rHitArray,
+ const vector<CbmRichRingLight*>& rings)
+{
+
for (UInt_t iRing = 0; iRing < rings.size(); iRing++) {
if (rings[iRing]->GetRecFlag() == -1) continue;
CbmRichRing* r = new CbmRichRing();
@@ -143,7 +138,8 @@ void CbmRichRingFinderHough::AddRingsToOutputArray(
}
}
r->SetTime(ringTime / (double) ringCounter);
- new ((*rRingArray)[fRingCount]) CbmRichRing(*r);
- fRingCount++;
+ int nofRings = rRingArray->GetEntriesFast();
+ new ((*rRingArray)[nofRings]) CbmRichRing(*r);
+ if (event != nullptr) event->AddData(ECbmDataType::kRichRing, nofRings);
}
}
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHough.h b/reco/detectors/rich/finder/CbmRichRingFinderHough.h
index f36343439cbedeb0a3a374a60e892493ad03d807..d85c319aca1c2efdade8199b43054a56e583192e 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHough.h
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHough.h
@@ -11,6 +11,7 @@
#define CBM_RICH_RING_FINDER_HOUGH
#include "CbmRichRingFinder.h"
+
#include <vector>
class CbmRichRingFinderHoughImpl;
@@ -32,18 +33,6 @@ using std::vector;
* \date 2008
**/
class CbmRichRingFinderHough : public CbmRichRingFinder {
-protected:
- Int_t fNEvent; // event number
- Int_t fRingCount; // number of found rings
-
-// choose between serial and SIMD implementation of the ring finder
-#ifdef HOUGH_SERIAL
- CbmRichRingFinderHoughImpl* fHTImpl;
-#endif
-
-#ifdef HOUGH_SIMD
- CbmRichRingFinderHoughSimd* fHTImpl;
-#endif
public:
/**
@@ -64,22 +53,30 @@ public:
/**
* \brief Inherited from CbmRichRingFinder.
*/
- virtual Int_t DoFind(TClonesArray* rHitArray,
- TClonesArray* rProjArray,
- TClonesArray* rRingArray);
+ virtual Int_t DoFind(CbmEvent* event, TClonesArray* rHitArray, TClonesArray* rProjArray, TClonesArray* rRingArray);
void SetUseAnnSelect(bool use) { fUseAnnSelect = use; }
private:
- bool fUseAnnSelect;
+ Int_t fEventNum = 0;
+ Bool_t fUseAnnSelect = true;
+
+// choose between serial and SIMD implementation of the ring finder
+#ifdef HOUGH_SERIAL
+ CbmRichRingFinderHoughImpl* fHTImpl = nullptr;
+#endif
+
+#ifdef HOUGH_SIMD
+ CbmRichRingFinderHoughSimd* fHTImpl = nullptr;
+#endif
+
/**
* \brief Add found rings to the output TClonesArray.
* \param[out] rRingArray Output array of CbmRichRing.
* \param[in] rHitArray Array of CbmRichHit.
* \param[in] rings Found rings.
*/
- void AddRingsToOutputArray(TClonesArray* rRingArray,
- TClonesArray* rHitArray,
+ void AddRingsToOutputArray(CbmEvent* event, TClonesArray* rRingArray, TClonesArray* rHitArray,
const vector<CbmRichRingLight*>& rings);
/**
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.cxx b/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.cxx
index eddad22e6be21d687be5472fe9229a67efdfef3f..e80c4e969112012f2c0cdc7cd11dc05e19f93d9f 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.cxx
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.cxx
@@ -6,19 +6,20 @@
**/
#include "CbmRichRingFinderHoughImpl.h"
-#include "CbmRichRingLight.h"
#include "../../littrack/std/utils/CbmLitMemoryManagment.h"
#include "CbmRichRingFitterCOP.h"
+#include "CbmRichRingLight.h"
#include "CbmRichRingSelectAnn.h"
#include "TSystem.h"
#include <algorithm>
-#include <cmath>
#include <iostream>
#include <set>
+#include <cmath>
+
using std::cout;
using std::endl;
using std::set;
@@ -88,14 +89,17 @@ CbmRichRingFinderHoughImpl::CbmRichRingFinderHoughImpl()
fCurTime(0.)
-{}
+{
+}
-CbmRichRingFinderHoughImpl::~CbmRichRingFinderHoughImpl() {
+CbmRichRingFinderHoughImpl::~CbmRichRingFinderHoughImpl()
+{
if (NULL != fFitCOP) delete fFitCOP;
if (NULL != fANNSelect) delete fANNSelect;
}
-void CbmRichRingFinderHoughImpl::Init() {
+void CbmRichRingFinderHoughImpl::Init()
+{
SetParameters();
fHist.resize(fNofBinsXY);
@@ -109,7 +113,8 @@ void CbmRichRingFinderHoughImpl::Init() {
}
}
-void CbmRichRingFinderHoughImpl::DoFind() {
+void CbmRichRingFinderHoughImpl::DoFind()
+{
// if (fData.size() > MAX_NOF_HITS) {
// cout << endl << endl << "-E- CbmRichRingFinderHoughImpl::DoFind" << ". Number of hits is more than " << MAX_NOF_HITS << endl << endl;
// return;
@@ -125,7 +130,8 @@ void CbmRichRingFinderHoughImpl::DoFind() {
fData.clear();
}
-void CbmRichRingFinderHoughImpl::SetParameters() {
+void CbmRichRingFinderHoughImpl::SetParameters()
+{
fMaxDistance = 11.5;
fMinDistance = 2.5;
fMinDistanceSq = fMinDistance * fMinDistance;
@@ -159,7 +165,8 @@ void CbmRichRingFinderHoughImpl::SetParameters() {
fNofBinsXY = fNofBinsX * fNofBinsY;
}
-void CbmRichRingFinderHoughImpl::HoughTransformReconstruction() {
+void CbmRichRingFinderHoughImpl::HoughTransformReconstruction()
+{
int indmin, indmax;
unsigned int size = fData.size();
for (unsigned int iHit = 0; iHit < size; iHit++) {
@@ -169,29 +176,23 @@ void CbmRichRingFinderHoughImpl::HoughTransformReconstruction() {
fCurMinY = fData[iHit].fHit.fY - fMaxDistance;
fCurTime = fData[iHit].fTime;
- DefineLocalAreaAndHits(
- fData[iHit].fHit.fX, fData[iHit].fHit.fY, &indmin, &indmax);
+ DefineLocalAreaAndHits(fData[iHit].fHit.fX, fData[iHit].fHit.fY, &indmin, &indmax);
HoughTransform(indmin, indmax);
FindPeak(indmin, indmax);
}
}
-void CbmRichRingFinderHoughImpl::DefineLocalAreaAndHits(float x0,
- float y0,
- int* indmin,
- int* indmax) {
+void CbmRichRingFinderHoughImpl::DefineLocalAreaAndHits(float x0, float y0, int* indmin, int* indmax)
+{
CbmRichHoughHit mpnt;
vector<CbmRichHoughHit>::iterator itmin, itmax;
//find all hits which are in the corridor
mpnt.fHit.fX = x0 - fMaxDistance;
- itmin =
- std::lower_bound(fData.begin(), fData.end(), mpnt, CbmRichHoughHitCmpUp());
+ itmin = std::lower_bound(fData.begin(), fData.end(), mpnt, CbmRichHoughHitCmpUp());
mpnt.fHit.fX = x0 + fMaxDistance;
- itmax =
- std::lower_bound(fData.begin(), fData.end(), mpnt, CbmRichHoughHitCmpUp())
- - 1;
+ itmax = std::lower_bound(fData.begin(), fData.end(), mpnt, CbmRichHoughHitCmpUp()) - 1;
*indmin = itmin - fData.begin();
*indmax = itmax - fData.begin();
@@ -225,16 +226,15 @@ void CbmRichRingFinderHoughImpl::DefineLocalAreaAndHits(float x0,
}
}
-void CbmRichRingFinderHoughImpl::HoughTransform(unsigned int indmin,
- unsigned int indmax) {
+void CbmRichRingFinderHoughImpl::HoughTransform(unsigned int indmin, unsigned int indmax)
+{
for (int iPart = 0; iPart < fNofParts; iPart++) {
HoughTransformGroup(indmin, indmax, iPart);
} //iPart
}
-void CbmRichRingFinderHoughImpl::HoughTransformGroup(unsigned int /*indmin*/,
- unsigned int /*indmax*/,
- int iPart) {
+void CbmRichRingFinderHoughImpl::HoughTransformGroup(unsigned int /*indmin*/, unsigned int /*indmax*/, int iPart)
+{
unsigned int nofHits = fHitInd[iPart].size();
float xcs, ycs; // xcs = xc - fCurMinX
float dx = 1.0f / fDx, dy = 1.0f / fDy, dr = 1.0f / fDr;
@@ -315,7 +315,8 @@ void CbmRichRingFinderHoughImpl::HoughTransformGroup(unsigned int /*indmin*/,
//hitIndPart.clear();
}
-void CbmRichRingFinderHoughImpl::FindPeak(int indmin, int indmax) {
+void CbmRichRingFinderHoughImpl::FindPeak(int indmin, int indmax)
+{
// Find MAX bin R and compare it with CUT
int maxBinR = -1, maxR = -1;
unsigned int size = fHistR.size();
@@ -408,9 +409,8 @@ void CbmRichRingFinderHoughImpl::FindPeak(int indmin, int indmax) {
delete ring2;
}
-void CbmRichRingFinderHoughImpl::RemoveHitsAroundRing(int indmin,
- int indmax,
- CbmRichRingLight* ring) {
+void CbmRichRingFinderHoughImpl::RemoveHitsAroundRing(int indmin, int indmax, CbmRichRingLight* ring)
+{
float rms = sqrt(ring->GetChi2() / ring->GetNofHits());
float dCut = fRmsCoeffEl * rms;
if (dCut > fMaxCutEl) dCut = fMaxCutEl;
@@ -425,7 +425,8 @@ void CbmRichRingFinderHoughImpl::RemoveHitsAroundRing(int indmin,
}
}
-void CbmRichRingFinderHoughImpl::RingSelection() {
+void CbmRichRingFinderHoughImpl::RingSelection()
+{
int nofRings = fFoundRings.size();
sort(fFoundRings.begin(), fFoundRings.end(), CbmRichRingComparatorMore());
set<unsigned int> usedHitsAll;
@@ -443,9 +444,7 @@ void CbmRichRingFinderHoughImpl::RingSelection() {
set<unsigned int>::iterator it = usedHitsAll.find(ring->GetHitId(iHit));
if (it != usedHitsAll.end()) { nofUsedHitsAll++; }
}
- if ((float) nofUsedHitsAll / (float) nofHits > fUsedHitsAllCut) {
- isGoodRingAll = false;
- }
+ if ((float) nofUsedHitsAll / (float) nofHits > fUsedHitsAllCut) { isGoodRingAll = false; }
if (isGoodRingAll) {
fFoundRings[iRing]->SetRecFlag(1);
@@ -464,8 +463,8 @@ void CbmRichRingFinderHoughImpl::RingSelection() {
goodRingIndex.clear();
}
-void CbmRichRingFinderHoughImpl::ReAssignSharedHits(int ringInd1,
- int ringInd2) {
+void CbmRichRingFinderHoughImpl::ReAssignSharedHits(int ringInd1, int ringInd2)
+{
CbmRichRingLight* ring1 = fFoundRings[ringInd1];
CbmRichRingLight* ring2 = fFoundRings[ringInd2];
int nofHits1 = ring1->GetNofHits();
@@ -487,16 +486,16 @@ void CbmRichRingFinderHoughImpl::ReAssignSharedHits(int ringInd1,
float ry2 = hitY - ring2->GetCenterY();
float dr2 = fabs(sqrt(rx2 * rx2 + ry2 * ry2) - ring2->GetRadius());
- if (dr1 > dr2) {
- ring1->RemoveHit(hitId1);
- } else {
+ if (dr1 > dr2) { ring1->RemoveHit(hitId1); }
+ else {
ring2->RemoveHit(hitId2);
}
} //iHit2
} //iHit1
}
-int CbmRichRingFinderHoughImpl::GetHitIndexById(unsigned int hitId) {
+int CbmRichRingFinderHoughImpl::GetHitIndexById(unsigned int hitId)
+{
unsigned int size = fData.size();
for (unsigned int i = 0; i < size; i++) {
if (fData[i].fHit.fId == hitId) return i;
@@ -504,11 +503,8 @@ int CbmRichRingFinderHoughImpl::GetHitIndexById(unsigned int hitId) {
return -1;
}
-void CbmRichRingFinderHoughImpl::CalculateRingParameters(float x[],
- float y[],
- float* xc,
- float* yc,
- float* r) {
+void CbmRichRingFinderHoughImpl::CalculateRingParameters(float x[], float y[], float* xc, float* yc, float* r)
+{
float t1, t2, t3, t4, t5, t6, t8, t9, t10, t11, t14, t16, t19, t21, t41;
t1 = x[1] * x[1];
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.h b/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.h
index 6a3cc250092a7007cbf0248c10aa28f75eca6069..c6d5078134f30a38e10fff4718e18bb91c9315eb 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.h
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHoughImpl.h
@@ -33,8 +33,7 @@ class CbmRichRingSelectAnn;
class CbmRichRingFinderHoughImpl {
protected:
- static const unsigned short MAX_NOF_HITS =
- 65000; // maximum number of hits in RICH detector
+ static const unsigned short MAX_NOF_HITS = 65000; // maximum number of hits in RICH detector
// parameters of the Hough Transform algorithm
unsigned short fNofParts; // number of groups of hits for HT
@@ -54,12 +53,10 @@ protected:
unsigned short fNofBinsY; // number of bins in Y direction
unsigned short fNofBinsXY; // fNofBinsX*fNofBinsY
- unsigned short
- fHTCut; // cut number of entries in maximum bin of XY histogram
+ unsigned short fHTCut; // cut number of entries in maximum bin of XY histogram
unsigned short fNofBinsR; // number of bins in radius histogram
- unsigned short
- fHTCutR; // cut number of entries in maximum bin of Radius histogram
+ unsigned short fHTCutR; // cut number of entries in maximum bin of Radius histogram
unsigned short fMinNofHitsInArea; // minimum number of hits in the local area
@@ -78,11 +75,10 @@ protected:
bool fUseAnnSelect;
- vector<CbmRichHoughHit> fData; // Rich hits
- vector<unsigned short> fHist; // XY histogram
- vector<unsigned short> fHistR; // Radius histogram
- vector<vector<unsigned int>>
- fHitInd; // store hit indexes for different group of hits
+ vector<CbmRichHoughHit> fData; // Rich hits
+ vector<unsigned short> fHist; // XY histogram
+ vector<unsigned short> fHistR; // Radius histogram
+ vector<vector<unsigned int>> fHitInd; // store hit indexes for different group of hits
vector<CbmRichRingLight*> fFoundRings; // collect found rings
CbmRichRingFitterCOP* fFitCOP; // COP ring fitter
CbmRichRingSelectAnn* fANNSelect; // ANN selection criteria
@@ -113,8 +109,7 @@ public:
* @param[out] yc Y coordinate of the ring center.
* @param[out] r Ring radius.
*/
- void
- CalculateRingParameters(float x[], float y[], float* xc, float* yc, float* r);
+ void CalculateRingParameters(float x[], float y[], float* xc, float* yc, float* r);
/**
* \brief Run HT for each hit.
@@ -128,8 +123,7 @@ public:
* \param[out] indmin Minimum index of the hit in local area.
* \param[out] indmax Maximum index of the hit in local area.
*/
- virtual void
- DefineLocalAreaAndHits(float x0, float y0, int* indmin, int* indmax);
+ virtual void DefineLocalAreaAndHits(float x0, float y0, int* indmin, int* indmax);
/**
* \brief Run HoughTransformGroup for each group of hits.
@@ -144,8 +138,7 @@ public:
* \param[in] indmax Maximum index of the hit in local area.
* \param[in] iPart Index of the hit group.
*/
- virtual void
- HoughTransformGroup(unsigned int indmin, unsigned int indmax, int iPart);
+ virtual void HoughTransformGroup(unsigned int indmin, unsigned int indmax, int iPart);
/**
* \brief Find peak in the HT histograms.
@@ -194,7 +187,8 @@ public:
* \brief Set array of hits.
* \param[in] data Array of hits.
*/
- void SetData(const vector<CbmRichHoughHit>& data) {
+ void SetData(const vector<CbmRichHoughHit>& data)
+ {
fData.clear();
fData = data;
}
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.cxx b/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.cxx
index d38f20de01dcbd885fadca7fc5a7f69cf81de2a0..cf429e4531ae6ce4c6c1ec8f250b59915b43721f 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.cxx
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.cxx
@@ -7,14 +7,15 @@
//#include "../L1/L1Algo/L1Types.h"
#include "CbmRichRingFinderHoughSimd.h"
-#include "../L1/L1Algo/vectors/P4_F32vec4.h"
+
#include <emmintrin.h>
+#include "../L1/L1Algo/vectors/P4_F32vec4.h"
+
CbmRichRingFinderHoughSimd::CbmRichRingFinderHoughSimd() {}
-void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
- unsigned short int indmax,
- Int_t iPart) {
+void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin, unsigned short int indmax, Int_t iPart)
+{
// Float_t r12, r13, r23;
// Float_t rx0, rx1, rx2, ry0, ry1,ry2; //rx[3], ry[3];//, x[3], y[3];
//Float_t xc, yc, r;
@@ -66,10 +67,8 @@ void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
iH13 = hitIndPart[iHit1 + 2];
iH14 = hitIndPart[iHit1 + 3];
- iH1X =
- fvec(fData[iH11]->fX, fData[iH12]->fX, fData[iH13]->fX, fData[iH14]->fX);
- iH1Y =
- fvec(fData[iH11]->fY, fData[iH12]->fY, fData[iH13]->fY, fData[iH14]->fY);
+ iH1X = fvec(fData[iH11]->fX, fData[iH12]->fX, fData[iH13]->fX, fData[iH14]->fX);
+ iH1Y = fvec(fData[iH11]->fY, fData[iH12]->fY, fData[iH13]->fY, fData[iH14]->fY);
for (unsigned short int iHit2 = iHit1 + 1; iHit2 < nofHits; iHit2++) {
iH2 = hitIndPart[iHit2];
@@ -80,16 +79,12 @@ void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
rx0 = iH1X - iH2X; //rx12
ry0 = iH1Y - iH2Y; //ry12
r12 = rx0 * rx0 + ry0 * ry0;
- if (_mm_comineq_ss(_mm_cmpgt_ss(r12, fMinDistanceSqv), fvec0) == 1)
- continue;
- if (_mm_comineq_ss(_mm_cmplt_ss(r12, fMaxDistanceSqv), fvec0) == 1)
- continue;
+ if (_mm_comineq_ss(_mm_cmpgt_ss(r12, fMinDistanceSqv), fvec0) == 1) continue;
+ if (_mm_comineq_ss(_mm_cmplt_ss(r12, fMaxDistanceSqv), fvec0) == 1) continue;
- t10 = fvec(fData[iH11]->fX2plusY2 - fData[iH2]->fX2plusY2,
- fData[iH12]->fX2plusY2 - fData[iH2]->fX2plusY2,
- fData[iH13]->fX2plusY2 - fData[iH2]->fX2plusY2,
- fData[iH14]->fX2plusY2 - fData[iH2]->fX2plusY2);
+ t10 = fvec(fData[iH11]->fX2plusY2 - fData[iH2]->fX2plusY2, fData[iH12]->fX2plusY2 - fData[iH2]->fX2plusY2,
+ fData[iH13]->fX2plusY2 - fData[iH2]->fX2plusY2, fData[iH14]->fX2plusY2 - fData[iH2]->fX2plusY2);
for (unsigned short int iHit3 = iHit2 + 1; iHit3 < nofHits; iHit3++) {
//iH3 = hitIndPart[iHit3];
h3 = fDataV[hitIndPart[iHit3]];
@@ -100,18 +95,14 @@ void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
rx1 = iH1X - iH3X; //rx13
ry1 = iH1Y - iH3Y; //ry13
r13 = rx1 * rx1 + ry1 * ry1;
- if (_mm_comineq_ss(_mm_cmpgt_ss(r13, fMinDistanceSqv), fvec0) == 1)
- continue;
- if (_mm_comineq_ss(_mm_cmplt_ss(r13, fMaxDistanceSqv), fvec0) == 1)
- continue;
+ if (_mm_comineq_ss(_mm_cmpgt_ss(r13, fMinDistanceSqv), fvec0) == 1) continue;
+ if (_mm_comineq_ss(_mm_cmplt_ss(r13, fMaxDistanceSqv), fvec0) == 1) continue;
rx2 = iH2X - h3.fX; //rx23
ry2 = iH2Y - h3.fY; //ry23
r23 = rx2 * rx2 + ry2 * ry2;
- if (_mm_comineq_ss(_mm_cmpgt_ss(r23, fMinDistanceSqv), fvec0) == 1)
- continue;
- if (_mm_comineq_ss(_mm_cmplt_ss(r23, fMaxDistanceSqv), fvec0) == 1)
- continue;
+ if (_mm_comineq_ss(_mm_cmpgt_ss(r23, fMinDistanceSqv), fvec0) == 1) continue;
+ if (_mm_comineq_ss(_mm_cmplt_ss(r23, fMaxDistanceSqv), fvec0) == 1) continue;
t19 = fvec05 / (rx2 * ry0 - rx0 * ry2);
@@ -151,25 +142,21 @@ void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
indXY = (int*) &indXYv;
intR = (int*) &intRv;
- if (intR[0] > 9 && intR[0] < fNofBinsR && indXY[0] > -1
- && indXY[0] < fNofBinsXY) {
+ if (intR[0] > 9 && intR[0] < fNofBinsR && indXY[0] > -1 && indXY[0] < fNofBinsXY) {
fHist[indXY[0]]++;
fHistR[intR[0]]++;
}
- if (intR[1] > 9 && intR[1] < fNofBinsR && indXY[1] > -1
- && indXY[1] < fNofBinsXY) {
+ if (intR[1] > 9 && intR[1] < fNofBinsR && indXY[1] > -1 && indXY[1] < fNofBinsXY) {
fHist[indXY[1]]++;
fHistR[intR[1]]++;
}
- if (intR[2] > 9 && intR[2] < fNofBinsR && indXY[2] > -1
- && indXY[2] < fNofBinsXY) {
+ if (intR[2] > 9 && intR[2] < fNofBinsR && indXY[2] > -1 && indXY[2] < fNofBinsXY) {
fHist[indXY[2]]++;
fHistR[intR[2]]++;
}
- if (intR[3] > 9 && intR[3] < fNofBinsR && indXY[3] > -1
- && indXY[3] < fNofBinsXY) {
+ if (intR[3] > 9 && intR[3] < fNofBinsR && indXY[3] > -1 && indXY[3] < fNofBinsXY) {
fHist[indXY[3]]++;
fHistR[intR[3]]++;
}
@@ -179,7 +166,8 @@ void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
}
-void CbmRichRingFinderHoughSimd::HoughTransformReconstruction() {
+void CbmRichRingFinderHoughSimd::HoughTransformReconstruction()
+{
Int_t indmin, indmax;
Float_t x0, y0;
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.h b/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.h
index d95125ed9eafc84bd34899d25d4eef9fdb76283d..9723cb9cee2de741e30c44595a2fe301f0ced37e 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.h
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHoughSimd.h
@@ -40,9 +40,7 @@ public:
virtual void HoughTransformReconstruction();
- virtual void HoughTransformGroup(unsigned short int indmin,
- unsigned short int indmax,
- Int_t iPart);
+ virtual void HoughTransformGroup(unsigned short int indmin, unsigned short int indmax, Int_t iPart);
std::vector<CbmRichHoughHitVec> fDataV;
};
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderIdeal.cxx b/reco/detectors/rich/finder/CbmRichRingFinderIdeal.cxx
index ae9c6d4c8e63c295be4b543bcd71f104a4fbe752..40dfe5d450eedbc83f12603544fbfd1953815b56 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderIdeal.cxx
+++ b/reco/detectors/rich/finder/CbmRichRingFinderIdeal.cxx
@@ -34,62 +34,50 @@
using namespace std;
-CbmRichRingFinderIdeal::CbmRichRingFinderIdeal()
- : CbmRichRingFinder()
- , fRichPoints(NULL)
- , fMcTracks(NULL)
- , fEventList(NULL)
- , fDigiMan(nullptr) {}
+CbmRichRingFinderIdeal::CbmRichRingFinderIdeal() {}
CbmRichRingFinderIdeal::~CbmRichRingFinderIdeal() {}
-void CbmRichRingFinderIdeal::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichRingFinderIdeal::Init", "RootManager is NULL!");
- }
+void CbmRichRingFinderIdeal::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (nullptr == manager) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): FairRootManager is nullptr.";
- CbmMCDataManager* mcManager =
- (CbmMCDataManager*) ioman->GetObject("MCDataManager");
- if (mcManager == nullptr)
- LOG(fatal) << "CbmRichRingFinderIdeal::Init() NULL MCDataManager.";
+ CbmMCDataManager* mcManager = (CbmMCDataManager*) manager->GetObject("MCDataManager");
+ if (mcManager == nullptr) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): MCDataManager is nullptr.";
fDigiMan = CbmDigiManager::Instance();
fDigiMan->Init();
- if (!fDigiMan->IsPresent(ECbmModuleId::kRich)) {
- Fatal("CbmRichHitProducer::Init", "No RichDigi array!");
- }
- if (!fDigiMan->IsMatchPresent(ECbmModuleId::kRich)) {
- Fatal("CbmRichHitProducer::Init", "No RichMatchDigi array!");
- }
+ if (!fDigiMan->IsPresent(ECbmModuleId::kRich)) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): No RichDigi.";
+
+ if (!fDigiMan->IsMatchPresent(ECbmModuleId::kRich)) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): No RichMatchDigi.";
fMcTracks = mcManager->InitBranch("MCTrack");
- if (NULL == fMcTracks) {
- LOG(fatal) << "CbmRichRingFinderIdeal::Init No MCTrack!";
- }
+ if (fMcTracks == nullptr) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): No MCTrack.";
fRichPoints = mcManager->InitBranch("RichPoint");
- if (NULL == fRichPoints) {
- LOG(fatal) << "CbmRichRingFinderIdeal::Init No RichPoint!";
- }
+ if (fRichPoints == nullptr) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): No RichPoint!";
- fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList.");
- if (NULL == fEventList) {
- LOG(fatal) << "CbmRichRingFinderIdeal::Init No MCEventList!";
- }
+ fEventList = (CbmMCEventList*) manager->GetObject("MCEventList.");
+ if (fEventList == nullptr) LOG(fatal) << "CbmRichRingFinderIdeal::Init(): No MCEventList.";
}
-Int_t CbmRichRingFinderIdeal::DoFind(TClonesArray* hitArray,
- TClonesArray* /*projArray*/,
- TClonesArray* ringArray) {
- if (NULL == hitArray) {
- cout << "-E- CbmRichRingFinderIdeal::DoFind, RichHit array missing!"
- << endl;
+Int_t CbmRichRingFinderIdeal::DoFind(CbmEvent* event, TClonesArray* hitArray, TClonesArray* /*projArray*/,
+ TClonesArray* ringArray)
+{
+
+ if (event != nullptr) {
+ LOG(fatal) << "CbmRichRingFinderIdeal::DoFind(): CbmEvent is not nullptr. "
+ "This class does not support time-based mode. Please switch to event-by-event mode.";
+ }
+
+ if (hitArray == nullptr) {
+ LOG(error) << "CbmRichRingFinderIdeal::DoFind(), hitArray is nullptr.";
return -1;
}
- if (NULL == ringArray) {
- cout << "-E- CbmRichRingFinderIdeal::DoFind, Ring array missing!" << endl;
+ if (ringArray == nullptr) {
+ LOG(error) << "CbmRichRingFinderIdeal::DoFind(): ringArray is nullptr.";
return -1;
}
@@ -98,11 +86,10 @@ Int_t CbmRichRingFinderIdeal::DoFind(TClonesArray* hitArray,
Int_t nofRichHits = hitArray->GetEntriesFast();
for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
const CbmRichHit* richHit = static_cast<CbmRichHit*>(hitArray->At(iHit));
- if (NULL == richHit) continue;
+ if (richHit == nullptr) continue;
Int_t eventId = GetEventIdForRichHit(richHit);
vector<pair<Int_t, Int_t>> motherIds =
- CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(
- fDigiMan, richHit, fRichPoints, fMcTracks, eventId);
+ CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(fDigiMan, richHit, fRichPoints, fMcTracks, eventId);
for (UInt_t i = 0; i < motherIds.size(); i++) {
hitMap[motherIds[i]]++;
}
@@ -119,9 +106,8 @@ Int_t CbmRichRingFinderIdeal::DoFind(TClonesArray* hitArray,
// Create RichRings for McTracks
Int_t nofMcTracks = fMcTracks->Size(fileId, eventId);
for (Int_t iT = 0; iT < nofMcTracks; iT++) {
- const CbmMCTrack* mcTrack =
- static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iT));
- if (NULL == mcTrack) continue;
+ const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iT));
+ if (mcTrack == nullptr) continue;
pair<Int_t, Int_t> val = std::make_pair(eventId, iT);
if (hitMap[val] <= 0) continue;
new ((*ringArray)[nofRings]) CbmRichRing();
@@ -132,35 +118,32 @@ Int_t CbmRichRingFinderIdeal::DoFind(TClonesArray* hitArray,
// Loop over RichHits. Get corresponding MCPoint and MCTrack index
for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
const CbmRichHit* richHit = static_cast<CbmRichHit*>(hitArray->At(iHit));
- if (NULL == richHit) continue;
+ if (richHit == nullptr) continue;
Int_t eventId = GetEventIdForRichHit(richHit);
vector<pair<Int_t, Int_t>> motherIds =
- CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(
- fDigiMan, richHit, fRichPoints, fMcTracks, eventId);
+ CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(fDigiMan, richHit, fRichPoints, fMcTracks, eventId);
for (UInt_t i = 0; i < motherIds.size(); i++) {
if (ringMap.find(motherIds[i]) == ringMap.end()) continue;
Int_t ringIndex = ringMap[motherIds[i]];
CbmRichRing* ring = (CbmRichRing*) ringArray->At(ringIndex);
- if (NULL == ring) continue;
+ if (ring == nullptr) continue;
ring->AddHit(iHit); // attach the hit to the ring
}
}
- LOG(info) << "-I- CbmRichRingFinderIdeal nofRings:" << nofRings;
-
return nofRings;
}
-Int_t CbmRichRingFinderIdeal::GetEventIdForRichHit(const CbmRichHit* richHit) {
- if (richHit == NULL) return -1;
+Int_t CbmRichRingFinderIdeal::GetEventIdForRichHit(const CbmRichHit* richHit)
+{
+ if (richHit == nullptr) return -1;
Int_t digiIndex = richHit->GetRefId();
if (digiIndex < 0) return -1;
- const CbmMatch* digiMatch =
- fDigiMan->GetMatch(ECbmModuleId::kRich, digiIndex);
+ const CbmMatch* digiMatch = fDigiMan->GetMatch(ECbmModuleId::kRich, digiIndex);
if (NULL == digiMatch) return -1;
return digiMatch->GetMatchedLink().GetEntry();
}
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderIdeal.h b/reco/detectors/rich/finder/CbmRichRingFinderIdeal.h
index 47042bfc7cbf8eac8eeb17c35e754a9a140a49f6..b142ab4e45b95711f2c66b0e1bb1b99ca6f23d49 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderIdeal.h
+++ b/reco/detectors/rich/finder/CbmRichRingFinderIdeal.h
@@ -12,6 +12,7 @@
#define CBM_RICH_RING_FINDER_IDEAL
#include "CbmRichRingFinder.h"
+
#include <vector>
using namespace std;
@@ -20,14 +21,9 @@ class CbmRichHit;
class CbmMCDataArray;
class CbmMCEventList;
class CbmDigiManager;
+class CbmEvent;
class CbmRichRingFinderIdeal : public CbmRichRingFinder {
-private:
- CbmMCDataArray* fRichPoints;
- CbmMCDataArray* fMcTracks;
- CbmMCEventList* fEventList;
- CbmDigiManager* fDigiMan;
-
public:
/**
@@ -48,11 +44,14 @@ public:
/**
* Inherited from CbmRichRingFinder.
*/
- virtual int DoFind(TClonesArray* hitArray,
- TClonesArray* projArray,
- TClonesArray* ringArray);
+ virtual int DoFind(CbmEvent* event, TClonesArray* hitArray, TClonesArray* projArray, TClonesArray* ringArray);
private:
+ CbmMCDataArray* fRichPoints = nullptr;
+ CbmMCDataArray* fMcTracks = nullptr;
+ CbmMCEventList* fEventList = nullptr;
+ CbmDigiManager* fDigiMan = nullptr;
+
/**
* \ brief Return evnetId from digiMatch corresponding to rich hit.
*/
diff --git a/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.cxx b/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.cxx
index 10d4e0c63de76997be45c70e852d0a7b4a7c3ce7..029740418ebb9618f9817a74875c0da6c5709891 100644
--- a/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.cxx
+++ b/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.cxx
@@ -7,6 +7,7 @@
#include "CbmRichProjectionProducerAnalytical.h"
+#include "CbmEvent.h"
#include "CbmMCTrack.h"
#include "CbmRichGeoManager.h"
@@ -33,78 +34,73 @@ using std::cout;
using std::endl;
-CbmRichProjectionProducerAnalytical::CbmRichProjectionProducerAnalytical()
- : fTrackParams(NULL), fNHits(0), fEventNum(0) {}
+CbmRichProjectionProducerAnalytical::CbmRichProjectionProducerAnalytical() {}
-CbmRichProjectionProducerAnalytical::~CbmRichProjectionProducerAnalytical() {
- FairRootManager* fManager = FairRootManager::Instance();
- fManager->Write();
-}
+CbmRichProjectionProducerAnalytical::~CbmRichProjectionProducerAnalytical() {}
-void CbmRichProjectionProducerAnalytical::Init() {
- LOG(info) << "CbmRichProjectionProducerAnalytical::Init()";
+void CbmRichProjectionProducerAnalytical::Init()
+{
FairRootManager* manager = FairRootManager::Instance();
+ if (nullptr == manager) LOG(fatal) << "CbmRichProjectionProducerAnalytical::Init(): FairRootManager is nullptr.";
fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
- if (NULL == fTrackParams) { LOG(fatal) << "No RichTrackParamZ array!"; }
+ if (fTrackParams == nullptr) LOG(fatal) << "CbmRichProjectionProducerAnalytical::Init(): No RichTrackParamZ array.";
}
-void CbmRichProjectionProducerAnalytical::DoProjection(TClonesArray* richProj) {
+void CbmRichProjectionProducerAnalytical::DoProjection(CbmEvent* event, TClonesArray* richProj)
+{
+
+ if (richProj == nullptr) {
+ LOG(error) << "CbmRichProjectionProducerAnalytical::DoExtrapolation(): richProj is nullptr.";
+ return;
+ }
+
fEventNum++;
- cout << "CbmRichProjectionProducerAnalytical:: event " << fEventNum << endl;
CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
- Double_t mirrorX = gp->fMirrorX;
- Double_t mirrorY = gp->fMirrorY;
- Double_t mirrorZ = gp->fMirrorZ;
- Double_t mirrorR = gp->fMirrorR;
+ Double_t mirX = gp->fMirrorX;
+ Double_t mirY = gp->fMirrorY;
+ Double_t mirZ = gp->fMirrorZ;
+ Double_t mirR = gp->fMirrorR;
- richProj->Delete();
TMatrixFSym covMat(5);
for (Int_t i = 0; i < 5; i++) {
for (Int_t j = 0; j <= i; j++) {
covMat(i, j) = 0;
}
}
- covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
- 1.e-4;
+ covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) = 1.e-4;
- for (Int_t j = 0; j < fTrackParams->GetEntriesFast(); j++) {
- FairTrackParam* point = (FairTrackParam*) fTrackParams->At(j);
- new ((*richProj)[j]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+ Int_t nofTrackParams = event ? event->GetNofData(ECbmDataType::kRichTrackParamZ) : fTrackParams->GetEntriesFast();
+ for (Int_t iT0 = 0; iT0 < nofTrackParams; iT0++) {
+ Int_t iT = event ? event->GetIndex(ECbmDataType::kRichTrackParamZ, iT0) : iT0;
+ FairTrackParam* trPar = static_cast<FairTrackParam*>(fTrackParams->At(iT));
+ new ((*richProj)[iT]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+ if (event != nullptr) event->AddData(ECbmDataType::kRichTrackProjection, iT);
// check if Array was filled
- if (point->GetX() == 0 && point->GetY() == 0 && point->GetZ() == 0
- && point->GetTx() == 0 && point->GetTy() == 0)
+ if (trPar->GetX() == 0 && trPar->GetY() == 0 && trPar->GetZ() == 0 && trPar->GetTx() == 0 && trPar->GetTy() == 0)
continue;
- if (point->GetQp() == 0) continue;
+ if (trPar->GetQp() == 0) continue;
Double_t rho1 = 0.;
TVector3 startP, momP, crossP, centerP;
- Double_t p = 1. / TMath::Abs(point->GetQp());
+ Double_t p = 1. / TMath::Abs(trPar->GetQp());
Double_t pz;
- if ((1 + point->GetTx() * point->GetTx() + point->GetTy() * point->GetTy())
- > 0.)
- pz = p
- / TMath::Sqrt(1 + point->GetTx() * point->GetTx()
- + point->GetTy() * point->GetTy());
+ Double_t pz2 = 1 + trPar->GetTx() * trPar->GetTx() + trPar->GetTy() * trPar->GetTy();
+ if (pz2 > 0.) pz = p / TMath::Sqrt(pz2);
else {
- cout << " -E- RichProjectionProducer: strange value for calculating pz: "
- << (1 + point->GetTx() * point->GetTx()
- + point->GetTy() * point->GetTy())
- << endl;
+ LOG(error) << "CbmRichProjectionProducerAnalytical::DoProjection(): strange value for calculating pz: " << pz2;
pz = 0.;
}
- Double_t px = pz * point->GetTx();
- Double_t py = pz * point->GetTy();
+ Double_t px = pz * trPar->GetTx();
+ Double_t py = pz * trPar->GetTy();
momP.SetXYZ(px, py, pz);
- point->Position(startP);
- if ((mirrorY * startP.y()) < 0)
- mirrorY =
- -mirrorY; // check that mirror center and startP are in same hemisphere
+ trPar->Position(startP);
+ if ((mirY * startP.y()) < 0) mirY = -mirY; // check that mirror center and startP are in same hemisphere
// calculation of intersection of track with selected mirror
// corresponds to calculation of intersection between a straight line and a sphere:
@@ -114,23 +110,17 @@ void CbmRichProjectionProducerAnalytical::DoProjection(TClonesArray* richProj) {
// dist = r^2 - fR^2
// -> rho1 = (-RxP+sqrt(RxP^2-normP2*dist))/normP2 extrapolation factor for:
// intersection point crossP = startP + rho1 * momP
- Double_t RxP =
- (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
- + momP.z() * (startP.z() - mirrorZ));
- Double_t normP2 =
- (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
+ Double_t RxP = (momP.x() * (startP.x() - mirX) + momP.y() * (startP.y() - mirY) + momP.z() * (startP.z() - mirZ));
+ Double_t normP2 = (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
Double_t dist =
- (startP.x() * startP.x() + mirrorX * mirrorX + startP.y() * startP.y()
- + mirrorY * mirrorY + startP.z() * startP.z() + mirrorZ * mirrorZ
- - 2 * startP.x() * mirrorX - 2 * startP.y() * mirrorY
- - 2 * startP.z() * mirrorZ - mirrorR * mirrorR);
+ (startP.x() * startP.x() + mirX * mirX + startP.y() * startP.y() + mirY * mirY + startP.z() * startP.z()
+ + mirZ * mirZ - 2 * startP.x() * mirX - 2 * startP.y() * mirY - 2 * startP.z() * mirZ - mirR * mirR);
if ((RxP * RxP - normP2 * dist) > 0.) {
- if (normP2 != 0.)
- rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
- if (normP2 == 0)
- cout << " Error in track extrapolation: momentum = 0 " << endl;
- } else {
+ if (normP2 != 0.) rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
+ if (normP2 == 0) LOG(error) << "CbmRichProjectionProducerAnalytical::DoProjection(): normP2 == 0";
+ }
+ else {
//cout << " -E- RichProjectionProducer: RxP*RxP-normP2*dist = " << RxP*RxP-normP2*dist << endl;
}
@@ -141,25 +131,18 @@ void CbmRichProjectionProducerAnalytical::DoProjection(TClonesArray* richProj) {
// check if crosspoint with mirror and chosen mirrorcenter (y) are in same hemisphere
// if not recalculate crossing point
- if ((mirrorY * crossP.y()) < 0) {
- mirrorY = -mirrorY;
- RxP =
- (momP.x() * (startP.x() - mirrorX) + momP.y() * (startP.y() - mirrorY)
- + momP.z() * (startP.z() - mirrorZ));
- normP2 =
- (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
- dist =
- (startP.x() * startP.x() + mirrorX * mirrorX + startP.y() * startP.y()
- + mirrorY * mirrorY + startP.z() * startP.z() + mirrorZ * mirrorZ
- - 2 * startP.x() * mirrorX - 2 * startP.y() * mirrorY
- - 2 * startP.z() * mirrorZ - mirrorR * mirrorR);
+ if ((mirY * crossP.y()) < 0) {
+ mirY = -mirY;
+ RxP = (momP.x() * (startP.x() - mirX) + momP.y() * (startP.y() - mirY) + momP.z() * (startP.z() - mirZ));
+ normP2 = (momP.x() * momP.x() + momP.y() * momP.y() + momP.z() * momP.z());
+ dist = (startP.x() * startP.x() + mirX * mirX + startP.y() * startP.y() + mirY * mirY + startP.z() * startP.z()
+ + mirZ * mirZ - 2 * startP.x() * mirX - 2 * startP.y() * mirY - 2 * startP.z() * mirZ - mirR * mirR);
if ((RxP * RxP - normP2 * dist) > 0.) {
- if (normP2 != 0.)
- rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
- if (normP2 == 0)
- cout << " Error in track extrapolation: momentum = 0 " << endl;
- } else {
+ if (normP2 != 0.) rho1 = (-RxP + TMath::Sqrt(RxP * RxP - normP2 * dist)) / normP2;
+ if (normP2 == 0) LOG(error) << "CbmRichProjectionProducerAnalytical::DoProjection(): normP2 == 0";
+ }
+ else {
//cout << " -E- RichProjectionProducer: RxP*RxP-normP2*dist = " << RxP*RxP-normP2*dist << endl;
}
@@ -169,26 +152,20 @@ void CbmRichProjectionProducerAnalytical::DoProjection(TClonesArray* richProj) {
crossP.SetXYZ(crossPx, crossPy, crossPz);
}
- centerP.SetXYZ(mirrorX, mirrorY, mirrorZ); // mirror center
+ centerP.SetXYZ(mirX, mirY, mirZ); // mirror center
// calculate normal on crosspoint with mirror
- TVector3 normP(crossP.x() - centerP.x(),
- crossP.y() - centerP.y(),
- crossP.z() - centerP.z());
+ TVector3 normP(crossP.x() - centerP.x(), crossP.y() - centerP.y(), crossP.z() - centerP.z());
normP = normP.Unit();
// check that normal has same z-direction as momentum
- if ((normP.z() * momP.z()) < 0.)
- normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
+ if ((normP.z() * momP.z()) < 0.) normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
// reflect track
- Double_t np =
- normP.x() * momP.x() + normP.y() * momP.y() + normP.z() * momP.z();
+ Double_t np = normP.x() * momP.x() + normP.y() * momP.y() + normP.z() * momP.z();
TVector3 ref;
- ref.SetXYZ(2 * np * normP.x() - momP.x(),
- 2 * np * normP.y() - momP.y(),
- 2 * np * normP.z() - momP.z());
+ ref.SetXYZ(2 * np * normP.x() - momP.x(), 2 * np * normP.y() - momP.y(), 2 * np * normP.z() - momP.z());
if (ref.z() != 0.) {
@@ -202,33 +179,30 @@ void CbmRichProjectionProducerAnalytical::DoProjection(TClonesArray* richProj) {
CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
//cout << "inPoint:" << inPos.X() << " " << inPos.Y() << " " << inPos.Z() << " outPoint:" << outPos.X() << " " << outPos.Y() << " " << outPos.Z() << endl;
-
-
- } else if (gp->fGeometryType == CbmRichGeometryTypeTwoWings) {
+ }
+ else if (gp->fGeometryType == CbmRichGeometryTypeTwoWings) {
GetPmtIntersectionPointTwoWings(¢erP, &crossP, &ref, &inPos);
// Transform intersection point in same way as MCPoints were transformed in HitProducer before stored as Hit
CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
- } else {
- Fatal("CbmRichProjectionProducerAnalytical ", "unnown geometry type");
+ }
+ else {
+ LOG(fatal) << "CbmRichProjectionProducerAnalytical::DoProjection(): unknown geometry type";
}
- Bool_t isInsidePmt =
- CbmRichGeoManager::GetInstance().IsPointInsidePmt(&outPos);
+ Bool_t isInsidePmt = CbmRichGeoManager::GetInstance().IsPointInsidePmt(&outPos);
if (isInsidePmt) {
- FairTrackParam richtrack(
- outPos.x(), outPos.y(), outPos.z(), 0., 0., 0., covMat);
- *(FairTrackParam*) (richProj->At(j)) = richtrack;
+ FairTrackParam richtrack(outPos.x(), outPos.y(), outPos.z(), 0., 0., 0., covMat);
+ *(FairTrackParam*) (richProj->At(iT)) = richtrack;
}
- } // if (refZ!=0.)
- } // j
+ }
+ }
}
-void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(
- const TVector3* centerP,
- const TVector3* crossP,
- const TVector3* ref,
- TVector3* outPoint) {
+void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(const TVector3* centerP,
+ const TVector3* crossP, const TVector3* ref,
+ TVector3* outPoint)
+{
// crosspoint whith photodetector plane:
// calculate intersection between straight line and (tilted) plane:
// normal on plane tilted by theta around x-axis: (0,-sin(theta),cos(theta)) = n
@@ -244,10 +218,8 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(
CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
if (gp->fGeometryType != CbmRichGeometryTypeTwoWings) {
- Fatal(
- "CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings ",
- "Wrong geometry, this method could be used only for "
- "CbmRichGeometryTypeTwoWings");
+ LOG(fatal) << "CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(): Wrong geometry, this method "
+ "could be used only for CbmRichGeometryTypeTwoWings";
}
Double_t pmtPhi = gp->fPmt.fPhi;
@@ -258,22 +230,18 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(
Double_t rho2 = 0.;
if (centerP->y() > 0) {
- rho2 =
- (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
- - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
- + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(pmtPhi) * ref->x()
- - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
- + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
+ rho2 = (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
+ - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
+ + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(pmtPhi) * ref->x() - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
+ + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
}
if (centerP->y() < 0) {
- rho2 =
- (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
- - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * (-pmtPlaneY - crossP->y())
- + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(pmtPhi) * ref->x()
- - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
- + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
+ rho2 = (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
+ - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * (-pmtPlaneY - crossP->y())
+ + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(pmtPhi) * ref->x() - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
+ + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
}
//rho2 = -1*(crossP.z() - fDetZ)/refZ; // only for theta = 0, phi=0
@@ -284,22 +252,16 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(
if (xX < 0) {
if (centerP->y() > 0) {
rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
- - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi)
- * (pmtPlaneY - crossP->y())
- + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi)
- * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(-pmtPhi) * ref->x()
- - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * (pmtPlaneY - crossP->y())
+ + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(-pmtPhi) * ref->x() - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
}
if (centerP->y() < 0) {
rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
- - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi)
- * (-pmtPlaneY - crossP->y())
- + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi)
- * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(-pmtPhi) * ref->x()
- - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * (-pmtPlaneY - crossP->y())
+ + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(-pmtPhi) * ref->x() - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
}
@@ -312,25 +274,20 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointTwoWings(
}
-void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(
- const TVector3* centerP,
- const TVector3* crossP,
- const TVector3* ref,
- TVector3* outPoint) {
+void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(const TVector3* centerP, const TVector3* crossP,
+ const TVector3* ref, TVector3* outPoint)
+{
CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
if (gp->fGeometryType != CbmRichGeometryTypeCylindrical) {
- Fatal("CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl ",
- "Wrong geometry, this method could be used only for "
- "CbmRichGeometryTypeCylindrical");
+ LOG(fatal) << "CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(): Wrong geometry, this method "
+ "could be used only for CbmRichGeometryTypeCylindrical";
}
Double_t xX = 0.;
Double_t yY = 0.;
Double_t zZ = 0.;
Double_t maxDist = 0.;
- for (map<string, CbmRichRecGeoParPmt>::iterator it = gp->fPmtMap.begin();
- it != gp->fPmtMap.end();
- it++) {
+ for (map<string, CbmRichRecGeoParPmt>::iterator it = gp->fPmtMap.begin(); it != gp->fPmtMap.end(); it++) {
Double_t pmtPlaneX = it->second.fPlaneX;
Double_t pmtPlaneY = it->second.fPlaneY;
Double_t pmtPlaneZ = it->second.fPlaneZ;
@@ -342,23 +299,17 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(
double rho2 = 0.;
if (centerP->y() > 0) {
- rho2 =
- (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
- - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
- + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi)
- * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(pmtPhi) * ref->x()
- - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
- + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
+ rho2 = (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
+ - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneY - crossP->y())
+ + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(pmtPhi) * ref->x() - TMath::Sin(pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
+ + TMath::Cos(pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
}
if (centerP->y() < 0) {
rho2 = (-TMath::Sin(pmtPhi) * (pmtPlaneX - crossP->x())
- - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi)
- * (-pmtPlaneY - crossP->y())
- + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi)
- * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(pmtPhi) * ref->x()
- - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
+ - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * (-pmtPlaneY - crossP->y())
+ + TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(pmtPhi) * ref->x() - TMath::Sin(-pmtTheta) * TMath::Cos(pmtPhi) * ref->y()
+ TMath::Cos(-pmtTheta) * TMath::Cos(pmtPhi) * ref->z());
}
@@ -369,22 +320,16 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(
if (cxX < 0) {
if (centerP->y() > 0) {
rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
- - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi)
- * (pmtPlaneY - crossP->y())
- + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi)
- * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(-pmtPhi) * ref->x()
- - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * (pmtPlaneY - crossP->y())
+ + TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(-pmtPhi) * ref->x() - TMath::Sin(pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ TMath::Cos(pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
}
if (centerP->y() < 0) {
rho2 = (-TMath::Sin(-pmtPhi) * (-pmtPlaneX - crossP->x())
- - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi)
- * (-pmtPlaneY - crossP->y())
- + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi)
- * (pmtPlaneZ - crossP->z()))
- / (-TMath::Sin(-pmtPhi) * ref->x()
- - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * (-pmtPlaneY - crossP->y())
+ + TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * (pmtPlaneZ - crossP->z()))
+ / (-TMath::Sin(-pmtPhi) * ref->x() - TMath::Sin(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->y()
+ TMath::Cos(-pmtTheta) * TMath::Cos(-pmtPhi) * ref->z());
}
@@ -393,8 +338,7 @@ void CbmRichProjectionProducerAnalytical::GetPmtIntersectionPointCyl(
czZ = crossP->z() + ref->z() * rho2;
}
- Double_t dP = TMath::Sqrt((crossP->x() - cxX) * (crossP->X() - cxX)
- + (crossP->y() - cyY) * (crossP->y() - cyY)
+ Double_t dP = TMath::Sqrt((crossP->x() - cxX) * (crossP->X() - cxX) + (crossP->y() - cyY) * (crossP->y() - cyY)
+ (crossP->z() - czZ) * (crossP->Z() - czZ));
if (dP >= maxDist) {
diff --git a/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.h b/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.h
index 27918ce756111cfd94da0eb1a51be6318aa79fdb..e089d3d6d16ef5a9dc8d01e9b13be50ee64dc088 100644
--- a/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.h
+++ b/reco/detectors/rich/tracks/CbmRichProjectionProducerAnalytical.h
@@ -11,6 +11,7 @@
#ifndef CBM_RICH_PROJECTION_PRODUCER_ANALYTICAL
#define CBM_RICH_PROJECTION_PRODUCER_ANALYTICAL
#include "CbmRichProjectionProducerBase.h"
+
#include "TObject.h"
#include "TVector3.h"
@@ -31,8 +32,7 @@ class FairTrackParam;
* \author S.Lebedev (initial version by P.Stolpovsky in 2005)
* \date 2016
**/
-class CbmRichProjectionProducerAnalytical :
- public CbmRichProjectionProducerBase {
+class CbmRichProjectionProducerAnalytical : public CbmRichProjectionProducerBase {
public:
/**
* \brief Standard constructor.
@@ -53,40 +53,32 @@ public:
* \brief Execute task.
* \param[out] richProj Output array of created projections.
*/
- virtual void DoProjection(TClonesArray* richProj);
+ virtual void DoProjection(CbmEvent* event, TClonesArray* richProj);
- void GetPmtIntersectionPointTwoWings(const TVector3* centerP,
- const TVector3* crossP,
- const TVector3* ref,
+ void GetPmtIntersectionPointTwoWings(const TVector3* centerP, const TVector3* crossP, const TVector3* ref,
TVector3* outPoint);
/*
* Find the intersection point with cylindrical PMT plane.
*
*/
- void GetPmtIntersectionPointCyl(const TVector3* centerP,
- const TVector3* crossP,
- const TVector3* ref,
+ void GetPmtIntersectionPointCyl(const TVector3* centerP, const TVector3* crossP, const TVector3* ref,
TVector3* outPoint);
private:
- TClonesArray* fTrackParams; // Starting points&directions
-
- int fNHits; // Number of hits
- int fEventNum; // number of events
+ TClonesArray* fTrackParams = nullptr;
+ int fEventNum = 0;
/**
* \brief Copy constructor.
*/
- CbmRichProjectionProducerAnalytical(
- const CbmRichProjectionProducerAnalytical&);
+ CbmRichProjectionProducerAnalytical(const CbmRichProjectionProducerAnalytical&);
/**
* \brief Assignment operator.
*/
- CbmRichProjectionProducerAnalytical&
- operator=(const CbmRichProjectionProducerAnalytical&);
+ CbmRichProjectionProducerAnalytical& operator=(const CbmRichProjectionProducerAnalytical&);
};
#endif
diff --git a/reco/detectors/rich/tracks/CbmRichProjectionProducerBase.h b/reco/detectors/rich/tracks/CbmRichProjectionProducerBase.h
index 48f1cd16b005683796859869ee57b040b14eeb74..0137cd317d57dd6cf43db949d92c7027051bbf27 100644
--- a/reco/detectors/rich/tracks/CbmRichProjectionProducerBase.h
+++ b/reco/detectors/rich/tracks/CbmRichProjectionProducerBase.h
@@ -11,6 +11,7 @@
#define CBM_RICH_PROJECTION_PRODUCER_BASE
class TClonesArray;
+class CbmEvent;
/**
* \class CbmRichProjectionProducerBase
@@ -42,7 +43,7 @@ public:
* Creates track projections onto the photodetector plane.
* \param[out] richProj Array of track projections onto the photodetector plane.
**/
- virtual void DoProjection(TClonesArray* richProj) = 0;
+ virtual void DoProjection(CbmEvent* event, TClonesArray* richProj) = 0;
private:
/**
@@ -53,8 +54,7 @@ private:
/**
* \brief Assignment operator.
*/
- CbmRichProjectionProducerBase&
- operator=(const CbmRichProjectionProducerBase&);
+ CbmRichProjectionProducerBase& operator=(const CbmRichProjectionProducerBase&);
};
#endif
diff --git a/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.cxx b/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.cxx
index 467f4302ffdcb890b5dfd68d5eb97d4e1b7a4593..1684f1bd9881391e10f9e94fee66dbbae9d35979 100644
--- a/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.cxx
+++ b/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.cxx
@@ -7,8 +7,10 @@
#include "CbmRichProjectionProducerTGeo.h"
+#include "CbmEvent.h"
#include "CbmMCTrack.h"
#include "CbmRichGeoManager.h"
+#include "utils/CbmRichNavigationUtil.h"
#include "FairGeoNode.h"
#include "FairGeoTransform.h"
@@ -29,54 +31,55 @@
#include <cmath>
-#include "utils/CbmRichNavigationUtil.h"
-
using std::cout;
using std::endl;
-CbmRichProjectionProducerTGeo::CbmRichProjectionProducerTGeo()
- : fTrackParams(NULL), fNHits(0), fEventNum(0) {}
+CbmRichProjectionProducerTGeo::CbmRichProjectionProducerTGeo() {}
-CbmRichProjectionProducerTGeo::~CbmRichProjectionProducerTGeo() {
- FairRootManager* fManager = FairRootManager::Instance();
- fManager->Write();
-}
+CbmRichProjectionProducerTGeo::~CbmRichProjectionProducerTGeo() {}
-void CbmRichProjectionProducerTGeo::Init() {
- LOG(info) << "CbmRichProjectionProducerTGeo::Init()";
+void CbmRichProjectionProducerTGeo::Init()
+{
FairRootManager* manager = FairRootManager::Instance();
+ if (nullptr == manager) LOG(fatal) << "CbmRichProjectionProducerTGeo::Init(): FairRootManager is nullptr.";
fTrackParams = (TClonesArray*) manager->GetObject("RichTrackParamZ");
- if (NULL == fTrackParams) { LOG(fatal) << "No RichTrackParamZ array!"; }
+ if (fTrackParams == nullptr) LOG(fatal) << "CbmRichProjectionProducerTGeo::Init(): No RichTrackParamZ array.";
}
-void CbmRichProjectionProducerTGeo::DoProjection(TClonesArray* richProj) {
+void CbmRichProjectionProducerTGeo::DoProjection(CbmEvent* event, TClonesArray* richProj)
+{
+
+ if (richProj == nullptr) {
+ LOG(error) << "CbmRichProjectionProducerTGeo::DoExtrapolation(): richProj is nullptr.";
+ return;
+ }
+
fEventNum++;
- LOG(info) << "CbmRichProjectionProducer: event " << fEventNum;
CbmRichRecGeoPar* gp = CbmRichGeoManager::GetInstance().fGP;
double mirrorX = gp->fMirrorX;
double mirrorY = gp->fMirrorY;
double mirrorZ = gp->fMirrorZ;
- richProj->Delete();
TMatrixFSym covMat(5);
for (Int_t i = 0; i < 5; i++) {
for (Int_t j = 0; j <= i; j++) {
covMat(i, j) = 0;
}
}
- covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
- 1.e-4;
+ covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) = 1.e-4;
- for (Int_t j = 0; j < fTrackParams->GetEntriesFast(); j++) {
- FairTrackParam* trackParam = (FairTrackParam*) fTrackParams->At(j);
- new ((*richProj)[j]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+ Int_t nofTrackParams = event ? event->GetNofData(ECbmDataType::kRichTrackParamZ) : fTrackParams->GetEntriesFast();
+ for (Int_t iT0 = 0; iT0 < nofTrackParams; iT0++) {
+ Int_t iT = event ? event->GetIndex(ECbmDataType::kRichTrackParamZ, iT0) : iT0;
+ FairTrackParam* trackParam = static_cast<FairTrackParam*>(fTrackParams->At(iT));
+ new ((*richProj)[iT]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+ if (event != nullptr) event->AddData(ECbmDataType::kRichTrackProjection, iT);
- if (trackParam->GetX() == 0 && trackParam->GetY() == 0
- && trackParam->GetZ() == 0 && trackParam->GetTx() == 0
+ if (trackParam->GetX() == 0 && trackParam->GetY() == 0 && trackParam->GetZ() == 0 && trackParam->GetTx() == 0
&& trackParam->GetTy() == 0)
continue;
if (trackParam->GetQp() == 0) continue;
@@ -87,30 +90,24 @@ void CbmRichProjectionProducerTGeo::DoProjection(TClonesArray* richProj) {
CbmRichNavigationUtil::GetDirCos(trackParam, nx, ny, nz);
dirCos.SetXYZ(nx, ny, nz);
- string volumeName = CbmRichNavigationUtil::FindIntersection(
- trackParam, crossP, "mirror_tile_type");
+ string volumeName = CbmRichNavigationUtil::FindIntersection(trackParam, crossP, "mirror_tile_type");
Bool_t mirrorIntersectionFound = (volumeName != string(""));
if (!mirrorIntersectionFound) continue;
// mirror center
- if (crossP.Y() > 0) {
- centerP.SetXYZ(mirrorX, mirrorY, mirrorZ);
- } else {
+ if (crossP.Y() > 0) { centerP.SetXYZ(mirrorX, mirrorY, mirrorZ); }
+ else {
centerP.SetXYZ(mirrorX, -mirrorY, mirrorZ);
}
// calculate normal on crosspoint with mirror
- TVector3 normP(crossP.x() - centerP.x(),
- crossP.y() - centerP.y(),
- crossP.z() - centerP.z());
+ TVector3 normP(crossP.x() - centerP.x(), crossP.y() - centerP.y(), crossP.z() - centerP.z());
normP = normP.Unit();
// check that normal has same z-direction as momentum
- if ((normP.z() * dirCos.z()) < 0.)
- normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
+ if ((normP.z() * dirCos.z()) < 0.) normP = TVector3(-1. * normP.x(), -1. * normP.y(), -1. * normP.z());
// reflect track
- Double_t np =
- normP.x() * dirCos.x() + normP.y() * dirCos.y() + normP.z() * dirCos.z();
+ Double_t np = normP.x() * dirCos.x() + normP.y() * dirCos.y() + normP.z() * dirCos.z();
Double_t refX = 2 * np * normP.x() - dirCos.x();
Double_t refY = 2 * np * normP.y() - dirCos.y();
Double_t refZ = 2 * np * normP.z() - dirCos.z();
@@ -119,8 +116,7 @@ void CbmRichProjectionProducerTGeo::DoProjection(TClonesArray* richProj) {
refl.Unit();
TVector3 pmtIntersectionPoint;
- volumeName = CbmRichNavigationUtil::FindIntersection(
- refl, crossP, pmtIntersectionPoint, "pmt");
+ volumeName = CbmRichNavigationUtil::FindIntersection(refl, crossP, pmtIntersectionPoint, "pmt");
Bool_t pmtIntersectionFound = (volumeName != string(""));
if (!pmtIntersectionFound) continue;
@@ -128,20 +124,17 @@ void CbmRichProjectionProducerTGeo::DoProjection(TClonesArray* richProj) {
// Transform intersection point in same way as MCPoints were
// transformed in HitProducer before stored as Hit:
TVector3 outPos;
- CbmRichGeoManager::GetInstance().RotatePoint(&pmtIntersectionPoint,
- &outPos);
+ CbmRichGeoManager::GetInstance().RotatePoint(&pmtIntersectionPoint, &outPos);
Double_t xDet = outPos.X();
Double_t yDet = outPos.Y();
Double_t zDet = outPos.Z();
//check that crosspoint inside the plane
- // if( xDet > (-fGP.fPmtXOrig-fGP.fPmtWidthX) && xDet < (fGP.fPmtXOrig+fGP.fPmtWidthX)){
- // if(TMath::Abs(yDet) > (fGP.fPmtY-fGP.fPmtWidthY) && TMath::Abs(yDet) < (fGP.fPmtY+fGP.fPmtWidthY)){
+ // if( xDet > (-fGP.fPmtXOrig-fGP.fPmtWidthX) && xDet < (fGP.fPmtXOrig+fGP.fPmtWidthX)){
+ // if(TMath::Abs(yDet) > (fGP.fPmtY-fGP.fPmtWidthY) && TMath::Abs(yDet) < (fGP.fPmtY+fGP.fPmtWidthY)){
FairTrackParam richtrack(xDet, yDet, zDet, 0., 0., 0., covMat);
- *(FairTrackParam*) (richProj->At(j)) = richtrack;
+ *(FairTrackParam*) (richProj->At(iT)) = richtrack;
// }
// }
} // j
-
- cout << "nofRichProjections:" << richProj->GetEntriesFast() << endl;
}
diff --git a/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.h b/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.h
index 72a29088eb4a92e61b7eb5f15a58dbe0b2b63fdc..38a3ed5b04f2264c0830bb9f72cc30740297c88f 100644
--- a/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.h
+++ b/reco/detectors/rich/tracks/CbmRichProjectionProducerTGeo.h
@@ -11,6 +11,7 @@
#ifndef CBM_RICH_PROJECTION_PRODUCER_TGEO
#define CBM_RICH_PROJECTION_PRODUCER_TGEO
#include "CbmRichProjectionProducerBase.h"
+
#include "TObject.h"
#include "TVector3.h"
@@ -52,14 +53,12 @@ public:
* \brief Execute task.
* \param[out] richProj Output array of created projections.
*/
- virtual void DoProjection(TClonesArray* richProj);
+ virtual void DoProjection(CbmEvent* event, TClonesArray* richProj);
private:
- TClonesArray* fTrackParams; // Starting points&directions
-
- int fNHits; // Number of hits
- int fEventNum; // number of events
+ TClonesArray* fTrackParams = nullptr;
+ int fEventNum = 0;
/**
* \brief Copy constructor.
@@ -69,8 +68,7 @@ private:
/**
* \brief Assignment operator.
*/
- CbmRichProjectionProducerTGeo&
- operator=(const CbmRichProjectionProducerTGeo&);
+ CbmRichProjectionProducerTGeo& operator=(const CbmRichProjectionProducerTGeo&);
};
#endif
diff --git a/reco/detectors/rich/tracks/CbmRichRingTrackAssignBase.h b/reco/detectors/rich/tracks/CbmRichRingTrackAssignBase.h
index 19d1e6ac97c2b616d9c8220c662da8013e2ed335..768d1298c28afe17b9a45a6d5d4ce471b3721dd3 100644
--- a/reco/detectors/rich/tracks/CbmRichRingTrackAssignBase.h
+++ b/reco/detectors/rich/tracks/CbmRichRingTrackAssignBase.h
@@ -11,6 +11,7 @@
#define CBM_RICH_RING_TRACK_ASSIGN_BASE
class TClonesArray;
+class CbmEvent;
/**
* \class CbmRichRingTrackAssignBase
@@ -25,7 +26,7 @@ public:
/**
* brief Standard constructor.
*/
- CbmRichRingTrackAssignBase() : fMaxDistance(999.), fMinNofHitsInRing(1) {}
+ CbmRichRingTrackAssignBase() {}
/**
* \brief Destructor.
@@ -43,12 +44,11 @@ public:
* \param[in] rings Array of RICH rings.
* \param[in] richProj Array of track projections onto the photodetector plane.
**/
- virtual void DoAssign(TClonesArray* rings, TClonesArray* richProj) = 0;
+ virtual void DoAssign(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj) = 0;
protected:
- double
- fMaxDistance; // max. distance between ring center and track extrapolation
- int fMinNofHitsInRing; // min number of hits per ring
+ double fMaxDistance = 999.; // max. distance between ring center and track extrapolation
+ int fMinNofHitsInRing = 1; // min number of hits per ring
private:
/**
diff --git a/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.cxx b/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.cxx
index 1100de85b02ddd32f8f1421b6e84fbcc71941a62..feee8de28bf1cab9ed676678dcb5a921e19802cf 100644
--- a/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.cxx
+++ b/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.cxx
@@ -7,13 +7,15 @@
#include "CbmRichRingTrackAssignClosestD.h"
-#include "CbmRichRing.h"
-
+#include "CbmEvent.h"
#include "CbmGlobalTrack.h"
#include "CbmMCTrack.h"
+#include "CbmRichRing.h"
#include "CbmTrdTrack.h"
+
#include "FairRootManager.h"
#include "FairTrackParam.h"
+#include <Logger.h>
#include "TClonesArray.h"
@@ -25,51 +27,50 @@ using std::cout;
using std::endl;
using std::vector;
-CbmRichRingTrackAssignClosestD::CbmRichRingTrackAssignClosestD()
- : fGlobalTracks(NULL)
- , fTrdTracks(NULL)
- , fTrdAnnCut(-0.5)
- , fUseTrd(false)
- , fAlgorithmType(RingTrack) {}
+CbmRichRingTrackAssignClosestD::CbmRichRingTrackAssignClosestD() {}
CbmRichRingTrackAssignClosestD::~CbmRichRingTrackAssignClosestD() {}
-void CbmRichRingTrackAssignClosestD::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichRingTrackAssignClosestD::Init",
- "RootManager not instantised!");
- }
+void CbmRichRingTrackAssignClosestD::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (nullptr == manager) LOG(fatal) << "CbmRichRingTrackAssignClosestD::Init(): FairRootManager is nullptr.";
- fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
- if (NULL == fGlobalTracks) {
- Fatal("CbmRichRingTrackAssignClosestD::Init", "No GlobalTrack array!");
- }
+ fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
+ if (fGlobalTracks == nullptr) LOG(fatal) << "CbmRichRingTrackAssignClosestD::Init(): No GlobalTrack.";
- fTrdTracks = (TClonesArray*) ioman->GetObject("TrdTrack");
+ fTrdTracks = (TClonesArray*) manager->GetObject("TrdTrack");
//if (NULL == fTrdTracks) {Fatal("CbmRichRingTrackAssignClosestD::Init", "No TrdTrack array!");}
}
-void CbmRichRingTrackAssignClosestD::DoAssign(TClonesArray* rings,
- TClonesArray* richProj) {
- if (fAlgorithmType == RingTrack) {
- DoAssignRingTrack(rings, richProj);
- } // RingTrack algorithm
+void CbmRichRingTrackAssignClosestD::DoAssign(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
+{
+ fEventNum++;
+ if (fAlgorithmType == RingTrack) { DoAssignRingTrack(event, rings, richProj); } // RingTrack algorithm
else if (fAlgorithmType == TrackRing) {
- DoAssignTrackRing(rings, richProj);
+ DoAssignTrackRing(event, rings, richProj);
} // TrackRing
else if (fAlgorithmType == Combined) {
- DoAssignRingTrack(rings, richProj);
- DoAssignTrackRing(rings, richProj);
+ DoAssignRingTrack(event, rings, richProj);
+ DoAssignTrackRing(event, rings, richProj);
} // Combined
+
+ Int_t nofTracks = event ? event->GetNofData(ECbmDataType::kRichTrackProjection) : richProj->GetEntriesFast();
+ Int_t nofRings = event ? event->GetNofData(ECbmDataType::kRichRing) : rings->GetEntriesFast();
+ LOG(info) << "CbmRichRingTrackAssignClosestD::DoAssign(): Event:" << fEventNum << " rings:" << nofRings
+ << " ringsInTS:" << rings->GetEntriesFast() << " tracks:" << nofTracks
+ << " tracksInTS:" << richProj->GetEntriesFast();
}
-void CbmRichRingTrackAssignClosestD::DoAssignRingTrack(TClonesArray* rings,
- TClonesArray* richProj) {
- Int_t nofTracks = richProj->GetEntriesFast();
- Int_t nofRings = rings->GetEntriesFast();
+void CbmRichRingTrackAssignClosestD::DoAssignRingTrack(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
+{
+
+ const Int_t nofTracks = event ? event->GetNofData(ECbmDataType::kRichTrackProjection) : richProj->GetEntriesFast();
+ const Int_t nofRings = event ? event->GetNofData(ECbmDataType::kRichRing) : rings->GetEntriesFast();
+ if (nofTracks <= 0 || nofRings <= 0) return;
+
vector<Int_t> trackIndex;
vector<Double_t> trackDist;
trackIndex.resize(nofRings);
@@ -80,40 +81,42 @@ void CbmRichRingTrackAssignClosestD::DoAssignRingTrack(TClonesArray* rings,
}
for (Int_t iIter = 0; iIter < 4; iIter++) {
- for (Int_t iRing = 0; iRing < nofRings; iRing++) {
- if (trackIndex[iRing] != -1) continue;
- CbmRichRing* pRing = (CbmRichRing*) rings->At(iRing);
- if (NULL == pRing) continue;
- if (pRing->GetNofHits() < fMinNofHitsInRing) continue;
-
- Double_t xRing = pRing->GetCenterX();
- Double_t yRing = pRing->GetCenterY();
+ for (Int_t iR0 = 0; iR0 < nofRings; iR0++) {
+ Int_t iR = event ? event->GetIndex(ECbmDataType::kRichRing, iR0) : iR0;
+
+ if (trackIndex[iR] != -1) continue;
+ CbmRichRing* ring = static_cast<CbmRichRing*>(rings->At(iR));
+ if (ring == nullptr) continue;
+ if (ring->GetNofHits() < fMinNofHitsInRing) continue;
+
+ Double_t xRing = ring->GetCenterX();
+ Double_t yRing = ring->GetCenterY();
Double_t rMin = 999.;
Int_t iTrackMin = -1;
- for (Int_t iTrack = 0; iTrack < nofTracks; iTrack++) {
- vector<Int_t>::iterator it =
- find(trackIndex.begin(), trackIndex.end(), iTrack);
+ for (Int_t iT0 = 0; iT0 < nofTracks; iT0++) {
+ Int_t iT = event ? event->GetIndex(ECbmDataType::kRichTrackProjection, iT0) : iT0;
+
+ vector<Int_t>::iterator it = find(trackIndex.begin(), trackIndex.end(), iT);
if (it != trackIndex.end()) continue;
- FairTrackParam* pTrack = (FairTrackParam*) richProj->At(iTrack);
- Double_t xTrack = pTrack->GetX();
- Double_t yTrack = pTrack->GetY();
+ FairTrackParam* track = static_cast<FairTrackParam*>(richProj->At(iT));
+ Double_t xTrack = track->GetX();
+ Double_t yTrack = track->GetY();
// no projection onto the photodetector plane
if (xTrack == 0 && yTrack == 0) continue;
- if (fUseTrd && fTrdTracks != NULL && !IsTrdElectron(iTrack)) continue;
+ if (fUseTrd && fTrdTracks != nullptr && !IsTrdElectron(iT)) continue;
- Double_t dist = TMath::Sqrt((xRing - xTrack) * (xRing - xTrack)
- + (yRing - yTrack) * (yRing - yTrack));
+ Double_t dist = TMath::Sqrt((xRing - xTrack) * (xRing - xTrack) + (yRing - yTrack) * (yRing - yTrack));
if (dist < rMin) {
rMin = dist;
- iTrackMin = iTrack;
+ iTrackMin = iT;
}
} // loop tracks
- trackIndex[iRing] = iTrackMin;
- trackDist[iRing] = rMin;
+ trackIndex[iR] = iTrackMin;
+ trackDist[iR] = rMin;
} //loop rings
for (UInt_t i1 = 0; i1 < trackIndex.size(); i1++) {
@@ -123,7 +126,8 @@ void CbmRichRingTrackAssignClosestD::DoAssignRingTrack(TClonesArray* rings,
if (trackDist[i1] >= trackDist[i2]) {
trackDist[i1] = 999.;
trackIndex[i1] = -1;
- } else {
+ }
+ else {
trackDist[i2] = 999.;
trackIndex[i2] = -1;
}
@@ -142,35 +146,36 @@ void CbmRichRingTrackAssignClosestD::DoAssignRingTrack(TClonesArray* rings,
}
}
-void CbmRichRingTrackAssignClosestD::DoAssignTrackRing(TClonesArray* rings,
- TClonesArray* richProj) {
- Int_t nofTracks = richProj->GetEntriesFast();
- Int_t nofRings = rings->GetEntriesFast();
-
- for (Int_t iTrack = 0; iTrack < nofTracks; iTrack++) {
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iTrack);
+void CbmRichRingTrackAssignClosestD::DoAssignTrackRing(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
+{
+ const Int_t nofTracks = event ? event->GetNofData(ECbmDataType::kRichTrackProjection) : richProj->GetEntriesFast();
+ const Int_t nofRings = event ? event->GetNofData(ECbmDataType::kRichRing) : rings->GetEntriesFast();
+ if (nofTracks <= 0 || nofRings <= 0) return;
+ for (Int_t iT0 = 0; iT0 < nofTracks; iT0++) {
+ Int_t iT = event ? event->GetIndex(ECbmDataType::kRichTrackProjection, iT0) : iT0;
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iT));
// track already has rich segment
- if (gTrack == NULL || gTrack->GetRichRingIndex() >= 0) continue;
+ if (gTrack == nullptr || gTrack->GetRichRingIndex() >= 0) continue;
- FairTrackParam* pTrack = (FairTrackParam*) richProj->At(iTrack);
- Double_t xTrack = pTrack->GetX();
- Double_t yTrack = pTrack->GetY();
+ FairTrackParam* track = static_cast<FairTrackParam*>(richProj->At(iT));
+ Double_t xTrack = track->GetX();
+ Double_t yTrack = track->GetY();
if (xTrack == 0 && yTrack == 0) continue;
- if (fUseTrd && fTrdTracks != NULL && !IsTrdElectron(iTrack)) continue;
+ if (fUseTrd && fTrdTracks != nullptr && !IsTrdElectron(iT)) continue;
Double_t rMin = 999.;
Int_t iRingMin = -1;
- for (Int_t iRing = 0; iRing < nofRings; iRing++) {
- CbmRichRing* pRing = (CbmRichRing*) rings->At(iRing);
- if (NULL == pRing) continue;
- if (pRing->GetNofHits() < fMinNofHitsInRing) continue;
-
- Double_t xRing = pRing->GetCenterX();
- Double_t yRing = pRing->GetCenterY();
- Double_t dist = TMath::Sqrt((xRing - xTrack) * (xRing - xTrack)
- + (yRing - yTrack) * (yRing - yTrack));
+ for (Int_t iR0 = 0; iR0 < nofRings; iR0++) {
+ Int_t iR = event ? event->GetIndex(ECbmDataType::kRichRing, iR0) : iR0;
+ CbmRichRing* ring = static_cast<CbmRichRing*>(rings->At(iR));
+ if (ring == nullptr) continue;
+ if (ring->GetNofHits() < fMinNofHitsInRing) continue;
+
+ Double_t xRing = ring->GetCenterX();
+ Double_t yRing = ring->GetCenterY();
+ Double_t dist = TMath::Sqrt((xRing - xTrack) * (xRing - xTrack) + (yRing - yTrack) * (yRing - yTrack));
if (dist < rMin) {
rMin = dist;
- iRingMin = iRing;
+ iRingMin = iR;
}
} // loop rings
if (iRingMin < 0) continue;
@@ -179,11 +184,12 @@ void CbmRichRingTrackAssignClosestD::DoAssignTrackRing(TClonesArray* rings,
} //loop tracks
}
-Bool_t CbmRichRingTrackAssignClosestD::IsTrdElectron(Int_t iTrack) {
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iTrack);
+Bool_t CbmRichRingTrackAssignClosestD::IsTrdElectron(Int_t iTrack)
+{
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
Int_t trdIndex = gTrack->GetTrdTrackIndex();
if (trdIndex == -1) return false;
- CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTracks->At(trdIndex);
+ CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(trdIndex));
if (NULL == trdTrack) return false;
if (trdTrack->GetPidANN() > fTrdAnnCut) { return true; }
diff --git a/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.h b/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.h
index d3280242cc5bcd96e701bf7f798f7d76e2c4adb0..126ec351611153da2681d916ed2334b80945cfa6 100644
--- a/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.h
+++ b/reco/detectors/rich/tracks/CbmRichRingTrackAssignClosestD.h
@@ -14,7 +14,8 @@
class TClonesArray;
-enum CbmRichRingTrackAssignClosestDAlgorithmEnum {
+enum CbmRichRingTrackAssignClosestDAlgorithmEnum
+{
TrackRing,
RingTrack,
Combined
@@ -48,25 +49,26 @@ public:
/**
* \brief Inherited from CbmRichRingTrackAssignBase.
*/
- void DoAssign(TClonesArray* rings, TClonesArray* richProj);
+ void DoAssign(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj);
/**
* \brief Implementation of the ring-track version of the algorithm.
*/
- void DoAssignRingTrack(TClonesArray* rings, TClonesArray* richProj);
+ void DoAssignRingTrack(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj);
/**
* \brief Implementation of the track-ring version of the algorithm.
*/
- void DoAssignTrackRing(TClonesArray* rings, TClonesArray* richProj);
+ void DoAssignTrackRing(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj);
private:
- TClonesArray* fGlobalTracks;
- TClonesArray* fTrdTracks;
+ TClonesArray* fGlobalTracks = nullptr;
+ TClonesArray* fTrdTracks = nullptr;
- double fTrdAnnCut; // ANN cut for electron identification in TRD
- bool fUseTrd; // if true electron identification in TRD will be performed
- CbmRichRingTrackAssignClosestDAlgorithmEnum fAlgorithmType;
+ double fTrdAnnCut = -0.5; // ANN cut for electron identification in TRD
+ bool fUseTrd = false; // if true electron identification in TRD will be used
+ CbmRichRingTrackAssignClosestDAlgorithmEnum fAlgorithmType = RingTrack;
+ int fEventNum = 0;
/**
* \brief Check if global track was identified as electron in the TRD detector.
@@ -83,8 +85,7 @@ private:
/**
* \brief Assignment operator.
*/
- CbmRichRingTrackAssignClosestD&
- operator=(const CbmRichRingTrackAssignClosestD&);
+ CbmRichRingTrackAssignClosestD& operator=(const CbmRichRingTrackAssignClosestD&);
};
#endif
diff --git a/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.cxx b/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.cxx
index c0e2a73d1531e5e2d5f1ed9cfc4eee077a6e39af..71ef0c1e879d00d7f00e7266d193605c59006bc5 100644
--- a/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.cxx
+++ b/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.cxx
@@ -7,91 +7,72 @@
#include "CbmRichRingTrackAssignIdeal.h"
-#include "TClonesArray.h"
-
-#include "CbmRichRing.h"
-
#include "CbmGlobalTrack.h"
#include "CbmMCTrack.h"
+#include "CbmRichRing.h"
#include "CbmTrackMatchNew.h"
+
#include "FairRootManager.h"
#include "FairTrackParam.h"
+#include <Logger.h>
-#include <iostream>
+#include "TClonesArray.h"
-using std::cout;
-using std::endl;
+#include <iostream>
-CbmRichRingTrackAssignIdeal::CbmRichRingTrackAssignIdeal()
- : fMcTracks(NULL)
- , fGlobalTracks(NULL)
- , fRingMatches(NULL)
- , fStsTrackMatches(NULL) {}
+CbmRichRingTrackAssignIdeal::CbmRichRingTrackAssignIdeal() {}
CbmRichRingTrackAssignIdeal::~CbmRichRingTrackAssignIdeal() {}
-void CbmRichRingTrackAssignIdeal::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichRingTrackAssignIdeal::Init", "RootManager not instantised!");
- }
+void CbmRichRingTrackAssignIdeal::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (nullptr == manager) LOG(fatal) << "CbmRichRingTrackAssignIdeal::Init(): FairRootManager is nullptr.";
- fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
- if (NULL == fMcTracks) {
- Fatal("CbmRichRingTrackAssignIdeal::Init", "No MCTrack array!");
- }
+ fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
+ if (fGlobalTracks == nullptr) LOG(fatal) << "CbmRichRingTrackAssignIdeal::Init():No GlobalTrack.";
- fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
- if (NULL == fGlobalTracks) {
- Fatal("CbmRichRingTrackAssignIdeal::Init", "No GlobalTrack array!");
- }
+ fRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
+ if (fRingMatches == nullptr) LOG(fatal) << "CbmRichRingTrackAssignIdeal::Init():No RichRingMatch.";
- fRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
- if (NULL == fRingMatches) {
- Fatal("CbmRichRingTrackAssignIdeal::Init", "No RichRingMatch array!");
- }
+ fStsTrackMatches = (TClonesArray*) manager->GetObject("StsTrackMatch");
+ if (fStsTrackMatches == nullptr) LOG(fatal) << "CbmRichRingTrackAssignIdeal::Init():No StsTrackMatch.";
+}
- fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
- if (NULL == fStsTrackMatches) {
- Fatal("CbmRichRingTrackAssignIdeal::Init", "No StsTrackMatch array!");
+void CbmRichRingTrackAssignIdeal::DoAssign(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
+{
+
+ if (event != nullptr) {
+ LOG(fatal) << "CbmRichRingTrackAssignIdeal::DoAssign(): CbmEvent is not nullptr. "
+ "This class does not support time-based mode. Please switch to event-by-event mode.";
}
-}
-void CbmRichRingTrackAssignIdeal::DoAssign(TClonesArray* rings,
- TClonesArray* richProj) {
Int_t nofTracks = richProj->GetEntriesFast();
Int_t nofRings = rings->GetEntriesFast();
- for (Int_t iRing = 0; iRing < nofRings; iRing++) {
- CbmRichRing* pRing = (CbmRichRing*) rings->At(iRing);
- if (NULL == pRing) continue;
- if (pRing->GetNofHits() < fMinNofHitsInRing) continue;
+ for (Int_t iR = 0; iR < nofRings; iR++) {
+ CbmRichRing* ring = static_cast<CbmRichRing*>(rings->At(iR));
+ if (ring == nullptr) continue;
+ if (ring->GetNofHits() < fMinNofHitsInRing) continue;
- CbmTrackMatchNew* pRingMatch = (CbmTrackMatchNew*) fRingMatches->At(iRing);
- if (NULL == pRingMatch) continue;
- Int_t ringID = pRingMatch->GetMatchedLink().GetIndex();
- // Double_t xRing = pRing->GetCenterX();
- // Double_t yRing = pRing->GetCenterY();
+ CbmTrackMatchNew* ringMatch = static_cast<CbmTrackMatchNew*>(fRingMatches->At(iR));
+ if (ringMatch == nullptr) continue;
+ Int_t ringId = ringMatch->GetMatchedLink().GetIndex();
- for (Int_t iTrack = 0; iTrack < nofTracks; iTrack++) {
- FairTrackParam* pTrack = (FairTrackParam*) richProj->At(iTrack);
- if (NULL == pTrack) continue;
- Double_t xTrack = pTrack->GetX();
- Double_t yTrack = pTrack->GetY();
+ for (Int_t iT = 0; iT < nofTracks; iT++) {
+ FairTrackParam* proj = static_cast<FairTrackParam*>(richProj->At(iT));
+ if (proj == nullptr) continue;
// no projection to photodetector plane
- if (xTrack == 0 && yTrack == 0) continue;
+ if (proj->GetX() == 0 && proj->GetY() == 0) continue;
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iTrack);
- if (NULL == gTrack) continue;
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iT));
+ if (gTrack == nullptr) continue;
if (gTrack->GetStsTrackIndex() == -1) continue;
- CbmTrackMatchNew* pTrackMatch =
- (CbmTrackMatchNew*) fStsTrackMatches->At(gTrack->GetStsTrackIndex());
- if (NULL == pTrackMatch) continue;
-
- if (pTrackMatch->GetMatchedLink().GetIndex() == ringID) {
- gTrack->SetRichRingIndex(iRing);
- } // ideal assignement
- } // loop tracks
- } // loop rings
+ CbmTrackMatchNew* trackMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(gTrack->GetStsTrackIndex()));
+ if (trackMatch == nullptr) continue;
+
+ if (trackMatch->GetMatchedLink().GetIndex() == ringId) { gTrack->SetRichRingIndex(iR); }
+ }
+ }
}
diff --git a/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.h b/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.h
index db3a69e14b3d92edf51b6c0b78a33b4f201ba2bd..f6072aff2b1bcb98d40503e22c8b083c6d68b064 100644
--- a/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.h
+++ b/reco/detectors/rich/tracks/CbmRichRingTrackAssignIdeal.h
@@ -44,13 +44,12 @@ public:
/**
* \brief Inherited from CbmRichRingTrackAssignBase.
*/
- virtual void DoAssign(TClonesArray* rings, TClonesArray* richProj);
+ virtual void DoAssign(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj);
private:
- TClonesArray* fMcTracks;
- TClonesArray* fGlobalTracks;
- TClonesArray* fRingMatches;
- TClonesArray* fStsTrackMatches;
+ TClonesArray* fGlobalTracks = nullptr;
+ TClonesArray* fRingMatches = nullptr;
+ TClonesArray* fStsTrackMatches = nullptr;
/**
* \brief Copy constructor.
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationBase.h b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationBase.h
index cc5424f31f04e89bb38a9a09249778ad1e080cdf..c5459765998a44e542ea3e3f6d088f5346aad5bb 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationBase.h
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationBase.h
@@ -11,6 +11,7 @@
#define CBM_RICH_TRACK_EXTRAPOLATION_BASE
class TClonesArray;
+class CbmEvent;
class CbmRichTrackExtrapolationBase {
public:
@@ -36,8 +37,7 @@ public:
* \param[out] extrapolatedTrackParams Output array of track parameters.
* \param[in] z Z coordinate to which track will be extrapolated.
*/
- virtual void DoExtrapolation(TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
+ virtual void DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks, TClonesArray* extrapolatedTrackParams,
double z) = 0;
private:
@@ -49,8 +49,7 @@ private:
/**
* \brief Assignment operator.
*/
- CbmRichTrackExtrapolationBase&
- operator=(const CbmRichTrackExtrapolationBase&);
+ CbmRichTrackExtrapolationBase& operator=(const CbmRichTrackExtrapolationBase&);
};
#endif
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.cxx b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.cxx
index acda6192030b006e522299ee1181a9ef0a01dc79..a4b5b5369de6e6d25f3e3b240954ea2f9df4be29 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.cxx
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.cxx
@@ -7,13 +7,16 @@
#include "CbmRichTrackExtrapolationIdeal.h"
+#include "CbmEvent.h"
#include "CbmGlobalTrack.h"
#include "CbmMCTrack.h"
#include "CbmRichPoint.h"
#include "CbmStsTrack.h"
#include "CbmTrackMatchNew.h"
+
#include "FairRootManager.h"
#include "FairTrackParam.h"
+#include <Logger.h>
#include "TClonesArray.h"
#include "TMatrixFSym.h"
@@ -23,95 +26,76 @@
using std::cout;
using std::endl;
-CbmRichTrackExtrapolationIdeal::CbmRichTrackExtrapolationIdeal()
- : fRefPlanePoints(NULL)
- , fMcTracks(NULL)
- , fStsTracks(NULL)
- , fStsTrackMatches(NULL) {}
+CbmRichTrackExtrapolationIdeal::CbmRichTrackExtrapolationIdeal() {}
CbmRichTrackExtrapolationIdeal::~CbmRichTrackExtrapolationIdeal() {}
-void CbmRichTrackExtrapolationIdeal::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichTrackExtrapolationIdeal::Init",
- "RootManager not instantised!");
- }
+void CbmRichTrackExtrapolationIdeal::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (manager = nullptr) LOG(fatal) << "CbmRichTrackExtrapolationIdeal::Init(): FairRootManager is nullptr.";
- fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
- if (NULL == fMcTracks) {
- Fatal("CbmRichTrackExtrapolationIdeal::Init", "No MCTrack array!");
- }
+ fRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
+ if (fRefPlanePoints = nullptr) LOG(fatal) << "CbmRichTrackExtrapolationIdeal::Init(): No RefPlanePoint array.";
- fRefPlanePoints = (TClonesArray*) ioman->GetObject("RefPlanePoint");
- if (NULL == fRefPlanePoints) {
- Fatal("CbmRichTrackExtrapolationIdeal::Init", "No RefPlanePoint array!");
- }
+ fStsTracks = (TClonesArray*) manager->GetObject("StsTrack");
+ if (fStsTracks = nullptr) LOG(fatal) << "CbmRichTrackExtrapolationIdeal::Init(): No StsTrack array.";
- fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
- if (NULL == fStsTracks) {
- Fatal("CbmRichTrackExtrapolationIdeal::Init", "No StsTrack array!");
- }
+ fStsTrackMatches = (TClonesArray*) manager->GetObject("StsTrackMatch");
+ if (fStsTrackMatches = nullptr) LOG(fatal) << "CbmRichTrackExtrapolationIdeal::Init(): No StsTrackMatch array.";
+}
+
+void CbmRichTrackExtrapolationIdeal::DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks,
+ TClonesArray* extrapolatedTrackParams, double /*z*/)
+{
- fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
- if (NULL == fStsTrackMatches) {
- Fatal("CbmRichTrackExtrapolationIdeal::Init", "No StsTrackMatch array!");
+ if (event != nullptr) {
+ LOG(fatal) << "CbmRichTrackExtrapolationIdeal::DoExtrapolation(): CbmEvent is not nullptr. "
+ "This class does not support time-based mode. Please switch to event-by-event mode.";
}
-}
-void CbmRichTrackExtrapolationIdeal::DoExtrapolation(
- TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
- double /*z*/) {
- if (NULL == extrapolatedTrackParams) {
- cout << "-E- CbmRichTrackExtrapolationIdeal::DoExtrapolate: TrackParam "
- "Array missing!"
- << endl;
+ if (extrapolatedTrackParams == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationIdeal::DoExtrapolation(): extrapolatedTrackParams missing!";
return;
}
- if (NULL == globalTracks) {
- cout << "-E- CbmRichTrackExtrapolationIdeal::DoExtrapolate: Global Track "
- "Array missing! "
- << endl;
+ if (globalTracks == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationIdeal::DoExtrapolation(): globalTracks missing!";
return;
}
- // some default variables
- Double_t tx, ty, qp;
Double_t charge = 1.;
TMatrixFSym covMat(5);
for (Int_t i = 0; i < 5; i++)
for (Int_t j = 0; j <= i; j++)
covMat(i, j) = 0;
- covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
- 1.e-4;
+ covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) = 1.e-4;
TVector3 pos, mom;
- Int_t nTracks = globalTracks->GetEntriesFast();
- for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) globalTracks->At(iTrack);
- new ((*extrapolatedTrackParams)[iTrack])
- FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
- Int_t idSTS = gTrack->GetStsTrackIndex();
- if (idSTS < 0) continue;
- CbmStsTrack* pSTStr = (CbmStsTrack*) fStsTracks->At(idSTS);
- if (NULL == pSTStr) continue;
- CbmTrackMatchNew* pTrackMatch =
- (CbmTrackMatchNew*) fStsTrackMatches->At(idSTS);
- if (NULL == pTrackMatch) continue;
- Int_t iMCmatch = pTrackMatch->GetMatchedLink().GetIndex();
- for (Int_t ii = 0; ii < fRefPlanePoints->GetEntriesFast(); ii++) {
- CbmRichPoint* pRefPlane = (CbmRichPoint*) fRefPlanePoints->At(ii);
- if (pRefPlane->GetTrackID() == iMCmatch) {
- pRefPlane->Momentum(mom);
- pRefPlane->Position(pos);
- tx = mom.Px() / mom.Pz();
- ty = mom.Py() / mom.Pz();
- qp = charge / mom.Mag();
+ Int_t nofGlobalTracks = globalTracks->GetEntriesFast();
+ for (Int_t iT = 0; iT < nofGlobalTracks; iT++) {
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(globalTracks->At(iT));
+ new ((*extrapolatedTrackParams)[iT]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+
+ Int_t stsInd = gTrack->GetStsTrackIndex();
+ if (stsInd < 0) continue;
+ CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
+ if (stsTrack == nullptr) continue;
+ CbmTrackMatchNew* stsTrackMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
+ if (stsTrackMatch == nullptr) continue;
+ Int_t stsMcInd = stsTrackMatch->GetMatchedLink().GetIndex();
+
+ for (Int_t iR = 0; iR < fRefPlanePoints->GetEntriesFast(); iR++) {
+ CbmRichPoint* refPlanePoint = static_cast<CbmRichPoint*>(fRefPlanePoints->At(iR));
+ if (refPlanePoint->GetTrackID() == stsMcInd) {
+ refPlanePoint->Momentum(mom);
+ refPlanePoint->Position(pos);
+ Double_t tx = mom.Px() / mom.Pz();
+ Double_t ty = mom.Py() / mom.Pz();
+ Double_t qp = charge / mom.Mag();
FairTrackParam richtrack(pos.X(), pos.Y(), pos.Z(), tx, ty, qp, covMat);
- *(FairTrackParam*) (extrapolatedTrackParams->At(iTrack)) = richtrack;
- } // select MCPoint to that track
- } // loop MCPoints
- } // loop tracks
+ *(FairTrackParam*) (extrapolatedTrackParams->At(iT)) = richtrack;
+ }
+ }
+ }
}
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.h b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.h
index 7081b7c4a81d8a822ab207f0316ec6ee01c41b53..abc3ef72e73725eac17cd26aca0750d8234f848c 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.h
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationIdeal.h
@@ -46,15 +46,13 @@ public:
/**
* \brief Inherited from CbmRichTrackExtrapolationBase.
*/
- virtual void DoExtrapolation(TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
+ virtual void DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks, TClonesArray* extrapolatedTrackParams,
double z);
private:
- TClonesArray* fRefPlanePoints;
- TClonesArray* fMcTracks;
- TClonesArray* fStsTracks;
- TClonesArray* fStsTrackMatches;
+ TClonesArray* fRefPlanePoints = nullptr;
+ TClonesArray* fStsTracks = nullptr;
+ TClonesArray* fStsTrackMatches = nullptr;
/**
* \brief Copy constructor.
@@ -64,8 +62,7 @@ private:
/**
* \brief Assignment operator.
*/
- CbmRichTrackExtrapolationIdeal&
- operator=(const CbmRichTrackExtrapolationIdeal&);
+ CbmRichTrackExtrapolationIdeal& operator=(const CbmRichTrackExtrapolationIdeal&);
};
#endif
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.cxx b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.cxx
index 08ba422493c6ec49c5473bc0e5856fed99d10caf..d89def7a4463922ff6d6bf8b062379a159287586 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.cxx
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.cxx
@@ -6,52 +6,47 @@
**/
#include "CbmRichTrackExtrapolationKF.h"
-#include "CbmStsKFTrackFitter.h"
-
+#include "CbmEvent.h"
#include "CbmGlobalTrack.h"
+#include "CbmStsKFTrackFitter.h"
#include "CbmStsTrack.h"
+
#include "FairRootManager.h"
#include "FairTrackParam.h"
+#include <Logger.h>
#include "TClonesArray.h"
#include "TMatrixFSym.h"
#include <iostream>
+
using std::cout;
using std::endl;
-CbmRichTrackExtrapolationKF::CbmRichTrackExtrapolationKF() : fStsTracks(0) {}
+CbmRichTrackExtrapolationKF::CbmRichTrackExtrapolationKF() {}
CbmRichTrackExtrapolationKF::~CbmRichTrackExtrapolationKF() {}
-void CbmRichTrackExtrapolationKF::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichTrackExtrapolationKF::Init", "RootManager not instantised!");
- }
+void CbmRichTrackExtrapolationKF::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (manager == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationKF::Init(): FairRootManager is nullptr.";
- fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
- if (NULL == fStsTracks) {
- Fatal("CbmRichTrackExtrapolationKF::Init", "No StsTrack array!");
- }
+ fStsTracks = (TClonesArray*) manager->GetObject("StsTrack");
+ if (fStsTracks == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationKF::Init(): No StsTrack array.";
}
-void CbmRichTrackExtrapolationKF::DoExtrapolation(
- TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
- double z) {
- if (NULL == extrapolatedTrackParams) {
- cout << "-E- CbmRichTrackExtrapolationKF::DoExtrapolate: TrackParamArray "
- "missing!"
- << endl;
+void CbmRichTrackExtrapolationKF::DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks,
+ TClonesArray* extrapolatedTrackParams, double z)
+{
+ if (extrapolatedTrackParams == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationKF::DoExtrapolation(): extrapolatedTrackParams is nullptr.";
return;
}
- if (NULL == globalTracks) {
- cout
- << "-E- CbmRichTrackExtrapolationKF::DoExtrapolate: Track Array missing!"
- << endl;
+ if (globalTracks == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationKF::DoExtrapolation(): globalTracks is nullptr.";
return;
}
@@ -59,26 +54,25 @@ void CbmRichTrackExtrapolationKF::DoExtrapolation(
for (Int_t i = 0; i < 5; i++)
for (Int_t j = 0; j <= i; j++)
covMat(i, j) = 0;
- covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
- 1.e-4;
+ covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) = 1.e-4;
TVector3 pos, mom;
-
- Int_t nTracks = globalTracks->GetEntriesFast();
- cout << "bmRichTrackExtrapolationKF nofGlobalTracks:" << nTracks << endl;
- for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) globalTracks->At(iTrack);
- new ((*extrapolatedTrackParams)[iTrack])
- FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
- Int_t idSTS = gTrack->GetStsTrackIndex();
- if (idSTS < 0) continue;
- CbmStsTrack* pSTStr = (CbmStsTrack*) fStsTracks->At(idSTS);
- if (NULL == pSTStr) continue;
+ Int_t nofGlobalTracks = event ? event->GetNofData(ECbmDataType::kGlobalTrack) : globalTracks->GetEntriesFast();
+ for (Int_t iT0 = 0; iT0 < nofGlobalTracks; iT0++) {
+ Int_t iT = event ? event->GetIndex(ECbmDataType::kGlobalTrack, iT0) : iT0;
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(globalTracks->At(iT));
+ new ((*extrapolatedTrackParams)[iT]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+ if (event != nullptr) event->AddData(ECbmDataType::kRichTrackParamZ, iT);
+
+ Int_t stsInd = gTrack->GetStsTrackIndex();
+ if (stsInd < 0) continue;
+ CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
+ if (stsTrack == nullptr) continue;
CbmStsKFTrackFitter KF;
FairTrackParam ExTrack;
- KF.Extrapolate(pSTStr, z, &ExTrack);
+ KF.Extrapolate(stsTrack, z, &ExTrack);
- *(FairTrackParam*) (extrapolatedTrackParams->At(iTrack)) = ExTrack;
+ *(FairTrackParam*) (extrapolatedTrackParams->At(iT)) = ExTrack;
}
}
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.h b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.h
index 856be739740f53ee0aef005d77cd162494dad4bf..9357566bf8aa975a9b2187601e917cef22ee3399 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.h
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationKF.h
@@ -46,12 +46,11 @@ public:
/**
* \brief Inherited from CbmRichTrackExtrapolationBase.
*/
- virtual void DoExtrapolation(TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
+ virtual void DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks, TClonesArray* extrapolatedTrackParams,
double z);
private:
- TClonesArray* fStsTracks;
+ TClonesArray* fStsTracks = nullptr;
private:
/**
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.cxx b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.cxx
index ed44d120c79465d085f557f379e53abd27764f5e..3d95493db2f19cbf95d87d3510d0c73ad4512da0 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.cxx
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.cxx
@@ -6,14 +6,16 @@
**/
#include "CbmRichTrackExtrapolationLittrack.h"
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmStsTrack.h"
#include "cbm/base/CbmLitToolFactory.h"
#include "cbm/utils/CbmLitConverterFairTrackParam.h"
#include "propagation/CbmLitTGeoTrackPropagator.h"
-#include "CbmGlobalTrack.h"
-#include "CbmStsTrack.h"
#include "FairRootManager.h"
#include "FairTrackParam.h"
+#include <Logger.h>
#include "TClonesArray.h"
@@ -22,66 +24,56 @@
using std::cout;
using std::endl;
-CbmRichTrackExtrapolationLittrack::CbmRichTrackExtrapolationLittrack()
- : CbmRichTrackExtrapolationBase(), fStsTracks(0), fLitPropagator() {}
+CbmRichTrackExtrapolationLittrack::CbmRichTrackExtrapolationLittrack() {}
CbmRichTrackExtrapolationLittrack::~CbmRichTrackExtrapolationLittrack() {}
-void CbmRichTrackExtrapolationLittrack::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichTrackExtrapolationLittrack::Init",
- "RootManager not instantised!");
- }
+void CbmRichTrackExtrapolationLittrack::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (manager == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationLittrack::Init(): FairRootManager is nullptr.";
- fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
- if (NULL == fStsTracks) {
- Fatal("CbmRichTrackExtrapolationLittrack::Init", "No StsTrack array!");
- }
+ fStsTracks = (TClonesArray*) manager->GetObject("StsTrack");
+ if (fStsTracks == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationLittrack::Init(): No StsTrack array.";
fLitPropagator = CbmLitToolFactory::CreateTrackPropagator("lit");
}
-void CbmRichTrackExtrapolationLittrack::DoExtrapolation(
- TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
- double z) {
- cout << "CbmRichTrackExtrapolationLittrack::DoExtrapolation" << endl;
- if (NULL == extrapolatedTrackParams) {
- cout << "-E- CbmRichTrackExtrapolationLittrack::DoExtrapolate: "
- "TrackParamArray missing!"
- << endl;
+void CbmRichTrackExtrapolationLittrack::DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks,
+ TClonesArray* extrapolatedTrackParams, double z)
+{
+ if (extrapolatedTrackParams == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationLittrack::DoExtrapolation(): extrapolatedTrackParams is nullptr.";
return;
}
- if (NULL == globalTracks) {
- cout << "-E- CbmRichTrackExtrapolationLittrack::DoExtrapolate: Track Array "
- "missing!"
- << endl;
+ if (globalTracks == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationLittrack::DoExtrapolation(): globalTracks is nullptr.";
return;
}
- Int_t nTracks = globalTracks->GetEntriesFast();
- for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) globalTracks->At(iTrack);
- new ((*extrapolatedTrackParams)[iTrack]) FairTrackParam();
- Int_t idSTS = gTrack->GetStsTrackIndex();
- if (idSTS < 0) continue;
- CbmStsTrack* pSTStr = (CbmStsTrack*) fStsTracks->At(idSTS);
- if (NULL == pSTStr) continue;
+ Int_t nofGlobalTracks = event ? event->GetNofData(ECbmDataType::kGlobalTrack) : globalTracks->GetEntriesFast();
+ for (Int_t iT0 = 0; iT0 < nofGlobalTracks; iT0++) {
+ Int_t iT = event ? event->GetIndex(ECbmDataType::kGlobalTrack, iT0) : iT0;
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(globalTracks->At(iT));
+ new ((*extrapolatedTrackParams)[iT]) FairTrackParam();
+ if (event != nullptr) event->AddData(ECbmDataType::kRichTrackParamZ, iT);
+
+ Int_t stsInd = gTrack->GetStsTrackIndex();
+ if (stsInd < 0) continue;
+ CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
+ if (stsTrack == nullptr) continue;
CbmLitTrackParam litInParam, litOutParam;
- CbmLitConverterFairTrackParam::FairTrackParamToCbmLitTrackParam(
- pSTStr->GetParamLast(), &litInParam);
+ CbmLitConverterFairTrackParam::FairTrackParamToCbmLitTrackParam(stsTrack->GetParamLast(), &litInParam);
std::vector<litfloat> F(25);
litfloat length = 0;
fLitPropagator->Propagate(&litInParam, &litOutParam, z, 11, &F, &length);
FairTrackParam outParam;
- CbmLitConverterFairTrackParam::CbmLitTrackParamToFairTrackParam(
- &litOutParam, &outParam);
+ CbmLitConverterFairTrackParam::CbmLitTrackParamToFairTrackParam(&litOutParam, &outParam);
- *(FairTrackParam*) (extrapolatedTrackParams->At(iTrack)) = outParam;
+ *(FairTrackParam*) (extrapolatedTrackParams->At(iT)) = outParam;
}
}
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.h b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.h
index 9efc3ee34b451307d47264ba702cd9e2e820c632..96d66ffda852120ea31aa24300f6c9c0372e0aae 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.h
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationLittrack.h
@@ -50,13 +50,12 @@ public:
/**
* \brief Inherited from CbmRichTrackExtrapolationBase.
*/
- virtual void DoExtrapolation(TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
+ virtual void DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks, TClonesArray* extrapolatedTrackParams,
double z);
private:
- TClonesArray* fStsTracks;
- TrackPropagatorPtr fLitPropagator;
+ TClonesArray* fStsTracks = nullptr;
+ TrackPropagatorPtr fLitPropagator = nullptr;
private:
/**
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.cxx b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.cxx
index 391b9ed77f23452a19d4b8204629661ce6434f20..a996e394e1d10c37ef3ac4e67e757170dd7c2ece 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.cxx
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.cxx
@@ -7,14 +7,16 @@
#include "CbmRichTrackExtrapolationMirrorIdeal.h"
-#include "CbmRichPoint.h"
-
+#include "CbmEvent.h"
#include "CbmGlobalTrack.h"
#include "CbmMCTrack.h"
+#include "CbmRichPoint.h"
#include "CbmStsTrack.h"
#include "CbmTrackMatchNew.h"
+
#include "FairRootManager.h"
#include "FairTrackParam.h"
+#include <Logger.h>
#include "TClonesArray.h"
#include "TMatrixFSym.h"
@@ -24,96 +26,79 @@
using std::cout;
using std::endl;
-CbmRichTrackExtrapolationMirrorIdeal::CbmRichTrackExtrapolationMirrorIdeal()
- : fRichMirrorPoints(NULL)
- , fMcTracks(NULL)
- , fSTSArray(NULL)
- , fTrackMatchArray(NULL) {}
+CbmRichTrackExtrapolationMirrorIdeal::CbmRichTrackExtrapolationMirrorIdeal() {}
CbmRichTrackExtrapolationMirrorIdeal::~CbmRichTrackExtrapolationMirrorIdeal() {}
-void CbmRichTrackExtrapolationMirrorIdeal::Init() {
- FairRootManager* ioman = FairRootManager::Instance();
- if (NULL == ioman) {
- Fatal("CbmRichTrackExtrapolationMirrorIdeal::Init",
- "RootManager not instantised!");
- }
+void CbmRichTrackExtrapolationMirrorIdeal::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (manager == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationMirrorIdeal::Init(): FairRootManager is nullptr.";
- fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
- if (NULL == fMcTracks) {
- Fatal("CbmRichTrackExtrapolationMirrorIdeal::Init", "No MCTrack array!");
- }
+ fMcTracks = (TClonesArray*) manager->GetObject("MCTrack");
+ if (fMcTracks == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationMirrorIdeal::Init(): No MCTrack array.";
- fRichMirrorPoints = (TClonesArray*) ioman->GetObject("RichMirrorPoint");
- if (NULL == fRichMirrorPoints) {
- Fatal("CbmRichTrackExtrapolationMirrorIdeal::Init",
- "No RichMirrorPoint array!");
- }
+ fRichMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
+ if (fRichMirrorPoints == nullptr)
+ LOG(fatal) << "CbmRichTrackExtrapolationMirrorIdeal::Init(): No RichMirrorPoint array.";
- fSTSArray = (TClonesArray*) ioman->GetObject("StsTrack");
- if (NULL == fSTSArray) {
- Fatal("CbmRichTrackExtrapolationMirrorIdeal::Init", "No StsTrack array!");
- }
+ fStsTracks = (TClonesArray*) manager->GetObject("StsTrack");
+ if (fStsTracks == nullptr) LOG(fatal) << "CbmRichTrackExtrapolationMirrorIdeal::Init(): No StsTrack array.";
- fTrackMatchArray = (TClonesArray*) ioman->GetObject("StsTrackMatch");
- if (NULL == fTrackMatchArray) {
- Fatal("CbmRichTrackExtrapolationMirrorIdeal::Init",
- "No StsTrackMatch array!");
- }
+ fStsTrackMatches = (TClonesArray*) manager->GetObject("StsTrackMatch");
+ if (fStsTrackMatches == nullptr)
+ LOG(fatal) << "CbmRichTrackExtrapolationMirrorIdeal::Init(): No StsTrackMatch array.";
}
-void CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolation(
- TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
- double /*z*/) {
- if (NULL == extrapolatedTrackParams) {
- cout << "-E- CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolate: "
- "TrackParam Array missing!"
- << endl;
+void CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks,
+ TClonesArray* extrapolatedTrackParams, double /*z*/)
+{
+
+ if (event != nullptr) {
+ LOG(fatal) << "CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolation(): CbmEvent is not nullptr. "
+ "This class does not support time-based mode. Please switch to event-by-event mode.";
+ }
+
+ if (extrapolatedTrackParams == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolation(): extrapolatedTrackParams missing!";
return;
}
- if (NULL == globalTracks) {
- cout << "-E- CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolate: Global "
- "Track Array missing!"
- << endl;
+ if (globalTracks == nullptr) {
+ LOG(error) << "CbmRichTrackExtrapolationMirrorIdeal::DoExtrapolation(): globalTracks missing!";
return;
}
- Double_t tx, ty, qp;
Double_t charge = 1.;
TMatrixFSym covMat(5);
for (Int_t i = 0; i < 5; i++)
for (Int_t j = 0; j <= i; j++)
covMat(i, j) = 0;
- covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) =
- 1.e-4;
+ covMat(0, 0) = covMat(1, 1) = covMat(2, 2) = covMat(3, 3) = covMat(4, 4) = 1.e-4;
TVector3 pos, mom;
- Int_t nTracks = globalTracks->GetEntriesFast();
- for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- CbmGlobalTrack* gTrack = (CbmGlobalTrack*) globalTracks->At(iTrack);
- new ((*extrapolatedTrackParams)[iTrack])
- FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
- Int_t idSTS = gTrack->GetStsTrackIndex();
-
- if (idSTS < 0) continue;
- CbmStsTrack* pSTStr = (CbmStsTrack*) fSTSArray->At(idSTS);
- if (NULL == pSTStr) continue;
- CbmTrackMatchNew* pTrackMatch =
- (CbmTrackMatchNew*) fTrackMatchArray->At(idSTS);
- if (NULL == pTrackMatch) continue;
- Int_t iMCmatch = pTrackMatch->GetMatchedLink().GetIndex();
- for (Int_t ii = 0; ii < fRichMirrorPoints->GetEntriesFast(); ii++) {
- CbmRichPoint* pMirror = (CbmRichPoint*) fRichMirrorPoints->At(ii);
- if (pMirror->GetTrackID() == iMCmatch) {
+ Int_t nofGlobalTracks = globalTracks->GetEntriesFast();
+ for (Int_t iT = 0; iT < nofGlobalTracks; iT++) {
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(globalTracks->At(iT));
+ new ((*extrapolatedTrackParams)[iT]) FairTrackParam(0., 0., 0., 0., 0., 0., covMat);
+
+ Int_t stsInd = gTrack->GetStsTrackIndex();
+ if (stsInd < 0) continue;
+ CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
+ if (stsTrack == nullptr) continue;
+ CbmTrackMatchNew* stsTrackMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
+ if (stsTrackMatch == nullptr) continue;
+ Int_t stsMcInd = stsTrackMatch->GetMatchedLink().GetIndex();
+ for (Int_t iM = 0; iM < fRichMirrorPoints->GetEntriesFast(); iM++) {
+ CbmRichPoint* pMirror = static_cast<CbmRichPoint*>(fRichMirrorPoints->At(iM));
+ if (pMirror->GetTrackID() == stsMcInd) {
pMirror->Momentum(mom);
pMirror->Position(pos);
- tx = mom.Px() / mom.Pz();
- ty = mom.Py() / mom.Pz();
- qp = charge / mom.Mag();
+ Double_t tx = mom.Px() / mom.Pz();
+ Double_t ty = mom.Py() / mom.Pz();
+ Double_t qp = charge / mom.Mag();
FairTrackParam richtrack(pos.X(), pos.Y(), pos.Z(), tx, ty, qp, covMat);
- *(FairTrackParam*) (extrapolatedTrackParams->At(iTrack)) = richtrack;
+ *(FairTrackParam*) (extrapolatedTrackParams->At(iT)) = richtrack;
}
}
}
diff --git a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.h b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.h
index 9a55b392375c1f6451825e0fd16b8c16acb99f36..bd5bb2e4612c6fee7d7514f6af9fc372381c4cfc 100644
--- a/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.h
+++ b/reco/detectors/rich/tracks/CbmRichTrackExtrapolationMirrorIdeal.h
@@ -32,8 +32,7 @@ class TClonesArray;
* \author Claudia Hoehne
* \date 2006
**/
-class CbmRichTrackExtrapolationMirrorIdeal :
- public CbmRichTrackExtrapolationBase {
+class CbmRichTrackExtrapolationMirrorIdeal : public CbmRichTrackExtrapolationBase {
public:
/**
* \brief Default constructor.
@@ -53,27 +52,24 @@ public:
/**
* \brief Inherited from CbmRichTrackExtrapolationBase.
*/
- virtual void DoExtrapolation(TClonesArray* globalTracks,
- TClonesArray* extrapolatedTrackParams,
+ virtual void DoExtrapolation(CbmEvent* event, TClonesArray* globalTracks, TClonesArray* extrapolatedTrackParams,
double z);
private:
- TClonesArray* fRichMirrorPoints;
- TClonesArray* fMcTracks;
- TClonesArray* fSTSArray;
- TClonesArray* fTrackMatchArray;
+ TClonesArray* fRichMirrorPoints = nullptr;
+ TClonesArray* fMcTracks = nullptr;
+ TClonesArray* fStsTracks = nullptr;
+ TClonesArray* fStsTrackMatches = nullptr;
/**
* \brief Copy constructor.
*/
- CbmRichTrackExtrapolationMirrorIdeal(
- const CbmRichTrackExtrapolationMirrorIdeal&);
+ CbmRichTrackExtrapolationMirrorIdeal(const CbmRichTrackExtrapolationMirrorIdeal&);
/**
* \brief Assignment operator.
*/
- CbmRichTrackExtrapolationMirrorIdeal&
- operator=(const CbmRichTrackExtrapolationMirrorIdeal&);
+ CbmRichTrackExtrapolationMirrorIdeal& operator=(const CbmRichTrackExtrapolationMirrorIdeal&);
};
#endif
diff --git a/reco/detectors/rich/utils/CbmRichNavigationUtil.h b/reco/detectors/rich/utils/CbmRichNavigationUtil.h
index 60dc0ea9fd0652e3131a0382b0aa8fa8d1187c75..afe3b1ce9a6d007ed395501e8777957fb8082299 100644
--- a/reco/detectors/rich/utils/CbmRichNavigationUtil.h
+++ b/reco/detectors/rich/utils/CbmRichNavigationUtil.h
@@ -5,6 +5,7 @@
#include "FairTrackParam.h"
#include <Logger.h>
+#include "TGeoManager.h"
#include "TObject.h"
#include "TVector3.h"
@@ -12,9 +13,8 @@
class CbmRichNavigationUtil {
public:
- static string FindIntersection(const FairTrackParam* par,
- TVector3& crossPoint,
- const string& volumeName) {
+ static string FindIntersection(const FairTrackParam* par, TVector3& crossPoint, const string& volumeName)
+ {
TVector3 dirCos, pos;
pos.SetXYZ(par->GetX(), par->GetY(), par->GetZ());
Double_t nx, ny, nz;
@@ -24,13 +24,11 @@ public:
return FindIntersection(dirCos, pos, crossPoint, volumeName);
}
- static string FindIntersection(const TVector3& dirCos,
- const TVector3& pos,
- TVector3& crossPoint,
- const string& volumeName) {
+ static string FindIntersection(const TVector3& dirCos, const TVector3& pos, TVector3& crossPoint,
+ const string& volumeName)
+ {
// if (volumeName == "pmt_pixel")cout << "InitTrack: " << pos.X() << " " << pos.Y() << " " << pos.Z() << " " << dirCos.X() << " " << dirCos.Y()<< " " << dirCos.Z() << endl;
- gGeoManager->InitTrack(
- pos.X(), pos.Y(), pos.Z(), dirCos.X(), dirCos.Y(), dirCos.Z());
+ gGeoManager->InitTrack(pos.X(), pos.Y(), pos.Z(), dirCos.X(), dirCos.Y(), dirCos.Z());
if (gGeoManager->IsOutside()) { return string(""); }
@@ -60,8 +58,7 @@ public:
return string("");
}
// Check for NaN values
- if (std::isnan(gGeoManager->GetCurrentPoint()[0])
- || std::isnan(gGeoManager->GetCurrentPoint()[1])
+ if (std::isnan(gGeoManager->GetCurrentPoint()[0]) || std::isnan(gGeoManager->GetCurrentPoint()[1])
|| std::isnan(gGeoManager->GetCurrentPoint()[2])) {
// if (volumeName == "pmt_pixel")std::cout << "Error! CbmRichNavigationUtil::FindIntersections: NaN values.\n";
gGeoManager->PopDummy();
@@ -75,14 +72,10 @@ public:
return string("");
}
- static void GetDirCos(const FairTrackParam* par,
- Double_t& nx,
- Double_t& ny,
- Double_t& nz) {
- Double_t p =
- (std::abs(par->GetQp()) != 0.) ? 1. / std::abs(par->GetQp()) : 1.e20;
- Double_t pz = std::sqrt(
- p * p / (par->GetTx() * par->GetTx() + par->GetTy() * par->GetTy() + 1));
+ static void GetDirCos(const FairTrackParam* par, Double_t& nx, Double_t& ny, Double_t& nz)
+ {
+ Double_t p = (std::abs(par->GetQp()) != 0.) ? 1. / std::abs(par->GetQp()) : 1.e20;
+ Double_t pz = std::sqrt(p * p / (par->GetTx() * par->GetTx() + par->GetTy() * par->GetTy() + 1));
Double_t px = par->GetTx() * pz;
Double_t py = par->GetTy() * pz;
TVector3 unit = TVector3(px, py, pz).Unit();