From 94b6940aae6dca18c5d313e1eeef4982d9c81f70 Mon Sep 17 00:00:00 2001
From: Vikas Singhal <vsinghal@gsi.de>
Date: Tue, 4 Jul 2023 18:02:45 +0530
Subject: [PATCH] mMuCh Hit reconstruction and including mMuCh in mcbm L1 macro
---
core/base/CbmDigiManager.cxx | 6 +-
core/base/CbmDigiManager.h | 3 +-
core/base/utils/CbmMcbmUtils.cxx | 11 +-
macro/beamtime/mcbm2022/mcbm_event_reco_L1.C | 48 ++++++++-
reco/detectors/much/CbmMuchFindHitsGem.cxx | 102 ++++++++++++++++---
reco/detectors/much/CbmMuchFindHitsGem.h | 1 +
6 files changed, 147 insertions(+), 24 deletions(-)
diff --git a/core/base/CbmDigiManager.cxx b/core/base/CbmDigiManager.cxx
index d47ba11247..c6161434a4 100644
--- a/core/base/CbmDigiManager.cxx
+++ b/core/base/CbmDigiManager.cxx
@@ -134,9 +134,11 @@ void CbmDigiManager::SetBranch()
// Get system ID and class name from digi class.
ECbmModuleId systemId = Digi::GetSystem();
string className = Digi::GetClassName();
-
+ std::cout << "systemId= " << systemId << " className= " << className << std::endl;
// TODO: Remove ugly fix for CbmMuchBeamTimeDigi once class has disappeared.
- if (systemId == ECbmModuleId::kMuch && fUseMuchBeamTimeDigi) className = "CbmMuchBeamTimeDigi";
+ // VS: CbmMuchBeamTimeDigi may be removed now. 26/06/23
+ //if (systemId == ECbmModuleId::kMuch && fUseMuchBeamTimeDigi) className = "CbmMuchBeamTimeDigi";
+ if (systemId == ECbmModuleId::kMuch) className = "CbmMuchDigi";
// --- Catch branch being already set
if (fBranches.find(systemId) != fBranches.end()) {
diff --git a/core/base/CbmDigiManager.h b/core/base/CbmDigiManager.h
index 6a88c15bf9..a803e09f15 100644
--- a/core/base/CbmDigiManager.h
+++ b/core/base/CbmDigiManager.h
@@ -160,7 +160,8 @@ public:
**
** Temporary solution until the classes are unified.
**/
- void UseMuchBeamTimeDigi(Bool_t choice = kTRUE) { fUseMuchBeamTimeDigi = choice; }
+ //void UseMuchBeamTimeDigi(Bool_t choice = kTRUE) { fUseMuchBeamTimeDigi = choice; }
+ void UseMuchBeamTimeDigi(Bool_t choice = kTRUE) { fUseMuchBeamTimeDigi = kFALSE; }
private:
diff --git a/core/base/utils/CbmMcbmUtils.cxx b/core/base/utils/CbmMcbmUtils.cxx
index 8cd1a700c3..291f600470 100644
--- a/core/base/utils/CbmMcbmUtils.cxx
+++ b/core/base/utils/CbmMcbmUtils.cxx
@@ -39,14 +39,17 @@ namespace cbm
/// Nickel runs: 2350 - 2397 = 23/05/2022 - 25/05/2022 (Lambda Benchmark but mTOF troubles)
sSetupName = "mcbm_beam_2022_05_23_nickel";
}
- else if (2454 <= ulRunId && ulRunId <= 2610) {
+ else if (2454 <= ulRunId && ulRunId <= 2497) {
/// Lambda Benchmark Gold runs: 2454 - 2497 = 16/06/2022 - 18/06/2022
sSetupName = "mcbm_beam_2022_06_16_gold";
}
- /// High Rate Gold runs with GEMs in Acceptance: 2498 - 2610 = 18/06/2022 - 20/06/2022
- else if (2610 < ulRunId) {
+ else if (2498 <= ulRunId && ulRunId <= 2610) {
+ /// High Rate Gold runs with GEMs in Acceptance: 2498 - 2610 = 18/06/2022 - 20/06/2022
+ sSetupName = "mcbm_beam_2022_06_18_gold";
+ }
+ else if (2611 < ulRunId) {
/// Future runs, exception there to force implementation and support from users side
- throw(std::invalid_argument("RunId bigger than latest run mapped (2610, mCBM 2022)! Please complete the map!"));
+ throw(std::invalid_argument("RunId bigger than latest run mapped (2611, mCBM 2022)! Please complete the map!"));
}
return sSetupName;
diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
index 49fd93511d..469aa44b8e 100644
--- a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
+++ b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
@@ -15,7 +15,7 @@
/// FIXME: Disable clang formatting to keep easy parameters overview
/* clang-format off */
-Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
+Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
Int_t nTimeslices = 10,
TString sInpDir = "./data/",
TString sOutDir = "./data/",
@@ -25,7 +25,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
Bool_t bTRD = kTRUE,
Bool_t bTRD2d = kTRUE,
Bool_t bRICH = kTRUE,
- Bool_t bMUCH = kFALSE,
+ Bool_t bMUCH = kTRUE,
Bool_t bTOF = kTRUE,
Bool_t bTOFtr = kTRUE,
Bool_t bPSD = kFALSE,
@@ -56,7 +56,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
TString sRunId = TString::Format("%04u", uRunId);
/// Initial pattern
- TString inFile = sInpFile; // + "/" + sRunId + ".digi";
+ TString inFile = sInpDir + sRunId + ".digi";
TString cFileId = sRunId;
//TString parFileIn = sInpDir + "/unp_mcbm_params_" + sRunId;
@@ -261,7 +261,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
const Int_t eb_TriggerMaxNumberT0 {2};
const Int_t eb_TriggerMinNumberSts {2};
const Int_t eb_TriggerMinNumberStsLayers {1};
- const Int_t eb_TriggerMinNumberMuch {0};
+ const Int_t eb_TriggerMinNumberMuch {1};
const Int_t eb_TriggerMinNumberTof {8};
const Int_t eb_TriggerMinNumberTofLayers {4};
const Int_t eb_TriggerMinNumberRich {0};
@@ -270,6 +270,8 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
const Int_t eb_TrigWinMaxT0 {50};
const Int_t eb_TrigWinMinSts {-60};
const Int_t eb_TrigWinMaxSts {60};
+ const Int_t eb_TrigWinMinMuch {-100};
+ const Int_t eb_TrigWinMaxMuch {100};
const Int_t eb_TrigWinMinTrd1d {-300};
const Int_t eb_TrigWinMaxTrd1d {300};
const Int_t eb_TrigWinMinTrd2d {-200};
@@ -284,7 +286,8 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
//Choose between NoOverlap, MergeOverlap, AllowOverlap
evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap);
-
+ //For time being it is needed. We will remove CbmMuchBeamTimeDigi.
+ evBuildRaw->ChangeMuchBeamtimeDigiFlag();
// Remove detectors where digis not found
evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
if (!bRICH || !geoSetup->IsActive(ECbmModuleId::kRich)) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
@@ -297,6 +300,8 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
// Set T0 as reference detector
evBuildRaw->SetReferenceDetector(kRawEventBuilderDetT0);
+ // For making MuCh as seed detector
+ // evBuildRaw->SetReferenceDetector(kRawEventBuilderDetMuch);
// Use sliding window seed builder with STS
//evBuildRaw->SetReferenceDetector(kRawEventBuilderDetUndef);
@@ -325,6 +330,10 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
evBuildRaw->SetTriggerMinLayersNumber(ECbmModuleId::kSts, eb_TriggerMinNumberStsLayers);
evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
}
+ if (bMUCH && geoSetup->IsActive(ECbmModuleId::kMuch)) {
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kMuch, eb_TriggerMinNumberMuch);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kMuch, -1);
+ }
if (geoSetup->IsActive(ECbmModuleId::kRich)) {
evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kRich, eb_TriggerMinNumberRich);
evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kRich, -1);
@@ -335,6 +344,8 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
evBuildRaw->SetTriggerWindow(ECbmModuleId::kTof, eb_TrigWinMinTof, eb_TrigWinMaxTof);
if (geoSetup->IsActive(ECbmModuleId::kSts))
evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, eb_TrigWinMinSts, eb_TrigWinMaxSts);
+ if (bMUCH && geoSetup->IsActive(ECbmModuleId::kMuch))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kMuch, eb_TrigWinMinMuch, eb_TrigWinMaxMuch);
if (geoSetup->IsActive(ECbmModuleId::kTrd))
evBuildRaw->SetTriggerWindow(ECbmModuleId::kTrd, eb_TrigWinMinTrd1d, eb_TrigWinMaxTrd1d);
if (geoSetup->IsActive(ECbmModuleId::kTrd))
@@ -444,6 +455,33 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
// ------------------------------------------------------------------------
}
+ // =========================================================================
+ // === local MUCH Reconstruction ===
+ // =========================================================================
+ if (bMUCH) {
+ TString muchGeoTag;
+ if (geoSetup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) {
+ // --- Parameter file name
+ TString geoTag;
+ geoSetup->GetGeoTag(ECbmModuleId::kMuch, geoTag);
+ Int_t muchFlag = 0;
+ if (geoTag.Contains("mcbm")) muchFlag = 1;
+ TString sectorFile = gSystem->Getenv("VMCWORKDIR");
+ sectorFile += "/parameters/much/much_" + geoTag(0, 4) + "_mcbm_digi_sector.root";
+ //sectorFile += "/parameters/much/much_v22j_mcbm_digi_sector.root";
+ std::cout << "Using parameter file " << sectorFile << std::endl;
+
+ // --- Initialization of the digi scheme
+ auto muchGeoScheme = CbmMuchGeoScheme::Instance();
+ if (!muchGeoScheme->IsInitialized()) { muchGeoScheme->Init(sectorFile.Data(), muchFlag); }
+ // --- Hit finder for GEMs
+ FairTask* muchReco = new CbmMuchFindHitsGem(sectorFile.Data(), muchFlag);
+ run->AddTask(muchReco);
+ std::cout << "-I- " << myName << ": Added task " << muchReco->GetName() << " with parameter file " << sectorFile
+ << std::endl;
+ }
+ }
+ // ------------------------------------------------------------------------
// =========================================================================
// === local TRD Reconstruction ===
diff --git a/reco/detectors/much/CbmMuchFindHitsGem.cxx b/reco/detectors/much/CbmMuchFindHitsGem.cxx
index 4f89b5ebf4..47e4fa8a04 100644
--- a/reco/detectors/much/CbmMuchFindHitsGem.cxx
+++ b/reco/detectors/much/CbmMuchFindHitsGem.cxx
@@ -123,9 +123,7 @@ InitStatus CbmMuchFindHitsGem::Init()
TDirectory* oldDir = gDirectory;
TFile* file = new TFile(fDigiFile);
- LOG_IF(fatal, !file) << "Could not open file " << fDigiFile;
- TObjArray* stations = file->Get<TObjArray>("stations");
- LOG_IF(fatal, !stations) << "TObjArray stations not found in file " << fDigiFile;
+ TObjArray* stations = (TObjArray*) file->Get("stations");
file->Close();
file->Delete();
/// Restore old global file and folder pointer to avoid messing with FairRoot
@@ -209,6 +207,11 @@ void CbmMuchFindHitsGem::Exec(Option_t*)
} //? event mode
fNofTimeslices++;
+ /* double *pos = IgnoredAddresses.data();
+ for (int i=0;i<IgnoredAddresses.size();i++)
+ LOG(info) << "Ignored Address " << *pos++ ; */
+ LOG(info) << "Total ignored address " << IgnoredAddresses.size();
+ IgnoredAddresses.clear();
}
// -------------------------------------------------------------------------
@@ -311,12 +314,47 @@ void CbmMuchFindHitsGem::FindClusters(CbmEvent* event)
// Double_t chanid = digi->GetChannelId();
UInt_t address = digi->GetAddress();
// UInt_t adc = digi->GetAdc();
+ // 22/02/23 for mCBM Data checking for GEM and RPC sector and channel limits in a module.
+ //
+ if (fFlag) {
+ if ((CbmMuchAddress::GetStationIndex(digi->GetAddress()) == 0)
+ || (CbmMuchAddress::GetStationIndex(digi->GetAddress()) == 1)) {
+ if (CbmMuchAddress::GetSectorIndex(digi->GetAddress()) < 0
+ || CbmMuchAddress::GetSectorIndex(digi->GetAddress()) > 95) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+ if (CbmMuchAddress::GetChannelIndex(digi->GetAddress()) < 0
+ || CbmMuchAddress::GetChannelIndex(digi->GetAddress()) > 22) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+ }
+ if ((CbmMuchAddress::GetStationIndex(digi->GetAddress()) == 2)
+ || (CbmMuchAddress::GetStationIndex(digi->GetAddress()) == 3)) {
+ if (CbmMuchAddress::GetSectorIndex(digi->GetAddress()) < 0
+ || CbmMuchAddress::GetSectorIndex(digi->GetAddress()) > 45) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+ if (CbmMuchAddress::GetChannelIndex(digi->GetAddress()) < 0
+ || CbmMuchAddress::GetChannelIndex(digi->GetAddress()) > 10) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+ }
+ if (fGeoScheme->GetModuleByDetId(address) == NULL) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+ }
fGeoScheme->GetModuleByDetId(address)->AddDigi(time, digiIndex);
}
// Find clusters module-by-module
for (UInt_t m = 0; m < modules.size(); m++) {
CbmMuchModuleGem* module = modules[m];
+ assert(module);
multimap<Double_t, Int_t> digis = modules[m]->GetDigis();
multimap<Double_t, Int_t>::iterator it, itmin, itmax;
@@ -341,7 +379,14 @@ void CbmMuchFindHitsGem::FindClusters(CbmEvent* event)
else
digi = static_cast<const CbmMuchDigi*>(fDigiManager->Get<CbmMuchBeamTimeDigi>(iDigi));
//const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(iDigi));
+
+ if (module->GetPad(digi->GetAddress()) == nullptr) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+
CbmMuchPad* pad = module->GetPad(digi->GetAddress());
+
pad->SetDigiIndex(iDigi);
fFiredPads.push_back(pad);
}
@@ -401,6 +446,7 @@ void CbmMuchFindHitsGem::ExecClusteringSimple(CbmMuchCluster* cluster, Int_t iCl
//CbmMuchDigi* digi = static_cast<CbmMuchDigi*>(fDigis->At(cluster->GetDigi(0)));
CbmMuchModule* m = fGeoScheme->GetModuleByDetId(digi->GetAddress());
CbmMuchModuleGem* module = (CbmMuchModuleGem*) m;
+ assert(module);
// Int_t iStation = CbmMuchAddress::GetStationIndex(digi->GetAddress());
Int_t maxCharge = 0;
@@ -427,6 +473,11 @@ void CbmMuchFindHitsGem::ExecClusteringSimple(CbmMuchCluster* cluster, Int_t iCl
digi = static_cast<const CbmMuchDigi*>(fDigiManager->Get<CbmMuchBeamTimeDigi>(digiIndex));
//digi = static_cast<CbmMuchDigi*> (fDigis->At(digiIndex));
if (digi->GetAdc() <= threshold) continue;
+ if (module->GetPad(digi->GetAddress()) == nullptr) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
+
CbmMuchPad* pad = module->GetPad(digi->GetAddress());
pad->SetDigiIndex(digiIndex);
fFiredPads.push_back(pad);
@@ -469,6 +520,11 @@ void CbmMuchFindHitsGem::ExecClusteringPeaks(CbmMuchCluster* cluster, Int_t iClu
//const CbmMuchDigi* digi = static_cast<const CbmMuchDigi*>(fDigis->At(iDigi));
UInt_t address = digi->GetAddress();
CbmMuchModuleGem* module = (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(address);
+ assert(module);
+ if (module->GetPad(digi->GetAddress()) == nullptr) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
CbmMuchPad* pad = module->GetPad(address);
Int_t adc = digi->GetAdc();
fClusterPads.push_back(pad);
@@ -526,8 +582,7 @@ void CbmMuchFindHitsGem::ExecClusteringPeaks(CbmMuchCluster* cluster, Int_t iClu
void CbmMuchFindHitsGem::CreateHits(CbmMuchCluster* cluster, Int_t iCluster, CbmEvent* event)
{
Int_t nDigis = cluster->GetNofDigis();
- Double_t sumq = 0, sumx = 0, sumy = 0, sumdx2 = 0, sumdy2 = 0, sumdxy2 = 0,
- sumdt2 = 0; // , sumt =0 // not used FU 22.03.23
+ Double_t sumq = 0, sumx = 0, sumy = 0, sumt = 0, sumdx2 = 0, sumdy2 = 0, sumdxy2 = 0, sumdt2 = 0;
Double_t q = 0, x = 0, y = 0, t = 0, z = 0, dx = 0, dy = 0, dxy = 0, dt = 0;
Double_t nX = 0, nY = 0, nZ = 0;
Int_t address = 0;
@@ -548,11 +603,18 @@ void CbmMuchFindHitsGem::CreateHits(CbmMuchCluster* cluster, Int_t iCluster, Cbm
if (i == 0) {
address = CbmMuchAddress::GetElementAddress(digi->GetAddress(), kMuchModule);
planeId = fGeoScheme->GetLayerSideNr(address);
- module = (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(address);
+ // Below address should be 32 bit CbmMuchAddress therefore changing it. Need to check for simulation. 22/02/23
+ module = (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(digi->GetAddress());
+ assert(module);
z = module->GetPosition()[2];
}
+ // One cluster must be from single module. No need of below line.
+ // CbmMuchModule* moduleDet = fGeoScheme->GetModuleByDetId(digi->GetAddress()); /// added
- CbmMuchModule* moduleDet = fGeoScheme->GetModuleByDetId(digi->GetAddress()); /// added
+ if (module->GetPad(digi->GetAddress()) == nullptr) {
+ IgnoredAddresses.push_back(digi->GetAddress());
+ continue;
+ }
CbmMuchPad* pad = module->GetPad(digi->GetAddress());
x = pad->GetX();
@@ -566,14 +628,14 @@ void CbmMuchFindHitsGem::CreateHits(CbmMuchCluster* cluster, Int_t iCluster, Cbm
2.3; // Hit time error for RPC = 2.3 ns, as residual dist width is 2.3. This is made to make pull dist width ~1
Double_t timeShift = 0.5; // this is added because residual time dist is shifted by -0.5 from 0.
- if (moduleDet->GetDetectorType() == 3) ///GEM
+ if (module->GetDetectorType() == 3) ///GEM
{
if (fFlag == 0) t = digi->GetTime() - gemDriftTimeCorrc + timeShift;
else
t = digi->GetTime(); // Not correcting Drift Time for mCBM data
dt = gemHitTimeError;
}
- if (moduleDet->GetDetectorType() == 4) ////RPC
+ if (module->GetDetectorType() == 4) ////RPC
{
t = digi->GetTime() - rpcDriftTimeCorrc + timeShift;
dt = rpcHitTimeError;
@@ -589,7 +651,7 @@ void CbmMuchFindHitsGem::CreateHits(CbmMuchCluster* cluster, Int_t iCluster, Cbm
sumq += q;
sumx += q * x;
sumy += q * y;
- // sumt += q * t; // not used FU 22.03.23
+ sumt += q * t;
sumdx2 += q * q * dx * dx;
sumdy2 += q * q * dy * dy;
sumdxy2 += q * q * dxy * dxy;
@@ -607,8 +669,24 @@ void CbmMuchFindHitsGem::CreateHits(CbmMuchCluster* cluster, Int_t iCluster, Cbm
dt = sqrt(sumdt2) / sumq;
Int_t iHit = fHits->GetEntriesFast();
- //------------------------------Added by O. Singh 11.12.2017 for mCbm ---------------------------
- Double_t tX = 18.5, tY = 80.0;
+ //------------------------------Added by A. Agarwal 29.09.2022 for mCbm ---------------------------
+ // Double_t tX = 8.5, tY = 81.0;
+ Double_t tX = 8.5, tY = 81.0; // commented for testing
+ //Double_t tX =18.5 , tY = 80.0 ;
+ if (module->GetDetectorType() == 3) // GEM
+ {
+ tX = 8.5;
+ tY = 81.0;
+ }
+ else if (module->GetDetectorType() == 4) // RPC
+ {
+ tX = -66.5;
+ tY = 72.0;
+ }
+ else {
+ LOG(error) << "Unknown detector type";
+ }
+
nX = x + tX; // Ajit + OS + Apar -> For miniMUCH setup in March 2019
nY = y + tY; // Ajit + OS + Apar -> For miniMUCH setup in March 2019
nZ = z;
diff --git a/reco/detectors/much/CbmMuchFindHitsGem.h b/reco/detectors/much/CbmMuchFindHitsGem.h
index 5fa323bd32..2ef31f995e 100644
--- a/reco/detectors/much/CbmMuchFindHitsGem.h
+++ b/reco/detectors/much/CbmMuchFindHitsGem.h
@@ -85,6 +85,7 @@ private:
std::vector<Bool_t> fLocalMax; //!
std::vector<CbmMuchPad*> fClusterPads; //!
std::vector<std::vector<Int_t>> fNeighbours; //!
+ std::vector<Double_t> IgnoredAddresses;
TClonesArray* fClusters; // Output array of CbmMuchCluster objects
TClonesArray* fHits; // Output array of CbmMuchHit
--
GitLab