Skip to content

[DEE] Update inverse depth calibration from pixel-based to strip-based #158

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
fhagemann wants to merge 1 commit into
base: develop/em
Choose a base branch
from
Draft

[DEE] Update inverse depth calibration from pixel-based to strip-based #158

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 7 additions & 3 deletions include/MSubModuleDepthReadout.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ class MSubModuleDepthReadout : public MSubModule
protected:

//! Load the splines file
bool LoadSplinesFile(MString FileName);
bool LoadStripCoeffsFile(MString FileName);


// private methods:
Expand All @@ -110,8 +110,12 @@ class MSubModuleDepthReadout : public MSubModule
//! Filename of the depth calibration coefficients (stretch, offset, timing noise, ...)
MString m_DepthCoefficientsFileName;

//! Map of the depth calibration coefficients
unordered_map<int, vector<double>> m_Coeffs;
//! Map: detector ID (int) -> mean stretch over all pixels / strips
unordered_map<int, double> m_MeanStretch;
//! Map: detector ID (int) -> mean offset over all pixels / strips
unordered_map<int, double> m_MeanOffset;
//! Map: detector ID (int) -> Side (LV=0, HV=1) -> Strip ID -> Depth calibration coefficients (per strip)
unordered_map<int, vector<unordered_map<int, vector<double>>>> m_StripCoeffs;

//! Reference energy of the depth calibration coefficients
double m_Coeffs_Energy;
Expand Down
185 changes: 94 additions & 91 deletions src/MSubModuleDepthReadout.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@

// ROOT libs:
#include "TRandom.h"
#include "TMath.h"

// MEGAlib libs:
#include "MModuleDepthCalibration.h"
Expand Down Expand Up @@ -75,7 +74,9 @@ bool MSubModuleDepthReadout::Initialize()
m_CTDMap.clear();
m_ElectronDriftTimes.clear();
m_HoleDriftTimes.clear();
m_Coeffs.clear();
m_StripCoeffs.clear();

m_Coeffs_Energy = 59.5;

// Load depth-related files using the parsers in MModuleDepthCalibration
MModuleDepthCalibration DepthCalibration;
Expand Down Expand Up @@ -111,21 +112,10 @@ bool MSubModuleDepthReadout::Initialize()
return false;
}

