Trapping correction branch by sophieehaight · Pull Request #181 · cositools/nuclearizer

sophieehaight · 2026-07-15T23:31:04Z

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

… file

…file

merging any upstream changes from the Cs-137 trapping app into the branch with the trapping correction module

fhagemann

I had a quick glimpse, and it seems like you're duplicating some function from the MModuleDepthCalibration to use here. Consider avoiding that code duplication by calling the original functions directly.

fhagemann · 2026-07-23T20:35:45Z

+int MModuleTrappingCorrection::GetHitGrade(MHit* H){
+  // Function for choosing which Depth-to-CTD relation to use for a given event.
+  // At time of writing, intention is to choose a CTD based on sub-pixel region determined via charge sharing (Event "grade").
+  // 5 possible grades, and one Error Grade, -1. GRADE 4 is as yet uncategorized complicated geometry. GRADE 5 means multiple, presumably separated strip hits.
+
+  //organize x and y strips into vectors
+  if (H == nullptr) {
+    return -1;
+  }
+  if (H->GetNStripHits() == 0) {
+    // Error if no strip hits listed. Bad grade is returned
+    if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleTrappingCorrection: HIT WITH NO STRIP HITS" << endl;
+    return -1;
+  }
+
+  // Take a Hit and separate its activated p and n strips into separate vectors.
+  std::vector<MStripHit*> LVStrips;
+  std::vector<MStripHit*> HVStrips;
+  vector<int> LVStripIDs;
+  vector<int> HVStripIDs;
+  for (unsigned int j = 0; j < H->GetNStripHits(); ++j) {
+    MStripHit* SH = H->GetStripHit(j);
+    if (SH == nullptr ) { 
+      if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleTrappingCorrection: Trapping Correction: got NULL strip hit :( " << endl;
+      return -1;
+    }
+    if (SH->GetEnergy() == 0 ) { 
+      if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleTrappingCorrection: Trapping Correction: got strip without energy :( " << endl; 
+      return -1;
+    }
+    if (SH->IsLowVoltageStrip()) {
+      LVStrips.push_back(SH); 
+      LVStripIDs.push_back(SH->GetStripID());
+    }
+    else {
+      HVStrips.push_back(SH);
+      HVStripIDs.push_back(SH->GetStripID());
+    }
+  }
+
+  // If the same strip has multiple hits, this is a bad grade.
+  bool MultiHitX = H->GetStripHitMultipleTimesLV();
+  bool MultiHitY = H->GetStripHitMultipleTimesHV();
+  if (MultiHitX || MultiHitY) {
+    return 5;  
+  }
+
+  if (LVStrips.size()>0 && HVStrips.size()>0) {
+    int HVmin = * std::min_element(HVStripIDs.begin(), HVStripIDs.end());
+    int HVmax = * std::max_element(HVStripIDs.begin(), HVStripIDs.end());
+
+    int LVmin = * std::min_element(LVStripIDs.begin(), LVStripIDs.end());
+    int LVmax = * std::max_element(LVStripIDs.begin(), LVStripIDs.end());
+
+    // If the strip hits are not all adjacent, it's a bad grade.
+    if ( ((HVmax - HVmin) >= (HVStrips.size())) || ((LVmax - LVmin) >= (LVStrips.size())) ) {
+      return 6;
+    }
+  }
+  else{
+    return -1;
+  }
+
+  int return_value;
+  // If 1 strip on each side, GRADE=0
+  // This represents the center of the pixel
+  if ( ((LVStrips.size() == 1) && (HVStrips.size() == 1)) || ((LVStrips.size() == 3) && (HVStrips.size() == 3)) ) {
+    return_value = 0;
+  } 
+  // If 2 hits on N side and 1 on P, GRADE=1
+  // This represents the middle of the edges of the pixel
+  else if ( (LVStrips.size() == 1) && (HVStrips.size() == 2) ) {
+    return_value = 1;
+  } 
+
+  // If 2 hits on P and 1 on N, GRADE=2
+  // This represents the middle of the edges of the pixel
+  else if ( (LVStrips.size() == 2) && (HVStrips.size() == 1) ) {
+    return_value = 2;
+  } 
+
+  // If 2 strip hits on both sides, GRADE=3
+  // This represents the corners the pixel
+  else if ( (LVStrips.size() == 2) && (HVStrips.size() == 2) ) {
+    return_value = 3;
+  } 
+
+  // If 3 hits on N side and 1 on P, GRADE=0
+  // This represents the middle of the pixel, near the p (LV) side of the detector.
+  else if ( (LVStrips.size() == 1) && (HVStrips.size() == 3) ) {
+    return_value = 0;
+  } 
+
+  // If 3 hits on P and 1 on N, GRADE=0
+  // This represents the middle of the pixel, near the n (HV) side of the detector.
+  else if ( (LVStrips.size() == 3) && (HVStrips.size() == 1) ) {
+    return_value = 0;
+  } 
+
+  // If 3 hits on N side and 2 on P, GRADE=0
+  // This represents the middle of the edge of the pixel, near the p (LV) side of the detector.
+  else if ( (LVStrips.size() == 2) && (HVStrips.size() == 3) ) {
+    return_value = 2;
+  } 
+
+  // If 3 hits on P and 2 on N, GRADE=0
+  // This represents the middle of the edge of the pixel, near the n (HV) side of the detector.
+  else if ( (LVStrips.size() == 3) && (HVStrips.size() == 2) ) {
+    return_value = 1;
+  } 
+
+  else {
+    // If more complicated than the above cases, return 4 for now.
+    // TODO: Handle more complicated charge distributions.
+    return_value = 4;
+  }
+
+  return return_value;
+}


Is this identical to MModuleDepthCalibration::GetHitGrade?
If so, we could avoid code duplication here by calling the original MModuleDepthCalibration::GetHitGrade instead.

fhagemann · 2026-07-23T20:37:56Z

+bool MModuleTrappingCorrection::Initialize()
+{
+
+  // The detectors need to be in the same order as DetIDs.
+  // ie DetID=0 should be the 0th detector in m_Detectors, DetID=1 should the 1st, etc.
+  vector<MDDetector*> DetList = m_Geometry->GetDetectorList();
+
+  // Look through the Geometry and get the names of all the detectors.
+  for (unsigned int i = 0; i < DetList.size(); ++i) {
+    // For now, DetID is in order of detectors, which puts contraints on how the geometry file should be written.
+    unsigned int DetID = i;
+
+    MDDetector* det = DetList[i];
+    vector<string> DetectorNames;
+    if (det->GetTypeName() == "Strip3D") {
+      if (det->GetNSensitiveVolumes() == 1) {
+        MDVolume* vol = det->GetSensitiveVolume(0);
+        string det_name = vol->GetName().GetString();
+
+        if (find(DetectorNames.begin(), DetectorNames.end(), det_name) == DetectorNames.end()) {
+          DetectorNames.push_back(det_name);
+
+          if (g_Verbosity >= c_Info) {
+            cout << "Found detector " << det_name << " corresponding to DetID=" << DetID << "." << endl;
+          }
+          m_DetectorIDs.push_back(DetID);
+          m_Detectors[DetID] = det;
+        } else {
+          if (g_Verbosity >= c_Error) {
+            cout<<"ERROR in MModuleTrappingCorrection::Initialize: Found a duplicate detector: "<<det_name<<endl;
+          }
+        }
+      }
+    }
+  } 


This is also done in MModuleDepthCalibration, so maybe we can export this to a function of MModuleDepthCalibration and call that new function here as well to avoid code duplication.
Might also be relevant for #182

fhagemann · 2026-07-23T20:40:39Z

+      if (Grade < 0){
+        H->SetNoDepth();
+      } else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent.
+        H->SetNoDepth();
+      } else { // If the Grade is 0-4, we can handle it.


Isn't this already done in MModuleDepthCalibration (setting H->SetNoDepth())?

nuclearizer/src/MModuleDepthCalibration.cxx

Lines 199 to 216 in 1a1b3ee

if (Grade < 0){

H->SetNoDepth();

Event->SetDepthCalibrationError("Error in depth calibration");

if (Grade == -1) {

++m_ErrorSH;

} else if (Grade == -2) {

++m_ErrorNullSH;

} else if (Grade == -3) {

++m_ErrorNoE;

}

} else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent.

H->SetNoDepth();

Event->SetDepthCalibrationError("Multiple hits on single strip");

if (Grade==5) {

++m_Error5;

} else if (Grade==6) {

++m_Error6;

}

If this charge trapping module is called AFTER the depth calibration, not sure if this is needed here (?)

fhagemann · 2026-07-23T20:43:47Z

+MStripHit* MModuleTrappingCorrection::GetDominantStrip(vector<MStripHit*>& Strips, double& EnergyFraction)
+{
+  double MaxEnergy = -numeric_limits<double>::max(); // AZ: When both energies are zero (which shouldn't happen) we still pick one
+  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;
+}


This also exists in MModuleDepthCalibration, so no need to duplicate ;)

fhagemann · 2026-07-23T20:47:07Z

+  //! Returns the strip with most energy from vector Strips, also gives back the energy fraction
+  MStripHit* GetDominantStrip(std::vector<MStripHit*>& Strips, double& EnergyFraction);
+
+  // //! Retrieve the appropriate Depth values given the DetID
+  // vector<double> GetDepth(int DetID);
+
+  //! Determine the Grade (geometry of charge sharing) of the Hit
+  int GetHitGrade(MHit* H);


This can also be removed (maybe check #144 for how we avoided code duplication in the DEE by reusing function from the energy calibration module).

fhagemann · 2026-07-24T17:38:59Z

This PR was merged by accident and we removed the commits from develop/em.

sophieehaight added 30 commits April 15, 2026 14:55

troubleshooting trapping app

3aee0c3

adapted trapping app to intake Cs-137 data and output lineshifts

92c5d45

changed fitting range around Cs line

c2bc846

added depth binning on pixel and detector level

c1f937a

added depth binning but need to remove pixel by pixel fitting

e2b1c0a

added CTD binning and photopeak fits to LV and HV data

260fb35

fixed fit plotting and reading out CTD bins dn centroids to param txt…

18e835a

… file

updated photopeak fit function w linear tail

53bb708

added output txt with spline fit for HV centroids

f9a12c0

now output only returns HV energy root fits and comprehensive txt fit…

29b13f9

… file

added centroid error to spline readout txt

b9e11f3

writes min and max CTD from data to splines txt file

8bc3f4f

adapted module for trapping correction, wrote function for sim input …

3d970d3

…file

Adapted header to trapping correction

ae68931

added params for sim cce input file

acf6b1f

added function for CCE corrected energy + interpolation helper function

17cfbd6

updating hit energies with correction factor

b194068

added LV CCE column to sim read in

124d92b

updated assembly files and removed GUI functionality

46eb4b4

correction compiles error free

50177b3

added fitting and param readout for HV strips

e1d03e2

removed spline fitting and consolidated to one param file output

a383733

fixed compilation error

0f13cb3

added GUI options for input param file

8187914

added expo with HV, LV photopeak plots

34381ba

implemented finish function to return line shift and corrected res

3f8f109

resolved conflicts with HDF loader

a9327e5

cleaned up commented code

b7c1cb7

cleaned up more commented code

b8f1556

updated HDF loader header file

df25e2b

sophieehaight added 3 commits July 14, 2026 17:04

I'm just trying to get the HDF loader module up to date

e1923d1

fixed compiliation error

5b3171e

Merge branch 'Cs_trapping_branch' into trapping_correction_branch

7bc061c

merging any upstream changes from the Cs-137 trapping app into the branch with the trapping correction module

sophieehaight requested a review from ckierans July 15, 2026 23:31

fhagemann self-requested a review July 23, 2026 19:25

fhagemann reviewed Jul 23, 2026

View reviewed changes

sophieehaight merged commit 7bc061c into cositools:develop/em Jul 24, 2026
1 check passed

Navigation Menu

Trapping correction branch - #181

Trapping correction branch#181
sophieehaight merged 33 commits into
cositools:develop/emfrom
sophieehaight:trapping_correction_branch

sophieehaight commented Jul 15, 2026

Uh oh!

fhagemann left a comment

Uh oh!

fhagemann Jul 23, 2026

Uh oh!

fhagemann Jul 23, 2026

Uh oh!

fhagemann Jul 23, 2026

Uh oh!

fhagemann Jul 23, 2026

Uh oh!

fhagemann Jul 23, 2026

Uh oh!

Uh oh!

fhagemann commented Jul 24, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

	if (Grade < 0){
	H->SetNoDepth();
	Event->SetDepthCalibrationError("Error in depth calibration");
	if (Grade == -1) {
	++m_ErrorSH;
	} else if (Grade == -2) {
	++m_ErrorNullSH;
	} else if (Grade == -3) {
	++m_ErrorNoE;
	}
	} else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent.
	H->SetNoDepth();
	Event->SetDepthCalibrationError("Multiple hits on single strip");
	if (Grade==5) {
	++m_Error5;
	} else if (Grade==6) {
	++m_Error6;
	}

Conversation

sophieehaight commented Jul 15, 2026

Uh oh!

fhagemann left a comment

Choose a reason for hiding this comment

Uh oh!

fhagemann Jul 23, 2026

Choose a reason for hiding this comment

Uh oh!

fhagemann Jul 23, 2026

Choose a reason for hiding this comment

Uh oh!

fhagemann Jul 23, 2026

Choose a reason for hiding this comment

Uh oh!

fhagemann Jul 23, 2026

Choose a reason for hiding this comment

Uh oh!

fhagemann Jul 23, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

fhagemann commented Jul 24, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants