diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C index 5fb13de91be96a717f5f458905d685d236b2343e..2035cac9ea773f3a72bc7a11b55b283afeb276ba 100644 --- a/macro/run/run_reco.C +++ b/macro/run/run_reco.C @@ -5,6 +5,7 @@ // --- Includes needed for IDE +#include <RtypesCore.h> #if !defined(__CLING__) #include "CbmBuildEventsFromTracksReal.h" #include "CbmBuildEventsIdeal.h" @@ -349,6 +350,11 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice = 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- " << myName << ": Added task " << trdCluster->GetName() << std::endl; diff --git a/reco/detectors/trd/CbmTrdClusterFinder.cxx b/reco/detectors/trd/CbmTrdClusterFinder.cxx index 5a25c00cfa5bdbd4486095977643ba77c0f18df6..40c5d2bb2e0d16918d2ba301b62e68e688fba1ad 100644 --- a/reco/detectors/trd/CbmTrdClusterFinder.cxx +++ b/reco/detectors/trd/CbmTrdClusterFinder.cxx @@ -1,6 +1,12 @@ #include "CbmTrdClusterFinder.h" +#include "CbmDefs.h" #include "CbmDigiManager.h" +#include "CbmTrdCluster.h" +#include "CbmTrdDigi.h" +#include "CbmTrdGeoHandler.h" +#include "CbmTrdModuleRecR.h" +#include "CbmTrdModuleRecT.h" #include "CbmTrdParAsic.h" #include "CbmTrdParModDigi.h" #include "CbmTrdParModGain.h" @@ -11,62 +17,34 @@ #include "CbmTrdParSetGas.h" #include "CbmTrdParSetGeo.h" -#include "CbmTrdCluster.h" -#include "CbmTrdDigi.h" -#include "CbmTrdGeoHandler.h" -#include "CbmTrdModuleRecR.h" -#include "CbmTrdModuleRecT.h" - #include <FairLogger.h> #include <FairRootManager.h> #include <FairRunAna.h> #include <FairRuntimeDb.h> +#include <RtypesCore.h> #include <TArray.h> -#include <TBits.h> #include <TClonesArray.h> #include <TGeoPhysicalNode.h> #include <TStopwatch.h> // #include "TCanvas.h" // #include "TImage.h" -#include <cmath> #include <iomanip> #include <iostream> -using std::cout; -using std::endl; -using std::fabs; -using std::setprecision; + +#include <cmath> + Int_t CbmTrdClusterFinder::fgConfig = 0; Float_t CbmTrdClusterFinder::fgMinimumChargeTH = .5e-06; //_____________________________________________________________________ -CbmTrdClusterFinder::CbmTrdClusterFinder() - : FairTask("TrdClusterFinder", 1) - , fClusters(NULL) - , fDigiMap() - , fModuleMap() - , fNeighbours() - , fModDigiMap() - , fDigiRow() - , fDigiCol() - , fDigiCharge() - , fClusterBuffer() - , fModClusterDigiMap() - , - //======================= - fModules() - , fAsicPar(NULL) - , fGasPar(NULL) - , fDigiPar(NULL) - , fGainPar(NULL) - , fGeoPar(NULL) - -{} +CbmTrdClusterFinder::CbmTrdClusterFinder() : FairTask("TrdClusterFinder", 1) { SetUseOnlyEventDigis(); } // -------------------------------------------------------------------- // ---- Destructor ---------------------------------------------------- -CbmTrdClusterFinder::~CbmTrdClusterFinder() { +CbmTrdClusterFinder::~CbmTrdClusterFinder() +{ if (fClusters) { fClusters->Clear("C"); @@ -80,87 +58,144 @@ CbmTrdClusterFinder::~CbmTrdClusterFinder() { } //_____________________________________________________________________ -Bool_t CbmTrdClusterFinder::AddCluster(CbmTrdCluster* c) { +Bool_t CbmTrdClusterFinder::AddCluster(CbmTrdCluster* c) +{ Int_t ncl(fClusters->GetEntriesFast()); new ((*fClusters)[ncl++]) CbmTrdCluster(*c); return kTRUE; } +// ---- addDigisToModules ---- +UInt_t CbmTrdClusterFinder::addDigisToModules() +{ + UInt_t ndigis = static_cast<UInt_t>(std::abs(CbmDigiManager::Instance()->GetNofDigis(ECbmModuleId::kTrd))); + for (size_t idigi = 0; idigi < ndigis; idigi++) { + addDigiToModule(idigi); + } + return ndigis; +} + +// ---- addDigisToModules ---- +UInt_t CbmTrdClusterFinder::addDigisToModules(CbmEvent* event) +{ + UInt_t ndigis = static_cast<UInt_t>(std::abs(event->GetNofData(ECbmDataType::kTrdDigi))); + for (size_t idigi = 0; idigi < ndigis; idigi++) { + auto digiindex = event->GetIndex(ECbmDataType::kTrdDigi, idigi); + addDigiToModule(digiindex); + } + return ndigis; +} + + +// ---- addDigiToModule ---- +void CbmTrdClusterFinder::addDigiToModule(UInt_t digiIdx) +{ + CbmTrdModuleRec* mod = nullptr; + + const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(digiIdx); + + Int_t moduleAddress = digi->GetAddressModule(); + + std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.find(moduleAddress); + if (imod == fModules.end()) mod = AddModule(digi); + else + mod = imod->second; + + mod->AddDigi(digi, digiIdx); +} + +// ---- processDigisInModules ---- +void CbmTrdClusterFinder::processDigisInModules(UInt_t ndigis, CbmEvent* event) +{ + CbmTrdModuleRec* mod(NULL); + Int_t digiCounter(0), clsCounter(0); + for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); imod != fModules.end(); imod++) { + mod = imod->second; + digiCounter += mod->GetOverThreshold(); + clsCounter += mod->FindClusters(); + AddClusters(mod->GetClusters(), event, kTRUE); + } + + // remove local data from all modules + for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); imod != fModules.end(); imod++) + imod->second->Clear("cls"); + + fNrDigis += ndigis; + fNrClusters += clsCounter; + + LOG(debug) << GetName() << "::Exec : Digis : " << ndigis << " / " << digiCounter << " above threshold (" + << 1e6 * fgMinimumChargeTH << " keV)"; + LOG(debug) << GetName() << "::Exec : Clusters : " << clsCounter; +} + //____________________________________________________________________________________ -CbmTrdModuleRec* CbmTrdClusterFinder::AddModule(const CbmTrdDigi* digi) { +CbmTrdModuleRec* CbmTrdClusterFinder::AddModule(const CbmTrdDigi* digi) +{ Int_t address = digi->GetAddressModule(); CbmTrdModuleRec* module(NULL); - if (digi->GetType() == CbmTrdDigi::kFASP) - module = fModules[address] = new CbmTrdModuleRecT(address); + if (digi->GetType() == CbmTrdDigi::kFASP) module = fModules[address] = new CbmTrdModuleRecT(address); else module = fModules[address] = new CbmTrdModuleRecR(address); // try to load Geometry parameters for module const CbmTrdParModGeo* pGeo(NULL); - if (!fGeoPar - || !(pGeo = (const CbmTrdParModGeo*) fGeoPar->GetModulePar(address))) { - LOG(fatal) << GetName() << "::AddModule : No Geo params for module " - << address << ". Using default."; - } else + if (!fGeoPar || !(pGeo = (const CbmTrdParModGeo*) fGeoPar->GetModulePar(address))) { + LOG(fatal) << GetName() << "::AddModule : No Geo params for module " << address << ". Using default."; + } + else module->SetGeoPar(pGeo); // try to load read-out parameters for module const CbmTrdParModDigi* pDigi(NULL); - if (!fDigiPar - || !(pDigi = (const CbmTrdParModDigi*) fDigiPar->GetModulePar(address))) { - LOG(warn) << GetName() << "::AddModule : No Read-Out params for modAddress " - << address << ". Using default."; - } else + if (!fDigiPar || !(pDigi = (const CbmTrdParModDigi*) fDigiPar->GetModulePar(address))) { + LOG(warn) << GetName() << "::AddModule : No Read-Out params for modAddress " << address << ". Using default."; + } + else module->SetDigiPar(pDigi); // try to load ASIC parameters for module CbmTrdParSetAsic* pAsic(NULL); - if (!fAsicPar - || !(pAsic = (CbmTrdParSetAsic*) fAsicPar->GetModuleSet(address))) { - LOG(warn) << GetName() << "::AddModule : No ASIC params for modAddress " - << address << ". Using default."; + if (!fAsicPar || !(pAsic = (CbmTrdParSetAsic*) fAsicPar->GetModuleSet(address))) { + LOG(warn) << GetName() << "::AddModule : No ASIC params for modAddress " << address << ". Using default."; // module->SetAsicPar(); // map ASIC channels to read-out channels - need ParModDigi already loaded - } else + } + else module->SetAsicPar(pAsic); // try to load Chamber parameters for module const CbmTrdParModGas* pChmb(NULL); - if (!fGasPar - || !(pChmb = (const CbmTrdParModGas*) fGasPar->GetModulePar(address))) { - LOG(warn) << GetName() << "::AddModule : No Gas params for modAddress " - << address << ". Using default."; - } else + if (!fGasPar || !(pChmb = (const CbmTrdParModGas*) fGasPar->GetModulePar(address))) { + LOG(warn) << GetName() << "::AddModule : No Gas params for modAddress " << address << ". Using default."; + } + else module->SetChmbPar(pChmb); // try to load Gain parameters for module if (digi->GetType() == CbmTrdDigi::kFASP) { const CbmTrdParModGain* pGain(NULL); - if (!fGainPar - || !(pGain = - (const CbmTrdParModGain*) fGainPar->GetModulePar(address))) { + if (!fGainPar || !(pGain = (const CbmTrdParModGain*) fGainPar->GetModulePar(address))) { //LOG(warn) << GetName() << "::AddModule : No Gain params for modAddress "<< address <<". Using default."; - } else + } + else module->SetGainPar(pGain); } return module; } //_____________________________________________________________________ -void CbmTrdClusterFinder::SetParContainers() { - fAsicPar = static_cast<CbmTrdParSetAsic*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetAsic")); - fGasPar = static_cast<CbmTrdParSetGas*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGas")); - fDigiPar = static_cast<CbmTrdParSetDigi*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetDigi")); - fGainPar = static_cast<CbmTrdParSetGain*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGain")); +void CbmTrdClusterFinder::SetParContainers() +{ + fAsicPar = static_cast<CbmTrdParSetAsic*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetAsic")); + fGasPar = static_cast<CbmTrdParSetGas*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGas")); + fDigiPar = static_cast<CbmTrdParSetDigi*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetDigi")); + fGainPar = static_cast<CbmTrdParSetGain*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGain")); // Get the full geometry information of the TRD modules fGeoPar = new CbmTrdParSetGeo(); } //_____________________________________________________________________ -InitStatus CbmTrdClusterFinder::Init() { +InitStatus CbmTrdClusterFinder::Init() +{ FairRootManager* ioman = FairRootManager::Instance(); @@ -169,16 +204,20 @@ InitStatus CbmTrdClusterFinder::Init() { if (!CbmDigiManager::Instance()->IsPresent(ECbmModuleId::kTrd)) LOG(fatal); fClusters = new TClonesArray("CbmTrdCluster", 100); - ioman->Register( - "TrdCluster", "TRD", fClusters, IsOutputBranchPersistent("TrdCluster")); + ioman->Register("TrdCluster", "TRD", fClusters, IsOutputBranchPersistent("TrdCluster")); - if (!IsTimeBased() - && nullptr == dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"))) { + // If not deactivated by the user, the clusterizer will look for the CbmEvent branch, to only use Digis connected to a CbmEvent. If no CbmEvent branch is found all digis in the TrdDigi branch are automatically used. + if (UseOnlyEventDigis()) { + fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent")); + if (fEvents == nullptr) { SetUseOnlyEventDigis(kFALSE); } + } + if (!IsTimeBased() && nullptr == dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"))) { LOG(warn) << GetName() << ": Event mode selected but no event array found! Run in " "time-based mode."; SetTimeBased(); } + // // Get the full geometry information of the detector gas layers and store // // them with the CbmTrdModuleRec. This information can then be used for // // transformation calculations @@ -200,18 +239,20 @@ InitStatus CbmTrdClusterFinder::Init() { // fDigiPar->Initialize(); LOG(info) << "================ TRD Cluster Finder ==============="; - LOG(info) << " Free streaming : " << (IsTimeBased() ? "yes" : "no"); - LOG(info) << " Multi hit detect : " << (HasMultiHit() ? "yes" : "no"); - LOG(info) << " Row merger : " << (HasRowMerger() ? "yes" : "no"); - LOG(info) << " c-Neighbour enable: " << (HasNeighbourCol() ? "yes" : "no"); - LOG(info) << " r-Neighbour enable: " << (HasNeighbourRow() ? "yes" : "no"); - LOG(info) << " Write clusters : " << (HasDumpClusters() ? "yes" : "no"); + LOG(info) << " Free streaming : " << (IsTimeBased() ? "yes" : "no"); + LOG(info) << " Multi hit detect : " << (HasMultiHit() ? "yes" : "no"); + LOG(info) << " Row merger : " << (HasRowMerger() ? "yes" : "no"); + LOG(info) << " c-Neighbour enable : " << (HasNeighbourCol() ? "yes" : "no"); + LOG(info) << " r-Neighbour enable : " << (HasNeighbourRow() ? "yes" : "no"); + LOG(info) << " Write clusters : " << (HasDumpClusters() ? "yes" : "no"); + LOG(info) << " Use only event digis : " << (UseOnlyEventDigis() ? "yes" : "no"); return kSUCCESS; } //_____________________________________________________________________ -void CbmTrdClusterFinder::Exec(Option_t* /*option*/) { +void CbmTrdClusterFinder::Exec(Option_t* /*option*/) +{ /** * Digis are sorted according to the moduleAddress. A combiId is calculted based * on the rowId and the colId to have a neighbouring criterion for digis within @@ -229,63 +270,34 @@ void CbmTrdClusterFinder::Exec(Option_t* /*option*/) { TStopwatch timer; timer.Start(); - - // Int_t nentries = fDigis->GetEntries(); - //printf("processing %d entries\n", nentries); - CbmTrdModuleRec* mod(NULL); - Int_t nDigis = CbmDigiManager::Instance()->GetNofDigis(ECbmModuleId::kTrd); - for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) { - - const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(iDigi); - - - //CbmTrdDigi *digi = (CbmTrdDigi*) fDigis->At(iDigi); - Int_t moduleAddress = digi->GetAddressModule(); - - std::map<Int_t, CbmTrdModuleRec*>::iterator imod = - fModules.find(moduleAddress); - if (imod == fModules.end()) - mod = AddModule(digi); - else - mod = imod->second; - // std::cout<<digi->GetTime()<<std::endl; - - mod->AddDigi(digi, iDigi); + UInt_t nDigis = 0; + + if (UseOnlyEventDigis()) { + for (auto eventobj : *fEvents) { + auto event = static_cast<CbmEvent*>(eventobj); + nDigis = addDigisToModules(event); + processDigisInModules(nDigis, event); + fNrEvents++; + } } - Int_t digiCounter(0), clsCounter(0); - for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); - imod != fModules.end(); - imod++) { - mod = imod->second; - digiCounter += mod->GetOverThreshold(); - clsCounter += mod->FindClusters(); - AddClusters(mod->GetClusters(), kTRUE); + if (!UseOnlyEventDigis()) { + nDigis = addDigisToModules(); + processDigisInModules(nDigis); + fNrEvents++; } - // remove local data from all modules - for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); - imod != fModules.end(); - imod++) - imod->second->Clear("cls"); timer.Stop(); - LOG(info) << GetName() << "::Exec : Digis : " << nDigis << " / " - << digiCounter << " above threshold (" << 1e6 * fgMinimumChargeTH - << " keV)"; - LOG(info) << GetName() << "::Exec : Clusters : " << clsCounter; - LOG(info) << GetName() << "::Exec : real time=" << timer.RealTime() - << " CPU time=" << timer.CpuTime(); + LOG(debug) << GetName() << "::Exec : real time=" << timer.RealTime() << " CPU time=" << timer.CpuTime(); - // //cout << " " << counterI << " (" << counterI*100.0/Float_t(counterJ) << "%)" << " are reconstructed with fRowClusterMerger" << endl; - // //printf(" %4d modules (%6.3f%%) are reconstructed with fRowClusterMerger\n",counterI,counterI*100/Float_t(counterJ)); - // LOG(info) << "CbmTrdClusterFinder::Exec : RowClusterMerger are used " << fRowMergerCounter << " times"; + fProcessTime += timer.RealTime(); } //_____________________________________________________________________ -Int_t CbmTrdClusterFinder::AddClusters(TClonesArray* clusters, - Bool_t /* move*/) { +Int_t CbmTrdClusterFinder::AddClusters(TClonesArray* clusters, CbmEvent* event, Bool_t /* move*/) +{ if (!clusters) return 0; CbmTrdCluster *cls(NULL), *clsSave(NULL); const CbmTrdDigi* digi(NULL); @@ -310,7 +322,7 @@ Int_t CbmTrdClusterFinder::AddClusters(TClonesArray* clusters, cols.Clear(); rows.Clear(); for (Int_t id = 0; id < cls->GetNofDigis(); id++) { - digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cls->GetDigi(id)); + digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cls->GetDigi(id)); Int_t digiChannel = digi->GetAddressChannel(); Int_t colId = digiChannel % ncols; Int_t globalRow = digiChannel / ncols; @@ -323,8 +335,9 @@ Int_t CbmTrdClusterFinder::AddClusters(TClonesArray* clusters, cls->SetNCols(cols.CountBits()); cls->SetNRows(rows.CountBits()); } - clsSave = new ((*fClusters)[ncl++]) - CbmTrdCluster(*cls); // TODO implement copy constructor + clsSave = new ((*fClusters)[ncl++]) CbmTrdCluster(*cls); // TODO implement copy constructor + // In case we have an event branch and we did only use digis from within the event, add the cluster to the event. This allows the hit producer to identify wether or not to add the corresponding hit to the event. + if (event) event->AddData(ECbmDataType::kTrdCluster, ncl); clsSave->SetTrianglePads(cls->HasTrianglePads()); if (cls->GetMatch() != NULL) delete cls; //only the matches have pointers to allocated memory, so otherwise the clear does the trick @@ -335,6 +348,22 @@ Int_t CbmTrdClusterFinder::AddClusters(TClonesArray* clusters, } //_____________________________________________________________________ -void CbmTrdClusterFinder::Finish() {} +void CbmTrdClusterFinder::Finish() +{ + std::cout << std::endl; + LOG(info) << "====================================="; + LOG(info) << GetName() << ": Finish run"; + LOG(info) << GetName() << ": Run summary "; + LOG(info) << GetName() << ": Processing time : " << std::fixed << std::setprecision(3) << fProcessTime; + LOG(info) << GetName() << ": Nr of events : " << fNrEvents; + LOG(info) << GetName() << ": Nr of input digis : " << fNrDigis; + LOG(info) << GetName() << ": Nr of produced clusters : " << fNrClusters; + LOG(info) << GetName() << ": Nr of clusters / event : " << std::fixed << std::setprecision(2) + << (fNrEvents > 0 ? fNrClusters / (Double_t) fNrEvents : 0); + LOG(info) << GetName() << ": Nr of digis / cluster : " << std::fixed << std::setprecision(2) + << (fNrClusters > 0 ? fNrDigis / (Double_t) fNrClusters : 0); + LOG(info) << "====================================="; + std::cout << std::endl; +} ClassImp(CbmTrdClusterFinder) diff --git a/reco/detectors/trd/CbmTrdClusterFinder.h b/reco/detectors/trd/CbmTrdClusterFinder.h index 545bdd249633767d5193d8ef7d78bf343ad47cd1..d7d88315575e9555a7a1fe7731dc389e3d8c0809 100644 --- a/reco/detectors/trd/CbmTrdClusterFinder.h +++ b/reco/detectors/trd/CbmTrdClusterFinder.h @@ -1,8 +1,23 @@ +/** + * @file CbmTrdClusterFinder.h + * @author Alexandru Bercuci + * @author Pascal Raisig (praisig@ikf.uni-frankfurt.de) + * @brief FairTask to produce TrdCluster objects from TrdHit objects + * @version 1.0 + * @date 2021-03-16 + * + * + */ + #ifndef CBMTRDCLUSTERFINDER_H #define CBMTRDCLUSTERFINDER_H +#include "CbmEvent.h" + #include "FairTask.h" +#include <RtypesCore.h> + #include <map> #include <set> #include <vector> @@ -39,14 +54,16 @@ class CbmTrdClusterFinder : public FairTask { friend class CbmTrdModuleRecT; public: - enum CbmTrdRecDef { - kTime = 0, ///< select Time based/Event by event reconstruction - kMultiHit, ///< multi hit detection - kRowMerger, ///< merge clusters over neighbour rows - kNeighbourCol, ///< use neighbour trigger; column wise - kNeighbourRow, ///< use neighbour trigger; row wise - kDumpClusters, ///< write clustered digis to output - kFASP ///< use FASP ASIC for triangular pad plane geometry + enum CbmTrdRecDef + { + kTime = 0, ///< select Time based/Event by event reconstruction + kMultiHit, ///< multi hit detection + kRowMerger, ///< merge clusters over neighbour rows + kNeighbourCol, ///< use neighbour trigger; column wise + kNeighbourRow, ///< use neighbour trigger; row wise + kDumpClusters, ///< write clustered digis to output + kFASP, ///< use FASP ASIC for triangular pad plane geometry + kOnlyEventDigis ///< use only digis connected to a CbmEvent }; /** @@ -67,6 +84,15 @@ public: static Bool_t HasRowMerger() { return TESTBIT(fgConfig, kRowMerger); } static Bool_t IsTimeBased() { return TESTBIT(fgConfig, kTime); } + /** + * @brief If true only digis connected ro a CbmEvent are processed + * Per default this is activated on construction, such that the user can + * turn it off before Init(), where it will also be autmatically + * deactivated if no CbmEvent branch is found. + * @return Bool_t + */ + static Bool_t UseOnlyEventDigis() { return TESTBIT(fgConfig, kOnlyEventDigis); } + /** Initialisation **/ //virtual InitStatus ReInit(); virtual InitStatus Init(); @@ -79,22 +105,30 @@ public: virtual void Finish(); - static void SetDumpClusters(Bool_t set = kTRUE) { + static void SetDumpClusters(Bool_t set = kTRUE) + { set ? SETBIT(fgConfig, kDumpClusters) : CLRBIT(fgConfig, kDumpClusters); } - static void SetRowMerger(Bool_t set = kTRUE) { - set ? SETBIT(fgConfig, kRowMerger) : CLRBIT(fgConfig, kRowMerger); - } - static void SetMultiHit(Bool_t set = kTRUE) { - set ? SETBIT(fgConfig, kMultiHit) : CLRBIT(fgConfig, kMultiHit); - } - static void SetNeighbourEnable(Bool_t col = kTRUE, Bool_t row = kFALSE) { + static void SetRowMerger(Bool_t set = kTRUE) { set ? SETBIT(fgConfig, kRowMerger) : CLRBIT(fgConfig, kRowMerger); } + static void SetMultiHit(Bool_t set = kTRUE) { set ? SETBIT(fgConfig, kMultiHit) : CLRBIT(fgConfig, kMultiHit); } + static void SetNeighbourEnable(Bool_t col = kTRUE, Bool_t row = kFALSE) + { col ? SETBIT(fgConfig, kNeighbourCol) : CLRBIT(fgConfig, kNeighbourCol); row ? SETBIT(fgConfig, kNeighbourRow) : CLRBIT(fgConfig, kNeighbourRow); } static void SetMinimumChargeTH(Float_t th) { fgMinimumChargeTH = th; } - static void SetTimeBased(Bool_t set = kTRUE) { - set ? SETBIT(fgConfig, kTime) : CLRBIT(fgConfig, kTime); + static void SetTimeBased(Bool_t set = kTRUE) { set ? SETBIT(fgConfig, kTime) : CLRBIT(fgConfig, kTime); } + + /** + * @brief Set the Use Only Event Digis + * Per default this is activated on construction, such that the user can + * turn it off before Init(), where it will also be autmatically + * deactivated if no CbmEvent branch is found. + * @param set + */ + static void SetUseOnlyEventDigis(Bool_t set = kTRUE) + { + set ? SETBIT(fgConfig, kOnlyEventDigis) : CLRBIT(fgConfig, kOnlyEventDigis); } protected: @@ -105,39 +139,76 @@ private: CbmTrdClusterFinder(const CbmTrdClusterFinder&); CbmTrdClusterFinder& operator=(const CbmTrdClusterFinder&); - Int_t AddClusters(TClonesArray* clusters, Bool_t moveOwner = kTRUE); - CbmTrdModuleRec* AddModule(const CbmTrdDigi* d); + Int_t AddClusters(TClonesArray* clusters, CbmEvent* event, Bool_t moveOwner = kTRUE); - static Int_t - fgConfig; ///< Configuration map for the clusterizer. See CbmTrdRecDef for details - static Float_t fgMinimumChargeTH; ///< + /** + * @brief Add all digis available from CbmDigiManager to the corresponding modules + * + * @return UInt_t + */ + UInt_t addDigisToModules(); + /** + * @brief Add all digis connected to the passed event to the corresponding modules + * + * @param event + * @return UInt_t + */ + UInt_t addDigisToModules(CbmEvent* event); - TClonesArray* fClusters; /** Output array of CbmTrdCluster **/ + /** + * @brief Add the digi in the TrdDigi branch at the passed digiIdx to its corresponding module + * + * @param digiIdx index in the std::vector + */ + void addDigiToModule(UInt_t digiIdx); + + /** + * @brief Call the clusterizer function of each module + * + * @param ndigis + * @param event + */ + void processDigisInModules(UInt_t ndigis, CbmEvent* event = nullptr); - std::map<Int_t, std::set<Int_t>> fDigiMap; //! /** sector digis **/ - std::map<Int_t, std::set<Int_t>> fModuleMap; //! /** sector id per module **/ + /** + * @brief Adds the module corresponding to the address of the passed digi to the ModuleMap (fModules) + * + * @param digi + * @return CbmTrdModuleRec* + */ + CbmTrdModuleRec* AddModule(const CbmTrdDigi* digi); - std::set<Int_t> fNeighbours; - std::map<Int_t, std::set<Int_t>> - fModDigiMap; //std::map<Int_t ModuleID, std::vector<Int_t DigiID> > - std::map<Int_t, Int_t> fDigiRow; - std::map<Int_t, Int_t> fDigiCol; - std::map<Int_t, Double_t> fDigiCharge; + static Int_t fgConfig; ///< Configuration map for the clusterizer. See CbmTrdRecDef for details + static Float_t fgMinimumChargeTH; ///< - std::vector<std::set<Int_t>> fClusterBuffer; - std::map<Int_t, std::vector<std::set<Int_t>>> fModClusterDigiMap; + + TClonesArray* fClusters = nullptr; /** Output array of CbmTrdCluster **/ + + /** @brief Array connected to the CbmEvent branch */ + TClonesArray* fEvents = nullptr; //================================================================== - std::map<Int_t, CbmTrdModuleRec*> - fModules; ///< list of modules being processed - CbmTrdParSetAsic* fAsicPar; ///< parameter list for ASIC characterization - CbmTrdParSetGas* fGasPar; ///< parameter list for HV status - CbmTrdParSetDigi* fDigiPar; ///< parameter list for read-out geometry - CbmTrdParSetGain* fGainPar; ///< parameter list for keV->ADC gain conversion - CbmTrdParSetGeo* fGeoPar; ///< parameter list for modules geometry - - ClassDef(CbmTrdClusterFinder, 1); + std::map<Int_t, CbmTrdModuleRec*> fModules = {}; ///< list of modules being processed + CbmTrdParSetAsic* fAsicPar = nullptr; ///< parameter list for ASIC characterization + CbmTrdParSetGas* fGasPar = nullptr; ///< parameter list for HV status + CbmTrdParSetDigi* fDigiPar = nullptr; ///< parameter list for read-out geometry + CbmTrdParSetGain* fGainPar = nullptr; ///< parameter list for keV->ADC gain conversion + CbmTrdParSetGeo* fGeoPar = nullptr; ///< parameter list for modules geometry + + /** @brief Number of processed events (without CbmEvent corresponds to nr of exec calls) */ + UInt_t fNrEvents = 0; + + /** @brief Number of digis as input for the hit production. */ + UInt_t fNrDigis = 0; + + /** @brief Number of produced clusters. */ + UInt_t fNrClusters = 0; + + /** @brief Total processing time [RealTime]. */ + Float_t fProcessTime = 0; + + ClassDef(CbmTrdClusterFinder, 2); }; #endif diff --git a/reco/detectors/trd/CbmTrdHitProducer.cxx b/reco/detectors/trd/CbmTrdHitProducer.cxx index 6ea6da0dd994c40c4a4d6c8251340976c967ef74..3606b64af03d74ccbb6ac20992fdd7492d523304 100644 --- a/reco/detectors/trd/CbmTrdHitProducer.cxx +++ b/reco/detectors/trd/CbmTrdHitProducer.cxx @@ -1,10 +1,13 @@ #include "CbmTrdHitProducer.h" + +#include "CbmDefs.h" #include "CbmDigiManager.h" #include "CbmTrdAddress.h" -#include "CbmTrdCluster.h" +#include "CbmTrdClusterFinder.h" #include "CbmTrdDigi.h" #include "CbmTrdGeoHandler.h" #include "CbmTrdHit.h" +#include "CbmTrdModuleRec.h" #include "CbmTrdModuleRecR.h" #include "CbmTrdModuleRecT.h" #include "CbmTrdParAsic.h" @@ -17,122 +20,191 @@ #include "CbmTrdParSetGas.h" #include "CbmTrdParSetGeo.h" -#include "TGeoPhysicalNode.h" -#include <TClonesArray.h> -#include <TGeoManager.h> -#include <TStopwatch.h> -#include <TVector3.h> - #include <FairLogger.h> #include <FairRootManager.h> #include <FairRunAna.h> #include <FairRuntimeDb.h> +#include "TGeoPhysicalNode.h" +#include <RtypesCore.h> +#include <TGeoManager.h> +#include <TStopwatch.h> +#include <TVector3.h> + #include <map> //____________________________________________________________________________________ -CbmTrdHitProducer::CbmTrdHitProducer() - : FairTask("TrdHitProducer") - , fClusters(NULL) - , fHits(NULL) - , fModules() - , fAsicPar(NULL) - , fGasPar(NULL) - , fDigiPar(NULL) - , fGainPar(NULL) - , fGeoPar(NULL) {} +CbmTrdHitProducer::CbmTrdHitProducer() : FairTask("TrdHitProducer") {} //____________________________________________________________________________________ -CbmTrdHitProducer::~CbmTrdHitProducer() { +CbmTrdHitProducer::~CbmTrdHitProducer() +{ fHits->Clear(); delete fHits; if (fGeoPar) delete fGeoPar; } //____________________________________________________________________________________ -Int_t CbmTrdHitProducer::AddHits(TClonesArray* hits, Bool_t /*moveOwner*/) { +void CbmTrdHitProducer::addModuleHits(CbmTrdModuleRec* mod, Int_t* hitCounter, CbmEvent* event) +{ /** Absorb the TClonesArrays of the individual modules into the global TClonesArray. */ - if (!hits) return 0; + if (!mod) return; - Int_t nhits {hits->GetEntriesFast()}; - fHits->AbsorbObjects(hits); + auto hits = mod->GetHits(); - return nhits; + if (!hits) return; + + while (!hits->IsEmpty()) { + fHits->AbsorbObjects(hits, 0, 0); + if (event) event->AddData(ECbmDataType::kTrdHit, *hitCounter); + (*hitCounter)++; + fNrHits++; + } } //____________________________________________________________________________________ -CbmTrdModuleRec* CbmTrdHitProducer::AddModule(Int_t address, - TGeoPhysicalNode* node) { +CbmTrdModuleRec* CbmTrdHitProducer::AddModule(Int_t address, TGeoPhysicalNode* node) +{ TString s(node->GetName()); Int_t typ = TString(s[s.Index("module") + 6]).Atoi(); CbmTrdModuleRec* module(NULL); - if (typ == 9) { - module = fModules[address] = new CbmTrdModuleRecT(address); - } else { + if (typ == 9) { module = fModules[address] = new CbmTrdModuleRecT(address); } + else { module = fModules[address] = new CbmTrdModuleRecR(address); } // Try to load geometry parameters for the module const CbmTrdParModGeo* pGeo = nullptr; - if (!fGeoPar) - LOG(warn) << GetName() << ": No geometry parameter container!"; + if (!fGeoPar) LOG(warn) << GetName() << ": No geometry parameter container!"; else pGeo = (const CbmTrdParModGeo*) fGeoPar->GetModulePar(address); - if (!pGeo) - LOG(warn) << GetName() << ": No geometry parameters for module " << address; + if (!pGeo) LOG(warn) << GetName() << ": No geometry parameters for module " << address; else module->SetGeoPar(pGeo); // try to load read-out parameters for module const CbmTrdParModDigi* pDigi(NULL); - if (!fDigiPar - || !(pDigi = (const CbmTrdParModDigi*) fDigiPar->GetModulePar(address))) { - LOG(warn) << GetName() << "::AddModule : No Read-Out params for modAddress " - << address << ". Using default."; - } else + if (!fDigiPar || !(pDigi = (const CbmTrdParModDigi*) fDigiPar->GetModulePar(address))) { + LOG(warn) << GetName() << "::AddModule : No Read-Out params for modAddress " << address << ". Using default."; + } + else module->SetDigiPar(pDigi); // try to load ASIC parameters for module CbmTrdParSetAsic* pAsic(NULL); - if (!fAsicPar - || !(pAsic = (CbmTrdParSetAsic*) fAsicPar->GetModuleSet(address))) { - LOG(warn) << GetName() << "::AddModule : No ASIC params for modAddress " - << address << ". Using default."; + if (!fAsicPar || !(pAsic = (CbmTrdParSetAsic*) fAsicPar->GetModuleSet(address))) { + LOG(warn) << GetName() << "::AddModule : No ASIC params for modAddress " << address << ". Using default."; // module->SetAsicPar(); // map ASIC channels to read-out channels - need ParModDigi already loaded - } else + } + else module->SetAsicPar(pAsic); // try to load Chamber parameters for module const CbmTrdParModGas* pChmb(NULL); - if (!fGasPar - || !(pChmb = (const CbmTrdParModGas*) fGasPar->GetModulePar(address))) { - LOG(warn) << GetName() << "::AddModule : No Gas params for modAddress " - << address << ". Using default."; - } else + if (!fGasPar || !(pChmb = (const CbmTrdParModGas*) fGasPar->GetModulePar(address))) { + LOG(warn) << GetName() << "::AddModule : No Gas params for modAddress " << address << ". Using default."; + } + else module->SetChmbPar(pChmb); // try to load Gain parameters for module if (typ == 9) { const CbmTrdParModGain* pGain(NULL); - if (!fGainPar - || !(pGain = - (const CbmTrdParModGain*) fGainPar->GetModulePar(address))) { + if (!fGainPar || !(pGain = (const CbmTrdParModGain*) fGainPar->GetModulePar(address))) { //LOG(warn) << GetName() << "::AddModule : No Gain params for modAddress "<< address <<". Using default."; - } else + } + else module->SetGainPar(pGain); } return module; } + +// ---- processCluster ---- +UInt_t CbmTrdHitProducer::processClusters() +{ + Int_t nclusters = fClusters->GetEntries(); + LOG(debug) << GetName() << "::Exec: " + << " Clusters : " << nclusters; + + for (Int_t icluster = 0; icluster < nclusters; icluster++) { + processCluster(icluster); + } + auto nhits = addHits(); + LOG(debug) << GetName() << "::Exec: " + << " Hits : " << nhits; + return nhits; +} + +// ---- processCluster ---- +UInt_t CbmTrdHitProducer::processClusters(CbmEvent* event) +{ + Int_t nclusters = event->GetNofData(ECbmDataType::kTrdCluster); + LOG(debug) << GetName() << "::Exec : " + << " Clusters : " << nclusters; + + for (Int_t icluster = 0; icluster < nclusters; icluster++) { + auto clusterIdx = event->GetIndex(ECbmDataType::kTrdCluster, icluster); + processCluster(clusterIdx); + } + auto nhits = addHits(event); + LOG(debug) << GetName() << "::Exec : " + << " Hits : " << nhits; + return nhits; +} + + +// ---- processCluster ---- +void CbmTrdHitProducer::processCluster(const Int_t clusterIdx) +{ + auto cluster = static_cast<const CbmTrdCluster*>(fClusters->At(clusterIdx)); + if (!cluster) return; + + // get/build module for current cluster + auto imod = fModules.find(cluster->GetAddress()); + auto mod = imod->second; + + std::vector<const CbmTrdDigi*> digivec = {}; + // get digis for current cluster + for (Int_t iDigi = 0; iDigi < cluster->GetNofDigis(); iDigi++) { + const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(iDigi)); + + if (digi->GetType() == CbmTrdDigi::kSPADIC && digi->GetCharge() <= 0) continue; + digivec.emplace_back(digi); + } + + mod->MakeHit(clusterIdx, cluster, &digivec); + + fNrClusters++; +} + +// ---- addHits ---- +Int_t CbmTrdHitProducer::addHits(CbmEvent* event) +{ + Int_t hitCounter(0); + for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); imod != fModules.end(); imod++) { + auto mod = imod->second; + + mod->Finalize(); + + addModuleHits(mod, &hitCounter, event); + } + + // AbsorberObjects takes care of cleaning up. + // Hence, remove local data from all modules is not needed + return hitCounter; +} + //____________________________________________________________________________________ -InitStatus CbmTrdHitProducer::Init() { - FairRootManager* rootMgr = FairRootManager::Instance(); - if (NULL == rootMgr) { +InitStatus CbmTrdHitProducer::Init() +{ + FairRootManager* ioman = FairRootManager::Instance(); + if (NULL == ioman) { LOG(error) << GetName() << "::Init: " << "ROOT manager is not instantiated!"; return kFATAL; @@ -141,7 +213,7 @@ InitStatus CbmTrdHitProducer::Init() { CbmDigiManager::Instance()->Init(); if (!CbmDigiManager::Instance()->IsPresent(ECbmModuleId::kTrd)) LOG(fatal); - fClusters = (TClonesArray*) rootMgr->GetObject("TrdCluster"); + fClusters = (TClonesArray*) ioman->GetObject("TrdCluster"); if (!fClusters) { LOG(error) << GetName() << "::Init: " << "no TrdCluster array!"; @@ -149,7 +221,13 @@ InitStatus CbmTrdHitProducer::Init() { } fHits = new TClonesArray("CbmTrdHit", 100); - rootMgr->Register("TrdHit", "TRD", fHits, IsOutputBranchPersistent("TrdHit")); + ioman->Register("TrdHit", "TRD", fHits, IsOutputBranchPersistent("TrdHit")); + + // If not deactivated by the user, the hitproducer will look for the CbmEvent branch, to only use Digis connected to a CbmEvent. If no CbmEvent branch is found all digis in the TrdDigi branch are automatically used. + if (CbmTrdClusterFinder::UseOnlyEventDigis()) { + fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent")); + if (fEvents == nullptr) { CbmTrdClusterFinder::SetUseOnlyEventDigis(kFALSE); } + } // Get the full geometry information of the detector gas layers and store // them with the CbmTrdModuleRec. This information can then be used for @@ -160,15 +238,13 @@ InitStatus CbmTrdHitProducer::Init() { Int_t nrModules = fDigiPar->GetNrOfModules(); Int_t nrNodes = moduleMap.size(); if (nrModules != nrNodes) - LOG(fatal) << "CbmTrdHitProducer::Init() - Geometry(" << nrNodes - << ") and parameter files(" << nrModules + LOG(fatal) << "CbmTrdHitProducer::Init() - Geometry(" << nrNodes << ") and parameter files(" << nrModules << ") have different number of modules."; for (Int_t loop = 0; loop < nrModules; ++loop) { - Int_t address = fDigiPar->GetModuleId(loop); + Int_t address = fDigiPar->GetModuleId(loop); std::map<Int_t, TGeoPhysicalNode*>::iterator it = moduleMap.find(address); if (it == moduleMap.end()) { - LOG(fatal) << "Expected module with address " << address - << " wasn't found in the map with TGeoNode information."; + LOG(fatal) << "Expected module with address " << address << " wasn't found in the map with TGeoNode information."; } AddModule(address, it->second); } @@ -176,80 +252,63 @@ InitStatus CbmTrdHitProducer::Init() { } //____________________________________________________________________________________ -void CbmTrdHitProducer::Exec(Option_t*) { +void CbmTrdHitProducer::Exec(Option_t*) +{ fHits->Delete(); TStopwatch timer; timer.Start(); - CbmTrdModuleRec* mod(NULL); - std::vector<const CbmTrdDigi*> digis; - Int_t nofCluster = fClusters->GetEntries(); - for (Int_t iCluster = 0; iCluster < nofCluster; iCluster++) { - - const CbmTrdCluster* cluster = - static_cast<const CbmTrdCluster*>(fClusters->At(iCluster)); - // if(!cluster) continue; - - // get/build module for current cluster - std::map<Int_t, CbmTrdModuleRec*>::iterator imod = - fModules.find(cluster->GetAddress()); - mod = imod->second; - - // get digi for current cluster - for (Int_t iDigi = 0; iDigi < cluster->GetNofDigis(); iDigi++) { - const CbmTrdDigi* digi = - CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(iDigi)); - //const CbmTrdDigi* digi = static_cast<const CbmTrdDigi*>(fDigis->At(cluster->GetDigi(iDigi))); - if (digi->GetType() == CbmTrdDigi::kSPADIC && digi->GetCharge() <= 0) - continue; - digis.push_back(digi); - } + UInt_t hitCounter = 0; - // run hit reconstruction - // std::cout<<" make hit"<<std::endl; - mod->MakeHit(iCluster, cluster, &digis); - digis.clear(); + if (CbmTrdClusterFinder::UseOnlyEventDigis()) { + for (auto eventobj : *fEvents) { + hitCounter = 0; + auto event = static_cast<CbmEvent*>(eventobj); + if (!event) continue; + hitCounter += processClusters(event); + fNrEvents++; + } } - Int_t hitCounter(0); - for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); - imod != fModules.end(); - imod++) { - mod = imod->second; - //hitCounter += mod->GetNhits(); - - mod->Finalize(); - // std::cout<<" add hit"<<std::endl; - hitCounter += AddHits(mod->GetHits(), kTRUE); + if (!CbmTrdClusterFinder::UseOnlyEventDigis()) { + hitCounter = processClusters(); + fNrEvents++; } - // remove local data from all modules - // Not needed any longer since AbsorbObjects take care - // for(std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); imod!=fModules.end(); imod++) imod->second->Clear("hit"); + timer.Stop(); - LOG(info) << GetName() << "::Exec: " - << " Clusters : " << fClusters->GetEntriesFast(); - LOG(info) << GetName() << "::Exec: " - << " Hits : " << hitCounter; - LOG(info) << GetName() << "::Exec: real time=" << timer.RealTime() - << " CPU time=" << timer.CpuTime(); + LOG(debug) << GetName() << "::Exec: real time=" << timer.RealTime() << " CPU time=" << timer.CpuTime(); + fProcessTime += timer.RealTime(); } //____________________________________________________________________________________ -void CbmTrdHitProducer::Finish() {} +void CbmTrdHitProducer::Finish() +{ + std::cout << std::endl; + LOG(info) << "====================================="; + LOG(info) << GetName() << ": Finish run"; + LOG(info) << GetName() << ": Run summary "; + LOG(info) << GetName() << ": Processing time : " << std::fixed << std::setprecision(3) << fProcessTime; + LOG(info) << GetName() << ": Nr of events : " << fNrEvents; + LOG(info) << GetName() << ": Nr of input clusters : " << fNrClusters; + LOG(info) << GetName() << ": Nr of produced hits : " << fNrHits; + LOG(info) << GetName() << ": Nr of hits / event : " << std::fixed << std::setprecision(2) + << (fNrEvents > 0 ? fNrHits / (Double_t) fNrEvents : 0); + LOG(info) << GetName() << ": Nr of hits / clusters: " << std::fixed << std::setprecision(2) + << (fNrClusters > 0 ? fNrHits / (Double_t) fNrClusters : 0); + LOG(info) << "====================================="; + std::cout << std::endl; +} //________________________________________________________________________________________ -void CbmTrdHitProducer::SetParContainers() { - fAsicPar = static_cast<CbmTrdParSetAsic*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetAsic")); - fGasPar = static_cast<CbmTrdParSetGas*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGas")); - fDigiPar = static_cast<CbmTrdParSetDigi*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetDigi")); - fGainPar = static_cast<CbmTrdParSetGain*>( - FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGain")); - fGeoPar = new CbmTrdParSetGeo(); +void CbmTrdHitProducer::SetParContainers() +{ + fAsicPar = static_cast<CbmTrdParSetAsic*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetAsic")); + fGasPar = static_cast<CbmTrdParSetGas*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGas")); + fDigiPar = static_cast<CbmTrdParSetDigi*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetDigi")); + fGainPar = static_cast<CbmTrdParSetGain*>(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdParSetGain")); + fGeoPar = new CbmTrdParSetGeo(); } diff --git a/reco/detectors/trd/CbmTrdHitProducer.h b/reco/detectors/trd/CbmTrdHitProducer.h index b25bff17a17edbededab67e62eea4f9e00e5af5a..2621328899d70d83db1e601866ee158c81dbc500 100644 --- a/reco/detectors/trd/CbmTrdHitProducer.h +++ b/reco/detectors/trd/CbmTrdHitProducer.h @@ -1,9 +1,27 @@ +/** + * @file CbmTrdHitProducer.h + * @author Alexandru Bercuci + * @author Pascal Raisig (praisig@ikf.uni-frankfurt.de) + * @brief FairTask to produce TrdHit objects from TrdCluster objects + * @version 1.0 + * @date 2021-03-16 + * + * + */ + #ifndef CBMTRDHITPRODUCER_H #define CBMTRDHITPRODUCER_H +#include "CbmEvent.h" +#include "CbmTrdCluster.h" + #include "FairTask.h" -class TClonesArray; +#include <RtypesCore.h> +#include <TClonesArray.h> + +#include <limits> + class TGeoPhysicalNode; class CbmTrdParSetAsic; class CbmTrdParSetGas; @@ -11,7 +29,6 @@ class CbmTrdParSetDigi; class CbmTrdParSetGain; class CbmTrdParSetGeo; class CbmTrdModuleRec; -class CbmTrdCluster; class CbmTrdHitProducer : public FairTask { public: /** @@ -44,22 +61,76 @@ private: CbmTrdHitProducer(const CbmTrdHitProducer&); CbmTrdHitProducer& operator=(const CbmTrdHitProducer&); - Int_t AddHits(TClonesArray* hits, Bool_t moveOwner = kTRUE); CbmTrdModuleRec* AddModule(Int_t address, TGeoPhysicalNode* node); + /** + * @brief Process all clusters found in the TrdClusters branch + * + * @return UInt_t + */ + UInt_t processClusters(); + + /** + * @brief Process all clusters attached to the given event + * + * @param event + * @return UInt_t + */ + UInt_t processClusters(CbmEvent* event); + + /** + * @brief Produce a hit from the cluster found at clusterIdx in fClusters + * + * @param clusterIdx + */ + void processCluster(const Int_t clusterIdx); + + /** + * @brief Pass all hits produced by the given module to the TrdHit branch + * + * @param mod + * @param hitCounter + * @param event + */ + void addModuleHits(CbmTrdModuleRec* mod, Int_t* hitCounter, CbmEvent* event); + + /** + * @brief Loop over all modules in the given geometry and call addModuleHits(imodule) + * + * @param event + * @return Int_t + */ + Int_t addHits(CbmEvent* event = nullptr); + + + /** @brief Input array of CbmTrdCluster */ + TClonesArray* fClusters = nullptr; + /** @brief Output array of CbmTrdHit */ + TClonesArray* fHits = nullptr; + /** @brief Array connected to the CbmEvent branch */ + TClonesArray* fEvents = nullptr; + //================================================================== = {} + std::map<Int_t, CbmTrdModuleRec*> fModules = {}; ///< list of modules being processed + CbmTrdParSetAsic* fAsicPar = nullptr; ///< parameter list for ASIC characterization + CbmTrdParSetGas* fGasPar = nullptr; ///< parameter list for HV status + CbmTrdParSetDigi* fDigiPar = nullptr; ///< parameter list for read-out geometry + CbmTrdParSetGain* fGainPar = nullptr; ///< parameter list for keV->ADC gain conversion + CbmTrdParSetGeo* fGeoPar = nullptr; ///< parameter list for modules geometry + + /** @brief Number of processed events (without CbmEvent corresponds to nr of exec calls) */ + UInt_t fNrEvents = 0; + + /** @brief Number of clusters as input for the hit production. */ + UInt_t fNrClusters = 0; + + /** @brief Number of produced hits. */ + UInt_t fNrHits = 0; + + /** @brief Total processing time [RealTime]. */ + Float_t fProcessTime = 0; - TClonesArray* fClusters; /** Input array of CbmTrdCluster **/ - TClonesArray* fHits; /** Output array of CbmTrdHit **/ - //================================================================== - std::map<Int_t, CbmTrdModuleRec*> - fModules; ///< list of modules being processed - CbmTrdParSetAsic* fAsicPar; ///< parameter list for ASIC characterization - CbmTrdParSetGas* fGasPar; ///< parameter list for HV status - CbmTrdParSetDigi* fDigiPar; ///< parameter list for read-out geometry - CbmTrdParSetGain* fGainPar; ///< parameter list for keV->ADC gain conversion - CbmTrdParSetGeo* fGeoPar; ///< parameter list for modules geometry - ClassDef(CbmTrdHitProducer, 1); + ClassDef(CbmTrdHitProducer, 2); }; #endif