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
31 changes: 28 additions & 3 deletions include/MModuleDepthCalibration.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,14 +91,31 @@ class MModuleDepthCalibration : public MModule
//! Load the detector and strip dimensions from the geometry object
bool LoadDetectorDimensions(MDGeometryQuest* Geometry);

//! Load in the specified coefficients file
//! Load the pixel-based coefficients file
bool LoadCoeffsFile(MString FName);
//! Load the strip-based coefficients file
bool LoadStripCoeffsFile(MString FName);

//! Set the depth calibration coefficients
//! Set the pixel-based depth calibration coefficients
void SetCoeffs( unordered_map<int, vector<double>> Coeffs ) { m_Coeffs = Coeffs; }
//! Get the depth calibration coefficients
//! Get the pixel-based depth calibration coefficients
unordered_map<int, vector<double>> GetCoeffs() { return m_Coeffs; }

//! Set the strip-based depth calibration coefficients
void SetStripCoeffs(unordered_map<int, vector<unordered_map<int, vector<double>>>> Coeffs) { m_StripCoeffs = Coeffs; }
//! Get the strip-based depth calibration coefficients
unordered_map<int, vector<unordered_map<int, vector<double>>>> GetStripCoeffs() { return m_StripCoeffs; }

//! Set the mean depth calibration offsets for all detectors
void SetMeanOffset(unordered_map<int, double> MeanOffset) { m_MeanOffset = MeanOffset; }
//! Get the mean depth calibration offsets for all detectors
unordered_map<int, double> GetMeanOffset() { return m_MeanOffset; }

//! Set the mean depth calibration stretches for all detectors
void SetMeanStretch(unordered_map<int, double> MeanStretch) { m_MeanStretch = MeanStretch; }
//! Get the mean depth calibration stretches for all detectors
unordered_map<int, double> GetMeanStretch() { return m_MeanStretch; }

//! Set the energy at which the depth calibration coefficients were determined
void SetCoeffsEnergy( double Coeffs_Energy ) { m_Coeffs_Energy = Coeffs_Energy; }
//! Get the energy at which the depth calibration coefficients were determined
Expand Down Expand Up @@ -188,6 +205,14 @@ class MModuleDepthCalibration : public MModule
MModuleEnergyCalibration* m_EnergyCalibration;
MGUIExpoDepthCalibration* m_ExpoDepthCalibration;


//! 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;

// The CTD Map maps each detector (int) to a 2D array of CTD values.
unordered_map<int, vector<vector<double>>> m_CTDMap;
unordered_map<int, vector<double>> m_DepthGrid;
Expand Down
10 changes: 7 additions & 3 deletions include/MSubModuleChargeTransport.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,11 +111,15 @@ class MSubModuleChargeTransport : public MSubModule

list<MDEEStripHit> m_ChargeTransportHits;

//! Filename of the depth calibration coefficients (stretch, offset, timing noise, ...)
//! Filename of the strip-based 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;

//! Filename of CTD->Depth splines
MString m_DepthSplinesFileName;
Expand Down
5 changes: 3 additions & 2 deletions include/MSubModuleDepthReadout.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,8 +96,9 @@ 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, vector<unordered_map<int, vector<double>>>> m_StripCoeffs;

//! Reference energy of the depth calibration coefficients
double m_Coeffs_Energy;

Expand Down
61 changes: 61 additions & 0 deletions src/MModuleDepthCalibration.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -538,6 +538,67 @@ bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName)
/////////////////////////////////////////////////////////////////////////////////


bool MModuleDepthCalibration::LoadStripCoeffsFile(MString FileName)

@fhagemann fhagemann Aug 5, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Up for debate:
I added the LoadStripCoeffsFile to read the strip-based depth coefficients file to MModuleDepthCalibration, even though it is never used in the forward pipeline.
The thought behind this was that both MSubModuleChargeTransport and MSubModuleDepthReadout require reading this file for z->drift time and timing smearing, respectively.

{

MFile CoeffsFile;
if (CoeffsFile.Open(FileName) == false) {
cout << "ERROR in MModuleDepthCalibration::LoadStripCoeffsFile: failed to open strip coefficients file." << endl;
return false;
}

MString Line;
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);
}

if (g_Verbosity >= c_Info) {
cout << "ERROR in MModuleDepthCalibration::LoadStripCoeffsFile: 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;
}
}
}

CoeffsFile.Close();

return true;
}


/////////////////////////////////////////////////////////////////////////////////



std::vector<double>* MModuleDepthCalibration::GetPixelCoeffs(int PixelCode)
{
// Check to see if the stretch and offset have been loaded. If so, try to get the coefficients for the specified pixel.
Expand Down
59 changes: 39 additions & 20 deletions src/MSubModuleChargeTransport.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ bool MSubModuleChargeTransport::Initialize()
return false;
}

m_Coeffs.clear();
m_StripCoeffs.clear();
m_DepthGrid.clear();
m_ElectronDriftTimes.clear();
m_HoleDriftTimes.clear();
Expand Down Expand Up @@ -182,10 +182,11 @@ bool MSubModuleChargeTransport::Initialize()
}

// Load depth calibration coefficients
DepthCalibration.SetCoeffsFileName(m_DepthCoefficientsFileName);
if (DepthCalibration.LoadCoeffsFile(m_DepthCoefficientsFileName) == true) {
if (DepthCalibration.LoadStripCoeffsFile(m_DepthCoefficientsFileName) == true) {
// Copy depth calibration coefficients
m_Coeffs = DepthCalibration.GetCoeffs();
m_StripCoeffs = DepthCalibration.GetStripCoeffs();
m_MeanStretch = DepthCalibration.GetMeanStretch();
m_MeanOffset = DepthCalibration.GetMeanOffset();

} else {
return false;
Expand Down Expand Up @@ -327,7 +328,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool

// Calculate strip ID by rounding down intentionally to avoid truncation towards zero
// TODO: Include mask metrology information when calculating the strip ID from the position.
int ID = static_cast<int>(std::floor((P + PWidth/2.0) / PPitch));
int StripID = static_cast<int>(std::floor((P + PWidth/2.0) / PPitch));

// Calculate the strip ID for the opposite side of the detector (and explicitly check for guard ring)
int OppositeStripID = static_cast<int>(std::floor((Q + QWidth/2.0) / QPitch));
Expand All @@ -337,26 +338,41 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool

// Check for strip ID and if the position is within the allowed strip length or on the guard ring
// TODO: Confirm the correct boundary of the guard ring based on SMEX detector models
if (ID >= 0 && ID < NStrips && std::abs(Q) <= QWidth/2.0 && std::hypot(P, Q) <= Radius) {
if (StripID >= 0 && StripID < NStrips && std::abs(Q) <= QWidth/2.0 && std::hypot(P, Q) <= Radius) {

// Determine the charge drift times in nanoseconds from simulations + stretch/offset from the depth calibration
// Set the default to a large number (here: 1e10 ns) in case no depth calibration coefficients exist
double FastPeakTime = 1e10;

TSpline3* DriftTimeSpline = isLV ? m_HoleDriftSplines[DetID] : m_ElectronDriftSplines[DetID];
int PixelCode = 10000*DetID + 100*(isLV ? ID : OppositeStripID) + (isLV ? OppositeStripID : ID);

double MeanStretch = 1.0;
double MeanOffset = 0.0;
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;
}
}
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;
}
}

// Apply stretch based on Eq. (3) in https://doi.org/10.1016/j.nima.2026.171332
// Apply no offset to the electron drift time --> add it fully to the hole (LV) signal
auto it = m_Coeffs.find(PixelCode);
if (it != m_Coeffs.end()) {
const vector<double>& Coeffs = it->second;
double Stretch = Coeffs[0];
double Offset = isLV ? Coeffs[1] : 0.0;
// Apply the mean offset fully to the hole (LV) signal
if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][isLV ? 0 : 1].count(StripID) == 1) {
vector<double> Coeffs = m_StripCoeffs[DetID][isLV ? 0 : 1][StripID];
double Stretch = MeanStretch * Coeffs[0];
double Offset = (isLV ? MeanOffset + Coeffs[1] : -Coeffs[1]);
FastPeakTime = (DriftTimeSpline->Eval(Z) + Offset) * Stretch;
} else {
if (g_Verbosity >= c_Warning) {
cout << "No depth calibration coefficients for pixel in DetID " << DetID << " HV " << (isLV ? OppositeStripID : ID) << " LV " << (isLV ? ID : OppositeStripID) << endl;
cout << "No depth calibration coefficients for " << (isLV ? "LV" : "HV") << " strip " << StripID << endl;
}
}

Expand Down Expand Up @@ -385,7 +401,7 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool

