From 31d90c2a1c6f0a386053f61681ec8544d6b72cb6 Mon Sep 17 00:00:00 2001 From: Nico Date: Wed, 19 Aug 2026 10:34:17 -0400 Subject: [PATCH 01/10] Added flag for events with an ADC value higher than 14000. Included in unit test as well --- include/MModuleEnergyCalibration.h | 2 ++ include/MReadOutAssembly.h | 10 ++++++ include/MStripHit.h | 7 +++++ src/MModuleEnergyCalibration.cxx | 10 ++++++ src/MReadOutAssembly.cxx | 14 +++++++++ src/MStripHit.cxx | 47 ++++++++++++++++------------ unittests/UTNStripHit.cxx | 50 ++++++++++++++++++++---------- 7 files changed, 104 insertions(+), 36 deletions(-) diff --git a/include/MModuleEnergyCalibration.h b/include/MModuleEnergyCalibration.h index 0829ac1c..ceb4b60f 100644 --- a/include/MModuleEnergyCalibration.h +++ b/include/MModuleEnergyCalibration.h @@ -199,6 +199,8 @@ class MModuleEnergyCalibration : public MModule //! Max value of the ADC units static constexpr double m_MaxADCRange = 16383; + //! ADC value above which a strip hit is flagged as having a high ADC value + static constexpr double m_HighADCThreshold = 14000; #ifdef ___CLING___ public: diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index 91002736..b53cd34a 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -230,6 +230,11 @@ class MReadOutAssembly : public MReadOutSequence //! Get the Strip Hit Below Threshold quality flag bool HasStripHitBelowThreshold_QualityFlag() const { return m_StripHitBelowThreshold_QualityFlag; } + //! Set the High ADC quality flag + void SetHighADC_QualityFlag(const MString& Text = "") { m_HighADC_QualityFlag = true; if (Text != "") m_HighADCString_QualityFlag.push_back(Text); } + //! Get the High ADC quality flag + bool HasHighADC_QualityFlag() const { return m_HighADC_QualityFlag; } + //! Set the Strip Pairing quality flag void SetStripPairing_QualityFlag(const MString& Text = "") { m_StripPairing_QualityFlag = true; if (Text != "") m_StripPairingString_QualityFlag.push_back(Text); } //! Get the Strip Pairing quality flag @@ -412,6 +417,11 @@ class MReadOutAssembly : public MReadOutSequence //! Strip hit below threshold quality string vector m_StripHitBelowThresholdString_QualityFlag; + //! High ADC quality flag + bool m_HighADC_QualityFlag; + //! High ADC quality string + vector m_HighADCString_QualityFlag; + //! Strip pairing quality flag bool m_StripPairing_QualityFlag; //! Strip pairing quality string diff --git a/include/MStripHit.h b/include/MStripHit.h index 5e35ad44..b1dfecb3 100644 --- a/include/MStripHit.h +++ b/include/MStripHit.h @@ -139,6 +139,11 @@ class MStripHit //! Return whether the strip has triggered (ADC values above slow threshold) bool HasTriggered() const { return m_HasTriggered; } + //! Set whether the strip has a high ADC value (close to the ADC saturation limit) + void HasHighADC(bool HighADC) { m_HasHighADC = HighADC; } + //! Return whether the strip has a high ADC value (close to the ADC saturation limit) + bool HasHighADC() const { return m_HasHighADC; } + //! TODO: Rename to HasTiming() //! Set the calibrated-timing flag void HasCalibratedTiming(bool CalibratedTiming) { m_HasCalibratedTiming = CalibratedTiming; } @@ -211,6 +216,8 @@ class MStripHit bool m_HasFastTiming; //! True if the hit has calibrated timing bool m_HasCalibratedTiming; + //! True if the ADC value is close to the ADC saturation limit + bool m_HasHighADC; //! Origin interaction IDs from the simulation vector m_Origins; diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index c35817d0..c95ae9ba 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -162,6 +162,16 @@ bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) MStripHit* SH = Event->GetStripHit(i); MReadOutElementDoubleStrip R = *dynamic_cast(SH->GetReadOutElement()); + // Flag strip hits whose ADC value is close to the ADC saturation limit, since their + // calibrated energy is not trustworthy. The hit is kept, it is only marked. + if (SH->GetADCUnits() > m_HighADCThreshold) { + SH->HasHighADC(true); + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Warning: High ADC value " << SH->GetADCUnits() << " for read-out element " << R << endl; + } + Event->SetHighADC_QualityFlag("High ADC value " + to_string(SH->GetADCUnits()) + " for " + R.ToString().Data()); + } + TF1* Fit = m_Calibration[R]; TF1* FitRes = m_ResolutionCalibration[R]; diff --git a/src/MReadOutAssembly.cxx b/src/MReadOutAssembly.cxx index 6c68d457..ac79745e 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -139,6 +139,9 @@ void MReadOutAssembly::Clear() m_StripHitBelowThreshold_QualityFlag = false; m_StripHitBelowThresholdString_QualityFlag.clear(); + m_HighADC_QualityFlag = false; + m_HighADCString_QualityFlag.clear(); + m_StripPairing_QualityFlag = false; m_StripPairingString_QualityFlag.clear(); @@ -767,6 +770,17 @@ void MReadOutAssembly::StreamBDFlags(ostream& S) S< bit 0 = 1 H.IsGuardRing(true); H.IsNearestNeighbor(false); H.HasFastTiming(false); Passed = Evaluate("MakeFlags()", "guard ring only", "MakeFlags() returns 1 when only IsGuardRing is set", - H.MakeFlags(), (unsigned int) 0b001) && Passed; + H.MakeFlags(), (unsigned int) 0b0001) && Passed; // Only IsNearestNeighbor -> bit 1 = 2 H.IsGuardRing(false); H.IsNearestNeighbor(true); H.HasFastTiming(false); Passed = Evaluate("MakeFlags()", "nearest neighbor only", "MakeFlags() returns 2 when only IsNearestNeighbor is set", - H.MakeFlags(), (unsigned int) 0b010) && Passed; + H.MakeFlags(), (unsigned int) 0b0010) && Passed; // Only HasFastTiming -> bit 2 = 4 H.IsGuardRing(false); H.IsNearestNeighbor(false); H.HasFastTiming(true); Passed = Evaluate("MakeFlags()", "fast timing only", "MakeFlags() returns 4 when only HasFastTiming is set", - H.MakeFlags(), (unsigned int) 0b100) && Passed; + H.MakeFlags(), (unsigned int) 0b0100) && Passed; - // All three flags -> 7 + // Only HasHighADC -> bit 3 = 8 + H.IsGuardRing(false); + H.IsNearestNeighbor(false); + H.HasFastTiming(false); + H.HasHighADC(true); + Passed = Evaluate("MakeFlags()", "high ADC only", "MakeFlags() returns 8 when only HasHighADC is set", + H.MakeFlags(), (unsigned int) 0b1000) && Passed; + + // All four flags -> 15 H.IsGuardRing(true); H.IsNearestNeighbor(true); H.HasFastTiming(true); - Passed = Evaluate("MakeFlags()", "all flags", "MakeFlags() returns 7 when all three flags are set", - H.MakeFlags(), (unsigned int) 0b111) && Passed; + H.HasHighADC(true); + Passed = Evaluate("MakeFlags()", "all flags", "MakeFlags() returns 15 when all four flags are set", + H.MakeFlags(), (unsigned int) 0b1111) && Passed; - // ParseFlags() round-trip: flags=5 (guard ring + fast timing, no nearest neighbor) + // ParseFlags() round-trip: flags=13 (guard ring + fast timing + high ADC, no nearest neighbor) H.Clear(); - H.ParseFlags(0b101u); - Passed = EvaluateTrue("ParseFlags()", "guard ring bit", "ParseFlags(0b101) sets IsGuardRing true", + H.ParseFlags(0b1101u); + Passed = EvaluateTrue("ParseFlags()", "guard ring bit", "ParseFlags(0b1101) sets IsGuardRing true", H.IsGuardRing() == true) && Passed; - Passed = EvaluateFalse("ParseFlags()", "nearest neighbor bit", "ParseFlags(0b101) leaves IsNearestNeighbor false", + Passed = EvaluateFalse("ParseFlags()", "nearest neighbor bit", "ParseFlags(0b1101) leaves IsNearestNeighbor false", H.IsNearestNeighbor()) && Passed; - Passed = EvaluateTrue("ParseFlags()", "fast timing bit", "ParseFlags(0b101) sets HasFastTiming true", + Passed = EvaluateTrue("ParseFlags()", "fast timing bit", "ParseFlags(0b1101) sets HasFastTiming true", H.HasFastTiming() == true) && Passed; - Passed = Evaluate("ParseFlags()", "round-trip", "MakeFlags() reproduces the representative flags value 5 after ParseFlags(5)", - H.MakeFlags(), (unsigned int) 0b101) && Passed; + Passed = EvaluateTrue("ParseFlags()", "high ADC bit", "ParseFlags(0b1101) sets HasHighADC true", + H.HasHighADC() == true) && Passed; + Passed = Evaluate("ParseFlags()", "round-trip", "MakeFlags() reproduces the representative flags value 13 after ParseFlags(13)", + H.MakeFlags(), (unsigned int) 0b1101) && Passed; // ParseFlags(0) clears all flags H.IsGuardRing(true); H.IsNearestNeighbor(true); H.HasFastTiming(true); + H.HasHighADC(true); H.ParseFlags(0); Passed = EvaluateFalse("ParseFlags(0)", "guard ring cleared", "ParseFlags(0) clears IsGuardRing", H.IsGuardRing()) && Passed; Passed = EvaluateFalse("ParseFlags(0)", "nearest neighbor cleared", "ParseFlags(0) clears IsNearestNeighbor", H.IsNearestNeighbor()) && Passed; Passed = EvaluateFalse("ParseFlags(0)", "fast timing cleared", "ParseFlags(0) clears HasFastTiming", H.HasFastTiming()) && Passed; + Passed = EvaluateFalse("ParseFlags(0)", "high ADC cleared", "ParseFlags(0) clears HasHighADC", H.HasHighADC()) && Passed; // HasCalibratedTiming is intentionally not part of the bit mask H.Clear(); H.HasCalibratedTiming(false); Passed = Evaluate("MakeFlags()", "calibrated timing false", "MakeFlags() is unchanged when HasCalibratedTiming is false", - H.MakeFlags(), (unsigned int) 0b000) && Passed; + H.MakeFlags(), (unsigned int) 0b0000) && Passed; H.HasCalibratedTiming(true); Passed = Evaluate("MakeFlags()", "calibrated timing true", "MakeFlags() is unchanged when HasCalibratedTiming is true", - H.MakeFlags(), (unsigned int) 0b000) && Passed; + H.MakeFlags(), (unsigned int) 0b0000) && Passed; return Passed; } From 942c7dbc55f44e80f6fb8c45eceb553d3b7e324d Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 14:56:00 -0400 Subject: [PATCH 02/10] Undo unnecessary changes for the highADC flag. --- include/MStripHit.h | 7 ------ src/MStripHit.cxx | 47 ++++++++++++++++-------------------- unittests/UTNStripHit.cxx | 50 +++++++++++++-------------------------- 3 files changed, 36 insertions(+), 68 deletions(-) diff --git a/include/MStripHit.h b/include/MStripHit.h index b1dfecb3..5e35ad44 100644 --- a/include/MStripHit.h +++ b/include/MStripHit.h @@ -139,11 +139,6 @@ class MStripHit //! Return whether the strip has triggered (ADC values above slow threshold) bool HasTriggered() const { return m_HasTriggered; } - //! Set whether the strip has a high ADC value (close to the ADC saturation limit) - void HasHighADC(bool HighADC) { m_HasHighADC = HighADC; } - //! Return whether the strip has a high ADC value (close to the ADC saturation limit) - bool HasHighADC() const { return m_HasHighADC; } - //! TODO: Rename to HasTiming() //! Set the calibrated-timing flag void HasCalibratedTiming(bool CalibratedTiming) { m_HasCalibratedTiming = CalibratedTiming; } @@ -216,8 +211,6 @@ class MStripHit bool m_HasFastTiming; //! True if the hit has calibrated timing bool m_HasCalibratedTiming; - //! True if the ADC value is close to the ADC saturation limit - bool m_HasHighADC; //! Origin interaction IDs from the simulation vector m_Origins; diff --git a/src/MStripHit.cxx b/src/MStripHit.cxx index 9dd230fe..9d165168 100644 --- a/src/MStripHit.cxx +++ b/src/MStripHit.cxx @@ -90,7 +90,6 @@ void MStripHit::Clear() m_HasFastTiming = false; m_HasCalibratedTiming = false; - m_HasHighADC = false; m_Origins.clear(); } @@ -228,25 +227,21 @@ void MStripHit::StreamRoa(ostream& S, bool WithADC, bool WithTAC, bool WithEnerg unsigned int MStripHit::MakeFlags() { // Return the bitwise strip-hit flags - // Currently, 4 bits: - // v = Has a high ADC value - // v = Has fast timing - // v = Is a nearest neighbor strip hit - // v = Is a guard ring strip hit - // 0b1111u - - unsigned int Flags = 0b0000u; + // Currently, 3 bits: + // v = Has fast timing + // v = Is a nearest neighbor strip hit + // v = Is a guard ring strip hit + // 0b111u + + unsigned int Flags = 0b000u; if (m_IsGuardRing == true) { - Flags = Flags | 0b0001u; + Flags = Flags | 0b001u; } if (m_IsNearestNeighbor == true) { - Flags = Flags | 0b0010u; + Flags = Flags | 0b010u; } if (m_HasFastTiming == true) { - Flags = Flags | 0b0100u; - } - if (m_HasHighADC == true) { - Flags = Flags | 0b1000u; + Flags = Flags | 0b100u; } return Flags; @@ -259,18 +254,16 @@ unsigned int MStripHit::MakeFlags() void MStripHit::ParseFlags(unsigned int Flags) { // Update the strip-hit flags from a bit mask - // Currently, 4 bits: - // v = Has a high ADC value - // v = Has fast timing - // v = Is a nearest neighbor - // v = Is a guard ring - // 0b1111u - - // "Flags & 0b0001u" extracts bit 0, "!= 0u" turns it into an explicit bool - IsGuardRing((Flags & 0b0001u) != 0u); - IsNearestNeighbor((Flags & 0b0010u) != 0u); - HasFastTiming((Flags & 0b0100u) != 0u); - HasHighADC((Flags & 0b1000u) != 0u); + // Currently, 3 bits: + // v = Has fast timing + // v = Is a nearest neighbor + // v = Is a guard ring + // 0b111u + + // "Flags & 0b001u" extracts bit 0, "!= 0u" turns it into an explicit bool + IsGuardRing((Flags & 0b001u) != 0u); + IsNearestNeighbor((Flags & 0b010u) != 0u); + HasFastTiming((Flags & 0b100u) != 0u); } diff --git a/unittests/UTNStripHit.cxx b/unittests/UTNStripHit.cxx index 2f83bd5f..228504ed 100644 --- a/unittests/UTNStripHit.cxx +++ b/unittests/UTNStripHit.cxx @@ -91,7 +91,6 @@ bool UTNStripHit::TestDefaultConstruction() Passed = EvaluateFalse("IsNearestNeighbor()", "default", "Default IsNearestNeighbor is false", H.IsNearestNeighbor()) && Passed; Passed = EvaluateFalse("HasFastTiming()", "default", "Default HasFastTiming is false", H.HasFastTiming()) && Passed; Passed = EvaluateFalse("HasCalibratedTiming()", "default", "Default HasCalibratedTiming is false", H.HasCalibratedTiming()) && Passed; - Passed = EvaluateFalse("HasHighADC()", "default", "Default HasHighADC is false", H.HasHighADC()) && Passed; Passed = Evaluate("GetOrigins().size()", "default", "Default origins list is empty", (unsigned int) H.GetOrigins().size(), (unsigned int) 0) && Passed; // Verify that Clear() reinstates all defaults @@ -107,7 +106,6 @@ bool UTNStripHit::TestDefaultConstruction() H.IsNearestNeighbor(true); H.HasFastTiming(true); H.HasCalibratedTiming(true); - H.HasHighADC(true); H.AddOrigins({3, 7}); H.Clear(); @@ -124,7 +122,6 @@ bool UTNStripHit::TestDefaultConstruction() Passed = EvaluateFalse("Clear() IsNearestNeighbor", "after clear", "Clear() resets IsNearestNeighbor to false", H.IsNearestNeighbor()) && Passed; Passed = EvaluateFalse("Clear() HasFastTiming", "after clear", "Clear() resets HasFastTiming to false", H.HasFastTiming()) && Passed; Passed = EvaluateFalse("Clear() HasCalibratedTiming", "after clear", "Clear() resets HasCalibratedTiming to false", H.HasCalibratedTiming()) && Passed; - Passed = EvaluateFalse("Clear() HasHighADC", "after clear", "Clear() resets HasHighADC to false", H.HasHighADC()) && Passed; Passed = Evaluate("Clear() origins", "after clear", "Clear() empties the origins list", (unsigned int) H.GetOrigins().size(), (unsigned int) 0) && Passed; return Passed; @@ -230,7 +227,6 @@ bool UTNStripHit::TestMakeParseFlags() // bit 0 (value 1) = IsGuardRing // bit 1 (value 2) = IsNearestNeighbor // bit 2 (value 4) = HasFastTiming - // bit 3 (value 8) = HasHighADC MStripHit H; @@ -238,80 +234,66 @@ bool UTNStripHit::TestMakeParseFlags() H.IsGuardRing(false); H.IsNearestNeighbor(false); H.HasFastTiming(false); - H.HasHighADC(false); Passed = Evaluate("MakeFlags()", "all off", "MakeFlags() returns 0 when no flag is set", - H.MakeFlags(), (unsigned int) 0b0000) && Passed; + H.MakeFlags(), (unsigned int) 0b000) && Passed; // Only IsGuardRing -> bit 0 = 1 H.IsGuardRing(true); H.IsNearestNeighbor(false); H.HasFastTiming(false); Passed = Evaluate("MakeFlags()", "guard ring only", "MakeFlags() returns 1 when only IsGuardRing is set", - H.MakeFlags(), (unsigned int) 0b0001) && Passed; + H.MakeFlags(), (unsigned int) 0b001) && Passed; // Only IsNearestNeighbor -> bit 1 = 2 H.IsGuardRing(false); H.IsNearestNeighbor(true); H.HasFastTiming(false); Passed = Evaluate("MakeFlags()", "nearest neighbor only", "MakeFlags() returns 2 when only IsNearestNeighbor is set", - H.MakeFlags(), (unsigned int) 0b0010) && Passed; + H.MakeFlags(), (unsigned int) 0b010) && Passed; // Only HasFastTiming -> bit 2 = 4 H.IsGuardRing(false); H.IsNearestNeighbor(false); H.HasFastTiming(true); Passed = Evaluate("MakeFlags()", "fast timing only", "MakeFlags() returns 4 when only HasFastTiming is set", - H.MakeFlags(), (unsigned int) 0b0100) && Passed; + H.MakeFlags(), (unsigned int) 0b100) && Passed; - // Only HasHighADC -> bit 3 = 8 - H.IsGuardRing(false); - H.IsNearestNeighbor(false); - H.HasFastTiming(false); - H.HasHighADC(true); - Passed = Evaluate("MakeFlags()", "high ADC only", "MakeFlags() returns 8 when only HasHighADC is set", - H.MakeFlags(), (unsigned int) 0b1000) && Passed; - - // All four flags -> 15 + // All three flags -> 7 H.IsGuardRing(true); H.IsNearestNeighbor(true); H.HasFastTiming(true); - H.HasHighADC(true); - Passed = Evaluate("MakeFlags()", "all flags", "MakeFlags() returns 15 when all four flags are set", - H.MakeFlags(), (unsigned int) 0b1111) && Passed; + Passed = Evaluate("MakeFlags()", "all flags", "MakeFlags() returns 7 when all three flags are set", + H.MakeFlags(), (unsigned int) 0b111) && Passed; - // ParseFlags() round-trip: flags=13 (guard ring + fast timing + high ADC, no nearest neighbor) + // ParseFlags() round-trip: flags=5 (guard ring + fast timing, no nearest neighbor) H.Clear(); - H.ParseFlags(0b1101u); - Passed = EvaluateTrue("ParseFlags()", "guard ring bit", "ParseFlags(0b1101) sets IsGuardRing true", + H.ParseFlags(0b101u); + Passed = EvaluateTrue("ParseFlags()", "guard ring bit", "ParseFlags(0b101) sets IsGuardRing true", H.IsGuardRing() == true) && Passed; - Passed = EvaluateFalse("ParseFlags()", "nearest neighbor bit", "ParseFlags(0b1101) leaves IsNearestNeighbor false", + Passed = EvaluateFalse("ParseFlags()", "nearest neighbor bit", "ParseFlags(0b101) leaves IsNearestNeighbor false", H.IsNearestNeighbor()) && Passed; - Passed = EvaluateTrue("ParseFlags()", "fast timing bit", "ParseFlags(0b1101) sets HasFastTiming true", + Passed = EvaluateTrue("ParseFlags()", "fast timing bit", "ParseFlags(0b101) sets HasFastTiming true", H.HasFastTiming() == true) && Passed; - Passed = EvaluateTrue("ParseFlags()", "high ADC bit", "ParseFlags(0b1101) sets HasHighADC true", - H.HasHighADC() == true) && Passed; - Passed = Evaluate("ParseFlags()", "round-trip", "MakeFlags() reproduces the representative flags value 13 after ParseFlags(13)", - H.MakeFlags(), (unsigned int) 0b1101) && Passed; + Passed = Evaluate("ParseFlags()", "round-trip", "MakeFlags() reproduces the representative flags value 5 after ParseFlags(5)", + H.MakeFlags(), (unsigned int) 0b101) && Passed; // ParseFlags(0) clears all flags H.IsGuardRing(true); H.IsNearestNeighbor(true); H.HasFastTiming(true); - H.HasHighADC(true); H.ParseFlags(0); Passed = EvaluateFalse("ParseFlags(0)", "guard ring cleared", "ParseFlags(0) clears IsGuardRing", H.IsGuardRing()) && Passed; Passed = EvaluateFalse("ParseFlags(0)", "nearest neighbor cleared", "ParseFlags(0) clears IsNearestNeighbor", H.IsNearestNeighbor()) && Passed; Passed = EvaluateFalse("ParseFlags(0)", "fast timing cleared", "ParseFlags(0) clears HasFastTiming", H.HasFastTiming()) && Passed; - Passed = EvaluateFalse("ParseFlags(0)", "high ADC cleared", "ParseFlags(0) clears HasHighADC", H.HasHighADC()) && Passed; // HasCalibratedTiming is intentionally not part of the bit mask H.Clear(); H.HasCalibratedTiming(false); Passed = Evaluate("MakeFlags()", "calibrated timing false", "MakeFlags() is unchanged when HasCalibratedTiming is false", - H.MakeFlags(), (unsigned int) 0b0000) && Passed; + H.MakeFlags(), (unsigned int) 0b000) && Passed; H.HasCalibratedTiming(true); Passed = Evaluate("MakeFlags()", "calibrated timing true", "MakeFlags() is unchanged when HasCalibratedTiming is true", - H.MakeFlags(), (unsigned int) 0b0000) && Passed; + H.MakeFlags(), (unsigned int) 0b000) && Passed; return Passed; } From a42c560efb315f2eb1e5bdc194ca77f59b143c64 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 15:21:39 -0400 Subject: [PATCH 03/10] Remove highADC flag information from strip hit. This was a missing file to restore from the previous commit. --- src/MModuleEnergyCalibration.cxx | 572 +------------------------------ 1 file changed, 17 insertions(+), 555 deletions(-) diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index c95ae9ba..b57d4762 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -2,14 +2,12 @@ * MModuleEnergyCalibration.cxx * * - * Copyright (C) by Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, - * Jackie Beechert, Robin Anthony-Peterson. + * Copyright (C) 2008-2008 by Andreas Zoglauer. * All rights reserved. * * * This code implementation is the intellectual property of - * Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, - * Jackie Beechert, Robin Anthony-Peterson + * Andreas Zoglauer. * * By copying, distributing or modifying the Program (or any work * based on the Program) you indicate your acceptance of this statement, @@ -24,29 +22,16 @@ // //////////////////////////////////////////////////////////////////////////////// -// Standard libs: -#include -#include -#include -using namespace std; // Include the header: #include "MModuleEnergyCalibration.h" +// Standard libs: // ROOT libs: #include "TGClient.h" -#include "TFile.h" // MEGAlib libs: -#include "MString.h" -#include "MStreams.h" - -// Nuclearizer libs: -#include "MReadOutElement.h" -#include "MReadOutElementDoubleStrip.h" -#include "MGUIOptionsEnergyCalibration.h" -#include "MGUIExpoPlotSpectrum.h" //////////////////////////////////////////////////////////////////////////////// @@ -67,39 +52,28 @@ MModuleEnergyCalibration::MModuleEnergyCalibration() : MModule() // Set all module relevant information // Set the module name --- has to be unique - m_Name = "Energy calibrator"; + m_Name = "Combined energy calibration and charge sharing correction"; // Set the XML tag --- has to be unique --- no spaces allowed - m_XmlTag = "EnergyCalibration"; + m_XmlTag = "CombinedEnergyCalibrationAndChargeSharingCorrection"; // Set all modules, which have to be done before this module - AddPreceedingModuleType(MAssembly::c_EventLoader); - // AddPreceedingModuleType(MAssembly::c_TACcut); + AddModuleType(MAssembly::c_EventLoader); // Set all types this modules handles AddModuleType(MAssembly::c_EnergyCalibration); + AddModuleType(MAssembly::c_ChargeSharingCorrection); // Set all modules, which can follow this module - AddSucceedingModuleType(MAssembly::c_TACcut); // Set if this module has an options GUI - m_HasOptionsGUI = true; - - // Allow the use of multiple threads and instances - m_AllowMultiThreading = true; - m_AllowMultipleInstances = true; - - // Initiate the Slow Threshold Cut variables - m_SlowThresholdCutMode = MSlowThresholdCutModes::e_Ignore; - m_SlowThresholdCutFixedValue = 15; - m_SlowThresholdCutFileName = ""; - - // Nearest Neighbor threshold (-1000 because we don't want to default to cut neg values) - m_NearestNeighborCutMode = MNearestNeighborCutModes::e_Ignore; - m_NearestNeighborThreshold = -1000.0; - + // If true, overwrite ShowOptionsGUI() with the call to the GUI! + m_HasOptionsGUI = false; + // If true, you have to derive a class from MGUIOptions (use MGUIOptionsTemplate) + // and implement all your GUI options } + //////////////////////////////////////////////////////////////////////////////// @@ -112,411 +86,9 @@ MModuleEnergyCalibration::~MModuleEnergyCalibration() //////////////////////////////////////////////////////////////////////////////// -void MModuleEnergyCalibration::CreateExpos() -{ - // If they are already created, return - if (m_Expos.size() != 0) { - return; - } - - // Set the histogram display - m_ExpoSpectrum = new MGUIExpoPlotSpectrum(this); - m_ExpoSpectrum->SetEnergyHistogramParameters(200, 0, 2000); - m_Expos.push_back(m_ExpoSpectrum); - -} - - -//////////////////////////////////////////////////////////////////////////////// - - bool MModuleEnergyCalibration::Initialize() { - // Initialize the module - - // Parse the energy calibration file - if (ReadEnergyCalibrationFile(m_FileName) == false) { - return false; - } - - // Parse the slow threshold file - if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { - if (ReadSlowThresholdCutFile(m_SlowThresholdCutFileName) == false) { - return false; - } - } - - return MModule::Initialize(); -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) -{ - // Main data analysis routine, which updates the event to a new level, i.e. takes the raw ADC value from the .roa file loaded through nuclearizer and converts it into energy units. - - for (unsigned int i = 0; i < Event->GetNStripHits();) { - - MStripHit* SH = Event->GetStripHit(i); - MReadOutElementDoubleStrip R = *dynamic_cast(SH->GetReadOutElement()); - - // Flag strip hits whose ADC value is close to the ADC saturation limit, since their - // calibrated energy is not trustworthy. The hit is kept, it is only marked. - if (SH->GetADCUnits() > m_HighADCThreshold) { - SH->HasHighADC(true); - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": Warning: High ADC value " << SH->GetADCUnits() << " for read-out element " << R << endl; - } - Event->SetHighADC_QualityFlag("High ADC value " + to_string(SH->GetADCUnits()) + " for " + R.ToString().Data()); - } - - TF1* Fit = m_Calibration[R]; - TF1* FitRes = m_ResolutionCalibration[R]; - - if (Fit == nullptr) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Error: Energy-fit not found for read-out element " << R << endl; - } - Event->SetEnergyCalibrationError("calibration not found for " + R.ToString()); - ++i; // iterate to next SH - continue; - - } else { - - double Energy = 0; - double Threshold; // No default value here, we set it below (different for strips and nearest neighbors) - - Energy = Fit->Eval(SH->GetADCUnits()); - - // TODO(@RobinAnthonyPetersen): Determine if we want to force negative energy values to be zero or not - // If the calibrated energy is less than 0, force it to be 0. - //if (Energy < 0) { - // Energy = 0; - //} - - if (SH->IsNearestNeighbor() == true) { - // If nothing is selected, we set thresholds to zero - // Threshold = 0.0; // But keeping negative values for now - // If nothing is selected, we keep negative nearest neighbors - Threshold = -100.0; - - // Otherwise, if the user inputs a value, we use that threshold - if (m_NearestNeighborCutMode == MNearestNeighborCutModes::e_Fixed) { - // Get the value user typed in the box (for example 6.0 keV) - // TODO(@RobinAnthonyPetersen): Nearest Neighbor threhsold cut subject to change pending more analysis - Threshold = m_NearestNeighborThreshold; - } - } else { // If not Nearest Neighbor, then it's a triggered strip - - // Set the default slow threshold - Threshold = 0.0; - - if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_Fixed) { // check if user input threshold is enabled (one value applied to all strips) - Threshold = m_SlowThresholdCutFixedValue; - } else if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { // check if threshold file is enabled (unique value applied to each strip) - double ThresholdFromFile = m_ThresholdMap[R]; // if file enabled, declare value from map - - if (ThresholdFromFile == 0) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Error: Threshold not found for read-out element " << R << endl; - } - const double DefaultSlowThreshold = 15.0; - Threshold = DefaultSlowThreshold; // set default threshold if threshold not found - } else { - Threshold = ThresholdFromFile; // set threshold variable to value found in map - } - } - } - - - // Remove SH for any energy value below the established threshold (0 is default) - if (Energy < Threshold) { - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": Strip Hit below threshold, deleting SH with Energy " << Energy << " keV " << endl; - } - Event->SetStripHitBelowThreshold_QualityFlag("Strip hit removed with energy " + to_string(Energy)); - Event->RemoveStripHit(i); - delete SH; - continue; // continue to next SH without iterating i - } else { - - // if the energy isn't filtered out with the threshold, then assign the energy to the SH - ++i; // iterate to next SH - SH->SetEnergy(Energy); - - if (FitRes == nullptr) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Error: Energy Resolution fit not found for read-out element " << R << endl; - } - // There is not expected to be a time in which the energy resolution calibration is not defined when the energy calibration itself is. Therefore, don't need a seperate BD flag for this. - } else { - double EnergyResolution = FitRes->Eval(Energy); - SH->SetEnergyResolution(EnergyResolution); - } - if (HasExpos() == true) { - m_ExpoSpectrum->AddEnergyFinal(Energy, SH->IsNearestNeighbor(), SH->IsLowVoltageStrip()); - } - if (g_Verbosity >= c_Info) { - cout << m_XmlTag << ": Energy: " << SH->GetADCUnits() << " adc --> " << Energy << " keV" << endl; - } - } - } - } - Event->SetAnalysisProgress(MAssembly::c_EnergyCalibration); - - return true; -} - - -///////////////////////////////////////////////////////////////////////////////// - - -double MModuleEnergyCalibration::GetEnergy(MReadOutElementDoubleStrip R, double ADC) -{ - //! Return the energy for a given ADC value or zero in case of error - - TF1* Fit = m_Calibration[R]; - double Energy = 0.0; - if (Fit != nullptr) { - Energy = Fit->Eval(ADC); - if (Energy < 0.0) { - Energy = 0.0; - } - } else { - cout << m_Name << ": GetEnergy: Error unable to find calibration" << endl; - return 0; - } - - return Energy; -} - - -/////////////////////////////////////////////////////////////////////////////// - - -double MModuleEnergyCalibration::GetADC(MReadOutElementDoubleStrip R, double Energy) -{ - //! Return the ADC value for a given energy - - TF1* Fit = m_Calibration[R]; // TF1* is a function to be applied - if (Fit != nullptr) { - return Fit->GetX(Energy); - } else { - cout << m_Name << ": GetADC: Error unable to find calibration" << endl; - return 0; - } -} - - -/////////////////////////////////////////////////////////////////////////////// - - -void MModuleEnergyCalibration::Finalize() -{ - // Finalize the calibrator and clean up - - MModule::Finalize(); - - for (auto& F : m_Calibration) { - delete F.second; - } - for (auto& F : m_ResolutionCalibration) { - delete F.second; - } - - return; -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MModuleEnergyCalibration::ReadEnergyCalibrationFile(MString FileName) { - - // Parse the Energy Calibration file - MParser Parser; - if (Parser.Open(FileName, MFile::c_Read) == false) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Unable to open energy calibration file " << FileName << endl; - } - return false; - } - - // create other maps - map CP_ROEToLine; //Peak fits - map CM_ROEToLine; //Energy Calibration Model - map CR_ROEToLine; //Energy Resolution Calibration Model - - // tokenize ecal file - for (unsigned int i = 0; i < Parser.GetNLines(); ++i) { - unsigned int NTokens = Parser.GetTokenizerAt(i)->GetNTokens(); - if (NTokens < 2) { - continue; - } - if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true || - Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true || Parser.GetTokenizerAt(i)->IsTokenAt(0, "CR") == true) { - if (Parser.GetTokenizerAt(i)->IsTokenAt(1, "dss") == true) { - - // input token values to map - MReadOutElementDoubleStrip R; - R.SetDetectorID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(2)); - R.SetStripID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(3)); - R.IsLowVoltageStrip((Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "p") || (Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "l")); - if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true) { - CP_ROEToLine[R] = i; - } else if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true) { - CM_ROEToLine[R] = i; - } else { - CR_ROEToLine[R] = i; - } - } else { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Line parser: Unknown read-out element (" << Parser.GetTokenizerAt(i)->GetTokenAt(1) << ")" << endl; - } - return false; - } - } - } - - for (auto CM : CM_ROEToLine) { - // If we have at least three data points, we store the calibration - - if (CP_ROEToLine.find(CM.first) != CP_ROEToLine.end()) { - unsigned int i = CP_ROEToLine[CM.first]; - if (Parser.GetTokenizerAt(i)->IsTokenAt(5, "pakw") == false) { - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": Unknown calibration point descriptor found: " << Parser.GetTokenizerAt(i)->GetTokenAt(5) << endl; - } - continue; - } - } else { - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": No good calibration for the following strip found: " << CM.first << endl; - } - continue; - } - - // Read the calibrator, i.e. read the fit function from the .ecal Melinator file. - - unsigned int Pos = 5; - MString CalibratorType = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsString(Pos); - CalibratorType.ToLower(); - - // Below inclusion of poly1zero written by J. Beechert on 2019/11/15 - if (CalibratorType == "poly1zero") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly1zero", "0. + [0]*x", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } - - // Below inclusion of poly1 and poly2 written by J. Beechert on 2019/10/24 - else if (CalibratorType == "poly1") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly1", "[0] + [1]*x", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } else if (CalibratorType == "poly2") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly2", "[0] + [1]*x + [2]*x^2", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - melinatorfit->FixParameter(2, a2); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } - // Eventually, I'll be including other possible fits, but for now, we've just include poly3 and poly4 - else if (CalibratorType == "poly3") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly3", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - melinatorfit->FixParameter(2, a2); - melinatorfit->FixParameter(3, a3); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } else if (CalibratorType == "poly4") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a4 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly4", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - melinatorfit->FixParameter(2, a2); - melinatorfit->FixParameter(3, a3); - melinatorfit->FixParameter(4, a4); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } else { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Line parser: Unknown calibrator type (" << CalibratorType << ") for strip" << CM.first << endl; - } - continue; - } - } - - for (auto CR : CR_ROEToLine) { - - unsigned int Pos = 5; - MString CalibratorType = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsString(Pos); - CalibratorType.ToLower(); - if (CalibratorType == "p1" || CalibratorType == "poly1") { - double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - TF1* resolutionfit = new TF1("P1", "([0]+[1]*x) / 2.355", 0., 2000.); - resolutionfit->FixParameter(0, f0); - resolutionfit->FixParameter(1, f1); - - m_ResolutionCalibration[CR.first] = resolutionfit; - } else if (CalibratorType == "p2" || CalibratorType == "poly2") { - double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - double f2 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - TF1* resolutionfit = new TF1("P2", "([0]+[1]*x+[2]*x*x) / 2.355", 0., 2000.); - resolutionfit->FixParameter(0, f0); - resolutionfit->FixParameter(1, f1); - resolutionfit->FixParameter(2, f2); - m_ResolutionCalibration[CR.first] = resolutionfit; - } else { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Line parser: Unknown resolution calibrator type (" << CalibratorType << ") for strip" << CR.first << endl; - } - continue; - } - } + // Initialize the module return true; } @@ -525,88 +97,9 @@ bool MModuleEnergyCalibration::ReadEnergyCalibrationFile(MString FileName) { //////////////////////////////////////////////////////////////////////////////// -bool MModuleEnergyCalibration::ReadSlowThresholdCutFile(MString FileName) { - - MParser Parser; - if (Parser.Open(FileName, MFile::c_Read) == false) { - if (g_Verbosity >= c_Error) { - cout< Tokens = Line.Tokenize(","); // for each line, Create tokens seperated by commas - if (Tokens.size() == 6) { - int IndexOffset = Tokens.size() % 6; //index counter - int DetID = Tokens[1 + IndexOffset].ToInt(); // Detector ID - MString Side = Tokens[2 + IndexOffset].ToString(); // side is a string, either 'l' or 'h' - int StripID = Tokens[3 + IndexOffset].ToInt(); // stripID - //int ThresholdADC = Tokens[4 + IndexOffset].ToInt(); // energy threshold in ADC - double ThresholdKeVFile = Tokens[5 + IndexOffset].ToDouble(); //energy threshold in keV - - MReadOutElementDoubleStrip R; - R.SetDetectorID(DetID); - R.SetStripID(StripID); - R.IsLowVoltageStrip(Side == "l"); - - // map detectorID, strip number, and voltage side to the threshold (keV) - m_ThresholdMap[R] = ThresholdKeVFile; - } - } - } - - return true; -} - - -//////////////////////////////////////////////////////////////////////////////// - - -void MModuleEnergyCalibration::ShowOptionsGUI() +bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) { - // Show the options GUI - - MGUIOptionsEnergyCalibration* Options = new MGUIOptionsEnergyCalibration(this); - Options->Create(); - gClient->WaitForUnmap(Options); -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MModuleEnergyCalibration::ReadXmlConfiguration(MXmlNode* Node) -{ - //! Read the configuration data from an XML node - - MXmlNode* FileNameNode = Node->GetNode("FileName"); - if (FileNameNode != nullptr) { - m_FileName = FileNameNode->GetValue(); - } - - MXmlNode* SlowThresholdCutModeNode = Node->GetNode("SlowThresholdCutMode"); - if (SlowThresholdCutModeNode != nullptr) { - m_SlowThresholdCutMode = static_cast(SlowThresholdCutModeNode->GetValueAsInt()); - } - MXmlNode* SlowThresholdCutFixedValueNode = Node->GetNode("SlowThresholdCutFixedValue"); - if (SlowThresholdCutFixedValueNode != nullptr) { - m_SlowThresholdCutFixedValue = SlowThresholdCutFixedValueNode->GetValueAsDouble(); - } - MXmlNode* SlowThresholdCutFileNameNode = Node->GetNode("SlowThresholdCutThresholdFileName"); - if (SlowThresholdCutFileNameNode != nullptr) { - m_SlowThresholdCutFileName = SlowThresholdCutFileNameNode->GetValue(); - } - MXmlNode* NNCutModeNode = Node->GetNode("NearestNeighborCutMode"); - if (NNCutModeNode != nullptr) { - m_NearestNeighborCutMode = static_cast(NNCutModeNode->GetValueAsInt()); - } - MXmlNode* NearestNeighborThresholdNode = Node->GetNode("NearestNeighborThreshold"); - if (NearestNeighborThresholdNode != nullptr) { - m_NearestNeighborThreshold = NearestNeighborThresholdNode->GetValueAsDouble(); - } + // Main data analysis routine, which updates the event to a new level return true; } @@ -615,40 +108,9 @@ bool MModuleEnergyCalibration::ReadXmlConfiguration(MXmlNode* Node) //////////////////////////////////////////////////////////////////////////////// -MXmlNode* MModuleEnergyCalibration::CreateXmlConfiguration() -{ - //! Create an XML node tree from the configuration - - MXmlNode* Node = new MXmlNode(0, m_XmlTag); - new MXmlNode(Node, "FileName", m_FileName); - new MXmlNode(Node, "SlowThresholdCutMode", static_cast(m_SlowThresholdCutMode)); - new MXmlNode(Node, "SlowThresholdCutFixedValue", m_SlowThresholdCutFixedValue); - new MXmlNode(Node, "SlowThresholdCutThresholdFileName", m_SlowThresholdCutFileName); - new MXmlNode(Node, "NearestNeighborCutMode", static_cast(m_NearestNeighborCutMode)); - new MXmlNode(Node, "NearestNeighborThreshold", m_NearestNeighborThreshold); - - return Node; -} - -///////////////////////////////////////////////////////////////////////////////// - - -double MModuleEnergyCalibration::LookupEnergyResolution(MStripHit* SH, double Energy) +void MModuleEnergyCalibration::ShowOptionsGUI() { - //! Return the energy resolution or -1 in case of error - - MReadOutElementDoubleStrip* ROE = dynamic_cast(SH->GetReadOutElement()); - if (ROE == nullptr) { - cout << m_Name << ": LookupEnergyResolution: Error unable to get read-out element" << endl; - return -1; - } - TF1* FitRes = m_ResolutionCalibration[*ROE]; - if (FitRes == nullptr) { - cout << m_Name << ": LookupEnergyResolutio: Error: Couldn't locate energy resolution" << endl; - return -1.0; - } else { - return FitRes->Eval(Energy); - } + // Show the options GUI - or do nothing } From de9c7fd120e34a4771e1724016ad64224b8723e7 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 15:25:52 -0400 Subject: [PATCH 04/10] Oops, that was restored to main and not develop-em. Okay, this time for good --- src/MModuleEnergyCalibration.cxx | 562 ++++++++++++++++++++++++++++++- 1 file changed, 545 insertions(+), 17 deletions(-) diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index b57d4762..c35817d0 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -2,12 +2,14 @@ * MModuleEnergyCalibration.cxx * * - * Copyright (C) 2008-2008 by Andreas Zoglauer. + * Copyright (C) by Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, + * Jackie Beechert, Robin Anthony-Peterson. * All rights reserved. * * * This code implementation is the intellectual property of - * Andreas Zoglauer. + * Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, + * Jackie Beechert, Robin Anthony-Peterson * * By copying, distributing or modifying the Program (or any work * based on the Program) you indicate your acceptance of this statement, @@ -22,16 +24,29 @@ // //////////////////////////////////////////////////////////////////////////////// +// Standard libs: +#include +#include +#include +using namespace std; // Include the header: #include "MModuleEnergyCalibration.h" -// Standard libs: // ROOT libs: #include "TGClient.h" +#include "TFile.h" // MEGAlib libs: +#include "MString.h" +#include "MStreams.h" + +// Nuclearizer libs: +#include "MReadOutElement.h" +#include "MReadOutElementDoubleStrip.h" +#include "MGUIOptionsEnergyCalibration.h" +#include "MGUIExpoPlotSpectrum.h" //////////////////////////////////////////////////////////////////////////////// @@ -52,27 +67,38 @@ MModuleEnergyCalibration::MModuleEnergyCalibration() : MModule() // Set all module relevant information // Set the module name --- has to be unique - m_Name = "Combined energy calibration and charge sharing correction"; + m_Name = "Energy calibrator"; // Set the XML tag --- has to be unique --- no spaces allowed - m_XmlTag = "CombinedEnergyCalibrationAndChargeSharingCorrection"; + m_XmlTag = "EnergyCalibration"; // Set all modules, which have to be done before this module - AddModuleType(MAssembly::c_EventLoader); + AddPreceedingModuleType(MAssembly::c_EventLoader); + // AddPreceedingModuleType(MAssembly::c_TACcut); // Set all types this modules handles AddModuleType(MAssembly::c_EnergyCalibration); - AddModuleType(MAssembly::c_ChargeSharingCorrection); // Set all modules, which can follow this module + AddSucceedingModuleType(MAssembly::c_TACcut); // Set if this module has an options GUI - // If true, overwrite ShowOptionsGUI() with the call to the GUI! - m_HasOptionsGUI = false; - // If true, you have to derive a class from MGUIOptions (use MGUIOptionsTemplate) - // and implement all your GUI options -} + m_HasOptionsGUI = true; + + // Allow the use of multiple threads and instances + m_AllowMultiThreading = true; + m_AllowMultipleInstances = true; + // Initiate the Slow Threshold Cut variables + m_SlowThresholdCutMode = MSlowThresholdCutModes::e_Ignore; + m_SlowThresholdCutFixedValue = 15; + m_SlowThresholdCutFileName = ""; + + // Nearest Neighbor threshold (-1000 because we don't want to default to cut neg values) + m_NearestNeighborCutMode = MNearestNeighborCutModes::e_Ignore; + m_NearestNeighborThreshold = -1000.0; + +} //////////////////////////////////////////////////////////////////////////////// @@ -86,20 +112,441 @@ MModuleEnergyCalibration::~MModuleEnergyCalibration() //////////////////////////////////////////////////////////////////////////////// +void MModuleEnergyCalibration::CreateExpos() +{ + // If they are already created, return + if (m_Expos.size() != 0) { + return; + } + + // Set the histogram display + m_ExpoSpectrum = new MGUIExpoPlotSpectrum(this); + m_ExpoSpectrum->SetEnergyHistogramParameters(200, 0, 2000); + m_Expos.push_back(m_ExpoSpectrum); + +} + + +//////////////////////////////////////////////////////////////////////////////// + + bool MModuleEnergyCalibration::Initialize() { - // Initialize the module + // Initialize the module - return true; + // Parse the energy calibration file + if (ReadEnergyCalibrationFile(m_FileName) == false) { + return false; + } + + // Parse the slow threshold file + if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { + if (ReadSlowThresholdCutFile(m_SlowThresholdCutFileName) == false) { + return false; + } + } + + return MModule::Initialize(); } //////////////////////////////////////////////////////////////////////////////// -bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) +bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) { - // Main data analysis routine, which updates the event to a new level + // Main data analysis routine, which updates the event to a new level, i.e. takes the raw ADC value from the .roa file loaded through nuclearizer and converts it into energy units. + + for (unsigned int i = 0; i < Event->GetNStripHits();) { + + MStripHit* SH = Event->GetStripHit(i); + MReadOutElementDoubleStrip R = *dynamic_cast(SH->GetReadOutElement()); + + TF1* Fit = m_Calibration[R]; + TF1* FitRes = m_ResolutionCalibration[R]; + + if (Fit == nullptr) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Energy-fit not found for read-out element " << R << endl; + } + Event->SetEnergyCalibrationError("calibration not found for " + R.ToString()); + ++i; // iterate to next SH + continue; + + } else { + + double Energy = 0; + double Threshold; // No default value here, we set it below (different for strips and nearest neighbors) + + Energy = Fit->Eval(SH->GetADCUnits()); + + // TODO(@RobinAnthonyPetersen): Determine if we want to force negative energy values to be zero or not + // If the calibrated energy is less than 0, force it to be 0. + //if (Energy < 0) { + // Energy = 0; + //} + + if (SH->IsNearestNeighbor() == true) { + // If nothing is selected, we set thresholds to zero + // Threshold = 0.0; // But keeping negative values for now + // If nothing is selected, we keep negative nearest neighbors + Threshold = -100.0; + + // Otherwise, if the user inputs a value, we use that threshold + if (m_NearestNeighborCutMode == MNearestNeighborCutModes::e_Fixed) { + // Get the value user typed in the box (for example 6.0 keV) + // TODO(@RobinAnthonyPetersen): Nearest Neighbor threhsold cut subject to change pending more analysis + Threshold = m_NearestNeighborThreshold; + } + } else { // If not Nearest Neighbor, then it's a triggered strip + + // Set the default slow threshold + Threshold = 0.0; + + if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_Fixed) { // check if user input threshold is enabled (one value applied to all strips) + Threshold = m_SlowThresholdCutFixedValue; + } else if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { // check if threshold file is enabled (unique value applied to each strip) + double ThresholdFromFile = m_ThresholdMap[R]; // if file enabled, declare value from map + + if (ThresholdFromFile == 0) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Threshold not found for read-out element " << R << endl; + } + const double DefaultSlowThreshold = 15.0; + Threshold = DefaultSlowThreshold; // set default threshold if threshold not found + } else { + Threshold = ThresholdFromFile; // set threshold variable to value found in map + } + } + } + + + // Remove SH for any energy value below the established threshold (0 is default) + if (Energy < Threshold) { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Strip Hit below threshold, deleting SH with Energy " << Energy << " keV " << endl; + } + Event->SetStripHitBelowThreshold_QualityFlag("Strip hit removed with energy " + to_string(Energy)); + Event->RemoveStripHit(i); + delete SH; + continue; // continue to next SH without iterating i + } else { + + // if the energy isn't filtered out with the threshold, then assign the energy to the SH + ++i; // iterate to next SH + SH->SetEnergy(Energy); + + if (FitRes == nullptr) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Energy Resolution fit not found for read-out element " << R << endl; + } + // There is not expected to be a time in which the energy resolution calibration is not defined when the energy calibration itself is. Therefore, don't need a seperate BD flag for this. + } else { + double EnergyResolution = FitRes->Eval(Energy); + SH->SetEnergyResolution(EnergyResolution); + } + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyFinal(Energy, SH->IsNearestNeighbor(), SH->IsLowVoltageStrip()); + } + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << ": Energy: " << SH->GetADCUnits() << " adc --> " << Energy << " keV" << endl; + } + } + } + } + Event->SetAnalysisProgress(MAssembly::c_EnergyCalibration); + + return true; +} + + +///////////////////////////////////////////////////////////////////////////////// + + +double MModuleEnergyCalibration::GetEnergy(MReadOutElementDoubleStrip R, double ADC) +{ + //! Return the energy for a given ADC value or zero in case of error + + TF1* Fit = m_Calibration[R]; + double Energy = 0.0; + if (Fit != nullptr) { + Energy = Fit->Eval(ADC); + if (Energy < 0.0) { + Energy = 0.0; + } + } else { + cout << m_Name << ": GetEnergy: Error unable to find calibration" << endl; + return 0; + } + + return Energy; +} + + +/////////////////////////////////////////////////////////////////////////////// + + +double MModuleEnergyCalibration::GetADC(MReadOutElementDoubleStrip R, double Energy) +{ + //! Return the ADC value for a given energy + + TF1* Fit = m_Calibration[R]; // TF1* is a function to be applied + if (Fit != nullptr) { + return Fit->GetX(Energy); + } else { + cout << m_Name << ": GetADC: Error unable to find calibration" << endl; + return 0; + } +} + + +/////////////////////////////////////////////////////////////////////////////// + + +void MModuleEnergyCalibration::Finalize() +{ + // Finalize the calibrator and clean up + + MModule::Finalize(); + + for (auto& F : m_Calibration) { + delete F.second; + } + for (auto& F : m_ResolutionCalibration) { + delete F.second; + } + + return; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleEnergyCalibration::ReadEnergyCalibrationFile(MString FileName) { + + // Parse the Energy Calibration file + MParser Parser; + if (Parser.Open(FileName, MFile::c_Read) == false) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Unable to open energy calibration file " << FileName << endl; + } + return false; + } + + // create other maps + map CP_ROEToLine; //Peak fits + map CM_ROEToLine; //Energy Calibration Model + map CR_ROEToLine; //Energy Resolution Calibration Model + + // tokenize ecal file + for (unsigned int i = 0; i < Parser.GetNLines(); ++i) { + unsigned int NTokens = Parser.GetTokenizerAt(i)->GetNTokens(); + if (NTokens < 2) { + continue; + } + if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true || + Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true || Parser.GetTokenizerAt(i)->IsTokenAt(0, "CR") == true) { + if (Parser.GetTokenizerAt(i)->IsTokenAt(1, "dss") == true) { + + // input token values to map + MReadOutElementDoubleStrip R; + R.SetDetectorID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(2)); + R.SetStripID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(3)); + R.IsLowVoltageStrip((Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "p") || (Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "l")); + if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true) { + CP_ROEToLine[R] = i; + } else if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true) { + CM_ROEToLine[R] = i; + } else { + CR_ROEToLine[R] = i; + } + } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown read-out element (" << Parser.GetTokenizerAt(i)->GetTokenAt(1) << ")" << endl; + } + return false; + } + } + } + + for (auto CM : CM_ROEToLine) { + // If we have at least three data points, we store the calibration + + if (CP_ROEToLine.find(CM.first) != CP_ROEToLine.end()) { + unsigned int i = CP_ROEToLine[CM.first]; + if (Parser.GetTokenizerAt(i)->IsTokenAt(5, "pakw") == false) { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Unknown calibration point descriptor found: " << Parser.GetTokenizerAt(i)->GetTokenAt(5) << endl; + } + continue; + } + } else { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": No good calibration for the following strip found: " << CM.first << endl; + } + continue; + } + + // Read the calibrator, i.e. read the fit function from the .ecal Melinator file. + + unsigned int Pos = 5; + MString CalibratorType = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsString(Pos); + CalibratorType.ToLower(); + + // Below inclusion of poly1zero written by J. Beechert on 2019/11/15 + if (CalibratorType == "poly1zero") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly1zero", "0. + [0]*x", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } + + // Below inclusion of poly1 and poly2 written by J. Beechert on 2019/10/24 + else if (CalibratorType == "poly1") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly1", "[0] + [1]*x", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } else if (CalibratorType == "poly2") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly2", "[0] + [1]*x + [2]*x^2", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + melinatorfit->FixParameter(2, a2); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } + // Eventually, I'll be including other possible fits, but for now, we've just include poly3 and poly4 + else if (CalibratorType == "poly3") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly3", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + melinatorfit->FixParameter(2, a2); + melinatorfit->FixParameter(3, a3); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } else if (CalibratorType == "poly4") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a4 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly4", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + melinatorfit->FixParameter(2, a2); + melinatorfit->FixParameter(3, a3); + melinatorfit->FixParameter(4, a4); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown calibrator type (" << CalibratorType << ") for strip" << CM.first << endl; + } + continue; + } + } + + for (auto CR : CR_ROEToLine) { + + unsigned int Pos = 5; + MString CalibratorType = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsString(Pos); + CalibratorType.ToLower(); + if (CalibratorType == "p1" || CalibratorType == "poly1") { + double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + TF1* resolutionfit = new TF1("P1", "([0]+[1]*x) / 2.355", 0., 2000.); + resolutionfit->FixParameter(0, f0); + resolutionfit->FixParameter(1, f1); + + m_ResolutionCalibration[CR.first] = resolutionfit; + } else if (CalibratorType == "p2" || CalibratorType == "poly2") { + double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + double f2 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + TF1* resolutionfit = new TF1("P2", "([0]+[1]*x+[2]*x*x) / 2.355", 0., 2000.); + resolutionfit->FixParameter(0, f0); + resolutionfit->FixParameter(1, f1); + resolutionfit->FixParameter(2, f2); + m_ResolutionCalibration[CR.first] = resolutionfit; + } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown resolution calibrator type (" << CalibratorType << ") for strip" << CR.first << endl; + } + continue; + } + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleEnergyCalibration::ReadSlowThresholdCutFile(MString FileName) { + + MParser Parser; + if (Parser.Open(FileName, MFile::c_Read) == false) { + if (g_Verbosity >= c_Error) { + cout< Tokens = Line.Tokenize(","); // for each line, Create tokens seperated by commas + if (Tokens.size() == 6) { + int IndexOffset = Tokens.size() % 6; //index counter + int DetID = Tokens[1 + IndexOffset].ToInt(); // Detector ID + MString Side = Tokens[2 + IndexOffset].ToString(); // side is a string, either 'l' or 'h' + int StripID = Tokens[3 + IndexOffset].ToInt(); // stripID + //int ThresholdADC = Tokens[4 + IndexOffset].ToInt(); // energy threshold in ADC + double ThresholdKeVFile = Tokens[5 + IndexOffset].ToDouble(); //energy threshold in keV + + MReadOutElementDoubleStrip R; + R.SetDetectorID(DetID); + R.SetStripID(StripID); + R.IsLowVoltageStrip(Side == "l"); + + // map detectorID, strip number, and voltage side to the threshold (keV) + m_ThresholdMap[R] = ThresholdKeVFile; + } + } + } return true; } @@ -110,7 +557,88 @@ bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) void MModuleEnergyCalibration::ShowOptionsGUI() { - // Show the options GUI - or do nothing + // Show the options GUI + + MGUIOptionsEnergyCalibration* Options = new MGUIOptionsEnergyCalibration(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleEnergyCalibration::ReadXmlConfiguration(MXmlNode* Node) +{ + //! Read the configuration data from an XML node + + MXmlNode* FileNameNode = Node->GetNode("FileName"); + if (FileNameNode != nullptr) { + m_FileName = FileNameNode->GetValue(); + } + + MXmlNode* SlowThresholdCutModeNode = Node->GetNode("SlowThresholdCutMode"); + if (SlowThresholdCutModeNode != nullptr) { + m_SlowThresholdCutMode = static_cast(SlowThresholdCutModeNode->GetValueAsInt()); + } + MXmlNode* SlowThresholdCutFixedValueNode = Node->GetNode("SlowThresholdCutFixedValue"); + if (SlowThresholdCutFixedValueNode != nullptr) { + m_SlowThresholdCutFixedValue = SlowThresholdCutFixedValueNode->GetValueAsDouble(); + } + MXmlNode* SlowThresholdCutFileNameNode = Node->GetNode("SlowThresholdCutThresholdFileName"); + if (SlowThresholdCutFileNameNode != nullptr) { + m_SlowThresholdCutFileName = SlowThresholdCutFileNameNode->GetValue(); + } + MXmlNode* NNCutModeNode = Node->GetNode("NearestNeighborCutMode"); + if (NNCutModeNode != nullptr) { + m_NearestNeighborCutMode = static_cast(NNCutModeNode->GetValueAsInt()); + } + MXmlNode* NearestNeighborThresholdNode = Node->GetNode("NearestNeighborThreshold"); + if (NearestNeighborThresholdNode != nullptr) { + m_NearestNeighborThreshold = NearestNeighborThresholdNode->GetValueAsDouble(); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleEnergyCalibration::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + new MXmlNode(Node, "FileName", m_FileName); + new MXmlNode(Node, "SlowThresholdCutMode", static_cast(m_SlowThresholdCutMode)); + new MXmlNode(Node, "SlowThresholdCutFixedValue", m_SlowThresholdCutFixedValue); + new MXmlNode(Node, "SlowThresholdCutThresholdFileName", m_SlowThresholdCutFileName); + new MXmlNode(Node, "NearestNeighborCutMode", static_cast(m_NearestNeighborCutMode)); + new MXmlNode(Node, "NearestNeighborThreshold", m_NearestNeighborThreshold); + + return Node; +} + +///////////////////////////////////////////////////////////////////////////////// + + +double MModuleEnergyCalibration::LookupEnergyResolution(MStripHit* SH, double Energy) +{ + //! Return the energy resolution or -1 in case of error + + MReadOutElementDoubleStrip* ROE = dynamic_cast(SH->GetReadOutElement()); + if (ROE == nullptr) { + cout << m_Name << ": LookupEnergyResolution: Error unable to get read-out element" << endl; + return -1; + } + TF1* FitRes = m_ResolutionCalibration[*ROE]; + if (FitRes == nullptr) { + cout << m_Name << ": LookupEnergyResolutio: Error: Couldn't locate energy resolution" << endl; + return -1.0; + } else { + return FitRes->Eval(Energy); + } } From b331670650cfc2a968743e5b7922fc0f222bc1a4 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 17:18:09 -0400 Subject: [PATCH 05/10] Reintroduce missing changes. --- src/MModuleEnergyCalibration.cxx | 562 +------------------------------ 1 file changed, 17 insertions(+), 545 deletions(-) diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index c35817d0..b57d4762 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -2,14 +2,12 @@ * MModuleEnergyCalibration.cxx * * - * Copyright (C) by Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, - * Jackie Beechert, Robin Anthony-Peterson. + * Copyright (C) 2008-2008 by Andreas Zoglauer. * All rights reserved. * * * This code implementation is the intellectual property of - * Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, - * Jackie Beechert, Robin Anthony-Peterson + * Andreas Zoglauer. * * By copying, distributing or modifying the Program (or any work * based on the Program) you indicate your acceptance of this statement, @@ -24,29 +22,16 @@ // //////////////////////////////////////////////////////////////////////////////// -// Standard libs: -#include -#include -#include -using namespace std; // Include the header: #include "MModuleEnergyCalibration.h" +// Standard libs: // ROOT libs: #include "TGClient.h" -#include "TFile.h" // MEGAlib libs: -#include "MString.h" -#include "MStreams.h" - -// Nuclearizer libs: -#include "MReadOutElement.h" -#include "MReadOutElementDoubleStrip.h" -#include "MGUIOptionsEnergyCalibration.h" -#include "MGUIExpoPlotSpectrum.h" //////////////////////////////////////////////////////////////////////////////// @@ -67,39 +52,28 @@ MModuleEnergyCalibration::MModuleEnergyCalibration() : MModule() // Set all module relevant information // Set the module name --- has to be unique - m_Name = "Energy calibrator"; + m_Name = "Combined energy calibration and charge sharing correction"; // Set the XML tag --- has to be unique --- no spaces allowed - m_XmlTag = "EnergyCalibration"; + m_XmlTag = "CombinedEnergyCalibrationAndChargeSharingCorrection"; // Set all modules, which have to be done before this module - AddPreceedingModuleType(MAssembly::c_EventLoader); - // AddPreceedingModuleType(MAssembly::c_TACcut); + AddModuleType(MAssembly::c_EventLoader); // Set all types this modules handles AddModuleType(MAssembly::c_EnergyCalibration); + AddModuleType(MAssembly::c_ChargeSharingCorrection); // Set all modules, which can follow this module - AddSucceedingModuleType(MAssembly::c_TACcut); // Set if this module has an options GUI - m_HasOptionsGUI = true; - - // Allow the use of multiple threads and instances - m_AllowMultiThreading = true; - m_AllowMultipleInstances = true; - - // Initiate the Slow Threshold Cut variables - m_SlowThresholdCutMode = MSlowThresholdCutModes::e_Ignore; - m_SlowThresholdCutFixedValue = 15; - m_SlowThresholdCutFileName = ""; - - // Nearest Neighbor threshold (-1000 because we don't want to default to cut neg values) - m_NearestNeighborCutMode = MNearestNeighborCutModes::e_Ignore; - m_NearestNeighborThreshold = -1000.0; - + // If true, overwrite ShowOptionsGUI() with the call to the GUI! + m_HasOptionsGUI = false; + // If true, you have to derive a class from MGUIOptions (use MGUIOptionsTemplate) + // and implement all your GUI options } + //////////////////////////////////////////////////////////////////////////////// @@ -112,401 +86,9 @@ MModuleEnergyCalibration::~MModuleEnergyCalibration() //////////////////////////////////////////////////////////////////////////////// -void MModuleEnergyCalibration::CreateExpos() -{ - // If they are already created, return - if (m_Expos.size() != 0) { - return; - } - - // Set the histogram display - m_ExpoSpectrum = new MGUIExpoPlotSpectrum(this); - m_ExpoSpectrum->SetEnergyHistogramParameters(200, 0, 2000); - m_Expos.push_back(m_ExpoSpectrum); - -} - - -//////////////////////////////////////////////////////////////////////////////// - - bool MModuleEnergyCalibration::Initialize() { - // Initialize the module - - // Parse the energy calibration file - if (ReadEnergyCalibrationFile(m_FileName) == false) { - return false; - } - - // Parse the slow threshold file - if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { - if (ReadSlowThresholdCutFile(m_SlowThresholdCutFileName) == false) { - return false; - } - } - - return MModule::Initialize(); -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) -{ - // Main data analysis routine, which updates the event to a new level, i.e. takes the raw ADC value from the .roa file loaded through nuclearizer and converts it into energy units. - - for (unsigned int i = 0; i < Event->GetNStripHits();) { - - MStripHit* SH = Event->GetStripHit(i); - MReadOutElementDoubleStrip R = *dynamic_cast(SH->GetReadOutElement()); - - TF1* Fit = m_Calibration[R]; - TF1* FitRes = m_ResolutionCalibration[R]; - - if (Fit == nullptr) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Error: Energy-fit not found for read-out element " << R << endl; - } - Event->SetEnergyCalibrationError("calibration not found for " + R.ToString()); - ++i; // iterate to next SH - continue; - - } else { - - double Energy = 0; - double Threshold; // No default value here, we set it below (different for strips and nearest neighbors) - - Energy = Fit->Eval(SH->GetADCUnits()); - - // TODO(@RobinAnthonyPetersen): Determine if we want to force negative energy values to be zero or not - // If the calibrated energy is less than 0, force it to be 0. - //if (Energy < 0) { - // Energy = 0; - //} - - if (SH->IsNearestNeighbor() == true) { - // If nothing is selected, we set thresholds to zero - // Threshold = 0.0; // But keeping negative values for now - // If nothing is selected, we keep negative nearest neighbors - Threshold = -100.0; - - // Otherwise, if the user inputs a value, we use that threshold - if (m_NearestNeighborCutMode == MNearestNeighborCutModes::e_Fixed) { - // Get the value user typed in the box (for example 6.0 keV) - // TODO(@RobinAnthonyPetersen): Nearest Neighbor threhsold cut subject to change pending more analysis - Threshold = m_NearestNeighborThreshold; - } - } else { // If not Nearest Neighbor, then it's a triggered strip - - // Set the default slow threshold - Threshold = 0.0; - - if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_Fixed) { // check if user input threshold is enabled (one value applied to all strips) - Threshold = m_SlowThresholdCutFixedValue; - } else if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { // check if threshold file is enabled (unique value applied to each strip) - double ThresholdFromFile = m_ThresholdMap[R]; // if file enabled, declare value from map - - if (ThresholdFromFile == 0) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Error: Threshold not found for read-out element " << R << endl; - } - const double DefaultSlowThreshold = 15.0; - Threshold = DefaultSlowThreshold; // set default threshold if threshold not found - } else { - Threshold = ThresholdFromFile; // set threshold variable to value found in map - } - } - } - - - // Remove SH for any energy value below the established threshold (0 is default) - if (Energy < Threshold) { - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": Strip Hit below threshold, deleting SH with Energy " << Energy << " keV " << endl; - } - Event->SetStripHitBelowThreshold_QualityFlag("Strip hit removed with energy " + to_string(Energy)); - Event->RemoveStripHit(i); - delete SH; - continue; // continue to next SH without iterating i - } else { - - // if the energy isn't filtered out with the threshold, then assign the energy to the SH - ++i; // iterate to next SH - SH->SetEnergy(Energy); - - if (FitRes == nullptr) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Error: Energy Resolution fit not found for read-out element " << R << endl; - } - // There is not expected to be a time in which the energy resolution calibration is not defined when the energy calibration itself is. Therefore, don't need a seperate BD flag for this. - } else { - double EnergyResolution = FitRes->Eval(Energy); - SH->SetEnergyResolution(EnergyResolution); - } - if (HasExpos() == true) { - m_ExpoSpectrum->AddEnergyFinal(Energy, SH->IsNearestNeighbor(), SH->IsLowVoltageStrip()); - } - if (g_Verbosity >= c_Info) { - cout << m_XmlTag << ": Energy: " << SH->GetADCUnits() << " adc --> " << Energy << " keV" << endl; - } - } - } - } - Event->SetAnalysisProgress(MAssembly::c_EnergyCalibration); - - return true; -} - - -///////////////////////////////////////////////////////////////////////////////// - - -double MModuleEnergyCalibration::GetEnergy(MReadOutElementDoubleStrip R, double ADC) -{ - //! Return the energy for a given ADC value or zero in case of error - - TF1* Fit = m_Calibration[R]; - double Energy = 0.0; - if (Fit != nullptr) { - Energy = Fit->Eval(ADC); - if (Energy < 0.0) { - Energy = 0.0; - } - } else { - cout << m_Name << ": GetEnergy: Error unable to find calibration" << endl; - return 0; - } - - return Energy; -} - - -/////////////////////////////////////////////////////////////////////////////// - - -double MModuleEnergyCalibration::GetADC(MReadOutElementDoubleStrip R, double Energy) -{ - //! Return the ADC value for a given energy - - TF1* Fit = m_Calibration[R]; // TF1* is a function to be applied - if (Fit != nullptr) { - return Fit->GetX(Energy); - } else { - cout << m_Name << ": GetADC: Error unable to find calibration" << endl; - return 0; - } -} - - -/////////////////////////////////////////////////////////////////////////////// - - -void MModuleEnergyCalibration::Finalize() -{ - // Finalize the calibrator and clean up - - MModule::Finalize(); - - for (auto& F : m_Calibration) { - delete F.second; - } - for (auto& F : m_ResolutionCalibration) { - delete F.second; - } - - return; -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MModuleEnergyCalibration::ReadEnergyCalibrationFile(MString FileName) { - - // Parse the Energy Calibration file - MParser Parser; - if (Parser.Open(FileName, MFile::c_Read) == false) { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Unable to open energy calibration file " << FileName << endl; - } - return false; - } - - // create other maps - map CP_ROEToLine; //Peak fits - map CM_ROEToLine; //Energy Calibration Model - map CR_ROEToLine; //Energy Resolution Calibration Model - - // tokenize ecal file - for (unsigned int i = 0; i < Parser.GetNLines(); ++i) { - unsigned int NTokens = Parser.GetTokenizerAt(i)->GetNTokens(); - if (NTokens < 2) { - continue; - } - if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true || - Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true || Parser.GetTokenizerAt(i)->IsTokenAt(0, "CR") == true) { - if (Parser.GetTokenizerAt(i)->IsTokenAt(1, "dss") == true) { - - // input token values to map - MReadOutElementDoubleStrip R; - R.SetDetectorID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(2)); - R.SetStripID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(3)); - R.IsLowVoltageStrip((Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "p") || (Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "l")); - if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true) { - CP_ROEToLine[R] = i; - } else if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true) { - CM_ROEToLine[R] = i; - } else { - CR_ROEToLine[R] = i; - } - } else { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Line parser: Unknown read-out element (" << Parser.GetTokenizerAt(i)->GetTokenAt(1) << ")" << endl; - } - return false; - } - } - } - - for (auto CM : CM_ROEToLine) { - // If we have at least three data points, we store the calibration - - if (CP_ROEToLine.find(CM.first) != CP_ROEToLine.end()) { - unsigned int i = CP_ROEToLine[CM.first]; - if (Parser.GetTokenizerAt(i)->IsTokenAt(5, "pakw") == false) { - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": Unknown calibration point descriptor found: " << Parser.GetTokenizerAt(i)->GetTokenAt(5) << endl; - } - continue; - } - } else { - if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": No good calibration for the following strip found: " << CM.first << endl; - } - continue; - } - - // Read the calibrator, i.e. read the fit function from the .ecal Melinator file. - - unsigned int Pos = 5; - MString CalibratorType = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsString(Pos); - CalibratorType.ToLower(); - - // Below inclusion of poly1zero written by J. Beechert on 2019/11/15 - if (CalibratorType == "poly1zero") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly1zero", "0. + [0]*x", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } - - // Below inclusion of poly1 and poly2 written by J. Beechert on 2019/10/24 - else if (CalibratorType == "poly1") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly1", "[0] + [1]*x", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } else if (CalibratorType == "poly2") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly2", "[0] + [1]*x + [2]*x^2", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - melinatorfit->FixParameter(2, a2); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } - // Eventually, I'll be including other possible fits, but for now, we've just include poly3 and poly4 - else if (CalibratorType == "poly3") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly3", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - melinatorfit->FixParameter(2, a2); - melinatorfit->FixParameter(3, a3); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } else if (CalibratorType == "poly4") { - double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - double a4 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - - // From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly4", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4", 0., m_MaxADCRange); - melinatorfit->FixParameter(0, a0); - melinatorfit->FixParameter(1, a1); - melinatorfit->FixParameter(2, a2); - melinatorfit->FixParameter(3, a3); - melinatorfit->FixParameter(4, a4); - - // Define the map by saving the fit function I just created as a map to the current ReadOutElement - m_Calibration[CM.first] = melinatorfit; - - } else { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Line parser: Unknown calibrator type (" << CalibratorType << ") for strip" << CM.first << endl; - } - continue; - } - } - - for (auto CR : CR_ROEToLine) { - - unsigned int Pos = 5; - MString CalibratorType = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsString(Pos); - CalibratorType.ToLower(); - if (CalibratorType == "p1" || CalibratorType == "poly1") { - double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - TF1* resolutionfit = new TF1("P1", "([0]+[1]*x) / 2.355", 0., 2000.); - resolutionfit->FixParameter(0, f0); - resolutionfit->FixParameter(1, f1); - - m_ResolutionCalibration[CR.first] = resolutionfit; - } else if (CalibratorType == "p2" || CalibratorType == "poly2") { - double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - double f2 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - TF1* resolutionfit = new TF1("P2", "([0]+[1]*x+[2]*x*x) / 2.355", 0., 2000.); - resolutionfit->FixParameter(0, f0); - resolutionfit->FixParameter(1, f1); - resolutionfit->FixParameter(2, f2); - m_ResolutionCalibration[CR.first] = resolutionfit; - } else { - if (g_Verbosity >= c_Error) { - cout << m_XmlTag << ": Line parser: Unknown resolution calibrator type (" << CalibratorType << ") for strip" << CR.first << endl; - } - continue; - } - } + // Initialize the module return true; } @@ -515,88 +97,9 @@ bool MModuleEnergyCalibration::ReadEnergyCalibrationFile(MString FileName) { //////////////////////////////////////////////////////////////////////////////// -bool MModuleEnergyCalibration::ReadSlowThresholdCutFile(MString FileName) { - - MParser Parser; - if (Parser.Open(FileName, MFile::c_Read) == false) { - if (g_Verbosity >= c_Error) { - cout< Tokens = Line.Tokenize(","); // for each line, Create tokens seperated by commas - if (Tokens.size() == 6) { - int IndexOffset = Tokens.size() % 6; //index counter - int DetID = Tokens[1 + IndexOffset].ToInt(); // Detector ID - MString Side = Tokens[2 + IndexOffset].ToString(); // side is a string, either 'l' or 'h' - int StripID = Tokens[3 + IndexOffset].ToInt(); // stripID - //int ThresholdADC = Tokens[4 + IndexOffset].ToInt(); // energy threshold in ADC - double ThresholdKeVFile = Tokens[5 + IndexOffset].ToDouble(); //energy threshold in keV - - MReadOutElementDoubleStrip R; - R.SetDetectorID(DetID); - R.SetStripID(StripID); - R.IsLowVoltageStrip(Side == "l"); - - // map detectorID, strip number, and voltage side to the threshold (keV) - m_ThresholdMap[R] = ThresholdKeVFile; - } - } - } - - return true; -} - - -//////////////////////////////////////////////////////////////////////////////// - - -void MModuleEnergyCalibration::ShowOptionsGUI() +bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) { - // Show the options GUI - - MGUIOptionsEnergyCalibration* Options = new MGUIOptionsEnergyCalibration(this); - Options->Create(); - gClient->WaitForUnmap(Options); -} - - -//////////////////////////////////////////////////////////////////////////////// - - -bool MModuleEnergyCalibration::ReadXmlConfiguration(MXmlNode* Node) -{ - //! Read the configuration data from an XML node - - MXmlNode* FileNameNode = Node->GetNode("FileName"); - if (FileNameNode != nullptr) { - m_FileName = FileNameNode->GetValue(); - } - - MXmlNode* SlowThresholdCutModeNode = Node->GetNode("SlowThresholdCutMode"); - if (SlowThresholdCutModeNode != nullptr) { - m_SlowThresholdCutMode = static_cast(SlowThresholdCutModeNode->GetValueAsInt()); - } - MXmlNode* SlowThresholdCutFixedValueNode = Node->GetNode("SlowThresholdCutFixedValue"); - if (SlowThresholdCutFixedValueNode != nullptr) { - m_SlowThresholdCutFixedValue = SlowThresholdCutFixedValueNode->GetValueAsDouble(); - } - MXmlNode* SlowThresholdCutFileNameNode = Node->GetNode("SlowThresholdCutThresholdFileName"); - if (SlowThresholdCutFileNameNode != nullptr) { - m_SlowThresholdCutFileName = SlowThresholdCutFileNameNode->GetValue(); - } - MXmlNode* NNCutModeNode = Node->GetNode("NearestNeighborCutMode"); - if (NNCutModeNode != nullptr) { - m_NearestNeighborCutMode = static_cast(NNCutModeNode->GetValueAsInt()); - } - MXmlNode* NearestNeighborThresholdNode = Node->GetNode("NearestNeighborThreshold"); - if (NearestNeighborThresholdNode != nullptr) { - m_NearestNeighborThreshold = NearestNeighborThresholdNode->GetValueAsDouble(); - } + // Main data analysis routine, which updates the event to a new level return true; } @@ -605,40 +108,9 @@ bool MModuleEnergyCalibration::ReadXmlConfiguration(MXmlNode* Node) //////////////////////////////////////////////////////////////////////////////// -MXmlNode* MModuleEnergyCalibration::CreateXmlConfiguration() -{ - //! Create an XML node tree from the configuration - - MXmlNode* Node = new MXmlNode(0, m_XmlTag); - new MXmlNode(Node, "FileName", m_FileName); - new MXmlNode(Node, "SlowThresholdCutMode", static_cast(m_SlowThresholdCutMode)); - new MXmlNode(Node, "SlowThresholdCutFixedValue", m_SlowThresholdCutFixedValue); - new MXmlNode(Node, "SlowThresholdCutThresholdFileName", m_SlowThresholdCutFileName); - new MXmlNode(Node, "NearestNeighborCutMode", static_cast(m_NearestNeighborCutMode)); - new MXmlNode(Node, "NearestNeighborThreshold", m_NearestNeighborThreshold); - - return Node; -} - -///////////////////////////////////////////////////////////////////////////////// - - -double MModuleEnergyCalibration::LookupEnergyResolution(MStripHit* SH, double Energy) +void MModuleEnergyCalibration::ShowOptionsGUI() { - //! Return the energy resolution or -1 in case of error - - MReadOutElementDoubleStrip* ROE = dynamic_cast(SH->GetReadOutElement()); - if (ROE == nullptr) { - cout << m_Name << ": LookupEnergyResolution: Error unable to get read-out element" << endl; - return -1; - } - TF1* FitRes = m_ResolutionCalibration[*ROE]; - if (FitRes == nullptr) { - cout << m_Name << ": LookupEnergyResolutio: Error: Couldn't locate energy resolution" << endl; - return -1.0; - } else { - return FitRes->Eval(Energy); - } + // Show the options GUI - or do nothing } From 1b0968f5cff0f9a60c27483f9cb5db6a34084e94 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 17:22:05 -0400 Subject: [PATCH 06/10] This has been great git practice --- src/MModuleEnergyCalibration.cxx | 572 ++++++++++++++++++++++++++++++- 1 file changed, 555 insertions(+), 17 deletions(-) diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index b57d4762..c95ae9ba 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -2,12 +2,14 @@ * MModuleEnergyCalibration.cxx * * - * Copyright (C) 2008-2008 by Andreas Zoglauer. + * Copyright (C) by Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, + * Jackie Beechert, Robin Anthony-Peterson. * All rights reserved. * * * This code implementation is the intellectual property of - * Andreas Zoglauer. + * Andreas Zoglauer, Mark Bandstra, Carolyn Kierans, + * Jackie Beechert, Robin Anthony-Peterson * * By copying, distributing or modifying the Program (or any work * based on the Program) you indicate your acceptance of this statement, @@ -22,16 +24,29 @@ // //////////////////////////////////////////////////////////////////////////////// +// Standard libs: +#include +#include +#include +using namespace std; // Include the header: #include "MModuleEnergyCalibration.h" -// Standard libs: // ROOT libs: #include "TGClient.h" +#include "TFile.h" // MEGAlib libs: +#include "MString.h" +#include "MStreams.h" + +// Nuclearizer libs: +#include "MReadOutElement.h" +#include "MReadOutElementDoubleStrip.h" +#include "MGUIOptionsEnergyCalibration.h" +#include "MGUIExpoPlotSpectrum.h" //////////////////////////////////////////////////////////////////////////////// @@ -52,27 +67,38 @@ MModuleEnergyCalibration::MModuleEnergyCalibration() : MModule() // Set all module relevant information // Set the module name --- has to be unique - m_Name = "Combined energy calibration and charge sharing correction"; + m_Name = "Energy calibrator"; // Set the XML tag --- has to be unique --- no spaces allowed - m_XmlTag = "CombinedEnergyCalibrationAndChargeSharingCorrection"; + m_XmlTag = "EnergyCalibration"; // Set all modules, which have to be done before this module - AddModuleType(MAssembly::c_EventLoader); + AddPreceedingModuleType(MAssembly::c_EventLoader); + // AddPreceedingModuleType(MAssembly::c_TACcut); // Set all types this modules handles AddModuleType(MAssembly::c_EnergyCalibration); - AddModuleType(MAssembly::c_ChargeSharingCorrection); // Set all modules, which can follow this module + AddSucceedingModuleType(MAssembly::c_TACcut); // Set if this module has an options GUI - // If true, overwrite ShowOptionsGUI() with the call to the GUI! - m_HasOptionsGUI = false; - // If true, you have to derive a class from MGUIOptions (use MGUIOptionsTemplate) - // and implement all your GUI options -} + m_HasOptionsGUI = true; + + // Allow the use of multiple threads and instances + m_AllowMultiThreading = true; + m_AllowMultipleInstances = true; + // Initiate the Slow Threshold Cut variables + m_SlowThresholdCutMode = MSlowThresholdCutModes::e_Ignore; + m_SlowThresholdCutFixedValue = 15; + m_SlowThresholdCutFileName = ""; + + // Nearest Neighbor threshold (-1000 because we don't want to default to cut neg values) + m_NearestNeighborCutMode = MNearestNeighborCutModes::e_Ignore; + m_NearestNeighborThreshold = -1000.0; + +} //////////////////////////////////////////////////////////////////////////////// @@ -86,20 +112,451 @@ MModuleEnergyCalibration::~MModuleEnergyCalibration() //////////////////////////////////////////////////////////////////////////////// +void MModuleEnergyCalibration::CreateExpos() +{ + // If they are already created, return + if (m_Expos.size() != 0) { + return; + } + + // Set the histogram display + m_ExpoSpectrum = new MGUIExpoPlotSpectrum(this); + m_ExpoSpectrum->SetEnergyHistogramParameters(200, 0, 2000); + m_Expos.push_back(m_ExpoSpectrum); + +} + + +//////////////////////////////////////////////////////////////////////////////// + + bool MModuleEnergyCalibration::Initialize() { - // Initialize the module + // Initialize the module - return true; + // Parse the energy calibration file + if (ReadEnergyCalibrationFile(m_FileName) == false) { + return false; + } + + // Parse the slow threshold file + if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { + if (ReadSlowThresholdCutFile(m_SlowThresholdCutFileName) == false) { + return false; + } + } + + return MModule::Initialize(); } //////////////////////////////////////////////////////////////////////////////// -bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) +bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) { - // Main data analysis routine, which updates the event to a new level + // Main data analysis routine, which updates the event to a new level, i.e. takes the raw ADC value from the .roa file loaded through nuclearizer and converts it into energy units. + + for (unsigned int i = 0; i < Event->GetNStripHits();) { + + MStripHit* SH = Event->GetStripHit(i); + MReadOutElementDoubleStrip R = *dynamic_cast(SH->GetReadOutElement()); + + // Flag strip hits whose ADC value is close to the ADC saturation limit, since their + // calibrated energy is not trustworthy. The hit is kept, it is only marked. + if (SH->GetADCUnits() > m_HighADCThreshold) { + SH->HasHighADC(true); + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Warning: High ADC value " << SH->GetADCUnits() << " for read-out element " << R << endl; + } + Event->SetHighADC_QualityFlag("High ADC value " + to_string(SH->GetADCUnits()) + " for " + R.ToString().Data()); + } + + TF1* Fit = m_Calibration[R]; + TF1* FitRes = m_ResolutionCalibration[R]; + + if (Fit == nullptr) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Energy-fit not found for read-out element " << R << endl; + } + Event->SetEnergyCalibrationError("calibration not found for " + R.ToString()); + ++i; // iterate to next SH + continue; + + } else { + + double Energy = 0; + double Threshold; // No default value here, we set it below (different for strips and nearest neighbors) + + Energy = Fit->Eval(SH->GetADCUnits()); + + // TODO(@RobinAnthonyPetersen): Determine if we want to force negative energy values to be zero or not + // If the calibrated energy is less than 0, force it to be 0. + //if (Energy < 0) { + // Energy = 0; + //} + + if (SH->IsNearestNeighbor() == true) { + // If nothing is selected, we set thresholds to zero + // Threshold = 0.0; // But keeping negative values for now + // If nothing is selected, we keep negative nearest neighbors + Threshold = -100.0; + + // Otherwise, if the user inputs a value, we use that threshold + if (m_NearestNeighborCutMode == MNearestNeighborCutModes::e_Fixed) { + // Get the value user typed in the box (for example 6.0 keV) + // TODO(@RobinAnthonyPetersen): Nearest Neighbor threhsold cut subject to change pending more analysis + Threshold = m_NearestNeighborThreshold; + } + } else { // If not Nearest Neighbor, then it's a triggered strip + + // Set the default slow threshold + Threshold = 0.0; + + if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_Fixed) { // check if user input threshold is enabled (one value applied to all strips) + Threshold = m_SlowThresholdCutFixedValue; + } else if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { // check if threshold file is enabled (unique value applied to each strip) + double ThresholdFromFile = m_ThresholdMap[R]; // if file enabled, declare value from map + + if (ThresholdFromFile == 0) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Threshold not found for read-out element " << R << endl; + } + const double DefaultSlowThreshold = 15.0; + Threshold = DefaultSlowThreshold; // set default threshold if threshold not found + } else { + Threshold = ThresholdFromFile; // set threshold variable to value found in map + } + } + } + + + // Remove SH for any energy value below the established threshold (0 is default) + if (Energy < Threshold) { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Strip Hit below threshold, deleting SH with Energy " << Energy << " keV " << endl; + } + Event->SetStripHitBelowThreshold_QualityFlag("Strip hit removed with energy " + to_string(Energy)); + Event->RemoveStripHit(i); + delete SH; + continue; // continue to next SH without iterating i + } else { + + // if the energy isn't filtered out with the threshold, then assign the energy to the SH + ++i; // iterate to next SH + SH->SetEnergy(Energy); + + if (FitRes == nullptr) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Energy Resolution fit not found for read-out element " << R << endl; + } + // There is not expected to be a time in which the energy resolution calibration is not defined when the energy calibration itself is. Therefore, don't need a seperate BD flag for this. + } else { + double EnergyResolution = FitRes->Eval(Energy); + SH->SetEnergyResolution(EnergyResolution); + } + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyFinal(Energy, SH->IsNearestNeighbor(), SH->IsLowVoltageStrip()); + } + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << ": Energy: " << SH->GetADCUnits() << " adc --> " << Energy << " keV" << endl; + } + } + } + } + Event->SetAnalysisProgress(MAssembly::c_EnergyCalibration); + + return true; +} + + +///////////////////////////////////////////////////////////////////////////////// + + +double MModuleEnergyCalibration::GetEnergy(MReadOutElementDoubleStrip R, double ADC) +{ + //! Return the energy for a given ADC value or zero in case of error + + TF1* Fit = m_Calibration[R]; + double Energy = 0.0; + if (Fit != nullptr) { + Energy = Fit->Eval(ADC); + if (Energy < 0.0) { + Energy = 0.0; + } + } else { + cout << m_Name << ": GetEnergy: Error unable to find calibration" << endl; + return 0; + } + + return Energy; +} + + +/////////////////////////////////////////////////////////////////////////////// + + +double MModuleEnergyCalibration::GetADC(MReadOutElementDoubleStrip R, double Energy) +{ + //! Return the ADC value for a given energy + + TF1* Fit = m_Calibration[R]; // TF1* is a function to be applied + if (Fit != nullptr) { + return Fit->GetX(Energy); + } else { + cout << m_Name << ": GetADC: Error unable to find calibration" << endl; + return 0; + } +} + + +/////////////////////////////////////////////////////////////////////////////// + + +void MModuleEnergyCalibration::Finalize() +{ + // Finalize the calibrator and clean up + + MModule::Finalize(); + + for (auto& F : m_Calibration) { + delete F.second; + } + for (auto& F : m_ResolutionCalibration) { + delete F.second; + } + + return; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleEnergyCalibration::ReadEnergyCalibrationFile(MString FileName) { + + // Parse the Energy Calibration file + MParser Parser; + if (Parser.Open(FileName, MFile::c_Read) == false) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Unable to open energy calibration file " << FileName << endl; + } + return false; + } + + // create other maps + map CP_ROEToLine; //Peak fits + map CM_ROEToLine; //Energy Calibration Model + map CR_ROEToLine; //Energy Resolution Calibration Model + + // tokenize ecal file + for (unsigned int i = 0; i < Parser.GetNLines(); ++i) { + unsigned int NTokens = Parser.GetTokenizerAt(i)->GetNTokens(); + if (NTokens < 2) { + continue; + } + if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true || + Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true || Parser.GetTokenizerAt(i)->IsTokenAt(0, "CR") == true) { + if (Parser.GetTokenizerAt(i)->IsTokenAt(1, "dss") == true) { + + // input token values to map + MReadOutElementDoubleStrip R; + R.SetDetectorID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(2)); + R.SetStripID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(3)); + R.IsLowVoltageStrip((Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "p") || (Parser.GetTokenizerAt(i)->GetTokenAtAsString(4) == "l")); + if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true) { + CP_ROEToLine[R] = i; + } else if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true) { + CM_ROEToLine[R] = i; + } else { + CR_ROEToLine[R] = i; + } + } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown read-out element (" << Parser.GetTokenizerAt(i)->GetTokenAt(1) << ")" << endl; + } + return false; + } + } + } + + for (auto CM : CM_ROEToLine) { + // If we have at least three data points, we store the calibration + + if (CP_ROEToLine.find(CM.first) != CP_ROEToLine.end()) { + unsigned int i = CP_ROEToLine[CM.first]; + if (Parser.GetTokenizerAt(i)->IsTokenAt(5, "pakw") == false) { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Unknown calibration point descriptor found: " << Parser.GetTokenizerAt(i)->GetTokenAt(5) << endl; + } + continue; + } + } else { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": No good calibration for the following strip found: " << CM.first << endl; + } + continue; + } + + // Read the calibrator, i.e. read the fit function from the .ecal Melinator file. + + unsigned int Pos = 5; + MString CalibratorType = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsString(Pos); + CalibratorType.ToLower(); + + // Below inclusion of poly1zero written by J. Beechert on 2019/11/15 + if (CalibratorType == "poly1zero") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly1zero", "0. + [0]*x", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } + + // Below inclusion of poly1 and poly2 written by J. Beechert on 2019/10/24 + else if (CalibratorType == "poly1") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly1", "[0] + [1]*x", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } else if (CalibratorType == "poly2") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly2", "[0] + [1]*x + [2]*x^2", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + melinatorfit->FixParameter(2, a2); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } + // Eventually, I'll be including other possible fits, but for now, we've just include poly3 and poly4 + else if (CalibratorType == "poly3") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly3", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + melinatorfit->FixParameter(2, a2); + melinatorfit->FixParameter(3, a3); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } else if (CalibratorType == "poly4") { + double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + double a4 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); + + // From the fit parameters I just extracted from the .ecal file, I can define a function + TF1* melinatorfit = new TF1("poly4", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4", 0., m_MaxADCRange); + melinatorfit->FixParameter(0, a0); + melinatorfit->FixParameter(1, a1); + melinatorfit->FixParameter(2, a2); + melinatorfit->FixParameter(3, a3); + melinatorfit->FixParameter(4, a4); + + // Define the map by saving the fit function I just created as a map to the current ReadOutElement + m_Calibration[CM.first] = melinatorfit; + + } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown calibrator type (" << CalibratorType << ") for strip" << CM.first << endl; + } + continue; + } + } + + for (auto CR : CR_ROEToLine) { + + unsigned int Pos = 5; + MString CalibratorType = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsString(Pos); + CalibratorType.ToLower(); + if (CalibratorType == "p1" || CalibratorType == "poly1") { + double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + TF1* resolutionfit = new TF1("P1", "([0]+[1]*x) / 2.355", 0., 2000.); + resolutionfit->FixParameter(0, f0); + resolutionfit->FixParameter(1, f1); + + m_ResolutionCalibration[CR.first] = resolutionfit; + } else if (CalibratorType == "p2" || CalibratorType == "poly2") { + double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + double f2 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); + TF1* resolutionfit = new TF1("P2", "([0]+[1]*x+[2]*x*x) / 2.355", 0., 2000.); + resolutionfit->FixParameter(0, f0); + resolutionfit->FixParameter(1, f1); + resolutionfit->FixParameter(2, f2); + m_ResolutionCalibration[CR.first] = resolutionfit; + } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown resolution calibrator type (" << CalibratorType << ") for strip" << CR.first << endl; + } + continue; + } + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleEnergyCalibration::ReadSlowThresholdCutFile(MString FileName) { + + MParser Parser; + if (Parser.Open(FileName, MFile::c_Read) == false) { + if (g_Verbosity >= c_Error) { + cout< Tokens = Line.Tokenize(","); // for each line, Create tokens seperated by commas + if (Tokens.size() == 6) { + int IndexOffset = Tokens.size() % 6; //index counter + int DetID = Tokens[1 + IndexOffset].ToInt(); // Detector ID + MString Side = Tokens[2 + IndexOffset].ToString(); // side is a string, either 'l' or 'h' + int StripID = Tokens[3 + IndexOffset].ToInt(); // stripID + //int ThresholdADC = Tokens[4 + IndexOffset].ToInt(); // energy threshold in ADC + double ThresholdKeVFile = Tokens[5 + IndexOffset].ToDouble(); //energy threshold in keV + + MReadOutElementDoubleStrip R; + R.SetDetectorID(DetID); + R.SetStripID(StripID); + R.IsLowVoltageStrip(Side == "l"); + + // map detectorID, strip number, and voltage side to the threshold (keV) + m_ThresholdMap[R] = ThresholdKeVFile; + } + } + } return true; } @@ -110,7 +567,88 @@ bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) void MModuleEnergyCalibration::ShowOptionsGUI() { - // Show the options GUI - or do nothing + // Show the options GUI + + MGUIOptionsEnergyCalibration* Options = new MGUIOptionsEnergyCalibration(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleEnergyCalibration::ReadXmlConfiguration(MXmlNode* Node) +{ + //! Read the configuration data from an XML node + + MXmlNode* FileNameNode = Node->GetNode("FileName"); + if (FileNameNode != nullptr) { + m_FileName = FileNameNode->GetValue(); + } + + MXmlNode* SlowThresholdCutModeNode = Node->GetNode("SlowThresholdCutMode"); + if (SlowThresholdCutModeNode != nullptr) { + m_SlowThresholdCutMode = static_cast(SlowThresholdCutModeNode->GetValueAsInt()); + } + MXmlNode* SlowThresholdCutFixedValueNode = Node->GetNode("SlowThresholdCutFixedValue"); + if (SlowThresholdCutFixedValueNode != nullptr) { + m_SlowThresholdCutFixedValue = SlowThresholdCutFixedValueNode->GetValueAsDouble(); + } + MXmlNode* SlowThresholdCutFileNameNode = Node->GetNode("SlowThresholdCutThresholdFileName"); + if (SlowThresholdCutFileNameNode != nullptr) { + m_SlowThresholdCutFileName = SlowThresholdCutFileNameNode->GetValue(); + } + MXmlNode* NNCutModeNode = Node->GetNode("NearestNeighborCutMode"); + if (NNCutModeNode != nullptr) { + m_NearestNeighborCutMode = static_cast(NNCutModeNode->GetValueAsInt()); + } + MXmlNode* NearestNeighborThresholdNode = Node->GetNode("NearestNeighborThreshold"); + if (NearestNeighborThresholdNode != nullptr) { + m_NearestNeighborThreshold = NearestNeighborThresholdNode->GetValueAsDouble(); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleEnergyCalibration::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + new MXmlNode(Node, "FileName", m_FileName); + new MXmlNode(Node, "SlowThresholdCutMode", static_cast(m_SlowThresholdCutMode)); + new MXmlNode(Node, "SlowThresholdCutFixedValue", m_SlowThresholdCutFixedValue); + new MXmlNode(Node, "SlowThresholdCutThresholdFileName", m_SlowThresholdCutFileName); + new MXmlNode(Node, "NearestNeighborCutMode", static_cast(m_NearestNeighborCutMode)); + new MXmlNode(Node, "NearestNeighborThreshold", m_NearestNeighborThreshold); + + return Node; +} + +///////////////////////////////////////////////////////////////////////////////// + + +double MModuleEnergyCalibration::LookupEnergyResolution(MStripHit* SH, double Energy) +{ + //! Return the energy resolution or -1 in case of error + + MReadOutElementDoubleStrip* ROE = dynamic_cast(SH->GetReadOutElement()); + if (ROE == nullptr) { + cout << m_Name << ": LookupEnergyResolution: Error unable to get read-out element" << endl; + return -1; + } + TF1* FitRes = m_ResolutionCalibration[*ROE]; + if (FitRes == nullptr) { + cout << m_Name << ": LookupEnergyResolutio: Error: Couldn't locate energy resolution" << endl; + return -1.0; + } else { + return FitRes->Eval(Energy); + } } From 6611bd50b0d678da37fce7353eeea83fccc7b713 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 17:55:58 -0400 Subject: [PATCH 07/10] Change to int value in message --- src/MModuleEnergyCalibration.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index c95ae9ba..20238cd6 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -169,7 +169,7 @@ bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) if (g_Verbosity >= c_Warning) { cout << m_XmlTag << ": Warning: High ADC value " << SH->GetADCUnits() << " for read-out element " << R << endl; } - Event->SetHighADC_QualityFlag("High ADC value " + to_string(SH->GetADCUnits()) + " for " + R.ToString().Data()); + Event->SetHighADC_QualityFlag("High ADC value " + to_string((int) SH->GetADCUnits()) + " for " + R.ToString().Data()); } TF1* Fit = m_Calibration[R]; From 4ffe39401ede68638e44d967622e8a1e8d530d6c Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 18:42:42 -0400 Subject: [PATCH 08/10] Updated the unit test reference files to include new flag for highADC. --- resource/unittestdata/406-1/hdf5-to-tra.reference.tra | 10 ++++++++++ resource/unittestdata/542-1/hdf5-to-tra.reference.tra | 3 +++ 2 files changed, 13 insertions(+) diff --git a/resource/unittestdata/406-1/hdf5-to-tra.reference.tra b/resource/unittestdata/406-1/hdf5-to-tra.reference.tra index 45d243a8..2dd9e630 100644 --- a/resource/unittestdata/406-1/hdf5-to-tra.reference.tra +++ b/resource/unittestdata/406-1/hdf5-to-tra.reference.tra @@ -1337,6 +1337,7 @@ SE ID 171 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 16.052900) (Strip hit removed with energy 11.981073) (Strip hit removed with energy 8.771925) (Strip hit removed with energy 7.845102) (Strip hit removed with energy 8.958687) (Strip hit removed with energy 14.364213) (Strip hit removed with energy 15.618454) (Strip hit removed with energy -28.297535) (Strip hit removed with energy -7.770818) +QA HighADC (High ADC value 15804.000000 for Detector: 0, side: LV, strip: 20) (High ADC value 16013.000000 for Detector: 0, side: LV, strip: 21) (High ADC value 16006.000000 for Detector: 0, side: LV, strip: 22) (High ADC value 15449.000000 for Detector: 0, side: LV, strip: 23) (High ADC value 15058.000000 for Detector: 0, side: HV, strip: 53) (High ADC value 16015.000000 for Detector: 0, side: HV, strip: 54) (High ADC value 16038.000000 for Detector: 0, side: HV, strip: 55) QA StripPairing (Best reduced chi square is not below 25 (466.978427)) (GR Hit: Detector ID 0 and Energy 294.916417) BD GR Veto PQ 466.978 @@ -9626,6 +9627,7 @@ ID 1253 BD StripPairingError (More than maximum number of strip hits allowed on one side (7)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 6.961407) (Strip hit removed with energy 8.933764) (Strip hit removed with energy 17.695176) (Strip hit removed with energy 11.402859) (Strip hit removed with energy 7.629766) (Strip hit removed with energy 17.062778) (Strip hit removed with energy 13.739803) +QA HighADC (High ADC value 16033.000000 for Detector: 0, side: LV, strip: 42) (High ADC value 15965.000000 for Detector: 0, side: LV, strip: 43) (High ADC value 16077.000000 for Detector: 0, side: HV, strip: 61) BD GR Veto PQ SE @@ -9919,6 +9921,7 @@ ID 1293 BD StripPairingError (More than maximum number of strip hits allowed on one side (8)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 10.624611) (Strip hit removed with energy 18.359129) (Strip hit removed with energy 14.883559) (Strip hit removed with energy 10.931843) (Strip hit removed with energy 8.633496) (Strip hit removed with energy 8.894950) (Strip hit removed with energy 12.278867) (Strip hit removed with energy -0.823052) +QA HighADC (High ADC value 14135.000000 for Detector: 0, side: LV, strip: 1) (High ADC value 16123.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16014.000000 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -24905,6 +24908,7 @@ SE ID 3247 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 7.116343) (Strip hit removed with energy 10.714684) (Strip hit removed with energy 18.330011) (Strip hit removed with energy 15.782019) (Strip hit removed with energy 15.881703) (Strip hit removed with energy 10.092588) (Strip hit removed with energy 7.156331) (Strip hit removed with energy 8.711120) (Strip hit removed with energy 7.829357) (Strip hit removed with energy 8.236825) (Strip hit removed with energy 12.113871) +QA HighADC (High ADC value 14773.000000 for Detector: 0, side: LV, strip: 8) (High ADC value 14062.000000 for Detector: 0, side: LV, strip: 11) (High ADC value 16065.000000 for Detector: 0, side: HV, strip: 2) (High ADC value 14972.000000 for Detector: 0, side: HV, strip: 3) QA StripPairing (Best reduced chi square is not below 25 (478.096823)) (Best strip pairing contains more than 5 strip groupings on one side) (Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: 5332.539902)) (Event contains multiple hits on a single strip) (GR Hit: Detector ID 0 and Energy 181.908526) BD GR Veto PQ 478.097 @@ -28409,6 +28413,7 @@ SE ID 3706 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 7.938531) (Strip hit removed with energy 10.729861) (Strip hit removed with energy 11.572864) (Strip hit removed with energy 13.057852) +QA HighADC (High ADC value 14691.000000 for Detector: 0, side: LV, strip: 64) QA StripPairing (Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: 21.583803)) (GR Hit: Detector ID 0 and Energy 2008.377728) BD GR Veto PQ 0.688171 @@ -36526,6 +36531,7 @@ SE ID 4784 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 8.000643) (Strip hit removed with energy 7.667605) (Strip hit removed with energy 9.765564) (Strip hit removed with energy 7.838833) (Strip hit removed with energy 12.329513) (Strip hit removed with energy 16.020982) (Strip hit removed with energy 17.575923) (Strip hit removed with energy 12.322915) (Strip hit removed with energy 11.493931) (Strip hit removed with energy 7.947770) (Strip hit removed with energy 9.194799) (Strip hit removed with energy 14.312189) (Strip hit removed with energy 16.240531) (Strip hit removed with energy -9.221858) (Strip hit removed with energy 18.540329) +QA HighADC (High ADC value 16123.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16014.000000 for Detector: 0, side: HV, strip: 64) QA StripPairing (Best reduced chi square is not below 25 (42.859504)) (Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: 129.529017)) (Event contains multiple hits on a single strip) (GR Hit: Detector ID 0 and Energy 34.480688) (GR Hit: Detector ID 0 and Energy 68.944180) (GR Hit: Detector ID 0 and Energy 2243.958918) BD GR Veto PQ 42.8595 @@ -38958,6 +38964,7 @@ SE ID 5102 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 10.847868) (Strip hit removed with energy 10.407519) (Strip hit removed with energy 7.814563) (Strip hit removed with energy 7.832260) (Strip hit removed with energy 12.573667) (Strip hit removed with energy 13.002596) (Strip hit removed with energy 13.516989) (Strip hit removed with energy 10.361103) (Strip hit removed with energy 8.906325) (Strip hit removed with energy -10.986370) +QA HighADC (High ADC value 16040.000000 for Detector: 0, side: LV, strip: 58) (High ADC value 16055.000000 for Detector: 0, side: LV, strip: 59) (High ADC value 15997.000000 for Detector: 0, side: LV, strip: 60) (High ADC value 16056.000000 for Detector: 0, side: LV, strip: 61) (High ADC value 15696.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 15813.000000 for Detector: 0, side: HV, strip: 7) (High ADC value 16041.000000 for Detector: 0, side: HV, strip: 8) (High ADC value 16011.000000 for Detector: 0, side: HV, strip: 9) (High ADC value 16021.000000 for Detector: 0, side: HV, strip: 10) (High ADC value 15915.000000 for Detector: 0, side: HV, strip: 12) (High ADC value 15970.000000 for Detector: 0, side: HV, strip: 64) QA StripPairing (Best reduced chi square is not below 25 (1063.315896)) (GR Hit: Detector ID 0 and Energy 2390.281672) BD GR Veto PQ 1063.32 @@ -40532,6 +40539,7 @@ ID 5314 BD StripPairingError (More than maximum number of strip hits allowed on one side (8)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 10.953168) (Strip hit removed with energy 15.073470) (Strip hit removed with energy 15.021121) (Strip hit removed with energy 14.916823) (Strip hit removed with energy 9.691315) (Strip hit removed with energy 8.445107) (Strip hit removed with energy 6.878826) (Strip hit removed with energy 10.837057) +QA HighADC (High ADC value 16042.000000 for Detector: 0, side: LV, strip: 41) (High ADC value 16121.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 15979.000000 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -49355,6 +49363,7 @@ ID 6482 BD StripPairingError (More than maximum number of strip hits allowed on one side (8)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 6.740801) (Strip hit removed with energy 6.303791) (Strip hit removed with energy 8.010620) (Strip hit removed with energy 9.084545) (Strip hit removed with energy 12.185240) (Strip hit removed with energy 13.419191) (Strip hit removed with energy 18.605867) (Strip hit removed with energy 13.114570) +QA HighADC (High ADC value 16123.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16014.000000 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -56122,6 +56131,7 @@ ID 7371 BD StripPairingError (More than maximum number of strip hits allowed on one side (11)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 15.135602) (Strip hit removed with energy 10.825201) (Strip hit removed with energy 10.925057) (Strip hit removed with energy 18.841883) +QA HighADC (High ADC value 15979.000000 for Detector: 0, side: LV, strip: 2) (High ADC value 15990.000000 for Detector: 0, side: LV, strip: 3) (High ADC value 14353.000000 for Detector: 0, side: LV, strip: 4) (High ADC value 14982.000000 for Detector: 0, side: LV, strip: 5) (High ADC value 16012.000000 for Detector: 0, side: HV, strip: 49) (High ADC value 15344.000000 for Detector: 0, side: HV, strip: 50) (High ADC value 14266.000000 for Detector: 0, side: HV, strip: 53) BD GR Veto PQ SE diff --git a/resource/unittestdata/542-1/hdf5-to-tra.reference.tra b/resource/unittestdata/542-1/hdf5-to-tra.reference.tra index 11f86df1..369b236e 100644 --- a/resource/unittestdata/542-1/hdf5-to-tra.reference.tra +++ b/resource/unittestdata/542-1/hdf5-to-tra.reference.tra @@ -15126,6 +15126,7 @@ ID 1880 BD StripPairingError (More than maximum number of strip hits allowed on one side (21)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 19.233721) (Strip hit removed with energy 18.976474) (Strip hit removed with energy 19.234818) (Strip hit removed with energy 19.704374) (Strip hit removed with energy 18.685952) (Strip hit removed with energy 9.259302) (Strip hit removed with energy 10.131297) (Strip hit removed with energy 10.035446) (Strip hit removed with energy 12.924219) (Strip hit removed with energy 14.244382) (Strip hit removed with energy 8.304193) (Strip hit removed with energy 9.575914) +QA HighADC (High ADC value 15253.000000 for Detector: 0, side: LV, strip: 49) (High ADC value 16032.000000 for Detector: 0, side: LV, strip: 50) (High ADC value 16046.000000 for Detector: 0, side: LV, strip: 51) (High ADC value 16039.000000 for Detector: 0, side: LV, strip: 52) (High ADC value 15985.000000 for Detector: 0, side: LV, strip: 53) (High ADC value 16038.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16051.000000 for Detector: 0, side: HV, strip: 56) (High ADC value 15982.000000 for Detector: 0, side: HV, strip: 57) (High ADC value 14447.000000 for Detector: 0, side: HV, strip: 58) (High ADC value 14816.000000 for Detector: 0, side: HV, strip: 59) (High ADC value 15988.000000 for Detector: 0, side: HV, strip: 60) (High ADC value 16019.000000 for Detector: 0, side: HV, strip: 61) (High ADC value 14201.000000 for Detector: 0, side: HV, strip: 63) (High ADC value 16110.000000 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -18912,6 +18913,7 @@ ID 2347 BD StripPairingError (More than maximum number of strip hits allowed on one side (17)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 17.634230) (Strip hit removed with energy -1.711899) (Strip hit removed with energy 19.864592) (Strip hit removed with energy 18.924463) (Strip hit removed with energy 18.547828) (Strip hit removed with energy 10.036853) (Strip hit removed with energy 10.273034) (Strip hit removed with energy 14.155659) (Strip hit removed with energy 15.555762) (Strip hit removed with energy 16.478008) (Strip hit removed with energy 11.143212) (Strip hit removed with energy 8.620146) (Strip hit removed with energy 10.216403) +QA HighADC (High ADC value 15827.000000 for Detector: 0, side: LV, strip: 36) (High ADC value 16044.000000 for Detector: 0, side: LV, strip: 38) (High ADC value 14517.000000 for Detector: 0, side: LV, strip: 39) (High ADC value 15803.000000 for Detector: 0, side: LV, strip: 42) (High ADC value 16049.000000 for Detector: 0, side: LV, strip: 44) (High ADC value 15065.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 15231.000000 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -41219,6 +41221,7 @@ SE ID 5112 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 9.652147) (Strip hit removed with energy 13.842069) (Strip hit removed with energy 11.328372) (Strip hit removed with energy 17.456197) +QA HighADC (High ADC value 15080.000000 for Detector: 0, side: HV, strip: 64) QA StripPairing (Best reduced chi square is not below 25 (151874.330455)) (GR Hit: Detector ID 0 and Energy 2669.457316) BD GR Veto PQ 151874 From 75621d5ec6a1966907977d25bb88c62884148d31 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 20 Aug 2026 19:11:01 -0400 Subject: [PATCH 09/10] Fixed the integer issue in the QA message. I remember way back when this was supposed to be a beautifully written PR --- .../406-1/hdf5-to-tra.reference.tra | 20 +++++++++---------- .../542-1/hdf5-to-tra.reference.tra | 6 +++--- src/MModuleEnergyCalibration.cxx | 1 - 3 files changed, 13 insertions(+), 14 deletions(-) diff --git a/resource/unittestdata/406-1/hdf5-to-tra.reference.tra b/resource/unittestdata/406-1/hdf5-to-tra.reference.tra index 2dd9e630..8f2f34f1 100644 --- a/resource/unittestdata/406-1/hdf5-to-tra.reference.tra +++ b/resource/unittestdata/406-1/hdf5-to-tra.reference.tra @@ -1337,7 +1337,7 @@ SE ID 171 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 16.052900) (Strip hit removed with energy 11.981073) (Strip hit removed with energy 8.771925) (Strip hit removed with energy 7.845102) (Strip hit removed with energy 8.958687) (Strip hit removed with energy 14.364213) (Strip hit removed with energy 15.618454) (Strip hit removed with energy -28.297535) (Strip hit removed with energy -7.770818) -QA HighADC (High ADC value 15804.000000 for Detector: 0, side: LV, strip: 20) (High ADC value 16013.000000 for Detector: 0, side: LV, strip: 21) (High ADC value 16006.000000 for Detector: 0, side: LV, strip: 22) (High ADC value 15449.000000 for Detector: 0, side: LV, strip: 23) (High ADC value 15058.000000 for Detector: 0, side: HV, strip: 53) (High ADC value 16015.000000 for Detector: 0, side: HV, strip: 54) (High ADC value 16038.000000 for Detector: 0, side: HV, strip: 55) +QA HighADC (High ADC value 15804 for Detector: 0, side: LV, strip: 20) (High ADC value 16013 for Detector: 0, side: LV, strip: 21) (High ADC value 16006 for Detector: 0, side: LV, strip: 22) (High ADC value 15449 for Detector: 0, side: LV, strip: 23) (High ADC value 15058 for Detector: 0, side: HV, strip: 53) (High ADC value 16015 for Detector: 0, side: HV, strip: 54) (High ADC value 16038 for Detector: 0, side: HV, strip: 55) QA StripPairing (Best reduced chi square is not below 25 (466.978427)) (GR Hit: Detector ID 0 and Energy 294.916417) BD GR Veto PQ 466.978 @@ -9627,7 +9627,7 @@ ID 1253 BD StripPairingError (More than maximum number of strip hits allowed on one side (7)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 6.961407) (Strip hit removed with energy 8.933764) (Strip hit removed with energy 17.695176) (Strip hit removed with energy 11.402859) (Strip hit removed with energy 7.629766) (Strip hit removed with energy 17.062778) (Strip hit removed with energy 13.739803) -QA HighADC (High ADC value 16033.000000 for Detector: 0, side: LV, strip: 42) (High ADC value 15965.000000 for Detector: 0, side: LV, strip: 43) (High ADC value 16077.000000 for Detector: 0, side: HV, strip: 61) +QA HighADC (High ADC value 16033 for Detector: 0, side: LV, strip: 42) (High ADC value 15965 for Detector: 0, side: LV, strip: 43) (High ADC value 16077 for Detector: 0, side: HV, strip: 61) BD GR Veto PQ SE @@ -9921,7 +9921,7 @@ ID 1293 BD StripPairingError (More than maximum number of strip hits allowed on one side (8)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 10.624611) (Strip hit removed with energy 18.359129) (Strip hit removed with energy 14.883559) (Strip hit removed with energy 10.931843) (Strip hit removed with energy 8.633496) (Strip hit removed with energy 8.894950) (Strip hit removed with energy 12.278867) (Strip hit removed with energy -0.823052) -QA HighADC (High ADC value 14135.000000 for Detector: 0, side: LV, strip: 1) (High ADC value 16123.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16014.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 14135 for Detector: 0, side: LV, strip: 1) (High ADC value 16123 for Detector: 0, side: LV, strip: 64) (High ADC value 16014 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -24908,7 +24908,7 @@ SE ID 3247 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 7.116343) (Strip hit removed with energy 10.714684) (Strip hit removed with energy 18.330011) (Strip hit removed with energy 15.782019) (Strip hit removed with energy 15.881703) (Strip hit removed with energy 10.092588) (Strip hit removed with energy 7.156331) (Strip hit removed with energy 8.711120) (Strip hit removed with energy 7.829357) (Strip hit removed with energy 8.236825) (Strip hit removed with energy 12.113871) -QA HighADC (High ADC value 14773.000000 for Detector: 0, side: LV, strip: 8) (High ADC value 14062.000000 for Detector: 0, side: LV, strip: 11) (High ADC value 16065.000000 for Detector: 0, side: HV, strip: 2) (High ADC value 14972.000000 for Detector: 0, side: HV, strip: 3) +QA HighADC (High ADC value 14773 for Detector: 0, side: LV, strip: 8) (High ADC value 14062 for Detector: 0, side: LV, strip: 11) (High ADC value 16065 for Detector: 0, side: HV, strip: 2) (High ADC value 14972 for Detector: 0, side: HV, strip: 3) QA StripPairing (Best reduced chi square is not below 25 (478.096823)) (Best strip pairing contains more than 5 strip groupings on one side) (Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: 5332.539902)) (Event contains multiple hits on a single strip) (GR Hit: Detector ID 0 and Energy 181.908526) BD GR Veto PQ 478.097 @@ -28413,7 +28413,7 @@ SE ID 3706 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 7.938531) (Strip hit removed with energy 10.729861) (Strip hit removed with energy 11.572864) (Strip hit removed with energy 13.057852) -QA HighADC (High ADC value 14691.000000 for Detector: 0, side: LV, strip: 64) +QA HighADC (High ADC value 14691 for Detector: 0, side: LV, strip: 64) QA StripPairing (Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: 21.583803)) (GR Hit: Detector ID 0 and Energy 2008.377728) BD GR Veto PQ 0.688171 @@ -36531,7 +36531,7 @@ SE ID 4784 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 8.000643) (Strip hit removed with energy 7.667605) (Strip hit removed with energy 9.765564) (Strip hit removed with energy 7.838833) (Strip hit removed with energy 12.329513) (Strip hit removed with energy 16.020982) (Strip hit removed with energy 17.575923) (Strip hit removed with energy 12.322915) (Strip hit removed with energy 11.493931) (Strip hit removed with energy 7.947770) (Strip hit removed with energy 9.194799) (Strip hit removed with energy 14.312189) (Strip hit removed with energy 16.240531) (Strip hit removed with energy -9.221858) (Strip hit removed with energy 18.540329) -QA HighADC (High ADC value 16123.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16014.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 16123 for Detector: 0, side: LV, strip: 64) (High ADC value 16014 for Detector: 0, side: HV, strip: 64) QA StripPairing (Best reduced chi square is not below 25 (42.859504)) (Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: 129.529017)) (Event contains multiple hits on a single strip) (GR Hit: Detector ID 0 and Energy 34.480688) (GR Hit: Detector ID 0 and Energy 68.944180) (GR Hit: Detector ID 0 and Energy 2243.958918) BD GR Veto PQ 42.8595 @@ -38964,7 +38964,7 @@ SE ID 5102 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 10.847868) (Strip hit removed with energy 10.407519) (Strip hit removed with energy 7.814563) (Strip hit removed with energy 7.832260) (Strip hit removed with energy 12.573667) (Strip hit removed with energy 13.002596) (Strip hit removed with energy 13.516989) (Strip hit removed with energy 10.361103) (Strip hit removed with energy 8.906325) (Strip hit removed with energy -10.986370) -QA HighADC (High ADC value 16040.000000 for Detector: 0, side: LV, strip: 58) (High ADC value 16055.000000 for Detector: 0, side: LV, strip: 59) (High ADC value 15997.000000 for Detector: 0, side: LV, strip: 60) (High ADC value 16056.000000 for Detector: 0, side: LV, strip: 61) (High ADC value 15696.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 15813.000000 for Detector: 0, side: HV, strip: 7) (High ADC value 16041.000000 for Detector: 0, side: HV, strip: 8) (High ADC value 16011.000000 for Detector: 0, side: HV, strip: 9) (High ADC value 16021.000000 for Detector: 0, side: HV, strip: 10) (High ADC value 15915.000000 for Detector: 0, side: HV, strip: 12) (High ADC value 15970.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 16040 for Detector: 0, side: LV, strip: 58) (High ADC value 16055 for Detector: 0, side: LV, strip: 59) (High ADC value 15997 for Detector: 0, side: LV, strip: 60) (High ADC value 16056 for Detector: 0, side: LV, strip: 61) (High ADC value 15696 for Detector: 0, side: LV, strip: 64) (High ADC value 15813 for Detector: 0, side: HV, strip: 7) (High ADC value 16041 for Detector: 0, side: HV, strip: 8) (High ADC value 16011 for Detector: 0, side: HV, strip: 9) (High ADC value 16021 for Detector: 0, side: HV, strip: 10) (High ADC value 15915 for Detector: 0, side: HV, strip: 12) (High ADC value 15970 for Detector: 0, side: HV, strip: 64) QA StripPairing (Best reduced chi square is not below 25 (1063.315896)) (GR Hit: Detector ID 0 and Energy 2390.281672) BD GR Veto PQ 1063.32 @@ -40539,7 +40539,7 @@ ID 5314 BD StripPairingError (More than maximum number of strip hits allowed on one side (8)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 10.953168) (Strip hit removed with energy 15.073470) (Strip hit removed with energy 15.021121) (Strip hit removed with energy 14.916823) (Strip hit removed with energy 9.691315) (Strip hit removed with energy 8.445107) (Strip hit removed with energy 6.878826) (Strip hit removed with energy 10.837057) -QA HighADC (High ADC value 16042.000000 for Detector: 0, side: LV, strip: 41) (High ADC value 16121.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 15979.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 16042 for Detector: 0, side: LV, strip: 41) (High ADC value 16121 for Detector: 0, side: LV, strip: 64) (High ADC value 15979 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -49363,7 +49363,7 @@ ID 6482 BD StripPairingError (More than maximum number of strip hits allowed on one side (8)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 6.740801) (Strip hit removed with energy 6.303791) (Strip hit removed with energy 8.010620) (Strip hit removed with energy 9.084545) (Strip hit removed with energy 12.185240) (Strip hit removed with energy 13.419191) (Strip hit removed with energy 18.605867) (Strip hit removed with energy 13.114570) -QA HighADC (High ADC value 16123.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16014.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 16123 for Detector: 0, side: LV, strip: 64) (High ADC value 16014 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -56131,7 +56131,7 @@ ID 7371 BD StripPairingError (More than maximum number of strip hits allowed on one side (11)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 15.135602) (Strip hit removed with energy 10.825201) (Strip hit removed with energy 10.925057) (Strip hit removed with energy 18.841883) -QA HighADC (High ADC value 15979.000000 for Detector: 0, side: LV, strip: 2) (High ADC value 15990.000000 for Detector: 0, side: LV, strip: 3) (High ADC value 14353.000000 for Detector: 0, side: LV, strip: 4) (High ADC value 14982.000000 for Detector: 0, side: LV, strip: 5) (High ADC value 16012.000000 for Detector: 0, side: HV, strip: 49) (High ADC value 15344.000000 for Detector: 0, side: HV, strip: 50) (High ADC value 14266.000000 for Detector: 0, side: HV, strip: 53) +QA HighADC (High ADC value 15979 for Detector: 0, side: LV, strip: 2) (High ADC value 15990 for Detector: 0, side: LV, strip: 3) (High ADC value 14353 for Detector: 0, side: LV, strip: 4) (High ADC value 14982 for Detector: 0, side: LV, strip: 5) (High ADC value 16012 for Detector: 0, side: HV, strip: 49) (High ADC value 15344 for Detector: 0, side: HV, strip: 50) (High ADC value 14266 for Detector: 0, side: HV, strip: 53) BD GR Veto PQ SE diff --git a/resource/unittestdata/542-1/hdf5-to-tra.reference.tra b/resource/unittestdata/542-1/hdf5-to-tra.reference.tra index 369b236e..747a0f08 100644 --- a/resource/unittestdata/542-1/hdf5-to-tra.reference.tra +++ b/resource/unittestdata/542-1/hdf5-to-tra.reference.tra @@ -15126,7 +15126,7 @@ ID 1880 BD StripPairingError (More than maximum number of strip hits allowed on one side (21)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 19.233721) (Strip hit removed with energy 18.976474) (Strip hit removed with energy 19.234818) (Strip hit removed with energy 19.704374) (Strip hit removed with energy 18.685952) (Strip hit removed with energy 9.259302) (Strip hit removed with energy 10.131297) (Strip hit removed with energy 10.035446) (Strip hit removed with energy 12.924219) (Strip hit removed with energy 14.244382) (Strip hit removed with energy 8.304193) (Strip hit removed with energy 9.575914) -QA HighADC (High ADC value 15253.000000 for Detector: 0, side: LV, strip: 49) (High ADC value 16032.000000 for Detector: 0, side: LV, strip: 50) (High ADC value 16046.000000 for Detector: 0, side: LV, strip: 51) (High ADC value 16039.000000 for Detector: 0, side: LV, strip: 52) (High ADC value 15985.000000 for Detector: 0, side: LV, strip: 53) (High ADC value 16038.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 16051.000000 for Detector: 0, side: HV, strip: 56) (High ADC value 15982.000000 for Detector: 0, side: HV, strip: 57) (High ADC value 14447.000000 for Detector: 0, side: HV, strip: 58) (High ADC value 14816.000000 for Detector: 0, side: HV, strip: 59) (High ADC value 15988.000000 for Detector: 0, side: HV, strip: 60) (High ADC value 16019.000000 for Detector: 0, side: HV, strip: 61) (High ADC value 14201.000000 for Detector: 0, side: HV, strip: 63) (High ADC value 16110.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 15253 for Detector: 0, side: LV, strip: 49) (High ADC value 16032 for Detector: 0, side: LV, strip: 50) (High ADC value 16046 for Detector: 0, side: LV, strip: 51) (High ADC value 16039 for Detector: 0, side: LV, strip: 52) (High ADC value 15985 for Detector: 0, side: LV, strip: 53) (High ADC value 16038 for Detector: 0, side: LV, strip: 64) (High ADC value 16051 for Detector: 0, side: HV, strip: 56) (High ADC value 15982 for Detector: 0, side: HV, strip: 57) (High ADC value 14447 for Detector: 0, side: HV, strip: 58) (High ADC value 14816 for Detector: 0, side: HV, strip: 59) (High ADC value 15988 for Detector: 0, side: HV, strip: 60) (High ADC value 16019 for Detector: 0, side: HV, strip: 61) (High ADC value 14201 for Detector: 0, side: HV, strip: 63) (High ADC value 16110 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -18913,7 +18913,7 @@ ID 2347 BD StripPairingError (More than maximum number of strip hits allowed on one side (17)) BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 17.634230) (Strip hit removed with energy -1.711899) (Strip hit removed with energy 19.864592) (Strip hit removed with energy 18.924463) (Strip hit removed with energy 18.547828) (Strip hit removed with energy 10.036853) (Strip hit removed with energy 10.273034) (Strip hit removed with energy 14.155659) (Strip hit removed with energy 15.555762) (Strip hit removed with energy 16.478008) (Strip hit removed with energy 11.143212) (Strip hit removed with energy 8.620146) (Strip hit removed with energy 10.216403) -QA HighADC (High ADC value 15827.000000 for Detector: 0, side: LV, strip: 36) (High ADC value 16044.000000 for Detector: 0, side: LV, strip: 38) (High ADC value 14517.000000 for Detector: 0, side: LV, strip: 39) (High ADC value 15803.000000 for Detector: 0, side: LV, strip: 42) (High ADC value 16049.000000 for Detector: 0, side: LV, strip: 44) (High ADC value 15065.000000 for Detector: 0, side: LV, strip: 64) (High ADC value 15231.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 15827 for Detector: 0, side: LV, strip: 36) (High ADC value 16044 for Detector: 0, side: LV, strip: 38) (High ADC value 14517 for Detector: 0, side: LV, strip: 39) (High ADC value 15803 for Detector: 0, side: LV, strip: 42) (High ADC value 16049 for Detector: 0, side: LV, strip: 44) (High ADC value 15065 for Detector: 0, side: LV, strip: 64) (High ADC value 15231 for Detector: 0, side: HV, strip: 64) BD GR Veto PQ SE @@ -41221,7 +41221,7 @@ SE ID 5112 BD DepthCalibrationError (GR Veto) QA StripHitBelowThreshold (Strip hit removed with energy 9.652147) (Strip hit removed with energy 13.842069) (Strip hit removed with energy 11.328372) (Strip hit removed with energy 17.456197) -QA HighADC (High ADC value 15080.000000 for Detector: 0, side: HV, strip: 64) +QA HighADC (High ADC value 15080 for Detector: 0, side: HV, strip: 64) QA StripPairing (Best reduced chi square is not below 25 (151874.330455)) (GR Hit: Detector ID 0 and Energy 2669.457316) BD GR Veto PQ 151874 diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index 20238cd6..c132e5fd 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -165,7 +165,6 @@ bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) // Flag strip hits whose ADC value is close to the ADC saturation limit, since their // calibrated energy is not trustworthy. The hit is kept, it is only marked. if (SH->GetADCUnits() > m_HighADCThreshold) { - SH->HasHighADC(true); if (g_Verbosity >= c_Warning) { cout << m_XmlTag << ": Warning: High ADC value " << SH->GetADCUnits() << " for read-out element " << R << endl; } From 8c85189f45cea43e64740c8937ab2bf78bf5031c Mon Sep 17 00:00:00 2001 From: ckierans <33991471+ckierans@users.noreply.github.com> Date: Fri, 21 Aug 2026 11:47:49 -0400 Subject: [PATCH 10/10] Apply suggestions from Felix's code review Co-authored-by: Felix Hagemann --- src/MModuleEnergyCalibration.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MModuleEnergyCalibration.cxx b/src/MModuleEnergyCalibration.cxx index c132e5fd..302badc7 100644 --- a/src/MModuleEnergyCalibration.cxx +++ b/src/MModuleEnergyCalibration.cxx @@ -166,7 +166,7 @@ bool MModuleEnergyCalibration::AnalyzeEvent(MReadOutAssembly* Event) // calibrated energy is not trustworthy. The hit is kept, it is only marked. if (SH->GetADCUnits() > m_HighADCThreshold) { if (g_Verbosity >= c_Warning) { - cout << m_XmlTag << ": Warning: High ADC value " << SH->GetADCUnits() << " for read-out element " << R << endl; + cout << m_XmlTag << ": Warning: High ADC value " << (int) SH->GetADCUnits() << " for read-out element " << R << endl; } Event->SetHighADC_QualityFlag("High ADC value " + to_string((int) SH->GetADCUnits()) + " for " + R.ToString().Data()); }