From 4dd1623c41909fe2c2f79b72abcc426e41b097bc Mon Sep 17 00:00:00 2001 From: vfioretti Date: Wed, 19 Aug 2026 14:37:51 +0200 Subject: [PATCH] DEE Shield Readout module - from energy to ADC for the shield hits --- include/MGUIOptionsDEESMEX.h | 4 +- include/MModuleDEESMEX.h | 11 ++ include/MReadOutElementVoxel3D.h | 27 ++- include/MSubModuleShieldReadout.h | 69 ++++++- resource/dee/ShieldEnergyCalibration_NRL.csv | 26 +++ src/MGUIOptionsDEESMEX.cxx | 10 +- src/MReadOutElementVoxel3D.cxx | 195 ++++++++++++++----- src/MSubModuleShieldReadout.cxx | 191 ++++++++++++++++-- 8 files changed, 453 insertions(+), 80 deletions(-) create mode 100644 resource/dee/ShieldEnergyCalibration_NRL.csv diff --git a/include/MGUIOptionsDEESMEX.h b/include/MGUIOptionsDEESMEX.h index 55cd6fa6..655ba271 100644 --- a/include/MGUIOptionsDEESMEX.h +++ b/include/MGUIOptionsDEESMEX.h @@ -95,8 +95,10 @@ class MGUIOptionsDEESMEX : public MGUIOptions TGCheckButton* m_StopAfter; //! Entry field for the maximum number of accepted events MGUIEEntry* m_MaximumAcceptedEvents; - //! ACS DEE energy correction file + //! Shield DEE energy correction file MGUIEFileSelector* m_ShieldEnergyCorrectionFileSelector; + //! Shield DEE Readout energy calibration file + MGUIEFileSelector* m_ShieldReadoutFileSelector; //! Select deadtime constants file MGUIEFileSelector* m_DeadtimeFileSelector; diff --git a/include/MModuleDEESMEX.h b/include/MModuleDEESMEX.h index ab78a88e..4a49ffe9 100644 --- a/include/MModuleDEESMEX.h +++ b/include/MModuleDEESMEX.h @@ -137,6 +137,17 @@ class MModuleDEESMEX : public MModule return m_ShieldEnergyCorrection.GetShieldEnergyCorrectionFileName(); } + //! Set shield energy calibration file name + void SetShieldReadoutFileName(const MString& FileName) + { + m_ShieldReadout.SetShieldReadoutFileName(FileName); + } + //! Get shield energy calibration file name + MString GetShieldReadoutFileName() const + { + return m_ShieldReadout.GetShieldReadoutFileName(); + } + //! Set dead time file name void SetDeadtimeFileName(const MString& FileName) { diff --git a/include/MReadOutElementVoxel3D.h b/include/MReadOutElementVoxel3D.h index de4a7c54..51707674 100644 --- a/include/MReadOutElementVoxel3D.h +++ b/include/MReadOutElementVoxel3D.h @@ -30,15 +30,18 @@ //////////////////////////////////////////////////////////////////////////////// -//! The read-out element of a BGO voxel. It reads the detector name and the voxel x, y, z, ID +//! Read-out element identified by detector/crystal ID and optionally by 3D voxel IDs class MReadOutElementVoxel3D : public MReadOutElement { // public interface: public: - //! default constructor - Read out element of a voxel 3D + //! Default constructor MReadOutElementVoxel3D(); - //! full constructor - Read out element of a voxel 3D + //! Constructor for a crystal-only read-out element + MReadOutElementVoxel3D(const MString& DetectorID, unsigned int CrystalID); + + //! Constructor for a fully specified 3D voxel read-out element MReadOutElementVoxel3D(const MString& DetectorID, unsigned int CrystalID, unsigned int VoxelXID, unsigned int VoxelYID, unsigned int VoxelZID); //! Simple default destructor @@ -53,6 +56,9 @@ class MReadOutElementVoxel3D : public MReadOutElement //! Compare two read-out elements virtual bool operator==(const MReadOutElement& R) const; + //! Compare two voxel read-out elements for ordered containers + bool operator<(const MReadOutElementVoxel3D& R) const; + //! Return true if this read-out element is of the given type virtual bool IsOfType(const MString& String) const; //! Return the type of this read-out element @@ -80,6 +86,18 @@ class MReadOutElementVoxel3D : public MReadOutElement return m_CrystalID; } + //! Return true if all voxel IDs are defined + bool HasVoxelIDs() const; + + //! Return true if only detector and crystal IDs are defined + bool IsCrystalOnly() const; + + //! Clear voxel IDs and keep only detector and crystal identification + void ClearVoxelIDs(); + + //! Set all voxel IDs + void SetVoxelIDs(unsigned int VoxelXID, unsigned int VoxelYID, unsigned int VoxelZID); + //! Set Voxel X ID as int void SetVoxelXID(unsigned int VoxelXID) { @@ -96,7 +114,7 @@ class MReadOutElementVoxel3D : public MReadOutElement { m_VoxelYID = VoxelYID; } - //! Set Voxel Y ID as int + //! Get Voxel Y ID as int unsigned int GetVoxelYID() const { return m_VoxelYID; @@ -157,3 +175,4 @@ ostream& operator<<(ostream& os, const MReadOutElementVoxel3D& R); //////////////////////////////////////////////////////////////////////////////// + diff --git a/include/MSubModuleShieldReadout.h b/include/MSubModuleShieldReadout.h index a0a31711..9e697bbc 100644 --- a/include/MSubModuleShieldReadout.h +++ b/include/MSubModuleShieldReadout.h @@ -1,7 +1,7 @@ /* * MSubModuleShieldReadout.h * - * Copyright (C) by Andreas Zoglauer. + * Copyright (C) by Andreas Zoglauer, Valentina Fioretti. * All rights reserved. * * Please see the source-file for the copyright-notice. @@ -17,12 +17,14 @@ // Standard libs: +#include // ROOT libs: // MEGAlib libs: #include "MGlobal.h" #include "MSubModule.h" +#include "MReadOutElementVoxel3D.h" // Forward declarations: @@ -41,21 +43,39 @@ class MSubModuleShieldReadout : public MSubModule MSubModuleShieldReadout(const MSubModuleShieldReadout&) = delete; //! No copy assignment MSubModuleShieldReadout& operator=(const MSubModuleShieldReadout&) = delete; - //! No move constructors + //! No move constructor MSubModuleShieldReadout(MSubModuleShieldReadout&&) = delete; - //! No move operators + //! No move assignment MSubModuleShieldReadout& operator=(MSubModuleShieldReadout&&) = delete; //! Default destructor virtual ~MSubModuleShieldReadout(); + //! Set the shield readout calibration file name + void SetShieldReadoutFileName(const MString& FileName) + { + m_ShieldReadoutFileName = FileName; + } + + //! Set the shield readout calibration file name - deprecated spelling kept for compatibility + void SetShieldRedoutFileName(const MString& FileName) + { + SetShieldReadoutFileName(FileName); + } + + //! Return the shield readout calibration file name + MString GetShieldReadoutFileName() const + { + return m_ShieldReadoutFileName; + } + //! Initialize the module virtual bool Initialize(); //! Clear event data from the module virtual void Clear(); - //! Main data analysis routine, which updates the event to a new level + //! Convert the shield hit energies into ADC values virtual bool AnalyzeEvent(MReadOutAssembly* Event); //! Finalize the module @@ -63,24 +83,60 @@ class MSubModuleShieldReadout : public MSubModule //! Read the configuration data from an XML node virtual bool ReadXmlConfiguration(MXmlNode* Node); + //! Create an XML node tree from the configuration virtual MXmlNode* CreateXmlConfiguration(MXmlNode* Node); // protected methods: protected: + //! Read and parse the shield readout calibration CSV file + bool ParseShieldReadoutFile(); + + //! Compute ADC from the calibration ADC = A*Energy^2 + B*Energy + C + double GetADC(double Energy, const MString& DetectorID, unsigned int CrystalID); // private methods: private: - // protected members: protected: - + //! Shield readout calibration file name + MString m_ShieldReadoutFileName; // private members: private: + //! Calibration data for ADC = A*Energy^2 + B*Energy + C + struct ShieldCalibration + { + //! Quadratic coefficient A + double m_QuadraticCoefficient = 0.0; + //! Linear coefficient B + double m_LinearCoefficient = 0.0; + //! Constant coefficient C + double m_ConstantCoefficient = 0.0; + }; + //! Calibration indexed by crystal-only MReadOutElementVoxel3D keys (detector ID + crystal ID) + std::map m_ADCCalibration; + + //! Number of columns expected in the NRL shield calibration CSV + static constexpr unsigned int m_NumberOfCSVColumns = 16; + + //! CSV column containing the detector/crystal identifier + static constexpr unsigned int m_ChannelIdentifierColumn = 2; + + //! CSV column containing the quadratic calibration coefficient + static constexpr unsigned int m_QuadraticCoefficientColumn = 8; + + //! CSV column containing the linear calibration coefficient + static constexpr unsigned int m_LinearCoefficientColumn = 10; + + //! CSV column containing the constant calibration coefficient + static constexpr unsigned int m_ConstantCoefficientColumn = 12; + + //! Maximum value of the 14-bit ADC + static constexpr double m_MaxADCRange = 16383.0; @@ -95,3 +151,4 @@ class MSubModuleShieldReadout : public MSubModule //////////////////////////////////////////////////////////////////////////////// + diff --git a/resource/dee/ShieldEnergyCalibration_NRL.csv b/resource/dee/ShieldEnergyCalibration_NRL.csv new file mode 100644 index 00000000..2d5182bb --- /dev/null +++ b/resource/dee/ShieldEnergyCalibration_NRL.csv @@ -0,0 +1,26 @@ +Detector Meta Data,,,,,,,,ADC to Energy Fit Parameters (ax^2 + bx + c),,,,,,Resolution, +scb,gamma_channel,sipm_crystal,dac_id,cathod,asic,pannel,measured_cathod_voltage,a,a_err,b,b_err,c,c_err,Resolution @ 662 keV,Resolution Err @ 662 keV +A,8,Z1-0,2,5,0,Z1,27.5,-8.19E-05,5.84E-05,8.50E+00,4.42E-02,1.60E+03,3.47E+00,14.82,0.09 +A,12,Z1-1,3,6,0,Z1,27.47,-3.68E-04,7.08E-04,9.53E+00,5.32E-01,1.52E+03,4.17E+01,14.8,0.67 +A,16,Z1-2,4,7,0,Z1,27.52,1.21E-04,6.21E-04,8.48E+00,4.66E-01,1.57E+03,3.65E+01,14.61,0.61 +A,20,Z1-3,5,8,0,Z1,27.43,-9.25E-05,3.58E-04,8.64E+00,2.69E-01,1.55E+03,2.10E+01,14.82,0.37 +A,8,Z0-0,7,1,1,Z0,27.46,-2.39E-04,3.07E-04,8.87E+00,2.31E-01,1.63E+03,1.80E+01,14.86,0.31 +A,12,Z0-1,6,2,1,Z0,27.47,-3.13E-04,8.99E-04,9.42E+00,6.75E-01,1.50E+03,5.28E+01,14.43,0.83 +A,16,Z0-2,0,3,1,Z0,27.5,-8.08E-04,6.26E-04,1.05E+01,4.72E-01,1.49E+03,3.69E+01,14.77,0.57 +A,20,Z0-3,1,4,1,Z0,27.48,-5.90E-04,8.38E-04,9.86E+00,6.31E-01,1.50E+03,4.94E+01,14.35,0.77 +B,8,Y1-0,2,5,0,Y1,27.5,3.80E-05,2.38E-04,7.87E+00,1.78E-01,1.55E+03,1.40E+01,15.34,0.27 +B,12,Y1-1,3,6,0,Y1,27.48,4.39E-04,3.85E-04,7.20E+00,2.88E-01,1.55E+03,2.27E+01,15.4,0.44 +B,16,Y1-2,4,7,0,Y1,27.48,-8.41E-05,1.60E-04,8.32E+00,1.20E-01,1.45E+03,9.80E+00,15.34,0.18 +B,20,Z1-4,5,8,0,Z1,27.51,-6.54E-04,4.75E-04,1.01E+01,3.60E-01,1.51E+03,2.96E+01,14.89,0.46 +B,8,Y0-0,7,1,1,Y0,27.5,3.00E-04,5.72E-04,7.69E+00,4.31E-01,1.49E+03,4.23E+01,14.5,0.6 +B,12,Y0-1,6,2,1,Y0,27.36,7.64E-04,7.81E-04,6.79E+00,5.87E-01,1.48E+03,5.26E+01,14.27,0.83 +B,16,Y0-2,0,3,1,Y0,27.45,1.13E-04,1.99E-04,7.88E+00,1.49E-01,1.58E+03,1.26E+01,15.42,0.23 +B,20,Z0-4,1,4,1,Z0,27.52,-6.31E-04,5.99E-04,9.85E+00,4.54E-01,1.64E+03,3.61E+01,15.15,0.59 +C,8,X1-0,2,5,0,X1,27.49,2.86E-05,8.83E-05,7.75E+00,6.62E-02,1.50E+03,5.48E+00,15.34,0.12 +C,12,X1-1,3,6,0,X1,27.49,2.27E-04,5.27E-04,7.56E+00,3.97E-01,1.46E+03,3.48E+01,15.1,0.58 +C,16,X1-2,4,7,0,X1,27.49,-4.08E-04,6.27E-04,8.70E+00,4.74E-01,1.34E+03,4.44E+01,15.77,0.7 +C,20,D1-0,5,8,0,D1,27.49,-2.94E-03,7.89E-04,1.45E+01,6.10E-01,1.27E+03,4.86E+01,17.01,0.78 +C,8,X0-0,7,1,1,X0,27.46,-2.47E-03,1.69E-03,1.40E+01,1.34E+00,7.85E+02,1.05E+02,23.7,2.11 +C,12,X0-1,6,2,1,X0,27.51,6.56E-05,1.05E-04,7.87E+00,7.84E-02,1.60E+03,6.11E+00,15.55,0.13 +C,16,X0-2,0,3,1,X0,27.43,-2.11E-03,1.37E-03,1.34E+01,1.07E+00,9.44E+02,8.37E+01,22.13,1.62 +C,20,D0-0,1,4,1,D0,27.42,-3.18E-03,3.46E-04,1.47E+01,3.10E-01,1.38E+03,2.74E+01,15.22,0.37 \ No newline at end of file diff --git a/src/MGUIOptionsDEESMEX.cxx b/src/MGUIOptionsDEESMEX.cxx index 37be2f79..696fb0eb 100644 --- a/src/MGUIOptionsDEESMEX.cxx +++ b/src/MGUIOptionsDEESMEX.cxx @@ -159,6 +159,13 @@ void MGUIOptionsDEESMEX::Create() m_ShieldEnergyCorrectionFileSelector->SetFileType("Shield DEE energy correction file", "*.csv"); m_OptionsFrame->AddFrame(m_ShieldEnergyCorrectionFileSelector, LabelLayout); + // shield energy calibration file + m_ShieldReadoutFileSelector = new MGUIEFileSelector(m_OptionsFrame, + "Please select an energy calibration file for the Shield DEE Readout:", + dynamic_cast(m_Module)->GetShieldReadoutFileName()); + m_ShieldReadoutFileSelector->SetFileType("Shield DEE Readout energy calibration file", "*.csv"); + m_OptionsFrame->AddFrame(m_ShieldReadoutFileSelector, LabelLayout); + /* m_ApplyFudgeFactorSelector = new TGCheckButton(m_OptionsFrame, "Apply fudge factor to better match fluxes", 1); m_ApplyFudgeFactorSelector->SetOn(dynamic_cast(m_Module)->GetApplyFudgeFactor()); @@ -211,7 +218,8 @@ bool MGUIOptionsDEESMEX::OnApply() // Shield options: dynamic_cast(m_Module)->SetShieldEnergyCorrectionFileName(m_ShieldEnergyCorrectionFileSelector->GetFileName()); - + dynamic_cast(m_Module)->SetShieldReadoutFileName(m_ShieldReadoutFileSelector->GetFileName()); + // GeD options: dynamic_cast(m_Module)->SetEnergyCalibrationFileName(m_EnergyCalibrationFileSelector->GetFileName()); dynamic_cast(m_Module)->SetDepthCoefficientsFileName(m_DepthCalibrationCoeffsFileSelector->GetFileName()); diff --git a/src/MReadOutElementVoxel3D.cxx b/src/MReadOutElementVoxel3D.cxx index 191bff77..b9068e55 100644 --- a/src/MReadOutElementVoxel3D.cxx +++ b/src/MReadOutElementVoxel3D.cxx @@ -55,6 +55,23 @@ MReadOutElementVoxel3D::MReadOutElementVoxel3D() { } + +//////////////////////////////////////////////////////////////////////////////// + + +//! Constructor for a crystal-only read-out element +MReadOutElementVoxel3D::MReadOutElementVoxel3D( + const MString& DetectorID, + unsigned int CrystalID) + : MReadOutElement(), + m_DetectorID(DetectorID), + m_CrystalID(CrystalID), + m_VoxelXID(g_UnsignedIntNotDefined), + m_VoxelYID(g_UnsignedIntNotDefined), + m_VoxelZID(g_UnsignedIntNotDefined) +{ +} + //////////////////////////////////////////////////////////////////////////////// //! Parameterized constructor @@ -88,14 +105,14 @@ MReadOutElementVoxel3D::~MReadOutElementVoxel3D() //! Clone this read-out element - the returned element must be deleted! MReadOutElementVoxel3D* MReadOutElementVoxel3D::Clone() const { - MReadOutElementVoxel3D* R = new MReadOutElementVoxel3D( - m_DetectorID, - m_CrystalID, - m_VoxelXID, - m_VoxelYID, - m_VoxelZID); - return R; + if (IsCrystalOnly() == true) { + return new MReadOutElementVoxel3D(m_DetectorID, m_CrystalID); + } + + return new MReadOutElementVoxel3D(m_DetectorID, m_CrystalID, m_VoxelXID, m_VoxelYID, m_VoxelZID); } + + //////////////////////////////////////////////////////////////////////////////// @@ -113,11 +130,61 @@ void MReadOutElementVoxel3D::Clear() //////////////////////////////////////////////////////////////////////////////// +//! Return true if all voxel IDs are defined +bool MReadOutElementVoxel3D::HasVoxelIDs() const +{ + return m_VoxelXID != g_UnsignedIntNotDefined && + m_VoxelYID != g_UnsignedIntNotDefined && + m_VoxelZID != g_UnsignedIntNotDefined; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Return true if only detector and crystal IDs are defined +bool MReadOutElementVoxel3D::IsCrystalOnly() const +{ + return m_VoxelXID == g_UnsignedIntNotDefined && + m_VoxelYID == g_UnsignedIntNotDefined && + m_VoxelZID == g_UnsignedIntNotDefined; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Clear voxel IDs and keep detector and crystal identification +void MReadOutElementVoxel3D::ClearVoxelIDs() +{ + m_VoxelXID = g_UnsignedIntNotDefined; + m_VoxelYID = g_UnsignedIntNotDefined; + m_VoxelZID = g_UnsignedIntNotDefined; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Set all voxel IDs +void MReadOutElementVoxel3D::SetVoxelIDs( + unsigned int VoxelXID, + unsigned int VoxelYID, + unsigned int VoxelZID) +{ + m_VoxelXID = VoxelXID; + m_VoxelYID = VoxelYID; + m_VoxelZID = VoxelZID; +} + + +//////////////////////////////////////////////////////////////////////////////// + + //! Return true if this read-out element is of the given type bool MReadOutElementVoxel3D::IsOfType(const MString& String) const { - if (String == "voxel3d") - return true; + if (String == "voxel3d") return true; return false; } @@ -133,28 +200,31 @@ MString MReadOutElementVoxel3D::GetType() const } +//////////////////////////////////////////////////////////////////////////////// + + //! Test for equality bool MReadOutElementVoxel3D::operator==(const MReadOutElement& R) const { const MReadOutElementVoxel3D* Other = dynamic_cast(&R); if (Other == nullptr) { - if (g_Verbosity >= c_Error) { - cout << "ERROR: Comparison with different read-out element type" << endl; - } + if (g_Verbosity >= c_Error) cout << "ERROR: Comparison with different read-out element type" << endl; return false; } - if (m_DetectorID != Other->m_DetectorID) - return false; - if (m_CrystalID != Other->m_CrystalID) - return false; + if (m_DetectorID != Other->m_DetectorID) return false; + if (m_CrystalID != Other->m_CrystalID) return false; - if (m_VoxelXID != Other->m_VoxelXID) - return false; - if (m_VoxelYID != Other->m_VoxelYID) - return false; - if (m_VoxelZID != Other->m_VoxelZID) - return false; + const bool ThisHasVoxelIDs = HasVoxelIDs(); + const bool OtherHasVoxelIDs = Other->HasVoxelIDs(); + + if (ThisHasVoxelIDs != OtherHasVoxelIDs) return false; + if (IsCrystalOnly() == true && Other->IsCrystalOnly() == true) return true; + if (ThisHasVoxelIDs == false) return false; + + if (m_VoxelXID != Other->m_VoxelXID) return false; + if (m_VoxelYID != Other->m_VoxelYID) return false; + if (m_VoxelZID != Other->m_VoxelZID) return false; return true; } @@ -163,6 +233,32 @@ bool MReadOutElementVoxel3D::operator==(const MReadOutElement& R) const //////////////////////////////////////////////////////////////////////////////// +//! Compare two voxel read-out elements for ordered containers +bool MReadOutElementVoxel3D::operator<(const MReadOutElementVoxel3D& R) const +{ + if (m_DetectorID != R.m_DetectorID) return m_DetectorID < R.m_DetectorID; + if (m_CrystalID != R.m_CrystalID) return m_CrystalID < R.m_CrystalID; + + const bool ThisHasVoxelIDs = HasVoxelIDs(); + const bool OtherHasVoxelIDs = R.HasVoxelIDs(); + + if (ThisHasVoxelIDs != OtherHasVoxelIDs) return ThisHasVoxelIDs < OtherHasVoxelIDs; + if (IsCrystalOnly() == true && R.IsCrystalOnly() == true) return false; + if (ThisHasVoxelIDs == false) { + return m_VoxelXID < R.m_VoxelXID || + (m_VoxelXID == R.m_VoxelXID && m_VoxelYID < R.m_VoxelYID) || + (m_VoxelXID == R.m_VoxelXID && m_VoxelYID == R.m_VoxelYID && m_VoxelZID < R.m_VoxelZID); + } + + if (m_VoxelXID != R.m_VoxelXID) return m_VoxelXID < R.m_VoxelXID; + if (m_VoxelYID != R.m_VoxelYID) return m_VoxelYID < R.m_VoxelYID; + return m_VoxelZID < R.m_VoxelZID; +} + + +//////////////////////////////////////////////////////////////////////////////// + + //! Return the number of parsable elements unsigned int MReadOutElementVoxel3D::GetNumberOfParsableElements() const { @@ -177,9 +273,8 @@ unsigned int MReadOutElementVoxel3D::GetNumberOfParsableElements() const bool MReadOutElementVoxel3D::Parse(const MTokenizer& T, unsigned int StartElement) { if (T.GetNTokens() < StartElement + GetNumberOfParsableElements()) { - if (g_Verbosity >= c_Error) { - cout << "ERROR: Not enough elements to parse. Number of tokens is " << T.GetNTokens() << " and less than 5" << endl; - } + if (g_Verbosity >= c_Error) cout << "ERROR: Not enough elements to parse. Number of tokens is " + << T.GetNTokens() << " and less than 5" << endl; return false; } @@ -190,9 +285,8 @@ bool MReadOutElementVoxel3D::Parse(const MTokenizer& T, unsigned int StartElemen m_VoxelZID = T.GetTokenAtAsUnsignedIntFast(StartElement + 4); if (m_DetectorID == "") { - if (g_Verbosity >= c_Warning) - cout << "WARNING: Parsed empty DetectorID (token index " - << StartElement << ")." << endl; + if (g_Verbosity >= c_Warning) cout << "WARNING: Parsed empty DetectorID (token index " + << StartElement << ")." << endl; } if (m_CrystalID == g_UnsignedIntNotDefined || @@ -200,12 +294,11 @@ bool MReadOutElementVoxel3D::Parse(const MTokenizer& T, unsigned int StartElemen m_VoxelYID == g_UnsignedIntNotDefined || m_VoxelZID == g_UnsignedIntNotDefined) { - if (g_Verbosity >= c_Warning) - cout << "WARNING: Parsed undefined ID(s): " - << "Crystal = " << m_CrystalID - << " Vx = " << m_VoxelXID - << " Vy = " << m_VoxelYID - << " Vz = " << m_VoxelZID << endl; + if (g_Verbosity >= c_Warning) cout << "WARNING: Parsed undefined ID(s): " + << "Crystal = " << m_CrystalID + << " Vx = " << m_VoxelXID + << " Vy = " << m_VoxelYID + << " Vz = " << m_VoxelZID << endl; } return true; @@ -240,33 +333,34 @@ MString MReadOutElementVoxel3D::ToParsableString(bool WithDescriptor) const //! Convert content to a string MString MReadOutElementVoxel3D::ToString() const { + if (m_DetectorID == "" && g_Verbosity >= c_Warning) { + cout << "WARNING: called an element with empty DetectorID" << endl; + } - if (m_DetectorID == "") { - if (g_Verbosity >= c_Warning) { - cout << "WARNING: called an element with empty DetectorID" << endl; - } + if (m_CrystalID == g_UnsignedIntNotDefined && g_Verbosity >= c_Warning) { + cout << "WARNING: called an element with undefined CrystalID" << endl; } - if (m_CrystalID == g_UnsignedIntNotDefined || - m_VoxelXID == g_UnsignedIntNotDefined || - m_VoxelYID == g_UnsignedIntNotDefined || - m_VoxelZID == g_UnsignedIntNotDefined) { + const bool AnyVoxelIDDefined = + m_VoxelXID != g_UnsignedIntNotDefined || + m_VoxelYID != g_UnsignedIntNotDefined || + m_VoxelZID != g_UnsignedIntNotDefined; - if (g_Verbosity >= c_Warning) { - cout << "WARNING: called undefined ID(s): " - << "Crystal = " << m_CrystalID - << " Vx = " << m_VoxelXID - << " Vy = " << m_VoxelYID - << " Vz = " << m_VoxelZID << endl; - } + if (AnyVoxelIDDefined == true && HasVoxelIDs() == false && g_Verbosity >= c_Warning) { + cout << "WARNING: called an element with partially defined voxel IDs" << endl; } ostringstream OS; - OS << "DetectorID: " << m_DetectorID << ", CrystalID: " << m_CrystalID << ", VoxelID: (" << m_VoxelXID << ", " << m_VoxelYID << ", " << m_VoxelZID << ")"; + OS << "DetectorID: " << m_DetectorID << ", CrystalID: " << m_CrystalID; + + if (HasVoxelIDs() == true) { + OS << ", VoxelID: (" << m_VoxelXID << ", " << m_VoxelYID << ", " << m_VoxelZID << ")"; + } return OS.str(); } + //////////////////////////////////////////////////////////////////////////////// @@ -280,3 +374,4 @@ ostream& operator<<(ostream& os, const MReadOutElementVoxel3D& R) // MReadOutElementVoxel3D.cxx: the end... //////////////////////////////////////////////////////////////////////////////// + diff --git a/src/MSubModuleShieldReadout.cxx b/src/MSubModuleShieldReadout.cxx index 765aa1c1..62fe6a13 100644 --- a/src/MSubModuleShieldReadout.cxx +++ b/src/MSubModuleShieldReadout.cxx @@ -2,7 +2,7 @@ * MSubModuleShieldReadout.cxx * * - * Copyright (C) by Andreas Zoglauer. + * Copyright (C) by Andreas Zoglauer, Valentina Fioretti. * All rights reserved. * * @@ -27,10 +27,13 @@ #include "MSubModuleShieldReadout.h" // Standard libs: +#include // ROOT libs: // MEGAlib libs: +#include "MDEECrystalHit.h" +#include "MParser.h" #include "MSubModule.h" @@ -42,6 +45,14 @@ ClassImp(MSubModuleShieldReadout) #endif +constexpr unsigned int MSubModuleShieldReadout::m_NumberOfCSVColumns; +constexpr unsigned int MSubModuleShieldReadout::m_ChannelIdentifierColumn; +constexpr unsigned int MSubModuleShieldReadout::m_QuadraticCoefficientColumn; +constexpr unsigned int MSubModuleShieldReadout::m_LinearCoefficientColumn; +constexpr unsigned int MSubModuleShieldReadout::m_ConstantCoefficientColumn; +constexpr double MSubModuleShieldReadout::m_MaxADCRange; + + //////////////////////////////////////////////////////////////////////////////// @@ -69,6 +80,12 @@ bool MSubModuleShieldReadout::Initialize() { // Initialize the module + if (!ParseShieldReadoutFile()) { + if (g_Verbosity >= c_Error) cout << "ERROR: Failed to parse shield readout calibration file " + << m_ShieldReadoutFileName << endl; + return false; + } + return MSubModule::Initialize(); } @@ -89,13 +106,23 @@ void MSubModuleShieldReadout::Clear() bool MSubModuleShieldReadout::AnalyzeEvent(MReadOutAssembly* Event) { - // Main data analysis routine, which updates the event to a new level + // Convert each shield hit energy into the corresponding ADC value + + std::list& Hits = Event->GetDEECrystalHitListReference(); - // Dummy code: - list& Hits = Event->GetDEECrystalHitListReference(); - for (MDEECrystalHit& CH: Hits) { - CH.m_ADC = 2000 + 4*CH.m_Energy; - if (CH.m_ADC > 16383) CH.m_ADC = 16383; + for (MDEECrystalHit& CrystalHit: Hits) { + CrystalHit.m_ADC = GetADC(CrystalHit.m_Energy, CrystalHit.m_DetectorID, CrystalHit.m_CrystalID); + + if (CrystalHit.m_ADC > m_MaxADCRange) { + if (g_Verbosity >= c_Warning) cout << "Shield ADC value above the maximum range defined as " + << m_MaxADCRange << endl; + CrystalHit.m_ADC = m_MaxADCRange; + } + + if (CrystalHit.m_ADC < 0.0) { + if (g_Verbosity >= c_Warning) cout << "Shield ADC value below 0, set to 0" << endl; + CrystalHit.m_ADC = 0.0; + } } return true; @@ -107,7 +134,7 @@ bool MSubModuleShieldReadout::AnalyzeEvent(MReadOutAssembly* Event) void MSubModuleShieldReadout::Finalize() { - // Finalize the analysis - do all cleanup, i.e., undo Initialize() + // Finalize the analysis - do all cleanup, i.e., undo Initialize() MSubModule::Finalize(); } @@ -116,16 +143,145 @@ void MSubModuleShieldReadout::Finalize() //////////////////////////////////////////////////////////////////////////////// +double MSubModuleShieldReadout::GetADC(double Energy, + const MString& DetectorID, + unsigned int CrystalID) +{ + // The two-argument constructor intentionally creates a detector/crystal-only map key. + const MReadOutElementVoxel3D ReadOutElement(DetectorID, CrystalID); + + const auto CalibrationIterator = m_ADCCalibration.find(ReadOutElement); + if (CalibrationIterator == m_ADCCalibration.end()) { + if (g_Verbosity >= c_Error) cout << "ERROR: Shield ADC calibration not found for detector " + << DetectorID << " and crystal " << CrystalID << endl; + return 0.0; + } + + const double QuadraticCoefficient = CalibrationIterator->second.m_QuadraticCoefficient; + const double LinearCoefficient = CalibrationIterator->second.m_LinearCoefficient; + const double ConstantCoefficient = CalibrationIterator->second.m_ConstantCoefficient; + + // The NRL CSV calibration directly provides ADC as a function of energy: + // ADC = A*Energy^2 + B*Energy + C + const double ADC = QuadraticCoefficient * Energy * Energy + + LinearCoefficient * Energy + + ConstantCoefficient; + + if (g_Verbosity >= c_Info) { + cout << "DEE shield energy-to-ADC conversion:" << endl; + cout << "Channel: " << DetectorID << "-" << CrystalID << endl; + cout << "Energy: " << Energy << endl; + cout << "ADC: " << ADC << endl; + } + + return ADC; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MSubModuleShieldReadout::ParseShieldReadoutFile() +{ + + m_ADCCalibration.clear(); + + if (m_ShieldReadoutFileName == "") { + if (g_Verbosity >= c_Error) cout << "ERROR: Shield readout calibration filename is empty." << endl; + return false; + } + + // The shield readout calibration file is a comma-separated CSV file + MParser Parser(','); + if (Parser.Open(m_ShieldReadoutFileName, MFile::c_Read) == false) { + if (g_Verbosity >= c_Error) cout << "Unable to open shield readout calibration file " << m_ShieldReadoutFileName << endl; + return false; + } + + unsigned int Parsed = 0; + unsigned int Skipped = 0; + + for (unsigned int i = 0; i < Parser.GetNLines(); i++) { + unsigned int NTokens = Parser.GetTokenizerAt(i)->GetNTokens(); + if (NTokens == 0) + continue; // skip empty lines + + // Skip the two header lines of the NRL calibration CSV + MString FirstToken = Parser.GetTokenizerAt(i)->GetTokenAtAsString(0); + if (FirstToken == "Detector Meta Data" || FirstToken == "scb") + continue; + + if (NTokens != m_NumberOfCSVColumns) { + if (g_Verbosity >= c_Warning) { + cout << "WARNING: Line " << i << ": expected " << m_NumberOfCSVColumns + << " tokens, got " << NTokens << " (skipping)" << endl; + } + ++Skipped; + continue; + } // this shouldn't happen but just in case + + // sipm_crystal contains the detector and crystal identifier, e.g. Z1-0 + MString ChannelIdentifier = Parser.GetTokenizerAt(i)->GetTokenAtAsString(m_ChannelIdentifierColumn); + MTokenizer ChannelTokenizer('-', false); + ChannelTokenizer.Analyze(ChannelIdentifier); + + if (ChannelTokenizer.GetNTokens() != 2) { + if (g_Verbosity >= c_Warning) { + cout << "WARNING: Line " << i << ": invalid sipm_crystal value " + << ChannelIdentifier << " (skipping)" << endl; + } + ++Skipped; + continue; + } + + // Detector ID and crystal ID + MString DetectorID = ChannelTokenizer.GetTokenAtAsString(0); + unsigned int CrystalID = ChannelTokenizer.GetTokenAtAsUnsignedInt(1); + + // ADC = a*Energy^2 + b*Energy + c + ShieldCalibration Calibration; + Calibration.m_QuadraticCoefficient = Parser.GetTokenizerAt(i)->GetTokenAtAsDouble(m_QuadraticCoefficientColumn); + Calibration.m_LinearCoefficient = Parser.GetTokenizerAt(i)->GetTokenAtAsDouble(m_LinearCoefficientColumn); + Calibration.m_ConstantCoefficient = Parser.GetTokenizerAt(i)->GetTokenAtAsDouble(m_ConstantCoefficientColumn); + + MReadOutElementVoxel3D ReadOutElement(DetectorID, CrystalID); + + if (m_ADCCalibration.find(ReadOutElement) != m_ADCCalibration.end()) { + if (g_Verbosity >= c_Warning) { + cout << "WARNING: Duplicate shield calibration for detector " << DetectorID + << " and crystal " << CrystalID << " on line " << i + << "; replacing previous entry" << endl; + } + } + + m_ADCCalibration[ReadOutElement] = Calibration; + + ++Parsed; + } + + if (Parsed == 0) { + if (g_Verbosity >= c_Error) cout << "ERROR: Parsed 0 valid calibration lines from " << m_ShieldReadoutFileName << endl; + return false; + } + + if (Skipped > 0) { + if (g_Verbosity >= c_Warning) cout << "WARNING: Skipped " << Skipped << " line(s) while parsing " << m_ShieldReadoutFileName << endl; + } + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// + + bool MSubModuleShieldReadout::ReadXmlConfiguration(MXmlNode* Node) { - //! Read the configuration data from an XML node + // Read the configuration data from an XML node - /* - MXmlNode* SomeTagNode = Node->GetNode("SomeTag"); - if (SomeTagNode != 0) { - m_SomeTagValue = SomeTagNode->GetValue(); + MXmlNode* ShieldReadoutFileName = Node->GetNode("ShieldReadoutFileName"); + if (ShieldReadoutFileName != nullptr) { + m_ShieldReadoutFileName = ShieldReadoutFileName->GetValue(); } - */ return true; } @@ -136,11 +292,9 @@ bool MSubModuleShieldReadout::ReadXmlConfiguration(MXmlNode* Node) MXmlNode* MSubModuleShieldReadout::CreateXmlConfiguration(MXmlNode* Node) { - //! Create an XML node tree from the configuration + // Create an XML node tree from the configuration - /* - MXmlNode* SomeTagNode = new MXmlNode(Node, "SomeTag", "SomeValue"); - */ + new MXmlNode(Node, "ShieldReadoutFileName", m_ShieldReadoutFileName); return Node; } @@ -148,3 +302,4 @@ MXmlNode* MSubModuleShieldReadout::CreateXmlConfiguration(MXmlNode* Node) // MSubModuleShieldReadout.cxx: the end... //////////////////////////////////////////////////////////////////////////////// +