From 9407f5bd8d45ed3c51ba74a8c025c8acf1317637 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 15 Jun 2026 12:51:03 -0700 Subject: [PATCH 1/3] 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/3] 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/3] 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"