diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 0926a61f..dbd9d378 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -91,14 +91,31 @@ class MModuleDepthCalibration : public MModule //! Load the detector and strip dimensions from the geometry object bool LoadDetectorDimensions(MDGeometryQuest* Geometry); - //! Load in the specified coefficients file + //! Load the pixel-based coefficients file bool LoadCoeffsFile(MString FName); + //! Load the strip-based coefficients file + bool LoadStripCoeffsFile(MString FName); - //! Set the depth calibration coefficients + //! Set the pixel-based depth calibration coefficients void SetCoeffs( unordered_map> Coeffs ) { m_Coeffs = Coeffs; } - //! Get the depth calibration coefficients + //! Get the pixel-based depth calibration coefficients unordered_map> GetCoeffs() { return m_Coeffs; } + //! Set the strip-based depth calibration coefficients + void SetStripCoeffs(unordered_map>>> Coeffs) { m_StripCoeffs = Coeffs; } + //! Get the strip-based depth calibration coefficients + unordered_map>>> GetStripCoeffs() { return m_StripCoeffs; } + + //! Set the mean depth calibration offsets for all detectors + void SetMeanOffset(unordered_map MeanOffset) { m_MeanOffset = MeanOffset; } + //! Get the mean depth calibration offsets for all detectors + unordered_map GetMeanOffset() { return m_MeanOffset; } + + //! Set the mean depth calibration stretches for all detectors + void SetMeanStretch(unordered_map MeanStretch) { m_MeanStretch = MeanStretch; } + //! Get the mean depth calibration stretches for all detectors + unordered_map GetMeanStretch() { return m_MeanStretch; } + //! Set the energy at which the depth calibration coefficients were determined void SetCoeffsEnergy( double Coeffs_Energy ) { m_Coeffs_Energy = Coeffs_Energy; } //! Get the energy at which the depth calibration coefficients were determined @@ -188,6 +205,14 @@ class MModuleDepthCalibration : public MModule MModuleEnergyCalibration* m_EnergyCalibration; MGUIExpoDepthCalibration* m_ExpoDepthCalibration; + + //! 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; + // The CTD Map maps each detector (int) to a 2D array of CTD values. unordered_map>> m_CTDMap; unordered_map> m_DepthGrid; diff --git a/include/MSubModuleChargeTransport.h b/include/MSubModuleChargeTransport.h index aaf1d54b..93c48922 100644 --- a/include/MSubModuleChargeTransport.h +++ b/include/MSubModuleChargeTransport.h @@ -111,11 +111,15 @@ class MSubModuleChargeTransport : public MSubModule list m_ChargeTransportHits; - //! Filename of the depth calibration coefficients (stretch, offset, timing noise, ...) + //! Filename of the strip-based 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; //! Filename of CTD->Depth splines MString m_DepthSplinesFileName; diff --git a/include/MSubModuleDepthReadout.h b/include/MSubModuleDepthReadout.h index 3f03b421..e3e9d5bf 100644 --- a/include/MSubModuleDepthReadout.h +++ b/include/MSubModuleDepthReadout.h @@ -96,8 +96,9 @@ 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_StripCoeffs; + //! Reference energy of the depth calibration coefficients double m_Coeffs_Energy; diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index adcb9044..bccfc628 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -538,6 +538,67 @@ bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName) ///////////////////////////////////////////////////////////////////////////////// +bool MModuleDepthCalibration::LoadStripCoeffsFile(MString FileName) +{ + + MFile CoeffsFile; + if (CoeffsFile.Open(FileName) == false) { + cout << "ERROR in MModuleDepthCalibration::LoadStripCoeffsFile: failed to open strip coefficients file." << endl; + return false; + } + + MString Line; + 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); + } + + if (g_Verbosity >= c_Info) { + cout << "ERROR in MModuleDepthCalibration::LoadStripCoeffsFile: 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; + } + } + } + + CoeffsFile.Close(); + + return true; +} + + +///////////////////////////////////////////////////////////////////////////////// + + + std::vector* MModuleDepthCalibration::GetPixelCoeffs(int PixelCode) { // Check to see if the stretch and offset have been loaded. If so, try to get the coefficients for the specified pixel. diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index 7b410ce2..7cb38673 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -139,7 +139,7 @@ bool MSubModuleChargeTransport::Initialize() return false; } - m_Coeffs.clear(); + m_StripCoeffs.clear(); m_DepthGrid.clear(); m_ElectronDriftTimes.clear(); m_HoleDriftTimes.clear(); @@ -182,10 +182,11 @@ bool MSubModuleChargeTransport::Initialize() } // Load depth calibration coefficients - DepthCalibration.SetCoeffsFileName(m_DepthCoefficientsFileName); - if (DepthCalibration.LoadCoeffsFile(m_DepthCoefficientsFileName) == true) { + if (DepthCalibration.LoadStripCoeffsFile(m_DepthCoefficientsFileName) == true) { // Copy depth calibration coefficients - m_Coeffs = DepthCalibration.GetCoeffs(); + m_StripCoeffs = DepthCalibration.GetStripCoeffs(); + m_MeanStretch = DepthCalibration.GetMeanStretch(); + m_MeanOffset = DepthCalibration.GetMeanOffset(); } else { return false; @@ -327,7 +328,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // Calculate strip ID by rounding down intentionally to avoid truncation towards zero // TODO: Include mask metrology information when calculating the strip ID from the position. - int ID = static_cast(std::floor((P + PWidth/2.0) / PPitch)); + int StripID = static_cast(std::floor((P + PWidth/2.0) / PPitch)); // Calculate the strip ID for the opposite side of the detector (and explicitly check for guard ring) int OppositeStripID = static_cast(std::floor((Q + QWidth/2.0) / QPitch)); @@ -337,26 +338,41 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // Check for strip ID and if the position is within the allowed strip length or on the guard ring // TODO: Confirm the correct boundary of the guard ring based on SMEX detector models - if (ID >= 0 && ID < NStrips && std::abs(Q) <= QWidth/2.0 && std::hypot(P, Q) <= Radius) { + if (StripID >= 0 && StripID < NStrips && std::abs(Q) <= QWidth/2.0 && std::hypot(P, Q) <= Radius) { // Determine the charge drift times in nanoseconds from simulations + stretch/offset from the depth calibration // Set the default to a large number (here: 1e10 ns) in case no depth calibration coefficients exist double FastPeakTime = 1e10; TSpline3* DriftTimeSpline = isLV ? m_HoleDriftSplines[DetID] : m_ElectronDriftSplines[DetID]; - int PixelCode = 10000*DetID + 100*(isLV ? ID : OppositeStripID) + (isLV ? OppositeStripID : ID); + + double MeanStretch = 1.0; + double MeanOffset = 0.0; + 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; + } + } + 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; + } + } // 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 (LV) signal - auto it = m_Coeffs.find(PixelCode); - if (it != m_Coeffs.end()) { - const vector& Coeffs = it->second; - double Stretch = Coeffs[0]; - double Offset = isLV ? Coeffs[1] : 0.0; + // Apply the mean offset fully to the hole (LV) signal + if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][isLV ? 0 : 1].count(StripID) == 1) { + vector Coeffs = m_StripCoeffs[DetID][isLV ? 0 : 1][StripID]; + double Stretch = MeanStretch * Coeffs[0]; + double Offset = (isLV ? MeanOffset + Coeffs[1] : -Coeffs[1]); FastPeakTime = (DriftTimeSpline->Eval(Z) + Offset) * Stretch; } else { if (g_Verbosity >= c_Warning) { - cout << "No depth calibration coefficients for pixel in DetID " << DetID << " HV " << (isLV ? OppositeStripID : ID) << " LV " << (isLV ? ID : OppositeStripID) << endl; + cout << "No depth calibration coefficients for " << (isLV ? "LV" : "HV") << " strip " << StripID << endl; } } @@ -385,7 +401,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // create entry for the main hit MDEEStripHit MainSH = SH; - MainSH.m_ROE.SetStripID(ID); + MainSH.m_ROE.SetStripID(StripID); MainSH.m_OppositeStripID = OppositeStripID; MainSH.m_Energy = MainStripEnergy; // TODO: Implement a more realistic parameterization to determine nearest-neighbor timing values @@ -398,8 +414,8 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool NNLeftSH.m_Energy = std::max(NNLeftStripEnergy, 0.0); NNLeftSH.m_FastPeakTime = FastPeakTime - 50 * (1 - NNLeftStripEnergy / SH.m_SimulatedEnergy); NNLeftSH.m_OppositeStripID = OppositeStripID; - if (ID > 0) { - NNLeftSH.m_ROE.SetStripID(ID - 1); + if (StripID > 0) { + NNLeftSH.m_ROE.SetStripID(StripID - 1); NNLeftSH.m_IsGuardRing = false; } else { NNLeftSH.m_ROE.SetStripID(NStrips); @@ -412,8 +428,8 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool NNRightSH.m_Energy = std::max(NNRightStripEnergy, 0.0); NNRightSH.m_FastPeakTime = FastPeakTime - 50 * (1 - NNRightStripEnergy / SH.m_SimulatedEnergy); NNRightSH.m_OppositeStripID = OppositeStripID; - if (ID < NStrips - 1) { - NNRightSH.m_ROE.SetStripID(ID + 1); + if (StripID < NStrips - 1) { + NNRightSH.m_ROE.SetStripID(StripID + 1); NNRightSH.m_IsGuardRing = false; } else { NNRightSH.m_ROE.SetStripID(NStrips); @@ -438,7 +454,10 @@ void MSubModuleChargeTransport::Finalize() { // Finalize the analysis - do all cleanup, i.e., undo Initialize() - m_Coeffs.clear(); + m_StripCoeffs.clear(); + m_MeanStretch.clear(); + m_MeanOffset.clear(); + m_DepthGrid.clear(); m_ElectronDriftTimes.clear(); m_HoleDriftTimes.clear(); diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index 9ee3efc5..16ef34d5 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -73,17 +73,16 @@ MSubModuleDepthReadout::~MSubModuleDepthReadout() bool MSubModuleDepthReadout::Initialize() { - m_Coeffs.clear(); + m_StripCoeffs.clear(); // Load depth-related files using the parsers in MModuleDepthCalibration MModuleDepthCalibration DepthCalibration; DepthCalibration.SetUCSDOverride(false); - // Load depth calibration coefficients - DepthCalibration.SetCoeffsFileName(m_DepthCoefficientsFileName); - if (DepthCalibration.LoadCoeffsFile(m_DepthCoefficientsFileName) == true) { - // Copy depth calibration coefficients - m_Coeffs = DepthCalibration.GetCoeffs(); + // Load the strip-based depth calibration coefficients + if (DepthCalibration.LoadStripCoeffsFile(m_DepthCoefficientsFileName) == true) { + // Copy strip-based depth calibration coefficients + m_StripCoeffs = DepthCalibration.GetStripCoeffs(); m_Coeffs_Energy = DepthCalibration.GetCoeffsEnergy(); // The reference energy for the timing noise should be in the file header of the depth coefficients file @@ -146,17 +145,13 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) SH.m_HasFastTiming = true; if (m_ApplyTimingResolutionCalibration == true){ - int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID; - if (m_Coeffs.count(PixelCode) == 1){ - vector Coeffs = m_Coeffs[PixelCode]; - double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; - double CTD_Sigma = CTD_FWHM / 2.355; - // Smear the timing value based on the given CTD resolution - // --> divide by √2 to obtain TAC resolution from CTD resolution - SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0)); + 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 TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; + SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma); } else { if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "< Coeffs = m_Coeffs[PixelCode]; - double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; - double CTD_Sigma = CTD_FWHM / 2.355; - // Smear the timing value based on the given CTD resolution - // --> divide by √2 to obtain TAC resolution from CTD resolution - SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0)); - + if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][1].count(StripID) == 1) { + vector Coeffs = m_StripCoeffs[DetID][1][StripID]; + double TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy; + SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma); } else { if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<