Skip to content

Trapping correction module and app - #184

Open
sophieehaight wants to merge 10 commits into
cositools:develop/emfrom
sophieehaight:trapping_correction_branch
Open

sophieehaight wants to merge 10 commits into
cositools:develop/emfrom
sophieehaight:trapping_correction_branch

Conversation

@sophieehaight

Copy link
Copy Markdown

Contains a new module for applying depth-based charge trapping correction to individual hit energies along with modules with GUI options and GUI expos. The branch also includes an app for characterizing trapping with Cs-137 data. The trapping correction module requires a csv parameter file as input:
detector_0_trapping_parameters.csv
This file can also be found in resource/dee

@sophieehaight
sophieehaight marked this pull request as ready for review July 24, 2026 17:56
Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread include/MModuleTrappingCorrection.h Outdated
Comment thread resource/dee/dummy_trapping_parameters_singledetector.csv
@fhagemann

fhagemann commented Jul 28, 2026 •

Copy link
Copy Markdown

This is what I'm getting if I'm plotting the content of the file (parsing the base CCE curve, and applying the calibration parameters from the first line):

image

@fhagemann fhagemann left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

It looks like the units of depth_val (mm) and m_Depths (cm) do not match..

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
@fhagemann

fhagemann commented Jul 30, 2026 •

Copy link
Copy Markdown

I looked at applying the inverse of the charge trapping correction to the DEE, and with the latest file, I would say that it gives reasonable results.

image image image

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
@fhagemann

Copy link
Copy Markdown

Can we also remove HasCalibratedTiming() == true from the trapping app and replace by something like HasFastTiming() == true && HasNearestNeighbor() == false (in an effort to completely get rid of HasCalibratedTiming)?

@cositools cositools deleted a comment from codecov-commenter Aug 17, 2026

@fhagemann fhagemann left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Here is my feedback from just looking at the code.
I will also try to compile and run both the module and the app to provide more feedback on that as well! :)

Comment thread apps/TrappingCorrectionAm241.cxx
Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Comment thread apps/TrappingCorrectionCs137.cxx
Comment thread apps/TrappingCorrectionCs137.cxx
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread include/MModuleTrappingCorrection.h Outdated
Comment thread include/MModuleTrappingCorrection.h Outdated
Comment thread include/MModuleTrappingCorrection.h
@fhagemann

Copy link
Copy Markdown

The example file in resource/dee and the results after applying the inverse charge trapping in the DEE looks reasonable:

image image image image image

@codecov-commenter

Copy link
Copy Markdown

Welcome to Codecov 🎉

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

ℹ️ You can also turn on project coverage checks and project coverage reporting on Pull Request comment

Thanks for integrating Codecov - We've got you covered ☂️

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MGUIOptionsLoaderMeasurementsHDF.cxx
@fhagemann fhagemann linked an issue Aug 18, 2026 that may be closed by this pull request

@fhagemann fhagemann left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Final code-style review, will give it a run to see that it compiles.
But then this should be good to go (tackling reading trapping parameters for multiple detectors in a follow-up PR)

Comment thread src/MModuleTrappingCorrection.cxx Outdated
m_Name = "Trapping Correction"; // - correcting energies for charge trapping (by Sophie);

// Set the XML tag --- has to be unique --- no spaces allowed
m_XmlTag = "TrappingCorrection";

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Let's be consistent with XML tag naming here (see #193)

Suggested change
m_XmlTag = "TrappingCorrection";
m_XmlTag = "XmlTagTrappingCorrection";

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.

I'm actually not seeing this convention in the depth or energy calibration modules? It looks like the XML tag is just the name of the module

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

This is tackled in #197 once it is merged

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.

Once I change the name of the xml tag, nuclearizer can no longer recognize the module. Is there somewhere else in assembly I need to update the xml tag? I was looking but I couldn't find anything

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

This will affect reading in a nuclearizer config file (*.cfg). In those, you might need to update all occurences of TrappingCorrection with XmlTrappingCorrection.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

So either update the .nuclearizer.cfg in your home directory, or the respective .cfg that you load when running nuclearizer.

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment on lines +277 to +280
cout << "\n" << m_XmlTag << " --- " << titleLabel << " ---" << endl;
cout << " Centroid (Mu) : " << mu << " keV" << endl;
cout << " Fitted Gaussian FWHM : " << fwhm << " keV" << endl;
cout << " Direct Histogram FWHM: " << directFWHM << " keV" << endl;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Put this in c_Info?

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.

