StripEnergyThresholdFinder for per-strip slow and fast threshold extraction with diagnostics - #166
Conversation
|
I still need to fix all of the code style issues and work on some optimizations to speed the code up a bit. |
| struct StripKey | ||
| { | ||
| int det; | ||
| char side; | ||
| int strip; | ||
|
|
||
| bool operator<(const StripKey& o) const | ||
| { | ||
| if(det!=o.det) return det<o.det; | ||
| if(side!=o.side) return side<o.side; | ||
| return strip<o.strip; | ||
| } | ||
| }; |
There was a problem hiding this comment.
Can you not reuse the code from MStripMap.cxx here?
There was a problem hiding this comment.
Good point. I checked MStripMap, but it looks to me like it's centered on ReadOutID lookups. I’m currently using a local StripKey for (det, side, strip) grouping. It made more sense to me to think of the strip mapping from a detector -> side -> strip direction, especially when reviewing the code. I'd be happy to adjust it if we want to prioritize tighter integration.
There was a problem hiding this comment.
Isn't that functionality exactly what MReadOutElementDoubleStrip does?
There was a problem hiding this comment.
I did take a look at MReadOutElementDoubleStrip, and it looks like it's doing something similar (detector, strip, side), but it also looks like it's part of the MReadOutElement hierarchy, and it includes additional things (parsing, cloning, comparing objects) that we don't need here. We only need a quick identifier to group histogram data as a key in a container. It seems like using MReadOutElementDoubleStrip could bring in additional dependencies that complicate things. If there's a preferred way to represent strip identifiers in this context, then I'd be more than happy to try to figure it out. Let me know!
There was a problem hiding this comment.
This can indeed be replaced by MReadOutElementDoubleStrip. I opened a PR onto your branch, where one of the commits replaces StripKey by MReadOutElementDoubleStrip.
| class EnergyCalHelper | ||
| { | ||
| public: | ||
|
|
||
| struct Coeff | ||
| { | ||
| double c0=0; | ||
| double c1=0; | ||
| double c2=0; | ||
| double c3=0; | ||
| }; | ||
|
|
||
| bool Load(const string& fileName,int nDet=64,int nStrip=65) | ||
| { | ||
| m_NDet=nDet; | ||
| m_NStrip=nStrip; | ||
|
|
||
| m_Coeffs.resize(m_NDet); | ||
|
|
||
| for(int d=0;d<m_NDet;d++) | ||
| { | ||
| m_Coeffs[d].resize(2); | ||
|
|
||
| for(int s=0;s<2;s++) | ||
| m_Coeffs[d][s].resize(m_NStrip); | ||
| } | ||
|
|
||
| ifstream in(fileName); | ||
|
|
||
| if(!in) | ||
| { | ||
| cout<<"Unable to open calibration file "<<fileName<<endl; | ||
| return false; | ||
| } | ||
|
|
||
| string line; | ||
|
|
||
| while(getline(in,line)) | ||
| { | ||
| if(line.empty()) continue; | ||
|
|
||
| stringstream ss(line); | ||
|
|
||
| string tag; | ||
| ss>>tag; | ||
|
|
||
| if(tag!="CM") continue; | ||
|
|
||
| string unused; | ||
| int det=-1; | ||
| int strip=-1; | ||
| string side; | ||
| string order; | ||
|
|
||
| ss>>unused>>det>>strip>>side>>order; | ||
|
|
||
| if(det<0||det>=m_NDet) continue; | ||
| if(strip<0||strip>=m_NStrip) continue; | ||
|
|
||
| int sideInt=(side=="l")?0:1; | ||
|
|
||
| Coeff c; | ||
|
|
||
| if(order=="poly1zero") | ||
| ss>>c.c1; | ||
| else if(order=="poly1") | ||
| ss>>c.c0>>c.c1; | ||
| else if(order=="poly2") | ||
| ss>>c.c0>>c.c1>>c.c2; | ||
| else | ||
| ss>>c.c0>>c.c1>>c.c2>>c.c3; | ||
|
|
||
| m_Coeffs[det][sideInt][strip]=c; | ||
| } | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| double ADCToEnergy(int det,char side,int strip,double adc) const | ||
| { | ||
| int sideInt=(side=='l')?0:1; | ||
| const Coeff& c=m_Coeffs[det][sideInt][strip]; | ||
|
|
||
| return c.c3*pow(adc,3)+c.c2*pow(adc,2)+c.c1*adc+c.c0; | ||
| } | ||
|
|
||
| private: | ||
|
|
||
| int m_NDet; | ||
| int m_NStrip; | ||
|
|
||
| vector<vector<vector<Coeff>>> m_Coeffs; | ||
| }; |
There was a problem hiding this comment.
This essentially does what MModuleEnergyCalibration does.
Can you create a MModuleEnergyCalibration variable in this app and use its functions here to avoid code duplication and keep the code for this app short?
There was a problem hiding this comment.
Yes. I wrote the threshold app to be standalone, not requiring the full MSupervisor pipeline for calibration, so I built a lightweight calibration helper for direct ADC-to-energy conversion to occur within the app. But, it definitely duplicates functionality found in MModuleEnergyCalibration. I'll refactor the app to use the standard module chain.
There was a problem hiding this comment.
Alright, so I worked on integrating MModuleEnergyCalibration into the event processing, and right now it's being used to calibrate the strip hits (AnalyzeEvent and SH->GetEnergy()). It worked just fine, there. I also took a look at replacing EnergyCalHelper entirely, but I ran into some issues that broke things. In this app we're mostly operating on histogrammed data (ADC counts), not individual hits within an event. At a basic level, we just need to convert arbitrary ADC values into keV (like when we need to determine thresholds or set histogram axes), and not just values associated with a specific hit. The MModule applies calibration within the event, but it doesn't seem to provide a way to convert the ADC values outside of the event context. So, I think it might be best to keep EnergyCalHelper for now to perform those conversions, while using the MModule for all event-level calibration. If there is a more "MEGAlib" way to handle the quick ADC to keV conversions outside of the event-level pipeline, I'd be happy to look into that!
There was a problem hiding this comment.
Same thing here: check out the PR JarredMRoberts#1 onto your branch, that replaces EnergyCalHelper with the MModuleEnergyCalibration functionality, outside of the event-level pipeline.
| class TACCalHelper | ||
| { | ||
| public: | ||
|
|
||
| struct Coeff | ||
| { | ||
| double slope = 0; | ||
| double offset = 0; | ||
| }; | ||
|
|
||
| bool Load(const string& fileName) | ||
| { | ||
| ifstream in(fileName); | ||
| if(!in) | ||
| { | ||
| cout<<"Failed to open TAC calibration file: "<<fileName<<endl; | ||
| return false; | ||
| } | ||
|
|
||
| string line; | ||
| getline(in,line); // skip header | ||
|
|
||
| while(getline(in,line)) | ||
| { | ||
| if(line.empty()) continue; | ||
|
|
||
| stringstream ss(line); | ||
|
|
||
| int strip_id, det, side, strip; | ||
| double slope, slope_err, offset, offset_err; | ||
|
|
||
| ss >> strip_id >> det >> side >> strip | ||
| >> slope >> slope_err >> offset >> offset_err; | ||
|
|
||
| char sideChar = (side==0) ? 'l' : 'h'; | ||
|
|
||
| StripKey key{det, sideChar, strip}; | ||
|
|
||
| m_Coeffs[key] = {slope, offset}; | ||
| } | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| double TACToEnergy(int det, char side, int strip, double tac) const | ||
| { | ||
| StripKey key{det, side, strip}; | ||
|
|
||
| auto it = m_Coeffs.find(key); | ||
| if(it == m_Coeffs.end()) return 0; | ||
|
|
||
| return it->second.slope * tac + it->second.offset; | ||
| } | ||
|
|
||
| private: | ||
| map<StripKey, Coeff> m_Coeffs; | ||
| }; |
There was a problem hiding this comment.
This already exists in MModuleTACcut --> use existing class to avoid code duplication here.
There was a problem hiding this comment.
I took another close look at MModuleTACCut. This issue I'm running into is that this requires both a TAC Calibration file and a TAC cut file to run smoothly. It's also performing event filtering as part of its operation. For the threshold finding app, I'd prefer that we only use the TAC calibration file as an input, and I'd like to avoid applying any TAC cuts to the data in order to find the fast shaper thresholds. I did attempt to see if I could figure out a workaround to trick MModuleTACCut into thinking it was receiving a TAC cut file with arbitrarily low values, but this ended up being a bit odd. I'd prefer to use the simple TACCalHelper I've written instead.
There was a problem hiding this comment.
It looks like TACCalHelper is never used. I believe for this app you don't need to calibrate TAC values, you just care if a strip had FastTiming or not. So, I believe this whole thing can go.
| private: | ||
|
|
||
| // --- Configuration --- | ||
| EnergyCalHelper m_EnergyCal; |
There was a problem hiding this comment.
| EnergyCalHelper m_EnergyCal; | |
| MModuleEnergyCalibration m_EnergyCal; |
There was a problem hiding this comment.
MModuleEnergyCalibration is now being used for the event-level calibration but since this code operates on histogram ADC data, we still need a way to convert arbitrary AC values to keV outside of the event-level data. I'd prefer to keep EnergyCalHelper for this reason.
There was a problem hiding this comment.
You can convert arbitrary ADC values to keV outside of the event-level data, by just using MModuleEnergyCalibration::GetEnergy instead of EnergyCalHelper::ADCToEnergy, which essentially does the same thing.
|
|
||
|
|
||
| //EnergyCalHelper m_EnergyCal; | ||
| TACCalHelper m_EnergyCal_TAC; |
There was a problem hiding this comment.
That's a weird name, is Energy intended in this name? I would just call it m_TACCal.. And if it's called m_..., consider making this a (private) member variable of the class -- just like m_EnergyCal
Also: where is this TACCalHelper ever used?
There was a problem hiding this comment.
Ugh, good point. It's a bad name, and it isn't actually used anywhere :|
It's removed now.
| //EnergyCalHelper m_EnergyCal; | ||
| TACCalHelper m_EnergyCal_TAC; | ||
|
|
||
| if(!m_EnergyCal.Load(m_CalibrationFile.Data())) |
There was a problem hiding this comment.
MModuleEnergyCalibration has a ReadEnergyCalibrationFile function that you could call here instead:
| if(!m_EnergyCal.Load(m_CalibrationFile.Data())) | |
| if (m_EnergyCal.ReadEnergyCalibrationFile(m_CalibrationFile) == false) |
There was a problem hiding this comment.
Based on the comments above about how I'd like to keep the EnergyCalHelper for ADC to keV conversions in histogrammed data outside of the event context, let's keep this one as it is if that's ok.
There was a problem hiding this comment.
See my PR JarredMRoberts#1 to see how this can still be achieved without the need of EnergyCalHelper.
| EnergyCalHelper m_EnergyCal; | ||
| vector<string> m_InputFiles; | ||
| MString m_CalibrationFile; | ||
| MString m_TACCalibrationFile; |
There was a problem hiding this comment.
Hm, I'm seeing that only m_TACCalibrationFile is no longer being used. I can't remember what it used to be a part of, but it can be removed now. It looks like the other ones are being used for histogram building and ADC conversion.
There was a problem hiding this comment.
See comment above: I don't think that you are applying any TAC calibration or TAC cuts, so this can be safely removed (or the code to calibrate TAC values from ADC units to ns still needs to be added).
Co-authored-by: Felix Hagemann <hagemann@berkeley.edu>
Co-authored-by: Felix Hagemann <hagemann@berkeley.edu>
…nventions and integrate energy calibration improvements
| #include <chrono> | ||
|
|
||
| /* YAML */ | ||
| #include <yaml-cpp/yaml.h> |
There was a problem hiding this comment.
I had to add -lyaml-cpp to ALLLIBS in the Makefile to compile this app:
Line 119 in 1a1b3ee
Is there an alternative to running this app without depending on YAML? 😇
There was a problem hiding this comment.
But even if it compiled, running StripEnergyThresholdFinder $NUCLEARIZER/apps/StripEnergyThresholdFinder_config.yaml resulted in the following error:
terminate called after throwing an instance of 'YAML::TypedBadConversion<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >'
what(): bad conversion
Aborted (core dumped)
There was a problem hiding this comment.
The error I posted went away once the mismatch in capitalization between strip_map in the config.yaml and parsing Strip_map in the app was resolved :)
There was a problem hiding this comment.
Hi Jarred,
I had an extensive look at this app and finally got this to run.
I replied to all comments in this PR, and opened a separate PR onto your fork/branch, addressing some code changes to
- get this app running
- remove the helper classes by replacing them with existing nuclearizer/megalib code
Here is the PR onto your branch with detailed code changes: JarredMRoberts#1
|
|
||
|
|
||
| m_CalibrationFile = config["input"]["calibration_file"].as<string>().c_str(); | ||
| m_StripMapFile = config["input"]["Strip_map"].as<string>().c_str(); |
There was a problem hiding this comment.
This is what broke my code:
| m_StripMapFile = config["input"]["Strip_map"].as<string>().c_str(); | |
| m_StripMapFile = config["input"]["strip_map"].as<string>().c_str(); |
You might need to add specific guards here to make sure that the fields in the YAML file exist, before trying to parse it.
Or avoid YAML completely and write your own custom parser here.
| #include <chrono> | ||
|
|
||
| /* YAML */ | ||
| #include <yaml-cpp/yaml.h> |
There was a problem hiding this comment.
The error I posted went away once the mismatch in capitalization between strip_map in the config.yaml and parsing Strip_map in the app was resolved :)
| struct StripKey | ||
| { | ||
| int det; | ||
| char side; | ||
| int strip; | ||
|
|
||
| bool operator<(const StripKey& o) const | ||
| { | ||
| if(det!=o.det) return det<o.det; | ||
| if(side!=o.side) return side<o.side; | ||
| return strip<o.strip; | ||
| } | ||
| }; |
There was a problem hiding this comment.
This can indeed be replaced by MReadOutElementDoubleStrip. I opened a PR onto your branch, where one of the commits replaces StripKey by MReadOutElementDoubleStrip.
| class EnergyCalHelper | ||
| { | ||
| public: | ||
|
|
||
| struct Coeff | ||
| { | ||
| double c0=0; | ||
| double c1=0; | ||
| double c2=0; | ||
| double c3=0; | ||
| }; | ||
|
|
||
| bool Load(const string& fileName,int nDet=64,int nStrip=65) | ||
| { | ||
| m_NDet=nDet; | ||
| m_NStrip=nStrip; | ||
|
|
||
| m_Coeffs.resize(m_NDet); | ||
|
|
||
| for(int d=0;d<m_NDet;d++) | ||
| { | ||
| m_Coeffs[d].resize(2); | ||
|
|
||
| for(int s=0;s<2;s++) | ||
| m_Coeffs[d][s].resize(m_NStrip); | ||
| } | ||
|
|
||
| ifstream in(fileName); | ||
|
|
||
| if(!in) | ||
| { | ||
| cout<<"Unable to open calibration file "<<fileName<<endl; | ||
| return false; | ||
| } | ||
|
|
||
| string line; | ||
|
|
||
| while(getline(in,line)) | ||
| { | ||
| if(line.empty()) continue; | ||
|
|
||
| stringstream ss(line); | ||
|
|
||
| string tag; | ||
| ss>>tag; | ||
|
|
||
| if(tag!="CM") continue; | ||
|
|
||
| string unused; | ||
| int det=-1; | ||
| int strip=-1; | ||
| string side; | ||
| string order; | ||
|
|
||
| ss>>unused>>det>>strip>>side>>order; | ||
|
|
||
| if(det<0||det>=m_NDet) continue; | ||
| if(strip<0||strip>=m_NStrip) continue; | ||
|
|
||
| int sideInt=(side=="l")?0:1; | ||
|
|
||
| Coeff c; | ||
|
|
||
| if(order=="poly1zero") | ||
| ss>>c.c1; | ||
| else if(order=="poly1") | ||
| ss>>c.c0>>c.c1; | ||
| else if(order=="poly2") | ||
| ss>>c.c0>>c.c1>>c.c2; | ||
| else | ||
| ss>>c.c0>>c.c1>>c.c2>>c.c3; | ||
|
|
||
| m_Coeffs[det][sideInt][strip]=c; | ||
| } | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| double ADCToEnergy(int det,char side,int strip,double adc) const | ||
| { | ||
| int sideInt=(side=='l')?0:1; | ||
| const Coeff& c=m_Coeffs[det][sideInt][strip]; | ||
|
|
||
| return c.c3*pow(adc,3)+c.c2*pow(adc,2)+c.c1*adc+c.c0; | ||
| } | ||
|
|
||
| private: | ||
|
|
||
| int m_NDet; | ||
| int m_NStrip; | ||
|
|
||
| vector<vector<vector<Coeff>>> m_Coeffs; | ||
| }; |
There was a problem hiding this comment.
Same thing here: check out the PR JarredMRoberts#1 onto your branch, that replaces EnergyCalHelper with the MModuleEnergyCalibration functionality, outside of the event-level pipeline.
| class TACCalHelper | ||
| { | ||
| public: | ||
|
|
||
| struct Coeff | ||
| { | ||
| double slope = 0; | ||
| double offset = 0; | ||
| }; | ||
|
|
||
| bool Load(const string& fileName) | ||
| { | ||
| ifstream in(fileName); | ||
| if(!in) | ||
| { | ||
| cout<<"Failed to open TAC calibration file: "<<fileName<<endl; | ||
| return false; | ||
| } | ||
|
|
||
| string line; | ||
| getline(in,line); // skip header | ||
|
|
||
| while(getline(in,line)) | ||
| { | ||
| if(line.empty()) continue; | ||
|
|
||
| stringstream ss(line); | ||
|
|
||
| int strip_id, det, side, strip; | ||
| double slope, slope_err, offset, offset_err; | ||
|
|
||
| ss >> strip_id >> det >> side >> strip | ||
| >> slope >> slope_err >> offset >> offset_err; | ||
|
|
||
| char sideChar = (side==0) ? 'l' : 'h'; | ||
|
|
||
| StripKey key{det, sideChar, strip}; | ||
|
|
||
| m_Coeffs[key] = {slope, offset}; | ||
| } | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| double TACToEnergy(int det, char side, int strip, double tac) const | ||
| { | ||
| StripKey key{det, side, strip}; | ||
|
|
||
| auto it = m_Coeffs.find(key); | ||
| if(it == m_Coeffs.end()) return 0; | ||
|
|
||
| return it->second.slope * tac + it->second.offset; | ||
| } | ||
|
|
||
| private: | ||
| map<StripKey, Coeff> m_Coeffs; | ||
| }; |
There was a problem hiding this comment.
It looks like TACCalHelper is never used. I believe for this app you don't need to calibrate TAC values, you just care if a strip had FastTiming or not. So, I believe this whole thing can go.
| private: | ||
|
|
||
| // --- Configuration --- | ||
| EnergyCalHelper m_EnergyCal; |
There was a problem hiding this comment.
You can convert arbitrary ADC values to keV outside of the event-level data, by just using MModuleEnergyCalibration::GetEnergy instead of EnergyCalHelper::ADCToEnergy, which essentially does the same thing.
| EnergyCalHelper m_EnergyCal; | ||
| vector<string> m_InputFiles; | ||
| MString m_CalibrationFile; | ||
| MString m_TACCalibrationFile; |
There was a problem hiding this comment.
See comment above: I don't think that you are applying any TAC calibration or TAC cuts, so this can be safely removed (or the code to calibrate TAC values from ADC units to ns still needs to be added).
| //EnergyCalHelper m_EnergyCal; | ||
| TACCalHelper m_EnergyCal_TAC; | ||
|
|
||
| if(!m_EnergyCal.Load(m_CalibrationFile.Data())) |
There was a problem hiding this comment.
See my PR JarredMRoberts#1 to see how this can still be achieved without the need of EnergyCalHelper.
| line->Draw("SAME"); | ||
|
|
||
| // Legend entry (fix from earlier) | ||
| leg->AddEntry(line, "Fast Thresholdeshold", "l"); |
There was a problem hiding this comment.
| leg->AddEntry(line, "Fast Thresholdeshold", "l"); | |
| leg->AddEntry(line, "Fast Threshold", "l"); |
| #include <chrono> | ||
|
|
||
| /* YAML */ | ||
| #include <yaml-cpp/yaml.h> |
There was a problem hiding this comment.
After reviewing #184: Please remove the YAML dependency.
I would strongly suggest that we pass the required parameters either via terminal command options (see #184 for example), or via a custom file format with a custom parser. But in my opinion, it is an overkill to add the YAML dependency to nuclearizer and requiring it in the make process just for this app (also for users who might not even need to run this app).
|
I second not to have an additional parser linked. |
fhagemann
left a comment
There was a problem hiding this comment.
One more remark: I got this app to run and got reasonable results for the fast and slow thresholds only when using a dataset taken with NN off.
We might want to always filter out NN events, for this app to also run using datasets taken with NN on.
| MModuleLoaderMeasurementsHDF* Loader = new MModuleLoaderMeasurementsHDF(); | ||
|
|
||
| Loader->SetFileName(m_InputFiles[0].c_str()); | ||
|
|
||
| cout << "Loading file: " << m_InputFiles[0] << endl; | ||
| cout << "Number of input files: " << m_InputFiles.size() << endl; | ||
|
|
||
| Loader->SetFileNameStripMap(m_StripMapFile.Data()); |
There was a problem hiding this comment.
This app only works when using a dataset where NN was off.
If NN was on, we might require filtering those NN events out for later analysis by adding
Loader->SetIncludeNearestNeighbor(false);:
| MModuleLoaderMeasurementsHDF* Loader = new MModuleLoaderMeasurementsHDF(); | |
| Loader->SetFileName(m_InputFiles[0].c_str()); | |
| cout << "Loading file: " << m_InputFiles[0] << endl; | |
| cout << "Number of input files: " << m_InputFiles.size() << endl; | |
| Loader->SetFileNameStripMap(m_StripMapFile.Data()); | |
| MModuleLoaderMeasurementsHDF* Loader = new MModuleLoaderMeasurementsHDF(); | |
| Loader->SetFileName(m_InputFiles[0].c_str()); | |
| cout << "Loading file: " << m_InputFiles[0] << endl; | |
| cout << "Number of input files: " << m_InputFiles.size() << endl; | |
| Loader->SetFileNameStripMap(m_StripMapFile.Data()); | |
| Loader->SetIncludeNearestNeighbor(false); |
Using gse_20260217T115220.hdf5:
Current state, not filtering NN events:
| Slow thresholds | Example energy spectrum (HV10) |
|---|---|
![]() |
![]() |
Filtering NN events:
| Slow thresholds | Example energy spectrum (HV10) |
|---|---|
![]() |
![]() |





Adds a standalone application, StripEnergyThresholdFinder, for computing per-strip slow and fast energy thresholds. Slow thresholds are determined from ADC spectra using a noise peak and trough method, while fast thresholds are determined from dt0/dt1 timing crossover. The tool reads calibrated HDF5 data via MModuleLoaderMeasurementsHDF, applies strip mapping and energy calibration from a YAML configuration, and produces ROOT diagnostic outputs (energy spectra with thresholds, dt0 vs dt1 per strip, and threshold distributions) along with CSV export files. The implementation is self-contained under apps/ and does not modify existing modules. Tested on COSI datasets with consistent threshold behavior and expected diagnostic results. Target branch is develop/em.