From 41a093d5942c6365aa5cb90f9670e8adbdfd30de Mon Sep 17 00:00:00 2001
From: Alexandru Bercuci <abercuci@niham.nipne.ro>
Date: Thu, 22 Jun 2023 16:36:48 +0300
Subject: [PATCH] add identification of EbyE/TB for simulation/data. Test TRD2D
reconstruction functionality on simulations and data with both EbyE and full
timeslice scenarios
---
reco/detectors/trd/CbmTrdClusterFinder.cxx | 62 +++++++++++++++-------
reco/detectors/trd/CbmTrdClusterFinder.h | 2 +-
reco/detectors/trd/CbmTrdHitProducer.cxx | 2 -
reco/detectors/trd/CbmTrdModuleRec2D.cxx | 21 ++++----
4 files changed, 55 insertions(+), 32 deletions(-)
diff --git a/reco/detectors/trd/CbmTrdClusterFinder.cxx b/reco/detectors/trd/CbmTrdClusterFinder.cxx
index abcdab7a8a..6fbd1cc512 100644
--- a/reco/detectors/trd/CbmTrdClusterFinder.cxx
+++ b/reco/detectors/trd/CbmTrdClusterFinder.cxx
@@ -6,6 +6,7 @@
#include "CbmDefs.h"
#include "CbmDigiManager.h"
+#include "CbmTimeSlice.h"
#include "CbmTrdCluster.h"
#include "CbmTrdDigi.h"
#include "CbmTrdGeoHandler.h"
@@ -50,7 +51,11 @@ using std::stringstream;
Int_t CbmTrdClusterFinder::fgConfig = 0;
Float_t CbmTrdClusterFinder::fgMinimumChargeTH = .5e-06;
//_____________________________________________________________________
-CbmTrdClusterFinder::CbmTrdClusterFinder() : FairTask("TrdClusterFinder", 1) { SetUseOnlyEventDigis(); }
+CbmTrdClusterFinder::CbmTrdClusterFinder() : FairTask("TrdClusterFinder", 1)
+{
+ SetUseOnlyEventDigis(true);
+ SetTimeBased(true);
+}
// --------------------------------------------------------------------
// ---- Destructor ----------------------------------------------------
@@ -79,7 +84,7 @@ Bool_t CbmTrdClusterFinder::AddCluster(CbmTrdCluster* c)
// ---- addDigisToModules ----
UInt_t CbmTrdClusterFinder::addDigisToModules()
{
- const int NDIGICHUNK = 100; // force flush of cluster buffer once every NDIGICHUNK digi to avoid memory exhaustion
+ const int NDIGICHUNK = 1000; // force flush of cluster buffer once every NDIGICHUNK digi to avoid memory exhaustion
UInt_t ndigis = static_cast<UInt_t>(std::abs(CbmDigiManager::Instance()->GetNofDigis(ECbmModuleId::kTrd)));
if (ndigis == 0) return 0;
@@ -87,16 +92,15 @@ UInt_t CbmTrdClusterFinder::addDigisToModules()
int jdigi(0);
for (size_t idigi(0); idigi < ndigis; idigi++) {
addDigiToModule(idigi);
+ jdigi++;
// once in a while dump finished clusters
// TODO ad hoc condition. Maybe a find a better one
- if (jdigi >= NDIGICHUNK) {
+ if (IsTimeBased() && jdigi >= NDIGICHUNK) {
processDigisInModules(jdigi, nullptr, false);
jdigi = 0;
}
- else
- jdigi++;
}
- return jdigi;
+ return ndigis;
}
// ---- addDigisToModules ----
@@ -132,8 +136,6 @@ void CbmTrdClusterFinder::addDigiToModule(UInt_t digiIdx)
// ---- processDigisInModules ----
void CbmTrdClusterFinder::processDigisInModules(UInt_t ndigis, CbmEvent* event, bool clr)
{
- // printf("AB :: processDigisInModules(%d, %p) clr[%c]\n", ndigis, (void*) event, (clr?'Y':'N'));
-
CbmTrdModuleRec* mod(NULL);
Int_t digiCounter(0), clsCounter(0);
for (std::map<Int_t, CbmTrdModuleRec*>::iterator imod = fModules.begin(); imod != fModules.end(); imod++) {
@@ -228,25 +230,44 @@ void CbmTrdClusterFinder::SetParContainers()
InitStatus CbmTrdClusterFinder::Init()
{
- FairRootManager* ioman = FairRootManager::Instance();
-
-
CbmDigiManager::Instance()->Init();
if (!CbmDigiManager::Instance()->IsPresent(ECbmModuleId::kTrd)) LOG(fatal) << GetName() << "Missing Trd digi branch.";
+ FairRootManager* ioman = FairRootManager::Instance();
+
fClusters = new TClonesArray("CbmTrdCluster", 100);
ioman->Register("TrdCluster", "TRD", fClusters, IsOutputBranchPersistent("TrdCluster"));
- // 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.
+ // Identify the time order of events based on the following conditions :
+ // 1. Existence of "CbmEvent" branch means that a form of DigiEventBuilder was run before (for
+ // both data and simulations) and it should be used UNLESS the user EXPLICITELY specify NOT TO.
+ // 2. Existence of "TimeSlice" branch "IsEvent() == true" for simulations
+ fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
+ bool digiEvent(false);
+ CbmTimeSlice* ts = dynamic_cast<CbmTimeSlice*>(ioman->GetObject("TimeSlice."));
+ if (ts != nullptr) { digiEvent = ts->IsEvent(); }
+ LOG(info) << GetName() << ": Event trigger " << (fEvents ? "found" : "miss") << "; digi organized "
+ << (digiEvent ? "EbyE" : "Timebased") << ".";
+
+ // If activated 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 (fEvents == nullptr) {
+ LOG(warn) << GetName()
+ << ": Event mode selected but no CbmEvent branch found ! Processing all digi from the list.";
+ SetUseOnlyEventDigis(false);
+ }
}
- 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();
+ else {
+ if (fEvents) {
+ LOG(warn) << GetName()
+ << ": CbmEvent branch found but full digi stream asked by user ! Processing all digi from the list.";
+ fEvents = nullptr;
+ }
+ }
+
+ if (IsTimeBased() && digiEvent) {
+ LOG(warn) << GetName() << ": Timebased mode selected but digi EbyE ! Processing digi from event (simulation EbyE).";
+ SetTimeBased(false);
}
// // Get the full geometry information of the detector gas layers and store
@@ -347,7 +368,8 @@ void CbmTrdClusterFinder::Exec(Option_t* /*option*/)
logOut << fixed << setw(8) << setprecision(1) << right << timerTs.RealTime() * 1000. << " ms] ";
logOut << "TS " << fNrTs;
if (UseOnlyEventDigis()) logOut << ", events " << nEvents;
- logOut << ", digis " << nDigisUsed << " / " << nDigisAll;
+ logOut << ", digis " << nDigisUsed << " / " << nDigisAll << " clusters "
+ << (fClusters ? fClusters->GetEntriesFast() : 0);
LOG(info) << logOut.str();
fNrTs++;
}
diff --git a/reco/detectors/trd/CbmTrdClusterFinder.h b/reco/detectors/trd/CbmTrdClusterFinder.h
index c01df96f91..fd558db750 100644
--- a/reco/detectors/trd/CbmTrdClusterFinder.h
+++ b/reco/detectors/trd/CbmTrdClusterFinder.h
@@ -60,7 +60,7 @@ class CbmTrdClusterFinder : public FairTask {
public:
enum ECbmTrdRecDef
{
- kTime = 0, ///< select Time based/Event by event reconstruction
+ kTime = 0, ///< select Time based/Event by event simulation (see CbmTrdDigitizer::kTime)
kMultiHit, ///< multi hit detection
kRowMerger, ///< merge clusters over neighbor rows
kNeighbourCol, ///< use neighbor trigger; column wise
diff --git a/reco/detectors/trd/CbmTrdHitProducer.cxx b/reco/detectors/trd/CbmTrdHitProducer.cxx
index c2179dc347..8ab1405b1a 100644
--- a/reco/detectors/trd/CbmTrdHitProducer.cxx
+++ b/reco/detectors/trd/CbmTrdHitProducer.cxx
@@ -218,10 +218,8 @@ 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->PreProcessHits();
mod->PostProcessHits();
-
hitCounter += addModuleHits(mod, event);
}
diff --git a/reco/detectors/trd/CbmTrdModuleRec2D.cxx b/reco/detectors/trd/CbmTrdModuleRec2D.cxx
index 27aecef2c3..ed8b845fa7 100644
--- a/reco/detectors/trd/CbmTrdModuleRec2D.cxx
+++ b/reco/detectors/trd/CbmTrdModuleRec2D.cxx
@@ -233,7 +233,7 @@ Int_t CbmTrdModuleRec2D::FindClusters(bool clr)
for (auto clRow : fBuffer)
ncl0 += clRow.second.size();
- if (CWRITE(0)) printf("AB :: FindClusters() cls_found = %d cls_write = %d cls_keep = %d\n", mcl, ncl, ncl0);
+ if (CWRITE(0)) printf("FindClusters() cls_found = %d cls_write = %d cls_keep = %d\n", mcl, ncl, ncl0);
return ncl;
}
@@ -274,7 +274,6 @@ Bool_t CbmTrdModuleRec2D::PostProcessHits()
/** Steering routine for classifying hits and apply further analysis
* -> hit merging for row-cross (see RowCross)
*/
- return true;
CbmTrdHit *h0(nullptr), *h1(nullptr);
Int_t a0, nhits = fHits->GetEntriesFast();
@@ -287,15 +286,17 @@ Bool_t CbmTrdModuleRec2D::PostProcessHits()
h1 = (CbmTrdHit*) ((*fHits)[jh]);
if (h1->IsUsed()) continue;
+ // basic check on Time
+ if (h1->GetTime() < 4000 - h0->GetTime()) continue; // TODO check with preoper time calibration
+ // skip next hits for being too far (10us) in the future
+ if (h1->GetTime() > 10000 + h0->GetTime()) break;
+
// basic check on Row
if (TMath::Abs(h1->GetY() - h0->GetY()) > Dy) continue;
// basic check on Col
if (TMath::Abs(h1->GetX() - h0->GetX()) > Dx) continue;
- // basic check on Time
- if (TMath::Abs(h1->GetTime() - h0->GetTime()) > 4000) continue; // TODO check with preoper time calibration
-
// go to topologic checks
if (!(a0 = CheckMerge(h0->GetRefId(), h1->GetRefId()))) continue;
@@ -795,7 +796,7 @@ Bool_t CbmTrdModuleRec2D::BuildHit(CbmTrdHit* h)
h->SetDz(0);
h->SetDxy(0.);
h->SetTime(CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kFASP) * (vt0 + time) - tdrift + 30.29 + fHitTimeOff, edt);
- h->SetELoss(e);
+ h->SetELoss(e * 1.e-9); // conversion to GeV
h->SetClassType();
h->SetMaxType(IsMaxTilt());
h->SetOverFlow(HasOvf());
@@ -913,7 +914,7 @@ Double_t CbmTrdModuleRec2D::GetXoffset(Int_t n0) const
dx += vs[ir] * vx[ir];
}
if (TMath::Abs(R) > 0) return dx / R;
- LOG(warn) << GetName() << "::GetXoffset : Null total charge for hit size " << n;
+ LOG(debug) << GetName() << "::GetXoffset : Null total charge for hit size " << n;
return 0.;
}
@@ -929,7 +930,7 @@ Double_t CbmTrdModuleRec2D::GetYoffset(Int_t n0) const
}
}
if (TMath::Abs(T) > 0) return dy / T;
- LOG(warn) << GetName() << "::GetYoffset : Null total charge for hit size " << n;
+ LOG(debug) << GetName() << "::GetYoffset : Null total charge for hit size " << n;
//if (CWRITE(1))
return 0.;
}
@@ -1066,6 +1067,7 @@ Int_t CbmTrdModuleRec2D::LoadDigis(vector<const CbmTrdDigi*>* din, Int_t cid)
/** Load RAW digis info in calibration aware strucuture CbmTrdDigiReco
* Do basic sanity checks; also incomplete adjacent digi and if found merge them.
*/
+
Int_t col(-1), /*row, */ colT(-1), colR(-1);
const CbmTrdDigi *dgT(nullptr), *dgR(nullptr);
for (vector<const CbmTrdDigi*>::iterator i = din->begin(), j = i + 1; i != din->end(); i++) {
@@ -1204,7 +1206,7 @@ Int_t CbmTrdModuleRec2D::ProjectDigis(Int_t cid, Int_t cjd)
if (CWRITE(1)) {
if (dg0) cout << "dgR :" << dg0->ToString();
if (dg1) cout << "dgT :" << dg1->ToString();
- cout << "-------------------------------------" << endl;
+ //cout << "-------------------------------------" << endl;
}
// process tilt signal/time
@@ -1376,6 +1378,7 @@ Int_t CbmTrdModuleRec2D::LoadDigis(vector<const CbmTrdDigi*>* digis, vector<CbmT
* (pad width [pw], DAQ time [clk], signal [ADC chs]) with rectangular signals at integer
* positions.
*/
+
vs.clear();
vse.clear();
vx.clear();
--
GitLab