I use this output to evaluate how well the trapping correction is working. If I change the verbosity to info, the comand line becomes overloaded (I believe w info from the depth cal) and nuclearizer crashes. I agree this info doesn't need to be repeated so I can change the verbosity for the following comments but I do need it read out somewhere.

I think we probably need to clean up the outputs for the info verbosity right now because right now it's outputting so much info that I can't receive the info I'm trying to get from the trapping correction.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

@zoglauer maybe we can introduce the chatty verbosity from megalib also to nuclearizer, and move some of the info messages to chatty?

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.

yeah that could be a good idea!

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated

@fhagemann fhagemann left a comment •

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

This would be good from my side. Just gathering here the open TODOs (for follow-up PRs):

  • Allow for a charge trapping parameters file for multiple detectors, at not just one
  • Think about flushing output to c_Info without overloading the output
  • Maybe find an alternative to Interpolate to avoid code duplication (maybe from an existing package or another module --> the depth calibration also used some sort of interpolation algorithm to interpolate the CTD-to-depth spline).

@fhagemann fhagemann left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Some code-style and some more major comments.

One general comment:
you are using a lot of snake-case when assigning variable names, we might want to use camel-case consistently here (e.g. SimCCESortedHoles instead of sim_cce_sorted_h)

Comment thread src/MModuleTrappingCorrection.cxx
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment on lines +357 to +359
// Evaluate full function FWHM over the fit window (645 keV to 675 keV)
double xMinFit = 645.0;
double xMaxFit = 675.0;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

This seems to be tailored to Cs137. Is this what we want? What if we pass data from other isotopes through MModuleTrappingCorrection?

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.

I set it up in commits 8f30f9b, b95dc1e, and f191838 (took a couple tries to get my axes to update correctly here). I let the user define the energy of the photopeak they are tryting to fit on the nuclearizer expo. The trapping correction module then takes that input and uses it to define the Mu and min and max energies for the photopeak fits measuring resolutions pre and post trapping correction. I kept the default energy for Cs-137 as a benchmark.

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment on lines +588 to +613
@@ -605,7 +600,7 @@ TF1* MModuleTrappingCorrection::GeneratePhotopeakFunction()

// Provide initial sensible guesses for a Cs137 photopeak
PhotopeakFunction->SetParameter("Amplitude", 1000);
PhotopeakFunction->SetParameter("x0 (Mu)", muInitial);
PhotopeakFunction->SetParameter("x0 (Mu)", 661.7);
PhotopeakFunction->SetParameter("Sigma Gauss", 2.0);
PhotopeakFunction->SetParameter("BoverA", 0.05);
PhotopeakFunction->SetParameter("Gamma", 0.5);
@@ -615,7 +610,7 @@ TF1* MModuleTrappingCorrection::GeneratePhotopeakFunction()

// Set boundary limits to stabilize convergence
PhotopeakFunction->SetParLimits(0, 1, 1e8);
PhotopeakFunction->SetParLimits(1, fitMin, fitMax); // Keeps peak centered around Mu
PhotopeakFunction->SetParLimits(1, 645, 675); // Keeps peak centered around 662 keV

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Your last pushes reverted back to the fixed window around 662keV.
Is this intended?

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.

my mistake! I was testing the fit with and without the dynamic fitting window. This is fixed in 6321870

@fhagemann

Copy link
Copy Markdown

I've lost track. Is this ready for a final review?

@fhagemann

Copy link
Copy Markdown

If the last changes conclude this PR, we might want to clean up the commit history and reduce the ~100 commits to a few major commits (adding the module, adding the GUI options, adding example files, ...), as this PR changes only a handful of files and adds a handful to the resource folder.

Happy to help here!

@sophieehaight
sophieehaight force-pushed the trapping_correction_branch branch from 33d2a05 to 4aba2dc Compare September 25, 2026 21:43

@fhagemann fhagemann left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Comments on the Cs137 trapping app only:

Comment on lines +68 to +69
double g_MinCTD = -250;
double g_MaxCTD = 250;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Is there ever a scenario where this range might no be enough?
For payload calibrations, we will have slightly higher temperatures than in unit-level calibrations, resulting in lower charge carrier mobilities, that might push the CTD out of the [-250,250] range (maybe?)

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.