// Load depth calibration coefficients
DepthCalibration.SetCoeffsFileName(m_DepthCoefficientsFileName);
if (DepthCalibration.LoadCoeffsFile(m_DepthCoefficientsFileName) == true) {
// Copy depth calibration coefficients
m_Coeffs = DepthCalibration.GetCoeffs();
m_Coeffs_Energy = DepthCalibration.GetCoeffsEnergy();

// The reference energy for the timing noise should be in the file header of the depth coefficients file
// Throw a warning if it was not retrieved and m_Coeffs_Energy is still at its default value of 0
if (m_ApplyTimingResolutionCalibration == true && m_Coeffs_Energy == 0) {
if (g_Verbosity >= c_Warning) {
cout << "Timing values will not be smeared as no reference energy found in depth spline file "<<m_DepthCoefficientsFileName<<endl;
}
if (LoadStripCoeffsFile(m_DepthCoefficientsFileName) == false) {
if (g_Verbosity >= c_Error) {
cout << "ERROR in MSubModuleDepthReadout::Initialize: Could not parse depth calibration coefficients file" << endl;
}
} else {
return false;
}

Expand Down Expand Up @@ -164,6 +154,24 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
list<MDEEStripHit>& LVHits = Event->GetDEEStripHitLVListReference();
for (MDEEStripHit& SH: LVHits) {
int DetID = SH.m_ROE.GetDetectorID();
double MeanStretch, MeanOffset;
if (m_MeanStretch.count(DetID) == 1){
MeanStretch = m_MeanStretch[DetID];
} else {
if (g_Verbosity >= c_Error) {
cout << "Detector " << DetID << " does not have a mean stretch defined" << endl;
}
return false;
}
if (m_MeanOffset.count(DetID) == 1){
MeanOffset = m_MeanOffset[DetID];
} else {
if (g_Verbosity >= c_Error) {
cout << "Detector " << DetID << " does not have a mean offset defined" << endl;
}
return false;
}

double Z = SH.m_SimulatedPositionInDetector.Z();
if (m_DepthGrid.count(DetID) == 1) {
if (SH.m_IsGuardRing == false) {
Expand All @@ -176,23 +184,21 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)

// Apply stretch and offset based on Eq. (3) in https://doi.org/10.1016/j.nima.2026.171332
unsigned int StripID = SH.m_ROE.GetStripID();
int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID;
if (m_Coeffs.count(PixelCode) == 1){
vector<double> Coeffs = m_Coeffs[PixelCode];
double Stretch = Coeffs[0];
double Offset = Coeffs[1];
double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
double CTD_Sigma = CTD_FWHM / 2.355;
if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][0].count(StripID) == 1){

vector<double> Coeffs = m_StripCoeffs[DetID][0][StripID];
double Stretch = MeanStretch * Coeffs[0];
double Offset = MeanOffset + Coeffs[1];
double TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
double HoleDriftTime = (HoleSpline.Eval(Z) + Offset) * Stretch;

// Convert drift time to timing by subtracting 4200ns (for now)
// TODO: Improve determining the timing from drift times
SH.m_Timing = 4200.0 - HoleDriftTime;

// Smear the timing value based on the given CTD resolution
// --> divide by √2 to obtain TAC resolution from CTD resolution
// Smear the timing value based on the given strip TAC resolution
if (m_ApplyTimingResolutionCalibration == true) {
SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0));
SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma);
}

// Apply the inverse TAC cal to obtain TAC in ADC units
Expand All @@ -219,7 +225,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
}
} else {
if (g_Verbosity >= c_Info) {
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<<PixelCode<<"."<<endl;
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficients found for LV strip "<<StripID<<"."<<endl;
}
SH.m_TAC = 0;
SH.m_HasTriggered = false;
Expand All @@ -234,6 +240,16 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
list<MDEEStripHit>& HVHits = Event->GetDEEStripHitHVListReference();
for (MDEEStripHit& SH: HVHits) {
int DetID = SH.m_ROE.GetDetectorID();
double MeanStretch;
if (m_MeanStretch.count(DetID) == 1){
MeanStretch = m_MeanStretch[DetID];
} else {
if (g_Verbosity >= c_Error) {
cout << "Detector " << DetID << " does not have a mean stretch defined" << endl;
}
return false;
}

double Z = SH.m_SimulatedPositionInDetector.Z();
if (m_DepthGrid.count(DetID) == 1) {
if (SH.m_IsGuardRing == false) {
Expand All @@ -247,23 +263,20 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
// Apply stretch based on Eq. (3) in https://doi.org/10.1016/j.nima.2026.171332
// Apply no offset to the electron drift time --> add it fully to the hole signal
unsigned int StripID = SH.m_ROE.GetStripID();
int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID;
if (m_Coeffs.count(PixelCode) == 1){
if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][1].count(StripID) == 1){

vector<double> Coeffs = m_Coeffs[PixelCode];
double Stretch = Coeffs[0];
double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
double CTD_Sigma = CTD_FWHM / 2.355;
vector<double> Coeffs = m_StripCoeffs[DetID][1][StripID];
double Stretch = MeanStretch * Coeffs[0];
double TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
double ElectronDriftTime = ElectronSpline.Eval(Z) * Stretch;

// Convert drift time to timing by subtracting 4200ns (for now)
// TODO: Improve determining the timing from drift times
SH.m_Timing = 4200.0 - ElectronDriftTime;

// Smear the timing value based on the given CTD resolution
// --> divide by √2 to obtain TAC resolution from CTD resolution
// Smear the timing value based on the given strip TAC resolution
if (m_ApplyTimingResolutionCalibration == true) {
SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0));
SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma);
}

// Apply the inverse TAC cal to obtain TAC in ADC units
Expand All @@ -283,7 +296,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
}
} else {
if (g_Verbosity >= c_Info) {
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<<PixelCode<<endl;
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficients found for HV strip "<<StripID<<endl;
}
SH.m_TAC = 0;
SH.m_HasTriggered = false;
Expand All @@ -303,72 +316,59 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
////////////////////////////////////////////////////////////////////////////////


bool MSubModuleDepthReadout::LoadSplinesFile(MString FileName)
bool MSubModuleDepthReadout::LoadStripCoeffsFile(MString FileName)
{
// Input spline files should have the following format:
// ### DetID, HV, Temperature, Photopeak Energy
// depth, ctd, electron_drift_time, hole_drift_time

MFile SplineFile;
if (SplineFile.Open(FileName) == false) {
if (g_Verbosity >= c_Error) {
cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: failed to open depth splines file." << endl;
}
MFile CoeffsFile;
if (CoeffsFile.Open(FileName) == false) {
cout << "ERROR in MModuleDepthCalibration::LoadStripCoeffsFile: failed to open strip coefficients file." << endl;
return false;
}

// Populate these vectors for each detector
int DetID = -1;
vector<double> DepthVector;
vector<double> CTDVector;
vector<double> ElectronDriftTimeVector;
vector<double> HoleDriftTimeVector;

MString Line;
while (SplineFile.ReadLine(Line) == true) {
if (Line.Length() != 0) {
if (Line.BeginsWith("#") == true) {

if (DetID != -1 && DepthVector.size() > 0) {
m_DepthGrid[DetID] = DepthVector;
m_CTDMap[DetID] = CTDVector;
m_ElectronDriftTimes[DetID] = ElectronDriftTimeVector;
m_HoleDriftTimes[DetID] = HoleDriftTimeVector;
}
while (CoeffsFile.ReadLine(Line) == true) {
if (Line.BeginsWith('#') == true) {
std::vector<MString> Tokens = Line.Tokenize(",");
if (Tokens.size() < 5 || Tokens[0] != "#detector_info" || Tokens[1] == "detector_id") continue;

int DetID = Tokens[1].ToInt();
// MString DetectorName = Tokens[2];
// MString DepthSplineFileName = Tokens[3];
m_MeanStretch[DetID] = Tokens[4].ToDouble();
m_MeanOffset[DetID] = Tokens[5].ToDouble();

if (m_StripCoeffs.count(DetID) == 0) {
vector<unordered_map<int, vector<double>>> TempVector;
unordered_map<int, vector<double>> TempMapLV;
unordered_map<int, vector<double>> TempMapHV;
m_StripCoeffs[DetID] = TempVector;
m_StripCoeffs[DetID].push_back(TempMapLV);
m_StripCoeffs[DetID].push_back(TempMapHV);
}

// Get the Detector ID from the commented lines
vector<MString> Tokens = Line.Tokenize(" ");
DetID = Tokens[1].ToInt();
DepthVector.clear();
CTDVector.clear();
ElectronDriftTimeVector.clear();
HoleDriftTimeVector.clear();

} else {
vector<MString> Tokens = Line.Tokenize(",");
if (Tokens.size() >= 4) {
m_DepthGrid[DetID].push_back(Tokens[0].ToDouble());
m_CTDMap[DetID].push_back(Tokens[1].ToDouble());
m_ElectronDriftTimes[DetID].push_back(Tokens[2].ToDouble());
m_HoleDriftTimes[DetID].push_back(Tokens[3].ToDouble());
} else {
if (g_Verbosity >= c_Error) {
cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: Empty line in depth splines file." << endl;
}
return false;
}
if (g_Verbosity >= c_Info) {
cout << m_Name << ": Detector " << DetID << " has a mean stretch of " << m_MeanStretch[DetID] << " and a mean offset of " << m_MeanOffset[DetID] << endl;
}
} else {
std::vector<MString> Tokens = Line.Tokenize(",");
if (Tokens.size() == 7) {
// unsigned int ReadOutID = Tokens[0].ToUnsignedInt();
unsigned int DetID = Tokens[1].ToUnsignedInt();
unsigned int Side = Tokens[2].ToUnsignedInt();
int StripID = Tokens[3].ToInt();
double Stretch = Tokens[4].ToDouble();
double Offset = Tokens[5].ToDouble();
double TimingNoiseSigma = Tokens[6].ToDouble();

vector<double> coeffs;
coeffs.push_back(Stretch); coeffs.push_back(Offset); coeffs.push_back(TimingNoiseSigma);
m_StripCoeffs[DetID][Side][StripID] = coeffs;
}
}
}

if (DetID == -1 || m_DepthGrid.size() == 0) {
if (g_Verbosity >= c_Error) {
cout << "ERROR in MSubModuleDepthReadout::LoadSplinesFile: No depth splines recovered from the file." << endl;
}
return false;
}
CoeffsFile.Close();

SplineFile.Close();
return true;
}

Expand All @@ -384,6 +384,9 @@ void MSubModuleDepthReadout::Finalize()
m_ElectronDriftTimes.clear();
m_HoleDriftTimes.clear();

m_StripCoeffs.clear();
m_TACCal.clear();

MSubModule::Finalize();
}

Expand Down