// create entry for the main hit
MDEEStripHit MainSH = SH;
MainSH.m_ROE.SetStripID(ID);
MainSH.m_ROE.SetStripID(StripID);
MainSH.m_OppositeStripID = OppositeStripID;
MainSH.m_Energy = MainStripEnergy;
// TODO: Implement a more realistic parameterization to determine nearest-neighbor timing values
Expand All @@ -398,8 +414,8 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool
NNLeftSH.m_Energy = std::max(NNLeftStripEnergy, 0.0);
NNLeftSH.m_FastPeakTime = FastPeakTime - 50 * (1 - NNLeftStripEnergy / SH.m_SimulatedEnergy);
NNLeftSH.m_OppositeStripID = OppositeStripID;
if (ID > 0) {
NNLeftSH.m_ROE.SetStripID(ID - 1);
if (StripID > 0) {
NNLeftSH.m_ROE.SetStripID(StripID - 1);
NNLeftSH.m_IsGuardRing = false;
} else {
NNLeftSH.m_ROE.SetStripID(NStrips);
Expand All @@ -412,8 +428,8 @@ void MSubModuleChargeTransport::RunChargeTransportForHit(MDEEStripHit& SH, bool
NNRightSH.m_Energy = std::max(NNRightStripEnergy, 0.0);
NNRightSH.m_FastPeakTime = FastPeakTime - 50 * (1 - NNRightStripEnergy / SH.m_SimulatedEnergy);
NNRightSH.m_OppositeStripID = OppositeStripID;
if (ID < NStrips - 1) {
NNRightSH.m_ROE.SetStripID(ID + 1);
if (StripID < NStrips - 1) {
NNRightSH.m_ROE.SetStripID(StripID + 1);
NNRightSH.m_IsGuardRing = false;
} else {
NNRightSH.m_ROE.SetStripID(NStrips);
Expand All @@ -438,7 +454,10 @@ void MSubModuleChargeTransport::Finalize()
{
// Finalize the analysis - do all cleanup, i.e., undo Initialize()

m_Coeffs.clear();
m_StripCoeffs.clear();
m_MeanStretch.clear();
m_MeanOffset.clear();

m_DepthGrid.clear();
m_ElectronDriftTimes.clear();
m_HoleDriftTimes.clear();
Expand Down
42 changes: 16 additions & 26 deletions src/MSubModuleDepthReadout.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,17 +73,16 @@ MSubModuleDepthReadout::~MSubModuleDepthReadout()
bool MSubModuleDepthReadout::Initialize()
{

m_Coeffs.clear();
m_StripCoeffs.clear();

// Load depth-related files using the parsers in MModuleDepthCalibration
MModuleDepthCalibration DepthCalibration;
DepthCalibration.SetUCSDOverride(false);

// Load depth calibration coefficients
DepthCalibration.SetCoeffsFileName(m_DepthCoefficientsFileName);
if (DepthCalibration.LoadCoeffsFile(m_DepthCoefficientsFileName) == true) {
// Copy depth calibration coefficients
m_Coeffs = DepthCalibration.GetCoeffs();
// Load the strip-based depth calibration coefficients
if (DepthCalibration.LoadStripCoeffsFile(m_DepthCoefficientsFileName) == true) {
// Copy strip-based depth calibration coefficients
m_StripCoeffs = DepthCalibration.GetStripCoeffs();
m_Coeffs_Energy = DepthCalibration.GetCoeffsEnergy();

// The reference energy for the timing noise should be in the file header of the depth coefficients file
Expand Down Expand Up @@ -146,17 +145,13 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
SH.m_HasFastTiming = true;

if (m_ApplyTimingResolutionCalibration == true){
int PixelCode = 10000*DetID + 100*StripID + SH.m_OppositeStripID;
if (m_Coeffs.count(PixelCode) == 1){
vector<double> Coeffs = m_Coeffs[PixelCode];
double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
double CTD_Sigma = CTD_FWHM / 2.355;
// Smear the timing value based on the given CTD resolution
// --> divide by √2 to obtain TAC resolution from CTD resolution
SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0));
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 TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma);
} else {
if (g_Verbosity >= c_Info) {
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<<PixelCode<<"."<<endl;
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for LV strip "<<StripID<<"."<<endl;
}
}
}
Expand Down Expand Up @@ -209,18 +204,13 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
SH.m_HasFastTiming = true;

if (m_ApplyTimingResolutionCalibration == true) {
int PixelCode = 10000*DetID + 100*SH.m_OppositeStripID + StripID;
if (m_Coeffs.count(PixelCode) == 1){
vector<double> Coeffs = m_Coeffs[PixelCode];
double CTD_FWHM = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
double CTD_Sigma = CTD_FWHM / 2.355;
// Smear the timing value based on the given CTD resolution
// --> divide by √2 to obtain TAC resolution from CTD resolution
SH.m_Timing = gRandom->Gaus(SH.m_Timing, CTD_Sigma / TMath::Sqrt(2.0));

if (m_StripCoeffs.count(DetID) == 1 && m_StripCoeffs[DetID].size() == 2 && m_StripCoeffs[DetID][1].count(StripID) == 1) {
vector<double> Coeffs = m_StripCoeffs[DetID][1][StripID];
double TAC_Sigma = Coeffs[2] * m_Coeffs_Energy / SH.m_Energy;
SH.m_Timing = gRandom->Gaus(SH.m_Timing, TAC_Sigma);
} else {
if (g_Verbosity >= c_Info) {
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for pixel with code "<<PixelCode<<"."<<endl;
cout<<"MSubModuleDepthReadout::AnalyzeEvent: No depth coefficient found for HV strip "<<StripID<<"."<<endl;
}
}
}
Expand Down Expand Up @@ -260,7 +250,7 @@ bool MSubModuleDepthReadout::AnalyzeEvent(MReadOutAssembly* Event)
void MSubModuleDepthReadout::Finalize()
{
// Finalize the analysis - do all cleanup, i.e., undo Initialize()
m_Coeffs.clear();
m_StripCoeffs.clear();
m_TACCal.clear();

MSubModule::Finalize();
Expand Down
Loading