Perhaps on the AC side we would see CTD >250. On the DC side most detectors can only get the required number of counts in a Cs-137 photopeak at CTD = -200 so it is not an issue. Right now the keeping CTD < +250 allows me to place an upper limit on the CTD value when map depth values to CTD values for the CCE parameterization.

Getting data at CTD >250 is something I want to take a look at in the future but right now we get plenty of counts on that side of the detector. I'll add a TO DO comment for this but for the purposes of this pull request I think we can leave these limits as is.

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Comment on lines +146 to +147
MString m_TACCalFile;
MString m_TACCalibrationFile;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

What is the difference between TACCalFile and TACCalibrationFile?
Edit: looks like m_TACCalibrationFile is used for what used to be the TAC cut file, so this can go:

Suggested change
MString m_TACCalFile;
MString m_TACCalibrationFile;
MString m_TACCalFile;

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Usage<<" --emax: maximum Event energy (default 700 kev)"<<endl;
Usage<<" -e: energy calibration file (.ecal)"<<endl;
Usage<<" --tcal: TAC calibration file"<<endl;
Usage<<" --tcut: TAC cut file"<<endl;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I don't think that we need TAC cut files anymore with the latest changes in #194

Comment thread apps/TrappingCorrectionCs137.cxx Outdated

// First check if each option has sufficient arguments:
// Single argument
if ((Option == "-i") || (Option == "-o") || (Option == "--emin") || (Option == "--emax") || (Option == "--tcal") || (Option == "--tcut") || (Option == "-m") || (Option == "--ctdmin") || (Option == "--ctdmax")) {

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Suggested change
if ((Option == "-i") || (Option == "-o") || (Option == "--emin") || (Option == "--emax") || (Option == "--tcal") || (Option == "--tcut") || (Option == "-m") || (Option == "--ctdmin") || (Option == "--ctdmax")) {
if ((Option == "-i") || (Option == "-o") || (Option == "--emin") || (Option == "--emax") || (Option == "--tcal") || (Option == "-m") || (Option == "--ctdmin") || (Option == "--ctdmax")) {

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Comment on lines +273 to +277
if (Option == "--tcut") {
m_TACCalibrationFile = argv[++i];
cout<<"Accepting file name: "<<m_TACCalibrationFile<<endl;
}

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Suggested change
if (Option == "--tcut") {
m_TACCalibrationFile = argv[++i];
cout<<"Accepting file name: "<<m_TACCalibrationFile<<endl;
}

Comment on lines +404 to +408
cout << "Creating energy calibrator" << endl;
EnergyCalibrator = new MModuleEnergyCalibration();
EnergyCalibrator->SetFileName(m_EcalFile);
S->SetModule(EnergyCalibrator, MNumber);
++MNumber;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Do you require any sort of energy (slow) thresholds to be set here?

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.

Yes, I set a slow threshold at 30 keV. I had originally had this set locally in the ecal module but I will add it to the app.

Comment thread apps/TrappingCorrectionCs137.cxx
Comment thread apps/TrappingCorrectionCs137.cxx
Comment thread apps/TrappingCorrectionCs137.cxx Outdated
Comment on lines +755 to +776
MStripHit* TrappingCorrectionCs137::GetDominantStrip(vector<MStripHit*>& Strips, double& EnergyFraction)
{
double MaxEnergy = -numeric_limits<double>::max();
double TotalEnergy = 0.0;
MStripHit* MaxStrip = nullptr;

// Iterate through strip hits and get the strip with highest energy
for (const auto SH : Strips) {
double Energy = SH->GetEnergy();
TotalEnergy += Energy;
if (Energy > MaxEnergy) {
MaxStrip = SH;
MaxEnergy = Energy;
}
}
if (TotalEnergy == 0) {
EnergyFraction = 0;
} else {
EnergyFraction = MaxEnergy/TotalEnergy;
}
return MaxStrip;
}

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

This seems like it's an identical copy of MModuleDepthCalibration::GetDominantStrip, which is a public function, so you can probably use that here instead:

  1. Load #include "MModuleDepthCalibration.h"
  2. In Analyze, define (with all other modules) MModuleDepthCalibration* Depth = new MModuleDepthCalibration, maybe with a comment that you're not running the module, just using its functions here.
  3. Call Depth->GetDominantStrip where you just have GetDominantStrip.
  4. Remove this function here (plus its declaration in the class).

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
signal(SIGINT, CatchSignal);

// Initialize global MEGALIB variables, especially mgui, etc.
MGlobal::Initialize("Standalone", "a standalone example program");

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

You could be more explicit here:

Suggested change
MGlobal::Initialize("Standalone", "a standalone example program");
MGlobal::Initialize("TrappingCorrectionCs137", "A standalone program to determine CTD-dependent shifts of the Cs137 peak");

or something like that.

Comment thread src/MModuleTrappingCorrection.cxx

@fhagemann fhagemann left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

And some general comment on the file format (header lines starting with # or not)

@@ -0,0 +1,154 @@
# A_HV,A_LV,B,C
1.0009264205609762,1.00084829470918,1.288689028046797,0.3010529785474152
# z_depth_cm,e_CCE_HV,h_CCE_HV,e_CCE_LV,h_CCE_LV

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

The example files here in resource/dee and in resource/calibration are a bit inconsistent when it comes to the line defining the headers z_depth (in terms of if they should start with a # or not).

I would prefer the line to always start with a #,

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I guess this file won't work because it doesn't start with ### 0 ?

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.

Fixed now!

@@ -0,0 +1,154 @@
### 0
1.0008508757266585,1.0009808031358887,1.1384578503113774,0.3538332986513531
z_depth_cm,e_CCE_HV,h_CCE_HV,e_CCE_LV,h_CCE_LV

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Here and elsewhere in resource/calibration:

Suggested change
z_depth_cm,e_CCE_HV,h_CCE_HV,e_CCE_LV,h_CCE_LV
# z_depth_cm,e_CCE_HV,h_CCE_HV,e_CCE_LV,h_CCE_LV

Comment thread src/MModuleTrappingCorrection.cxx Outdated
Comment on lines +448 to +449
// Skip text column headers
if (Line.BeginsWith("z_depth")) continue;

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

This could go if we force the header lines to start with #.
Or is there any reason you chose not to?

@sophieehaight
sophieehaight force-pushed the trapping_correction_branch branch from 5809849 to ba643dc Compare September 26, 2026 23:36
created app to return depth dep centroid energy shift for Cs-137 data
Created charge trapping correction module with GUI options
@fhagemann
fhagemann force-pushed the trapping_correction_branch branch from 890d158 to 8a136f6 Compare September 27, 2026 01:05
@fhagemann
fhagemann dismissed their stale review September 27, 2026 01:08

The missing MGUIOptionsTrappingCorrection files were added

@ckierans ckierans left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This is a ton of good work, @sophieehaight !

I've done a review of just the Nuclearizer App for now, I'll try to get through the rest tomorrow.

It looks like there are a few things that you need to update to account for recent merged PRs as identified in Felix's comments and the ones below.

I've run the app, but it took a bit of editing to make things work. Do you usually run with or without the "-n" option for nearest neighbors?

At first I was getting strange results, namely a double peaked spectra.

Image

The double peak goes away when I don't use the -n option, but the resulting photopeak is much broader and shifted higher in energy, similar to what Julian has shown when we add NN strip hits without an appropriate energy threshold in at higher energy).

I then tried with a SlowThresholdCut of 15 keV and the double peak structure goes away. How were you doing this without the SlowThresholdCut?

Image

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


//! The user interface for the universal energy calibration

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

//! The user interface for the trapping correction

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.

yes, I am using the -n option to exclude nearest neighbors and I am applying a slowthreshold cut at 30 keV. I had originally just updated my Ecal module locally to apply the slow threshold cut but I am going to code the cut into the trapping app itself when I define the instance of the energy calibration. Will push that version soon!

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Do we actually want to have the option to run this without excluding NN strip hits? Maybe we should force that in the analysis?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

It still says "The user interface for the universal energy calibration" here. Update it to "trapping correction"

Comment thread include/MGUIOptionsTrappingCorrection.h Outdated
/*
* MGUIOptionsTrappingCorrection.h
*
* Copyright (C) by Andreas Zoglauer.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Add you name :)

Comment thread apps/TrappingCorrectionCs137.cxx Outdated

if (m_Interrupt == true) return false;

// [CTD, HV Energy, LV Energy] for full detector values

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Missing the DetID map in the comment:

// [CTD][DetID][HV Energy, LV Energy] for full detector values

MString InputFile = FileNames[0];
cout << " input file stored as: " << InputFile << endl;
vector<MString> HDFNames;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Add comment:

"// Enable sum of weights squared for TH1 by default"

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
EventFilter->SetMinimumHits(0);
EventFilter->SetMaximumHits(100);
EventFilter->SetMinimumTotalEnergy(m_MinEnergy);
EventFilter->SetMaximumTotalEnergy(m_MaxEnergy * 2);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Why *2 for the MaxEnergy? Maybe add a comment in the code to explain this logic?

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.

This is limiting the event level energy (summed over hv and lv) so it need to be twice our desired max for a hit energy. Let me add that comment!

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction);
MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction);

if ((LVSH->HasCalibratedTiming() == true) && (HVSH->HasCalibratedTiming() == true) && (LVSH != nullptr) && (HVSH != nullptr)) {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

HasCalibratedTiming is depreciated after PR #171. See #171 (comment)

I had to change this line to get any data in the histograms.

Suggested change
if ((LVSH->HasCalibratedTiming() == true) && (HVSH->HasCalibratedTiming() == true) && (LVSH != nullptr) && (HVSH != nullptr)) {
if ((LVSH->HasFastTiming() == true) && (HVSH->HasFastTiming() == true) && (LVSH != nullptr) && (HVSH != nullptr)) {

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
}
}

// Place this outside/before your CTD bin and Detector loops!

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Change this comment to be descriptive of the following functionality. I.e. "// Define an output file to save the photopeak fitting results"

Comment thread apps/TrappingCorrectionCs137.cxx Outdated
}
}

// Place this at the absolute end of your Analyze() function

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

"// Close output fit parameter file"


cout << "Creating energy calibrator" << endl;
EnergyCalibrator = new MModuleEnergyCalibration();
EnergyCalibrator->SetFileName(m_EcalFile);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I needed to add the following lines to include a SlowThreshold to remove the double peak spectra I was getting

Suggested change
EnergyCalibrator->SetFileName(m_EcalFile);
EnergyCalibrator->SetFileName(m_EcalFile);
EnergyCalibrator->SetSlowThresholdCutMode(MSlowThresholdCutModes::e_Fixed);
EnergyCalibrator->SetSlowThresholdCutFixedValue(15);

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.

I had been make a cut at 30 keV originally to solve this but it has the same effect at 15 keV so I can set the value there : )

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I believe you still need the line
EnergyCalibrator->SetSlowThresholdCutMode(MSlowThresholdCutModes::e_Fixed);?

if ((InputFile.GetSubString(InputFile.Length() - 4)) == "hdf5") {
HDFNames.push_back(InputFile);
cout << "hdf names loaded correctly" << endl;
} else if ((InputFile.GetSubString(InputFile.Length() - 3)) == "txt") {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Are you ever expecting to combing data from different runs?

If this .txt file option is only meant to load multiple files from the same run, then the MModuleLoaderMeasurementsHDF already has the SetLoadContinuationFiles option that should handle this.

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.

Having the txt file as an input for loading the HDF files allows me to control how many files I'm using to characterize the data. I have more data (30+ files of Cs data) for some of the detectors so I've been using only 20 data files via the txt file option to characterize individual detectors.

cleaned up outputs and added a check for fast timing
Created a trapping correction expo
@sophieehaight
sophieehaight force-pushed the trapping_correction_branch branch 2 times, most recently from 7b071e1 to 3f32b95 Compare September 29, 2026 00:15

@ckierans ckierans left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Continuing on with comments on the Nuclearizer Module for the forward pipeline...

Some larger comments, and plenty of more below:

I'm a bit confused at where the .csv files that are uploaded in the resource folder come from. When running the TrappingCorrectionCs137 app, I get a parameters file that has the Det_ID, CTD_Bin CTD_BinMidpoint_ns, HV_Centroid_keV, HV_Centroid_error_keV, LV_Centroid_keV, LV_Centroid_error_keV. But the csv file has the parameters from the z_depth_cm, e_CCE_HV, h_CCE_HV, e_CCE_LV, h_CCE_LV from a fit to the simulated CCE curves. Where does that fitting happen? That should be within the Nuclearizer app so that it generates the file that is then read into the forward pipeline.

You use the name SimCCEFile to describe the trapping correction files that you have. I know the shape is based on simulations, but you modify those simulated curves based on measurements. So I find this naming confusing. Can you change this to CCECorrectionFile or something instead? That's a lot of variables to modify, but it might be more easy to understand what those files are.

In the Nuclearizer Module Expo you have an automatic photopeak fit to report the FWHM of the uncorrected and corrected spectra. While this is a very nice confirmation that the module is doing what it's supposed to, I don't think this validation belongs in the general pipeline so I think that piece of code can be removed.

Comment thread include/MHit.h Outdated

// Local Position:

//! Set the local/raw position of the hit

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Suggested change
//! Set the local/raw position of the hit
//! Set the position of the hit in local detector coordinates

Comment thread include/MHit.h Outdated

//! Set the local/raw position of the hit
void SetLocalPosition(const MVector& Position) { m_LocalPosition = Position; }
//! Return the local/raw position of the hit

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Suggested change
//! Return the local/raw position of the hit
//! Retrurn the position of the hit in local detector coordinates

Comment thread include/MModuleTrappingCorrection.h Outdated

//! Set the photopeak mu value
void SetPhotopeakMu(double mu) { m_PhotopeakMu = mu; }

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Suggested change

Remove the line between the Set/Get functions.

Comment thread include/MModuleTrappingCorrection.h Outdated

//! Read the XML configuration
bool ReadXmlConfiguration(MXmlNode* Node);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Suggested change

Remove the line between the related Read/Create XML node functions


//! Initialize the module
virtual bool Initialize();

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Move Finalize after Initialize to keep the similar functions close together

Suggested change
//! Finalize the module
virtual void Finalize();

Comment thread src/MModuleTrappingCorrection.cxx Outdated

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

double MModuleTrappingCorrection::GetSimBasedCorrectedEnergy(double depth_val, double uncorrected_energy, const std::vector<double>& depths, const std::vector<double>& sim_cce_sorted_e, const std::vector<double>& sim_cce_sorted_h, double paramB, double paramC)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This function should use upper camel case to follow the coding conventions, e.g. "sim_cce_sorted_e" -> "SimCCESortedE"
Can you modify the variables to match the convention?

Comment thread src/MModuleTrappingCorrection.cxx Outdated
return uncorrected_energy;
}
double cce_base_e = Interpolate(depth_val, depths, sim_cce_sorted_e);
double cce_base_h = Interpolate(depth_val, depths, sim_cce_sorted_h);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

You're calling Interpolate for each event. This might be slow, so a more efficient way to do this is to store a ROOT Interpolator object in your DetectorTrappingData struct and initialize the Interpolators in the LoadSimCCEFile function at the beginning of the analysis. Then you can then evaluate them quickly for each event and not interpolating for each and every event. But, unless we should that this module is too slow, I wouldn't worry about that. I just wanted to mention it as a potential future improvement.

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.

sg! The module usually runs in a matter of seconds so I don't think it's an issue right now. Happy to update in the future though if we decide it's necessary.

@@ -0,0 +1,154 @@
### 0
1.0008508757266585,1.0009808031358887,1.1384578503113774,0.3538332986513531

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

What are these numbers at the top of the file? Can you add a comment to define them?

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.

I was waiting to correct this until I could produce the param files with the Cs137 trapping app but it's an easy fix for me in python too so I can just upload new param files with added comments :)

Comment thread src/MModuleTrappingCorrection.cxx Outdated
// Check depth bounds using this specific detector's CCE grid limits
double min_param_depth = std::min(detData.m_Depths.front(), detData.m_Depths.back());
double max_param_depth = std::max(detData.m_Depths.front(), detData.m_Depths.back());
bool is_depth_out_of_bounds = (depth_val < min_param_depth || depth_val > max_param_depth);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

It would be preferred if you used upper camel case here:
"min_param_depth" -> "MinParamDepth", etc.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Same as for the HV loop below

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.

I also updated is_depth_out_of_bounds -> IsDepthOutofBounds to match

Comment thread src/MModuleTrappingCorrection.cxx Outdated
}

double correctedLVEnergy = rawLVEnergy;
if (is_depth_out_of_bounds == false) {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I'm wondering about this forced bounds. The DepthCalibration should push the position of each hit into the active detector volume. If you find a hit outside of that, what should you do? I think it might be better if you do not place any constraints and assume the closest position, i.e. if the depth value is < min_param_depth, then just use min_param_depth as the z for the correction. That way you'll always analyze each event that the DepthCalibration passes on to you without making another down select on the position of the hit.

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.

Yes, I think this makes sense.
I originally put this check in because my old CCE were not simulating all the way to the detector edges and it was messing up the corrected energy when depths were not associated with a CCE.
Since updating the simulated CCEs and trapping param files, I haven't seen any events outside of the depth boundaries.
Still, in case it we do get out of bounds depths, I will add a corrected energy assignment for the depth.

@sophieehaight

Copy link
Copy Markdown
Author

I'm a bit confused at where the .csv files that are uploaded in the resource folder come from. When running the TrappingCorrectionCs137 app, I get a parameters file that has the Det_ID, CTD_Bin CTD_BinMidpoint_ns, HV_Centroid_keV, HV_Centroid_error_keV, LV_Centroid_keV, LV_Centroid_error_keV. But the csv file has the parameters from the z_depth_cm, e_CCE_HV, h_CCE_HV, e_CCE_LV, h_CCE_LV from a fit to the simulated CCE curves. Where does that fitting happen? That should be within the Nuclearizer app so that it generates the file that is then read into the forward pipeline.

The Cs-137 trapping at app produces an output files with the experimental line shifts in each depth bin. Then the fitting to the simulated CCE is done in python so the param generation process. In the calibrations meeting in the past we discussed putting the entire parameter generation process in nuclearizer but decided against it because the we will only use Cs-137 to characterize electron trapping for pre-flight calibrations. To characterize in-flight trapping we will need to develop another method likely using the continuum.

I would argue that it's probably not necessary to put the CCE fitting process in nuclearizer since I may be the only one who produces these pre-flight parameter files characterizing electron trapping but happy to discuss more!

Added trapping parameter files for unit level as well as dummy payload and single detector files
Record a local position in depth cal to access in trapping correction

Updated comments on local position
changed depth cal functions to public for use in the trapping correction
Created MGUIOptionsTrappingCorrection src file

Create MGUIOptionsTrappingCorrection.h
Updated for clarity and to align with TAC cut / TAC cal separation

cleaned up the trapping app comments + functionality
Removed energy histogram fits and resolution outputs

Added comments + Updates to formatting + style

Added corrected energy assignment for out of bounds depths
@sophieehaight
sophieehaight force-pushed the trapping_correction_branch branch from 83709ac to 54705c1 Compare October 5, 2026 03:14

@ckierans ckierans left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Thanks for the changes, Sophie! A few small corrections.

Where are you at with getting the CCE fits in there? I'm wondering if it makes sense to have this complete before merging, or accept it as is with the expectation that there will be modifications coming.

@sophieehaight

Copy link
Copy Markdown
Author

Where are you at with getting the CCE fits in there? I'm wondering if it makes sense to have this complete before merging, or accept it as is with the expectation that there will be modifications coming.

I'm troubleshooting the CCE fitting now. I'm hoping to have it done either today or tomorrow but that will depend on when I find the fix.

I think it is okay to merge as is, since the parameter generation pipeline is working (just not fully contained to nuclearizer). And then we can merge the updated trapping app/ param generation pipeline once I have it ironed out.

@fhagemann

Copy link
Copy Markdown

Someone should check if all the comments in this PR ended up being resolved (and flag the ones that might not have been). Is someone looking into this, or want me to do it?

@fhagemann

Copy link
Copy Markdown

I'd also be interested to see the plots of the 16 fits to determine the charge trapping correction parameters 😇

@sophieehaight

Copy link
Copy Markdown
Author

I just went through this pr now and it looks like all the comments have been addressed except concerning adding comments the trapping param files.

I was hoping to fix this once I get the trapping app outputting accurate trapping parameter fits (I'm still having trouble with the fitting functions in ROOT).

I've uploaded a presentation describing the unit level charge trapping pipeline, showing the fitted trapping parameter plots, and summarizing the improvement to resolution for each detector to the COSI drive. It can be found here.
There may be a couple of slide describing the repairs to resolution for some of the more stubborn detectors (e.g M61617-FM2) that I need to add but otherwise it should be complete!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Implement Charge Trapping correction

4 participants