diff --git a/include/MSubModuleDepthReadout.h b/include/MSubModuleDepthReadout.h index d894eaba..a4fbda6b 100644 --- a/include/MSubModuleDepthReadout.h +++ b/include/MSubModuleDepthReadout.h @@ -93,7 +93,7 @@ class MSubModuleDepthReadout : public MSubModule protected: //! Load the splines file - bool LoadSplinesFile(MString FileName); + bool LoadStripCoeffsFile(MString FileName); // private methods: @@ -110,8 +110,12 @@ class MSubModuleDepthReadout : public MSubModule //! Filename of the depth calibration coefficients (stretch, offset, timing noise, ...) MString m_DepthCoefficientsFileName; - //! Map of the depth calibration coefficients - unordered_map> m_Coeffs; + //! Map: detector ID (int) -> mean stretch over all pixels / strips + unordered_map m_MeanStretch; + //! Map: detector ID (int) -> mean offset over all pixels / strips + unordered_map m_MeanOffset; + //! Map: detector ID (int) -> Side (LV=0, HV=1) -> Strip ID -> Depth calibration coefficients (per strip) + unordered_map>>> m_StripCoeffs; //! Reference energy of the depth calibration coefficients double m_Coeffs_Energy; diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index b79fda4b..bb704bb0 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -30,7 +30,6 @@ // ROOT libs: #include "TRandom.h" -#include "TMath.h" // MEGAlib libs: #include "MModuleDepthCalibration.h" @@ -75,7 +74,9 @@ bool MSubModuleDepthReadout::Initialize() m_CTDMap.clear(); m_ElectronDriftTimes.clear(); m_HoleDriftTimes.clear(); - m_Coeffs.clear(); + m_StripCoeffs.clear(); + + m_Coeffs_Energy = 59.5; // Load depth-related files using the parsers in MModuleDepthCalibration MModuleDepthCalibration DepthCalibration; @@ -111,21 +112,10 @@ bool MSubModuleDepthReadout::Initialize() return false; } - // Load depth calibration coefficients - DepthCalibration.SetCoeffsFileName(m_DepthCoefficientsFileName); - if (DepthCalibration.LoadCoeffsFile(m_DepthCoefficientsFileName) == true) { - // Copy depth calibration coefficients - m_Coeffs = DepthCalibration.GetCoeffs(); - m_Coeffs_Energy = DepthCalibration.GetCoeffsEnergy(); - - // The reference energy for the timing noise should be in the file header of the depth coefficients file - // Throw a warning if it was not retrieved and m_Coeffs_Energy is still at its default value of 0 - if (m_ApplyTimingResolutionCalibration == true && m_Coeffs_Energy == 0) { - if (g_Verbosity >= c_Warning) { - cout << "Timing values will not be smeared as no reference energy found in depth spline file "<= c_Error) { + cout << "ERROR in MSubModuleDepthReadout::Initialize: Could not parse depth calibration coefficients file" << endl; } - } else { return false; } @@ -164,6 +154,24 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) list& LVHits = Event->GetDEEStripHitLVListReference(); for (MDEEStripHit& SH: LVHits) { int DetID = SH.m_ROE.GetDetectorID(); + double MeanStretch, MeanOffset; + if (m_MeanStretch.count(DetID) == 1){ + MeanStretch = m_MeanStretch[DetID]; + } else { + if (g_Verbosity >= c_Error) { + cout << "Detector " << DetID << " does not have a mean stretch defined" << endl; + } + return false; + } + if (m_MeanOffset.count(DetID) == 1){ + MeanOffset = m_MeanOffset[DetID]; + } else { + if (g_Verbosity >= c_Error) { + cout << "Detector " << DetID << " does not have a mean offset defined" << endl; + } + return false; + } + double Z = SH.m_SimulatedPositionInDetector.Z(); if (m_DepthGrid.count(DetID) == 1) { if (SH.m_IsGuardRing == false) { @@ -176,23 +184,21 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) // Apply stretch and offset based on Eq. (3) in https://doi.org/10.1016/j.nima.2026.171332 unsigned int StripID = SH.m_ROE.GetStripID(); - int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID; - if (m_Coeffs.count(PixelCode) == 1){ - vector Coeffs = m_Coeffs[PixelCode]; - double Stretch = Coeffs[0]; - double Offset = Coeffs[1]; - double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; - double CTD_Sigma = CTD_FWHM / 2.355; + if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][0].count(StripID) == 1){ + + vector Coeffs = m_StripCoeffs[DetID][0][StripID]; + double Stretch = MeanStretch * Coeffs[0]; + double Offset = MeanOffset + Coeffs[1]; + double TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; double HoleDriftTime = (HoleSpline.Eval(Z) + Offset) * Stretch; // Convert drift time to timing by subtracting 4200ns (for now) // TODO: Improve determining the timing from drift times SH.m_Timing = 4200.0 - HoleDriftTime; - // Smear the timing value based on the given CTD resolution - // --> divide by √2 to obtain TAC resolution from CTD resolution + // Smear the timing value based on the given strip TAC resolution if (m_ApplyTimingResolutionCalibration == true) { - SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0)); + SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma); } // Apply the inverse TAC cal to obtain TAC in ADC units @@ -219,7 +225,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) } } else { if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<& HVHits = Event->GetDEEStripHitHVListReference(); for (MDEEStripHit& SH: HVHits) { int DetID = SH.m_ROE.GetDetectorID(); + double MeanStretch; + if (m_MeanStretch.count(DetID) == 1){ + MeanStretch = m_MeanStretch[DetID]; + } else { + if (g_Verbosity >= c_Error) { + cout << "Detector " << DetID << " does not have a mean stretch defined" << endl; + } + return false; + } + double Z = SH.m_SimulatedPositionInDetector.Z(); if (m_DepthGrid.count(DetID) == 1) { if (SH.m_IsGuardRing == false) { @@ -247,23 +263,20 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) // Apply stretch based on Eq. (3) in https://doi.org/10.1016/j.nima.2026.171332 // Apply no offset to the electron drift time --> add it fully to the hole signal unsigned int StripID = SH.m_ROE.GetStripID(); - int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID; - if (m_Coeffs.count(PixelCode) == 1){ + if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][1].count(StripID) == 1){ - vector Coeffs = m_Coeffs[PixelCode]; - double Stretch = Coeffs[0]; - double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; - double CTD_Sigma = CTD_FWHM / 2.355; + vector Coeffs = m_StripCoeffs[DetID][1][StripID]; + double Stretch = MeanStretch * Coeffs[0]; + double TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; double ElectronDriftTime = ElectronSpline.Eval(Z) * Stretch; // Convert drift time to timing by subtracting 4200ns (for now) // TODO: Improve determining the timing from drift times SH.m_Timing = 4200.0 - ElectronDriftTime; - // Smear the timing value based on the given CTD resolution - // --> divide by √2 to obtain TAC resolution from CTD resolution + // Smear the timing value based on the given strip TAC resolution if (m_ApplyTimingResolutionCalibration == true) { - SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0)); + SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma); } // Apply the inverse TAC cal to obtain TAC in ADC units @@ -283,7 +296,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) } } else { if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<= c_Error) { - cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: failed to open depth splines file." << endl; - } + MFile CoeffsFile; + if (CoeffsFile.Open(FileName) == false) { + cout << "ERROR in MModuleDepthCalibration::LoadStripCoeffsFile: failed to open strip coefficients file." << endl; return false; } - // Populate these vectors for each detector - int DetID = -1; - vector DepthVector; - vector CTDVector; - vector ElectronDriftTimeVector; - vector HoleDriftTimeVector; - MString Line; - while (SplineFile.ReadLine(Line) == true) { - if (Line.Length() != 0) { - if (Line.BeginsWith("#") == true) { - - if (DetID != -1 && DepthVector.size() > 0) { - m_DepthGrid[DetID] = DepthVector; - m_CTDMap[DetID] = CTDVector; - m_ElectronDriftTimes[DetID] = ElectronDriftTimeVector; - m_HoleDriftTimes[DetID] = HoleDriftTimeVector; - } + while (CoeffsFile.ReadLine(Line) == true) { + if (Line.BeginsWith('#') == true) { + std::vector Tokens = Line.Tokenize(","); + if (Tokens.size() < 5 || Tokens[0] != "#detector_info" || Tokens[1] == "detector_id") continue; + + int DetID = Tokens[1].ToInt(); + // MString DetectorName = Tokens[2]; + // MString DepthSplineFileName = Tokens[3]; + m_MeanStretch[DetID] = Tokens[4].ToDouble(); + m_MeanOffset[DetID] = Tokens[5].ToDouble(); + + if (m_StripCoeffs.count(DetID) == 0) { + vector>> TempVector; + unordered_map> TempMapLV; + unordered_map> TempMapHV; + m_StripCoeffs[DetID] = TempVector; + m_StripCoeffs[DetID].push_back(TempMapLV); + m_StripCoeffs[DetID].push_back(TempMapHV); + } - // Get the Detector ID from the commented lines - vector Tokens = Line.Tokenize(" "); - DetID = Tokens[1].ToInt(); - DepthVector.clear(); - CTDVector.clear(); - ElectronDriftTimeVector.clear(); - HoleDriftTimeVector.clear(); - - } else { - vector Tokens = Line.Tokenize(","); - if (Tokens.size() >= 4) { - m_DepthGrid[DetID].push_back(Tokens[0].ToDouble()); - m_CTDMap[DetID].push_back(Tokens[1].ToDouble()); - m_ElectronDriftTimes[DetID].push_back(Tokens[2].ToDouble()); - m_HoleDriftTimes[DetID].push_back(Tokens[3].ToDouble()); - } else { - if (g_Verbosity >= c_Error) { - cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: Empty line in depth splines file." << endl; - } - return false; - } + if (g_Verbosity >= c_Info) { + cout << m_Name << ": Detector " << DetID << " has a mean stretch of " << m_MeanStretch[DetID] << " and a mean offset of " << m_MeanOffset[DetID] << endl; + } + } else { + std::vector Tokens = Line.Tokenize(","); + if (Tokens.size() == 7) { + // unsigned int ReadOutID = Tokens[0].ToUnsignedInt(); + unsigned int DetID = Tokens[1].ToUnsignedInt(); + unsigned int Side = Tokens[2].ToUnsignedInt(); + int StripID = Tokens[3].ToInt(); + double Stretch = Tokens[4].ToDouble(); + double Offset = Tokens[5].ToDouble(); + double TimingNoiseSigma = Tokens[6].ToDouble(); + + vector coeffs; + coeffs.push_back(Stretch); coeffs.push_back(Offset); coeffs.push_back(TimingNoiseSigma); + m_StripCoeffs[DetID][Side][StripID] = coeffs; } } } - if (DetID == -1 || m_DepthGrid.size() == 0) { - if (g_Verbosity >= c_Error) { - cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: No depth splines recovered from the file." << endl; - } - return false; - } + CoeffsFile.Close(); - SplineFile.Close(); return true; } @@ -384,6 +384,9 @@ void MSubModuleDepthReadout::Finalize() m_ElectronDriftTimes.clear(); m_HoleDriftTimes.clear(); + m_StripCoeffs.clear(); + m_TACCal.clear(); + MSubModule::Finalize(); }