From 9407f5bd8d45ed3c51ba74a8c025c8acf1317637 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 15 Jun 2026 12:51:03 -0700 Subject: [PATCH 1/8] Add `m_DriftTime` to `MDEEStripHit` --- include/MDEEStripHit.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/MDEEStripHit.h b/include/MDEEStripHit.h index 5deea325..1b6feb32 100644 --- a/include/MDEEStripHit.h +++ b/include/MDEEStripHit.h @@ -88,6 +88,8 @@ struct MDEEStripHit unsigned int m_TAC; //! The calibrated timing in ns double m_Timing; + //! The drift time in ns since creation of event + double m_DriftTime; //! The measured temperature value double m_Temperature; From fc5903bfdf4ff810f59ec04792bdefc9eefe66a5 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 15 Jun 2026 13:07:54 -0700 Subject: [PATCH 2/8] Move calculating drift times from `MSubModuleDepthReadout` to `MSubModuleChargeTransport` --- include/MModuleDEESMEX.h | 26 +-- include/MSubModuleChargeTransport.h | 31 +++ include/MSubModuleDepthReadout.h | 28 --- src/MModuleDEESMEX.cxx | 20 +- src/MSubModuleChargeTransport.cxx | 118 +++++++++-- src/MSubModuleDepthReadout.cxx | 292 +++++++--------------------- 6 files changed, 227 insertions(+), 288 deletions(-) diff --git a/include/MModuleDEESMEX.h b/include/MModuleDEESMEX.h index 421c026c..17ab2792 100644 --- a/include/MModuleDEESMEX.h +++ b/include/MModuleDEESMEX.h @@ -105,26 +105,14 @@ class MModuleDEESMEX : public MModule } //! Set depth coefficients file name - void SetDepthCoefficientsFileName(const MString& FileName) - { - m_DepthReadout.SetDepthCoefficientsFileName(FileName); - } + void SetDepthCoefficientsFileName(const MString& FileName) { m_DepthCoefficientsFileName = FileName; } //! Get depth coefficients file name - MString GetDepthCoefficientsFileName() const - { - return m_DepthReadout.GetDepthCoefficientsFileName(); - } + MString GetDepthCoefficientsFileName() const { return m_DepthCoefficientsFileName; } //! Set depth splines file name - void SetDepthSplinesFileName(const MString& FileName) - { - m_DepthReadout.SetDepthSplinesFileName(FileName); - } + void SetDepthSplinesFileName(const MString& FileName) { m_DepthSplinesFileName = FileName; } //! Get depth splines file name - MString GetDepthSplinesFileName() const - { - return m_DepthReadout.GetDepthSplinesFileName(); - } + MString GetDepthSplinesFileName() const { return m_DepthSplinesFileName;} //! Set TAC calibration file name void SetTACCalFileName(const MString& FileName) @@ -223,6 +211,12 @@ class MModuleDEESMEX : public MModule //! The sub module handling the output of the DEE in to the standard nuclearizer classes MSubModuleDEEOutput m_Output; + //! The file name of the simulated CTD and drift time splines + MString m_DepthSplinesFileName; + + //! The file name of the depth-calibration coefficients + MString m_DepthCoefficientsFileName; + //! Option to add noise to the strip energies bool m_ApplyResolutionCalibration; //! Option to enable shield veto effects diff --git a/include/MSubModuleChargeTransport.h b/include/MSubModuleChargeTransport.h index 7d277a47..40001c7c 100644 --- a/include/MSubModuleChargeTransport.h +++ b/include/MSubModuleChargeTransport.h @@ -19,6 +19,7 @@ // Standard libs: // ROOT libs: +#include "TSpline.h" // MEGAlib libs: #include "MGlobal.h" @@ -48,6 +49,15 @@ class MSubModuleChargeTransport : public MSubModule //! No move operators MSubModuleChargeTransport& operator=(MSubModuleChargeTransport&&) = delete; + //! Set filename for CTD->Depth splines + void SetDepthSplinesFileName( const MString& FileName) { m_DepthSplinesFileName = FileName; } + //! Get filename for CTD->Depth splines + MString GetDepthSplinesFileName() const {return m_DepthSplinesFileName; } + //! Set filename for coefficients file + void SetDepthCoefficientsFileName( const MString& FileName) { m_DepthCoefficientsFileName = FileName; } + //! Get filename for coefficients file + MString GetDepthCoefficientsFileName() const { return m_DepthCoefficientsFileName; } + //! Default destructor virtual ~MSubModuleChargeTransport(); @@ -101,6 +111,27 @@ class MSubModuleChargeTransport : public MSubModule list m_ChargeTransportHits; + //! Filename of the depth calibration coefficients (stretch, offset, timing noise, ...) + MString m_DepthCoefficientsFileName; + + //! Map of the depth calibration coefficients + unordered_map> m_Coeffs; + + //! Filename of CTD->Depth splines + MString m_DepthSplinesFileName; + + // Analog of the CTD-to-depth splines in MModuleDepthCalibration: + //! Map: detector ID (int) -> vector containing depth values + unordered_map> m_DepthGrid; + //! Map: detector ID (int) -> simulated electron drift times (+ electronics) for the depth values in m_DepthGrid + unordered_map> m_ElectronDriftTimes; + //! Map: detector ID (int) -> simulated hole drift times (+ electronics) for the depth values in m_DepthGrid + unordered_map> m_HoleDriftTimes; + + //! The corresponding interpolation splines + std::unordered_map m_HoleDriftSplines; + std::unordered_map m_ElectronDriftSplines; + // private members: private: diff --git a/include/MSubModuleDepthReadout.h b/include/MSubModuleDepthReadout.h index d894eaba..3f03b421 100644 --- a/include/MSubModuleDepthReadout.h +++ b/include/MSubModuleDepthReadout.h @@ -49,9 +49,6 @@ class MSubModuleDepthReadout : public MSubModule //! Default destructor virtual ~MSubModuleDepthReadout(); - //! Set geometry - void SetGeometry(MDGeometryQuest* Geometry) { m_Geometry = Geometry; } - //! Initialize the module virtual bool Initialize(); @@ -66,11 +63,6 @@ class MSubModuleDepthReadout : public MSubModule //! Get filename for coefficients file MString GetDepthCoefficientsFileName() const { return m_DepthCoefficientsFileName; } - //! Set filename for CTD->Depth splines - void SetDepthSplinesFileName( const MString& FileName) { m_DepthSplinesFileName = FileName; } - //! Get filename for CTD->Depth splines - MString GetDepthSplinesFileName() const {return m_DepthSplinesFileName; } - //! Set filename for TAC calibration void SetTACCalFileName( const MString& FileName) { m_TACCalFileName = FileName; } //! Get filename for TAC calibration @@ -92,9 +84,6 @@ class MSubModuleDepthReadout : public MSubModule // protected methods: protected: - //! Load the splines file - bool LoadSplinesFile(MString FileName); - // private methods: private: @@ -104,31 +93,14 @@ class MSubModuleDepthReadout : public MSubModule // protected members: protected: - //! The geometry - MDGeometryQuest* m_Geometry; - //! Filename of the depth calibration coefficients (stretch, offset, timing noise, ...) MString m_DepthCoefficientsFileName; //! Map of the depth calibration coefficients unordered_map> m_Coeffs; - //! Reference energy of the depth calibration coefficients double m_Coeffs_Energy; - //! Filename of CTD->Depth splines - MString m_DepthSplinesFileName; - - // Analog of the CTD-to-depth splines in MModuleDepthCalibration: - //! Map: detector ID (int) -> vector containing depth values - unordered_map> m_DepthGrid; - //! Map: detector ID (int) -> simulated CTD values for the depth values in m_DepthGrid - unordered_map> m_CTDMap; - //! Map: detector ID (int) -> simulated electron drift times (+ electronics) for the depth values in m_DepthGrid - unordered_map> m_ElectronDriftTimes; - //! Map: detector ID (int) -> simulated hole drift times (+ electronics) for the depth values in m_DepthGrid - unordered_map> m_HoleDriftTimes; - //! Filename of the TAC calibration file MString m_TACCalFileName; diff --git a/src/MModuleDEESMEX.cxx b/src/MModuleDEESMEX.cxx index 394387eb..2c533039 100644 --- a/src/MModuleDEESMEX.cxx +++ b/src/MModuleDEESMEX.cxx @@ -111,12 +111,14 @@ bool MModuleDEESMEX::Initialize() // Set the geometry to the SubModules using it m_ChargeTransport.SetGeometry(m_Geometry); - m_DepthReadout.SetGeometry(m_Geometry); m_StripReadout.SetApplyResolutionCalibration(m_ApplyResolutionCalibration); m_DepthReadout.SetApplyTimingResolutionCalibration(m_ApplyTimingResolutionCalibration); - // Initialize the module + // Pass the depth-calibration-related files to the SubModules using it + m_ChargeTransport.SetDepthSplinesFileName(m_DepthSplinesFileName); + m_ChargeTransport.SetDepthCoefficientsFileName(m_DepthCoefficientsFileName); + m_DepthReadout.SetDepthCoefficientsFileName(m_DepthCoefficientsFileName); // Each Initialize() should handle its own error messaging if (m_Intake.Initialize() == false) return false; @@ -294,6 +296,16 @@ bool MModuleDEESMEX::ReadXmlConfiguration(MXmlNode* Node) m_StripTrigger.ReadXmlConfiguration(Node); m_DepthReadout.ReadXmlConfiguration(Node); m_Output.ReadXmlConfiguration(Node); + + // Add depth-calibration-related file names (used by several submodules) + MXmlNode* DepthSplineFile = Node->GetNode("DepthSplineFileName"); + if (DepthSplineFile != nullptr) { + m_DepthSplinesFileName = DepthSplineFile->GetValue(); + } + MXmlNode* DepthCoefficientsFileName = Node->GetNode("DepthCoefficientsFileName"); + if (DepthCoefficientsFileName != nullptr) { + m_DepthCoefficientsFileName = DepthCoefficientsFileName->GetValue(); + } // Add noise button for energies MXmlNode* ResolutionCalibrationNode = Node->GetNode("ApplyResolutionCalibration"); @@ -339,6 +351,10 @@ MXmlNode* MModuleDEESMEX::CreateXmlConfiguration() m_DepthReadout.CreateXmlConfiguration(Node); m_Output.CreateXmlConfiguration(Node); + // Add depth-calibration-related file names (used by several submodules) + new MXmlNode(Node, "DepthSplineFileName", m_DepthSplinesFileName); + new MXmlNode(Node, "DepthCoefficientsFileName", m_DepthCoefficientsFileName); + // Add noise button for energies new MXmlNode(Node, "ApplyResolutionCalibration", m_ApplyResolutionCalibration); // Add shield veto effects button diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index 6760fac5..ed10034d 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -34,6 +34,7 @@ // MEGAlib libs: #include "MSubModule.h" +#include "MModuleDepthCalibration.h" #include "MDShapeIntersection.h" #include "MDShapeTUBS.h" @@ -136,6 +137,58 @@ bool MSubModuleChargeTransport::Initialize() return false; } + m_Coeffs.clear(); + m_DepthGrid.clear(); + m_ElectronDriftTimes.clear(); + m_HoleDriftTimes.clear(); + + // Load depth-related files using the parsers in MModuleDepthCalibration + MModuleDepthCalibration DepthCalibration; + DepthCalibration.SetUCSDOverride(false); + + // Determine the thicknesses of the individual detectors from the geometry + if (DepthCalibration.LoadDetectorDimensions(m_Geometry) == false) { + return false; + } + + // Load CTD-to-depth splines + DepthCalibration.SetSplinesFileName(m_DepthSplinesFileName); + if (DepthCalibration.LoadSplinesFile(m_DepthSplinesFileName) == true) { + + // There should be at least three more columns with 1. the CTDmap, 2. the electron drift times and 3. the hole drift times + unordered_map>> Columns = DepthCalibration.GetCTDMap(); + + // Copy the depth grid, CTD splines and the charge carrier drift times + m_DepthGrid = DepthCalibration.GetDepthGrid(); + for (auto const& [DetID, Column] : Columns) { + if (Column.size() < 3) { + if (g_Verbosity >= c_Error) { + cout<<"ERROR in MSubModuleChargeTransport::Initialize: Expected (at least) 4 columns for detector "<m_ROE == IterLV2->m_ROE) { IterLV1->m_Energy += IterLV2->m_Energy; + IterLV1->m_DriftTime = IterLV1->m_Energy > IterLV2->m_Energy ? IterLV1->m_DriftTime : IterLV2->m_DriftTime; IterLV2 = LVHits.erase(IterLV2); } else { ++IterLV2; @@ -206,6 +260,7 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) while (IterHV2 != HVHits.end()) { if (IterHV1->m_ROE == IterHV2->m_ROE) { IterHV1->m_Energy += IterHV2->m_Energy; + IterHV1->m_DriftTime = IterHV1->m_Energy > IterHV2->m_Energy ? IterHV1->m_DriftTime : IterHV2->m_DriftTime; IterHV2 = HVHits.erase(IterHV2); } else { ++IterHV2; @@ -262,7 +317,8 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool MVector Pos = SH.m_SimulatedPositionInDetector; double P = isLV ? Pos.X() : Pos.Y(); double Q = isLV ? Pos.Y() : Pos.X(); - double ΔZ = isLV ? Pos.Z() + Thickness / 2.0 : Thickness / 2.0 - Pos.Z(); + double Z = Pos.Z(); + double DeltaZ = isLV ? Z + Thickness / 2.0 : Thickness / 2.0 - Z; double MeanElectricField = BiasVoltage / Thickness; // unit: V/cm @@ -280,13 +336,35 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // 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) { + // Determine the charge drift times in nanoseconds from simulations + stretch/offset from the depth calibration + // Set the default to 0ns in case no depth calibration coefficients exist + double DriftTime = 1e10; + + TSpline3* DriftTimeSpline = isLV ? m_HoleDriftSplines[DetID] : m_ElectronDriftSplines[DetID]; + int PixelCode = 10000*DetID + 100*(isLV ? ID : OppositeStripID) + (isLV ? OppositeStripID : ID); + + // 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; + DriftTime = (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; + } + DriftTime = 1e10; + } + // Apply charge sharing based on relative coordinate to the gap of that strip (0 <= X < XPitch) double FromGap = std::fmod(P + PWidth/2.0, PPitch); // Charge transport based on Eq. (7) in https://doi.org/10.1016/j.nima.2023.168310 // calculate σ and η, assuming t = z / v = z / (µ * E) - double Sigma = std::sqrt(2.0 * kB * Temperature * ΔZ / (ElementaryCharge * MeanElectricField)); // in cm - double Eta = std::cbrt(std::pow(InitialChargeCloudSize, 3) + 3.0 * N * ElementaryCharge * ΔZ / (4.0 * TMath::Pi() * Epsilon0 * EpsilonR * MeanElectricField)); // in cm + double Sigma = std::sqrt(2.0 * kB * Temperature * DeltaZ / (ElementaryCharge * MeanElectricField)); // in cm + double Eta = std::cbrt(std::pow(InitialChargeCloudSize, 3) + 3.0 * N * ElementaryCharge * DeltaZ / (4.0 * TMath::Pi() * Epsilon0 * EpsilonR * MeanElectricField)); // in cm auto Lambda = [&](double x) -> double { double a = (x - Eta) / (TMath::Sqrt2() * Sigma); double b = (x + Eta) / (TMath::Sqrt2() * Sigma); @@ -307,6 +385,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool MainSH.m_ROE.SetStripID(ID); MainSH.m_OppositeStripID = OppositeStripID; MainSH.m_Energy = MainStripEnergy; + MainSH.m_DriftTime = DriftTime - 50 * (1 - MainStripEnergy / SH.m_SimulatedEnergy); MainSH.m_IsGuardRing = false; m_ChargeTransportHits.push_back(MainSH); @@ -314,11 +393,11 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool if (NNLeftStripEnergy > IonizationEnergy) { MDEEStripHit NNLeftSH = SH; NNLeftSH.m_Energy = NNLeftStripEnergy; + NNLeftSH.m_DriftTime = DriftTime - 50 * (1 - NNLeftStripEnergy / SH.m_SimulatedEnergy); NNLeftSH.m_OppositeStripID = OppositeStripID; if (ID > 0) { NNLeftSH.m_ROE.SetStripID(ID - 1); NNLeftSH.m_IsGuardRing = false; - // NNLeftSH.m_IsNearestNeighbor = true; } else { NNLeftSH.m_ROE.SetStripID(NStrips); NNLeftSH.m_IsGuardRing = true; @@ -330,11 +409,11 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool if (NNRightStripEnergy > IonizationEnergy) { MDEEStripHit NNRightSH = SH; NNRightSH.m_Energy = NNRightStripEnergy; + NNRightSH.m_DriftTime = DriftTime - 50 * (1 - NNRightStripEnergy / SH.m_SimulatedEnergy); NNRightSH.m_OppositeStripID = OppositeStripID; if (ID < NStrips - 1) { NNRightSH.m_ROE.SetStripID(ID + 1); NNRightSH.m_IsGuardRing = false; - // NNRightSH.m_IsNearestNeighbor = true; } else { NNRightSH.m_ROE.SetStripID(NStrips); NNRightSH.m_IsGuardRing = true; @@ -359,6 +438,21 @@ void MSubModuleChargeTransport::Finalize() { // Finalize the analysis - do all cleanup, i.e., undo Initialize() + m_Coeffs.clear(); + m_DepthGrid.clear(); + m_ElectronDriftTimes.clear(); + m_HoleDriftTimes.clear(); + + // Delete the allocated spline objects + for (auto& [DetID, spline] : m_ElectronDriftSplines) { + delete spline; + } + for (auto& [DetID, spline] : m_HoleDriftSplines) { + delete spline; + } + m_ElectronDriftSplines.clear(); + m_HoleDriftSplines.clear(); + MSubModule::Finalize(); } @@ -369,14 +463,6 @@ void MSubModuleChargeTransport::Finalize() bool MSubModuleChargeTransport::ReadXmlConfiguration(MXmlNode* Node) { //! Read the configuration data from an XML node - - /* - MXmlNode* SomeTagNode = Node->GetNode("SomeTag"); - if (SomeTagNode != 0) { - m_SomeTagValue = SomeTagNode->GetValue(); - } - */ - return true; } @@ -387,14 +473,8 @@ bool MSubModuleChargeTransport::ReadXmlConfiguration(MXmlNode* Node) MXmlNode* MSubModuleChargeTransport::CreateXmlConfiguration(MXmlNode* Node) { //! Create an XML node tree from the configuration - - /* - MXmlNode* SomeTagNode = new MXmlNode(Node, "SomeTag", "SomeValue"); - */ - return Node; } - // MSubModuleChargeTransport.cxx: the end... //////////////////////////////////////////////////////////////////////////////// diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index b79fda4b..6860b003 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -71,45 +71,11 @@ MSubModuleDepthReadout::~MSubModuleDepthReadout() bool MSubModuleDepthReadout::Initialize() { - m_DepthGrid.clear(); - m_CTDMap.clear(); - m_ElectronDriftTimes.clear(); - m_HoleDriftTimes.clear(); m_Coeffs.clear(); // Load depth-related files using the parsers in MModuleDepthCalibration MModuleDepthCalibration DepthCalibration; - - // Determine the thicknesses of the individual detectors from the geometry DepthCalibration.SetUCSDOverride(false); - if (DepthCalibration.LoadDetectorDimensions(m_Geometry) == false) { - return false; - } - - // Load CTD-to-depth splines - DepthCalibration.SetSplinesFileName(m_DepthSplinesFileName); - if (DepthCalibration.LoadSplinesFile(m_DepthSplinesFileName) == true) { - - // There should be at least three more columns with 1. the CTDmap, 2. the electron drift times and 3. the hole drift times - unordered_map>> Columns = DepthCalibration.GetCTDMap(); - - // Copy the depth grid, CTD splines and the charge carrier drift times - m_DepthGrid = DepthCalibration.GetDepthGrid(); - for (auto const& [DetID, Column] : Columns) { - if (Column.size() < 3) { - if (g_Verbosity >= c_Error) { - cout<<"ERROR in MSubModuleDepthReadout::Initialize: Expected (at least) 4 columns for detector "<& LVHits = Event->GetDEEStripHitLVListReference(); for (MDEEStripHit& SH: LVHits) { int DetID = SH.m_ROE.GetDetectorID(); - double Z = SH.m_SimulatedPositionInDetector.Z(); - if (m_DepthGrid.count(DetID) == 1) { - if (SH.m_IsGuardRing == false) { - - // Determine the hole drift times (accounting for electronics) - // TODO: This only applies to main strips, find a better implementation for nearest neighbors - vector DepthGrid = m_DepthGrid[DetID]; - vector HoleDriftTimes = m_HoleDriftTimes[DetID]; - TSpline3 HoleSpline = TSpline3("", &DepthGrid[0], &HoleDriftTimes[0], DepthGrid.size()); - - // 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; - 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 - if (m_ApplyTimingResolutionCalibration == true) { + if (SH.m_IsGuardRing == false) { + unsigned int StripID = SH.m_ROE.GetStripID(); + if (SH.m_DriftTime > -200.0 ) { + SH.m_Timing = 4200.0 - SH.m_DriftTime; + 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)); + } else { + if (g_Verbosity >= c_Info) { + cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<= 1 && StripID < m_TACCal[DetID][0].size()) { - vector TACCal = m_TACCal[DetID][0][StripID]; - double TACCalSlope = TACCal[0]; - double TACCalOffset = TACCal[1]; - if ((SH.m_Timing - TACCalOffset) / TACCalSlope >= 0) { - SH.m_TAC = (SH.m_Timing - TACCalOffset) / TACCalSlope; - } else { - if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: Simulated TAC value would be negative, setting it to zero."<= 1 && StripID < m_TACCal[DetID][0].size()) { + vector TACCal = m_TACCal[DetID][0][StripID]; + double TACCalSlope = TACCal[0]; + double TACCalOffset = TACCal[1]; + if ((SH.m_Timing - TACCalOffset) / TACCalSlope >= 0) { + SH.m_TAC = (SH.m_Timing - TACCalOffset) / TACCalSlope; } else { - if (g_Verbosity >= c_Error) { - cout<<"ERROR in MSubModuleDepthReadout::AnalyzeEvent: No TAC calibration found for LV strip "<= c_Info) { + cout<<"MSubModuleDepthReadout::AnalyzeEvent: Simulated TAC value would be negative, setting it to zero."< 16383) { - SH.m_TAC = 16383; - } } else { - if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<= c_Error) { + cout<<"ERROR in MSubModuleDepthReadout::AnalyzeEvent: No TAC calibration found for LV strip "<= c_Warning) { - cout << "No Depth Spline for Event with DetID " << DetID << endl; + if (SH.m_TAC > 16383) { + SH.m_TAC = 16383; + } + } else { + if (g_Verbosity >= c_Info) { + cout<<"MSubModuleDepthReadout::AnalyzeEvent: Unphysical drift time."<& HVHits = Event->GetDEEStripHitHVListReference(); for (MDEEStripHit& SH: HVHits) { int DetID = SH.m_ROE.GetDetectorID(); - double Z = SH.m_SimulatedPositionInDetector.Z(); - if (m_DepthGrid.count(DetID) == 1) { - if (SH.m_IsGuardRing == false) { - - // Determine the electron drift times (accounting for electronics) - // TODO: This only applies to main strips, find a better implementation for nearest neighbors - vector DepthGrid = m_DepthGrid[DetID]; - vector ElectronDriftTimes = m_ElectronDriftTimes[DetID]; - TSpline3 ElectronSpline = TSpline3("", &DepthGrid[0], &ElectronDriftTimes[0], DepthGrid.size()); - - // 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){ - - 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; - 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 - if (m_ApplyTimingResolutionCalibration == true) { + if (SH.m_IsGuardRing == false) { + unsigned int StripID = SH.m_ROE.GetStripID(); + if (SH.m_DriftTime > -200.0){ + SH.m_Timing = 4200.0 - SH.m_DriftTime; + if (m_ApplyTimingResolutionCalibration == true) { + int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID; + 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)); - } - - // Apply the inverse TAC cal to obtain TAC in ADC units - if (m_TACCal.count(DetID) == 1 && m_TACCal[DetID].size() >= 2 && StripID < m_TACCal[DetID][1].size()) { - vector TACCal = m_TACCal[DetID][1][StripID]; - double TACCalSlope = TACCal[0]; - double TACCalOffset = TACCal[1]; - SH.m_TAC = (SH.m_Timing - TACCalOffset) / TACCalSlope; } else { - if (g_Verbosity >= c_Error) { - cout<<"ERROR in MSubModuleDepthReadout::AnalyzeEvent: No TAC calibration found for HV strip "<= c_Info) { + cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "< 16383) { - SH.m_TAC = 16383; } + } + + // Apply the inverse TAC cal to obtain TAC in ADC units + if (m_TACCal.count(DetID) == 1 && m_TACCal[DetID].size() >= 2 && StripID < m_TACCal[DetID][1].size()) { + vector TACCal = m_TACCal[DetID][1][StripID]; + double TACCalSlope = TACCal[0]; + double TACCalOffset = TACCal[1]; + SH.m_TAC = (SH.m_Timing - TACCalOffset) / TACCalSlope; } else { - if (g_Verbosity >= c_Info) { - cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<= c_Error) { + cout<<"ERROR in MSubModuleDepthReadout::AnalyzeEvent: No TAC calibration found for HV strip "<= c_Info) { - cout << "No Depth Spline for Event with DetID " << DetID << endl; - } - } - } - - return true; -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MSubModuleDepthReadout::LoadSplinesFile(MString FileName) -{ - // Input spline files should have the following format: - // ### DetID, HV, Temperature, Photopeak Energy - // depth, ctd, electron_drift_time, hole_drift_time - - MFile SplineFile; - if (SplineFile.Open(FileName) == false) { - if (g_Verbosity >= c_Error) { - cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: failed to open depth splines 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; + if (SH.m_TAC > 16383) { + SH.m_TAC = 16383; } - - // 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<<"MSubModuleDepthReadout::AnalyzeEvent: Unphysical drift time."<= c_Error) { - cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: No depth splines recovered from the file." << endl; - } - return false; - } - - SplineFile.Close(); return true; } @@ -379,10 +239,8 @@ bool MSubModuleDepthReadout::LoadSplinesFile(MString FileName) void MSubModuleDepthReadout::Finalize() { // Finalize the analysis - do all cleanup, i.e., undo Initialize() - m_DepthGrid.clear(); - m_CTDMap.clear(); - m_ElectronDriftTimes.clear(); - m_HoleDriftTimes.clear(); + m_Coeffs.clear(); + m_TACCal.clear(); MSubModule::Finalize(); } @@ -394,16 +252,6 @@ void MSubModuleDepthReadout::Finalize() bool MSubModuleDepthReadout::ReadXmlConfiguration(MXmlNode* Node) { //! Read the configuration data from an XML node - MXmlNode* DepthSplineFile = Node->GetNode("DepthSplineFileName"); - if (DepthSplineFile != nullptr) { - m_DepthSplinesFileName = DepthSplineFile->GetValue(); - } - - MXmlNode* DepthCoefficientsFileName = Node->GetNode("DepthCoefficientsFileName"); - if (DepthCoefficientsFileName != nullptr) { - m_DepthCoefficientsFileName = DepthCoefficientsFileName->GetValue(); - } - MXmlNode* TACCalFileName = Node->GetNode("TACCalFileName"); if (TACCalFileName != nullptr) { m_TACCalFileName = TACCalFileName->GetValue(); @@ -419,8 +267,6 @@ bool MSubModuleDepthReadout::ReadXmlConfiguration(MXmlNode* Node) MXmlNode* MSubModuleDepthReadout::CreateXmlConfiguration(MXmlNode* Node) { //! Create an XML node tree from the configuration - new MXmlNode(Node, "DepthSplineFileName", m_DepthSplinesFileName); - new MXmlNode(Node, "DepthCoefficientsFileName", m_DepthCoefficientsFileName); new MXmlNode(Node, "TACCalFileName", m_TACCalFileName); return Node; From bbd325ac35b957c7f6873c845c1c916bf4bdaf60 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 15 Jun 2026 13:30:07 -0700 Subject: [PATCH 3/8] Always record NN in the DEE (to keep information about NN timing) --- include/MModuleDEESMEX.h | 2 +- include/MSubModuleChargeTransport.h | 9 ++++----- src/MSubModuleChargeTransport.cxx | 16 +++++++++------- src/MSubModuleDepthReadout.cxx | 2 ++ src/MSubModuleStripReadout.cxx | 2 ++ 5 files changed, 18 insertions(+), 13 deletions(-) diff --git a/include/MModuleDEESMEX.h b/include/MModuleDEESMEX.h index 17ab2792..ab78a88e 100644 --- a/include/MModuleDEESMEX.h +++ b/include/MModuleDEESMEX.h @@ -23,7 +23,7 @@ // MEGAlib libs: #include "MGlobal.h" -// Nuclearizer libs +// Nuclearizer libs: #include "MModule.h" #include "MSubModuleDEEIntake.h" #include "MSubModuleRandomCoincidence.h" diff --git a/include/MSubModuleChargeTransport.h b/include/MSubModuleChargeTransport.h index 40001c7c..aaf1d54b 100644 --- a/include/MSubModuleChargeTransport.h +++ b/include/MSubModuleChargeTransport.h @@ -125,13 +125,12 @@ class MSubModuleChargeTransport : public MSubModule unordered_map> m_DepthGrid; //! Map: detector ID (int) -> simulated electron drift times (+ electronics) for the depth values in m_DepthGrid unordered_map> m_ElectronDriftTimes; + //! The corresponding electron-drift-time interpolation spline + unordered_map m_ElectronDriftSplines; //! Map: detector ID (int) -> simulated hole drift times (+ electronics) for the depth values in m_DepthGrid unordered_map> m_HoleDriftTimes; - - //! The corresponding interpolation splines - std::unordered_map m_HoleDriftSplines; - std::unordered_map m_ElectronDriftSplines; - + //! The corresponding hole-drift-time interpolation spline + unordered_map m_HoleDriftSplines; // private members: private: diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index ed10034d..e975e6f5 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -34,10 +34,12 @@ // MEGAlib libs: #include "MSubModule.h" -#include "MModuleDepthCalibration.h" #include "MDShapeIntersection.h" #include "MDShapeTUBS.h" +// Nuclearizer libs: +#include "MModuleDepthCalibration.h" + //////////////////////////////////////////////////////////////////////////////// @@ -390,9 +392,9 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool m_ChargeTransportHits.push_back(MainSH); // create MDEEStripHit for the left NN - if (NNLeftStripEnergy > IonizationEnergy) { + // if (NNLeftStripEnergy > IonizationEnergy) { MDEEStripHit NNLeftSH = SH; - NNLeftSH.m_Energy = NNLeftStripEnergy; + NNLeftSH.m_Energy = std::max(NNLeftStripEnergy, 0.0); NNLeftSH.m_DriftTime = DriftTime - 50 * (1 - NNLeftStripEnergy / SH.m_SimulatedEnergy); NNLeftSH.m_OppositeStripID = OppositeStripID; if (ID > 0) { @@ -403,12 +405,12 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool NNLeftSH.m_IsGuardRing = true; } m_ChargeTransportHits.push_back(NNLeftSH); - } + // } // create MDEEStripHit for the right NN - if (NNRightStripEnergy > IonizationEnergy) { + // if (NNRightStripEnergy > IonizationEnergy) { MDEEStripHit NNRightSH = SH; - NNRightSH.m_Energy = NNRightStripEnergy; + NNRightSH.m_Energy = std::max(NNRightStripEnergy, 0.0); NNRightSH.m_DriftTime = DriftTime - 50 * (1 - NNRightStripEnergy / SH.m_SimulatedEnergy); NNRightSH.m_OppositeStripID = OppositeStripID; if (ID < NStrips - 1) { @@ -419,7 +421,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool NNRightSH.m_IsGuardRing = true; } m_ChargeTransportHits.push_back(NNRightSH); - } + // } } else { // TODO: implement charge sharing also for GR events diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index 6860b003..52e02b55 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -33,6 +33,8 @@ #include "TMath.h" // MEGAlib libs: + +// Nuclearizer libs: #include "MModuleDepthCalibration.h" #include "MModuleTACcut.h" diff --git a/src/MSubModuleStripReadout.cxx b/src/MSubModuleStripReadout.cxx index a5e0aef1..6c8c089c 100644 --- a/src/MSubModuleStripReadout.cxx +++ b/src/MSubModuleStripReadout.cxx @@ -33,6 +33,8 @@ // MEGAlib libs: #include "MParser.h" + +// Nuclearizer libs: #include "MModuleEnergyCalibration.h" From 06e69195111977b315dd83c09b97ad022b90af5c Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 30 Jul 2026 11:38:37 -0700 Subject: [PATCH 4/8] Rename `m_DriftTime` to `m_FastPeakTime` --- include/MDEEStripHit.h | 4 ++-- src/MSubModuleChargeTransport.cxx | 17 ++++++++--------- src/MSubModuleDepthReadout.cxx | 8 ++++---- 3 files changed, 14 insertions(+), 15 deletions(-) diff --git a/include/MDEEStripHit.h b/include/MDEEStripHit.h index 1b6feb32..73d602f6 100644 --- a/include/MDEEStripHit.h +++ b/include/MDEEStripHit.h @@ -88,8 +88,8 @@ struct MDEEStripHit unsigned int m_TAC; //! The calibrated timing in ns double m_Timing; - //! The drift time in ns since creation of event - double m_DriftTime; + //! The fast shaper peaking time in ns since creation of event + double m_FastPeakTime; //! The measured temperature value double m_Temperature; diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index e975e6f5..83f10a23 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -250,7 +250,7 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) while (IterLV2 != LVHits.end()) { if (IterLV1->m_ROE == IterLV2->m_ROE) { IterLV1->m_Energy += IterLV2->m_Energy; - IterLV1->m_DriftTime = IterLV1->m_Energy > IterLV2->m_Energy ? IterLV1->m_DriftTime : IterLV2->m_DriftTime; + IterLV1->m_FastPeakTime = IterLV1->m_Energy > IterLV2->m_Energy ? IterLV1->m_FastPeakTime : IterLV2->m_FastPeakTime; IterLV2 = LVHits.erase(IterLV2); } else { ++IterLV2; @@ -262,7 +262,7 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) while (IterHV2 != HVHits.end()) { if (IterHV1->m_ROE == IterHV2->m_ROE) { IterHV1->m_Energy += IterHV2->m_Energy; - IterHV1->m_DriftTime = IterHV1->m_Energy > IterHV2->m_Energy ? IterHV1->m_DriftTime : IterHV2->m_DriftTime; + IterHV1->m_FastPeakTime = IterHV1->m_Energy > IterHV2->m_Energy ? IterHV1->m_FastPeakTime : IterHV2->m_FastPeakTime; IterHV2 = HVHits.erase(IterHV2); } else { ++IterHV2; @@ -339,8 +339,8 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool if (ID >= 0 && ID < 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 0ns in case no depth calibration coefficients exist - double DriftTime = 1e10; + // 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); @@ -352,12 +352,11 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool const vector& Coeffs = it->second; double Stretch = Coeffs[0]; double Offset = isLV ? Coeffs[1] : 0.0; - DriftTime = (DriftTimeSpline->Eval(Z) + Offset) * Stretch; + 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; } - DriftTime = 1e10; } // Apply charge sharing based on relative coordinate to the gap of that strip (0 <= X < XPitch) @@ -387,7 +386,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool MainSH.m_ROE.SetStripID(ID); MainSH.m_OppositeStripID = OppositeStripID; MainSH.m_Energy = MainStripEnergy; - MainSH.m_DriftTime = DriftTime - 50 * (1 - MainStripEnergy / SH.m_SimulatedEnergy); + MainSH.m_FastPeakTime = FastPeakTime - 50 * (1 - MainStripEnergy / SH.m_SimulatedEnergy); MainSH.m_IsGuardRing = false; m_ChargeTransportHits.push_back(MainSH); @@ -395,7 +394,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // if (NNLeftStripEnergy > IonizationEnergy) { MDEEStripHit NNLeftSH = SH; NNLeftSH.m_Energy = std::max(NNLeftStripEnergy, 0.0); - NNLeftSH.m_DriftTime = DriftTime - 50 * (1 - NNLeftStripEnergy / SH.m_SimulatedEnergy); + NNLeftSH.m_FastPeakTime = FastPeakTime - 50 * (1 - NNLeftStripEnergy / SH.m_SimulatedEnergy); NNLeftSH.m_OppositeStripID = OppositeStripID; if (ID > 0) { NNLeftSH.m_ROE.SetStripID(ID - 1); @@ -411,7 +410,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // if (NNRightStripEnergy > IonizationEnergy) { MDEEStripHit NNRightSH = SH; NNRightSH.m_Energy = std::max(NNRightStripEnergy, 0.0); - NNRightSH.m_DriftTime = DriftTime - 50 * (1 - NNRightStripEnergy / SH.m_SimulatedEnergy); + NNRightSH.m_FastPeakTime = FastPeakTime - 50 * (1 - NNRightStripEnergy / SH.m_SimulatedEnergy); NNRightSH.m_OppositeStripID = OppositeStripID; if (ID < NStrips - 1) { NNRightSH.m_ROE.SetStripID(ID + 1); diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index 52e02b55..0f1bdecf 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -134,8 +134,8 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) int DetID = SH.m_ROE.GetDetectorID(); if (SH.m_IsGuardRing == false) { unsigned int StripID = SH.m_ROE.GetStripID(); - if (SH.m_DriftTime > -200.0 ) { - SH.m_Timing = 4200.0 - SH.m_DriftTime; + if (SH.m_FastPeakTime > -200.0 ) { + SH.m_Timing = 4200.0 - SH.m_FastPeakTime; if (m_ApplyTimingResolutionCalibration == true){ int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID; if (m_Coeffs.count(PixelCode) == 1){ @@ -188,8 +188,8 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) int DetID = SH.m_ROE.GetDetectorID(); if (SH.m_IsGuardRing == false) { unsigned int StripID = SH.m_ROE.GetStripID(); - if (SH.m_DriftTime > -200.0){ - SH.m_Timing = 4200.0 - SH.m_DriftTime; + if (SH.m_FastPeakTime > -200.0){ + SH.m_Timing = 4200.0 - SH.m_FastPeakTime; if (m_ApplyTimingResolutionCalibration == true) { int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID; if (m_Coeffs.count(PixelCode) == 1){ From 075605c6b065e0c09d453242c6c55dfbc0675f37 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 30 Jul 2026 11:45:45 -0700 Subject: [PATCH 5/8] Add TODOs --- src/MSubModuleChargeTransport.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index 83f10a23..d5cf2faa 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -245,6 +245,7 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) // Merge hits: // TODO: how to deal with flags like "m_IsNearestNeighbor" etc. ? + // TODO: how to properly combine m_FastPeakTime, maybe merge strip hits only in MSubModuleDepthReadout ? for (auto IterLV1 = LVHits.begin(); IterLV1 != LVHits.end(); ++IterLV1) { auto IterLV2 = std::next(IterLV1); while (IterLV2 != LVHits.end()) { @@ -364,6 +365,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool // Charge transport based on Eq. (7) in https://doi.org/10.1016/j.nima.2023.168310 // calculate σ and η, assuming t = z / v = z / (µ * E) + // TODO: Reevaluate whether we want to use t = z / v = z / (µ * E), or use the simulated drift times here instead double Sigma = std::sqrt(2.0 * kB * Temperature * DeltaZ / (ElementaryCharge * MeanElectricField)); // in cm double Eta = std::cbrt(std::pow(InitialChargeCloudSize, 3) + 3.0 * N * ElementaryCharge * DeltaZ / (4.0 * TMath::Pi() * Epsilon0 * EpsilonR * MeanElectricField)); // in cm auto Lambda = [&](double x) -> double { @@ -386,6 +388,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool MainSH.m_ROE.SetStripID(ID); MainSH.m_OppositeStripID = OppositeStripID; MainSH.m_Energy = MainStripEnergy; + // TODO: Implement a more realistic parameterization to determine nearest-neighbor timing values MainSH.m_FastPeakTime = FastPeakTime - 50 * (1 - MainStripEnergy / SH.m_SimulatedEnergy); MainSH.m_IsGuardRing = false; m_ChargeTransportHits.push_back(MainSH); From 2e6fa6769de91840f3dbfc9dc626143cad669a64 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 30 Jul 2026 11:57:09 -0700 Subject: [PATCH 6/8] Add comment on `m_FastPeakTime` cut values in `MSubModuleDepthReadout` --- src/MSubModuleDepthReadout.cxx | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index 0f1bdecf..207d1b2a 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -134,7 +134,12 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) int DetID = SH.m_ROE.GetDetectorID(); if (SH.m_IsGuardRing == false) { unsigned int StripID = SH.m_ROE.GetStripID(); - if (SH.m_FastPeakTime > -200.0 ) { + // Only compute timing values for strip hits with reasonable timing values + // Note: m_FastPeakTime can be negative if the charge drift time is shorter + // than the timing offset between main strips and non-collecting adjacent strips + // Note: The upper limit of 10000.0 is chosen to remove strip hits with no depth calibration, + // where the default value is set to 1e10 + if (SH.m_FastPeakTime > -200.0 && SH.m_FastPeakTime < 10000.0) { SH.m_Timing = 4200.0 - SH.m_FastPeakTime; if (m_ApplyTimingResolutionCalibration == true){ int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID; @@ -188,7 +193,12 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) int DetID = SH.m_ROE.GetDetectorID(); if (SH.m_IsGuardRing == false) { unsigned int StripID = SH.m_ROE.GetStripID(); - if (SH.m_FastPeakTime > -200.0){ + // Only compute timing values for strip hits with reasonable timing values + // Note: m_FastPeakTime can be negative if the charge drift time is shorter + // than the timing offset between main strips and non-collecting adjacent strips + // Note: The upper limit of 10000.0 is chosen to remove strip hits with no depth calibration, + // where the default value is set to 1e10 + if (SH.m_FastPeakTime > -200.0 && SH.m_FastPeakTime < 10000.0){ SH.m_Timing = 4200.0 - SH.m_FastPeakTime; if (m_ApplyTimingResolutionCalibration == true) { int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID; From bb05fb247899ebe933ad6045cea743893e7331d3 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 30 Jul 2026 12:01:07 -0700 Subject: [PATCH 7/8] Add `m_HasFastTiming` to `MDEEStripHit` --- include/MDEEStripHit.h | 2 ++ src/MDEEStripHit.cxx | 3 ++- src/MSubModuleDepthReadout.cxx | 13 +++++++++++-- 3 files changed, 15 insertions(+), 3 deletions(-) diff --git a/include/MDEEStripHit.h b/include/MDEEStripHit.h index 73d602f6..60ae4d69 100644 --- a/include/MDEEStripHit.h +++ b/include/MDEEStripHit.h @@ -84,6 +84,8 @@ struct MDEEStripHit unsigned int m_ADC; //! If the strip exceeds thresholds bool m_HasTriggered; + //! If the strip has fast timing + bool m_HasFastTiming; //! The measured TAC value in ADC units unsigned int m_TAC; //! The calibrated timing in ns diff --git a/src/MDEEStripHit.cxx b/src/MDEEStripHit.cxx index 6ef7a8fb..0286b1b8 100644 --- a/src/MDEEStripHit.cxx +++ b/src/MDEEStripHit.cxx @@ -44,7 +44,7 @@ ClassImp(MDEEStripHit) //////////////////////////////////////////////////////////////////////////////// -MDEEStripHit::MDEEStripHit() : m_SimulatedEventID(0), m_SimulatedPosition(0,0,0), m_SimulatedPositionInDetector(0,0,0), m_SimulatedRelativeDepth(0), m_SimulatedEnergy(0), m_SimulatedIsGuardRing(false), m_SimulatedHitIndex(0), m_IsGuardRing(false), m_Energy(0), m_ADC(0), m_HasTriggered(false), m_TAC(0), m_Temperature(0) +MDEEStripHit::MDEEStripHit() : m_SimulatedEventID(0), m_SimulatedPosition(0,0,0), m_SimulatedPositionInDetector(0,0,0), m_SimulatedRelativeDepth(0), m_SimulatedEnergy(0), m_SimulatedIsGuardRing(false), m_SimulatedHitIndex(0), m_IsGuardRing(false), m_Energy(0), m_ADC(0), m_HasTriggered(false), m_HasFastTiming(false), m_TAC(0), m_Temperature(0) { // Construct an instance of MDEEStripHit } @@ -64,6 +64,7 @@ MStripHit* MDEEStripHit::Convert() SH->HasTriggered(m_HasTriggered); SH->SetADCUnits(m_ADC); SH->SetTAC(m_TAC); + SH->HasFastTiming(m_HasFastTiming); //SH->AddOrigins(); SH->IsGuardRing(m_IsGuardRing); diff --git a/src/MSubModuleDepthReadout.cxx b/src/MSubModuleDepthReadout.cxx index 207d1b2a..9ee3efc5 100644 --- a/src/MSubModuleDepthReadout.cxx +++ b/src/MSubModuleDepthReadout.cxx @@ -141,6 +141,10 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) // where the default value is set to 1e10 if (SH.m_FastPeakTime > -200.0 && SH.m_FastPeakTime < 10000.0) { SH.m_Timing = 4200.0 - SH.m_FastPeakTime; + + // TODO: apply fast threshold + SH.m_HasFastTiming = true; + if (m_ApplyTimingResolutionCalibration == true){ int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID; if (m_Coeffs.count(PixelCode) == 1){ @@ -184,7 +188,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) cout<<"MSubModuleDepthReadout::AnalyzeEvent: Unphysical drift time."< -200.0 && SH.m_FastPeakTime < 10000.0){ SH.m_Timing = 4200.0 - SH.m_FastPeakTime; + + // TODO: apply fast threshold + SH.m_HasFastTiming = true; + if (m_ApplyTimingResolutionCalibration == true) { int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID; if (m_Coeffs.count(PixelCode) == 1){ @@ -209,6 +217,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event) // 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)); + } else { if (g_Verbosity >= c_Info) { cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "< Date: Thu, 30 Jul 2026 12:08:59 -0700 Subject: [PATCH 8/8] Clean up --- src/MSubModuleChargeTransport.cxx | 52 ++++++++++++++----------------- 1 file changed, 24 insertions(+), 28 deletions(-) diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index d5cf2faa..7b410ce2 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -394,36 +394,32 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool m_ChargeTransportHits.push_back(MainSH); // create MDEEStripHit for the left NN - // if (NNLeftStripEnergy > IonizationEnergy) { - MDEEStripHit NNLeftSH = SH; - 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); - NNLeftSH.m_IsGuardRing = false; - } else { - NNLeftSH.m_ROE.SetStripID(NStrips); - NNLeftSH.m_IsGuardRing = true; - } - m_ChargeTransportHits.push_back(NNLeftSH); - // } + MDEEStripHit NNLeftSH = SH; + 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); + NNLeftSH.m_IsGuardRing = false; + } else { + NNLeftSH.m_ROE.SetStripID(NStrips); + NNLeftSH.m_IsGuardRing = true; + } + m_ChargeTransportHits.push_back(NNLeftSH); // create MDEEStripHit for the right NN - // if (NNRightStripEnergy > IonizationEnergy) { - MDEEStripHit NNRightSH = SH; - 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); - NNRightSH.m_IsGuardRing = false; - } else { - NNRightSH.m_ROE.SetStripID(NStrips); - NNRightSH.m_IsGuardRing = true; - } - m_ChargeTransportHits.push_back(NNRightSH); - // } + MDEEStripHit NNRightSH = SH; + 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); + NNRightSH.m_IsGuardRing = false; + } else { + NNRightSH.m_ROE.SetStripID(NStrips); + NNRightSH.m_IsGuardRing = true; + } + m_ChargeTransportHits.push_back(NNRightSH); } else { // TODO: implement charge sharing also for GR events