Merge branch 'release/v10_06_01'

Author: nathanielerowe

.github/pull_request_template.md

@@ -1,13 +1,20 @@
 ## Description 
 Please provide a detailed description of the changes this pull request introduces. 
 
+$${\color{red}!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!}$$
+
+$${\color{red}\bf{\textrm{IMPORTANT UPDATE June 22nd 2025:}}}$$  If you are making a PR which is intended as a patch for the CURRENT production (which started in Spring 2025), you must make two PRs: one for develop and one for the production/v10_06_00 branch.
+
+$${\color{red}!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!}$$
+
 ## Checklist
 - [ ] Added at least 1 label from [available labels](https://github.com/SBNSoftware/sbndcode/issues/labels?sort=name-asc).
 - [ ] Assigned at least 1 reviewer under `Reviewers`,
 - [ ] Assigned all contributers including yourself under `Assignees`
 - [ ] Linked any relevant issues under `Developement`
 - [ ] Does this PR affect CAF data format? If so, please assign a CAF maintainer ([PetrilloAtWork](https://github.com/PetrilloAtWork) or [JosiePaton](https://github.com/JosiePaton)) as additional reviewer.
 - [ ] Does this affect the standard workflow? 
+- [ ] Is this PR a patch for the ongoing production? If so, separate PR must also be made for production/v10_06_00 branch! 
 
 ### Relevant PR links (optional)
 Does this PR require merging another PR in a different repository (such as sbnanobj/sbnobj etc.)?

CMakeLists.txt

@@ -15,7 +15,7 @@
 
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
-set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.06.00.02)
+set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.06.01)
 find_package(cetmodules REQUIRED)
 project(sbndcode LANGUAGES CXX)
 

README.md

@@ -0,0 +1,4 @@
+# sbndcode
+
+## Getting Started
+[SBN Young Guide](https://sbnsoftware.github.io/SBNYoung/SBNYoungIndex.html)

sbndcode/JobConfigurations/base/cafmakerjob_sbnd_data_base.fcl

@@ -125,8 +125,11 @@ physics:
 # change tools in sbn to have 2D TPC sim/sp input
 physics.producers.pandoraShowerRazzle.SimChannelLabel: "sptpc2d:gauss"
 physics.producers.pandoraTrackDazzle.SimChannelLabel: "sptpc2d:gauss"
-
 physics.producers.cnnid.WireLabel: "sptpc2d:gauss"
+# uncomment below lines to use DNN ROI SP for Razzle/Dazzle/CNNID
+#physics.producers.pandoraShowerRazzle.SimChannelLabel: "sptpc2d:dnnsp"
+#physics.producers.pandoraTrackDazzle.SimChannelLabel: "sptpc2d:dnnsp"
+#physics.producers.cnnid.WireLabel: "sptpc2d:dnnsp"
 
 ### Calorimetry for data
 physics.producers.cnnid.PointIdAlg.CalorimetryAlg: @local::sbnd_calorimetryalgdata

sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl

@@ -127,6 +127,8 @@ physics:
 physics.producers.pandoraShowerRazzle.SimChannelLabel: "simtpc2d:simpleSC"
 physics.producers.pandoraTrackDazzle.SimChannelLabel: "simtpc2d:simpleSC"
 physics.producers.cnnid.WireLabel: "simtpc2d:gauss"
+# uncomment below line for to use DNN ROI SP for CNNID scores
+#physics.producers.cnnid.WireLabel: "simtpc2d:dnnsp"
 
 physics.producers.vertexCharge.CaloAlg: @local::sbnd_calorimetryalgmc
 physics.producers.vertexStub.CaloAlg: @local::sbnd_calorimetryalgmc

sbndcode/JobConfigurations/standard/reco/reco1_data.fcl

@@ -35,3 +35,12 @@ physics.ana: [superadata]
 outputs.out1.SelectEvents: [ "reco1" ]
 
 physics.producers.gaushit.CalDataModuleLabel: "sptpc2d:gauss"
+
+# uncomment below 4 lines to run DNN ROI finding SP
+#physics.producers.sptpc2d.wcls_main.outputers: ["wclsFrameSaver:spsaver" , "wclsFrameSaver:dnnsaver"]
+#physics.producers.sptpc2d.wcls_main.structs.use_dnnroi: true
+#physics.producers.sptpc2d.wcls_main.structs.nchunks: 2 # should match training config
+#physics.producers.sptpc2d.wcls_main.structs.tick_per_slice: 4 # should match training config
+
+# uncomment below line to run DNN ROI finding SP for reco
+#physics.producers.gaushit.CalDataModuleLabel: "sptpc2d:dnnsp"

sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl

@@ -92,3 +92,9 @@ outputs:
                     ]
   }
 }
+
+# uncomment below 4 lines to run DNN ROI finding SP
+# physics.producers.simtpc2d.wcls_main.outputers: ["wclsDepoFluxWriter:postdrift", "wclsFrameSaver:spsaver", "wclsFrameSaver:dnnsaver"] # "wclsFrameSaver:simdigits" <- by default, do not save RawDigits. Uncomment this line to save RawDigits and set
+# physics.producers.simtpc2d.wcls_main.structs.use_dnnroi: true
+# physics.producers.simtpc2d.wcls_main.structs.nchunks: 2 # should match training config
+# physics.producers.simtpc2d.wcls_main.structs.tick_per_slice: 4 # should match training config

sbndcode/JobConfigurations/standard/standard_reco1_sbnd.fcl

@@ -35,3 +35,6 @@ physics.end_paths:      [ @sequence::physics.end_paths, ana ]
 #outputs table overrides
 outputs.out1.dataTier:       "reconstructed"
 outputs.out1.outputCommands:  [@sequence::outputs.out1.outputCommands, @sequence::sbnd_reco1_drops]
+
+# uncomment below line to run DNN ROI finding SP for reco -- must have run DNN ROI finding in detsim stage
+#physics.producers.gaushit.CalDataModuleLabel: "simtpc2d:dnnsp"

sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlgWienerData_tool.cc

@@ -27,6 +27,8 @@
 
 #include "sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlg.hh"
 
+#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h"
+#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h"
 
 namespace opdet {
   class OpDeconvolutionAlgWiener;
@@ -37,7 +39,7 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
 public:
   explicit OpDeconvolutionAlgWiener(fhicl::ParameterSet const& p);
 
-  ~OpDeconvolutionAlgWiener() {}
+  ~OpDeconvolutionAlgWiener() {delete fFilterTF1;}
 
   // Required functions.
   std::vector<raw::OpDetWaveform> RunDeconvolution(std::vector<raw::OpDetWaveform> const& wfHandle) override;
@@ -65,7 +67,6 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
   double fPMTChargeToADC;
   double fDecoWaveformPrecision;
   short unsigned int fBaselineSample;
-  std::string fOpDetDataFile;
   std::string fFilter;
   std::string fElectronics;
   bool fScaleHypoSignal;
@@ -80,6 +81,7 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
   unsigned int NDecoWf;
 
   TF1 *fFilterTF1;
+
   std::vector<double> fSignalHypothesis;
   std::vector<double> fNoiseHypothesis;
 
@@ -88,6 +90,8 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
   short unsigned int fBaseSampleBins;
   double fBaseVarCut;
 
+  sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService;
+
   // Declare member functions
   void ApplyExpoAvSmoothing(std::vector<double>& wf);
   void ApplyUnAvSmoothing(std::vector<double>& wf);
@@ -125,7 +129,6 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co
   fPMTChargeToADC = p.get< double >("PMTChargeToADC");
   fDecoWaveformPrecision = p.get< double >("DecoWaveformPrecision");
   fBaselineSample = p.get< short unsigned int >("BaselineSample");
-  fOpDetDataFile = p.get< std::string >("OpDetDataFile");
   fSkipChannelList = p.get< std::vector<int>>("SkipChannelList");
   fFilter = p.get< std::string >("Filter");
   fElectronics = p.get< std::string >("Electronics");
@@ -138,6 +141,8 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co
   fBaseVarCut = p.get< double >("BaseVarCut");
   fFilterParams = p.get< std::vector<double> >("FilterParams");
 
+  fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str());
+
   fNormUnAvSmooth=1./(2*fUnAvNeighbours+1);
   NDecoWf=0;
   MaxBinsFFT=std::pow(2, fMaxFFTSizePow);
@@ -147,38 +152,10 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co
   if (fElectronics=="Daphne") fSamplingFreq=fDaphne_Freq/1000.;//in GHz
   auto const* lar_prop = lar::providerFrom<detinfo::LArPropertiesService>();
 
-  //Load SER
-  std::string fname;
-  cet::search_path sp("FW_SEARCH_PATH");
-  sp.find_file(fOpDetDataFile, fname);
-  TFile* file = TFile::Open(fname.c_str(), "READ");
-  std::vector<std::vector<double>>* SinglePEVec_p;
-  std::vector<int>* fSinglePEChannels_p;
-  std::vector<double>* fPeakAmplitude_p;
-  std::vector<std::vector<double>> * fFilterParamVector_p;
-
-  file->GetObject("SERChannels", fSinglePEChannels_p);
-  file->GetObject("SinglePEVec", SinglePEVec_p);
-  file->GetObject("PeakAmplitude",  fPeakAmplitude_p);
-  file->GetObject("FilterParams",  fFilterParamVector_p);
-
-  if (fElectronics=="Daphne") file->GetObject("SinglePEVec_40ftCable_Daphne", SinglePEVec_p);
-  fSinglePEWaveVector = *SinglePEVec_p;
-  fSinglePEChannels = *fSinglePEChannels_p;
-  fPeakAmplitude = *fPeakAmplitude_p;
-  fFilterParamVector = *fFilterParamVector_p;          
-  
-  mf::LogInfo("OpDeconvolutionAlg")<<"Loaded SER from "<<fOpDetDataFile<<"... size="<<fSinglePEWave.size()<<std::endl;
-  file->Close();
-
-  if(fUseParamFilter && !fUseParamFilterInidividualChannel){
-    //If use param filter but not one for each channel, build it here
-    fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str());
-    for(size_t k=0; k<fFilterParams.size(); k++)
-      fFilterTF1->SetParameter(k, fFilterParams[k]);
-    mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<<fFilter<<std::endl;
-  }
-  else if (!fUseParamFilter){
+  //Load PMT Calibration Database
+  fPMTCalibrationDatabaseService = lar::providerFrom<sbndDB::IPMTCalibrationDatabaseService const>();
+
+  if (!fUseParamFilter){
     //Create light signal hypothesis for "on-fly" Wiener filter
     fSignalHypothesis.resize(MaxBinsFFT, 0);
     if(fFilter=="Wiener")
@@ -200,38 +177,26 @@ std::vector<raw::OpDetWaveform> opdet::OpDeconvolutionAlgWiener::RunDeconvolutio
   {
     int channelNumber = wf.ChannelNumber();
     auto it = std::find(fSkipChannelList.begin(), fSkipChannelList.end(), channelNumber);
-    bool AnalyseChannel = false;
     if (it == fSkipChannelList.end()) {
-      //If it's not try to find its SER in the file 
-      for(size_t i=0; i<fSinglePEChannels.size(); i++)
+      fSinglePEWave = fPMTCalibrationDatabaseService->getSER(channelNumber);
+      double SPEAmplitude =  fPMTCalibrationDatabaseService->getSPEAmplitude(channelNumber);
+      double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);});
+      double SinglePENormalization = std::abs(SPEAmplitude/SPEPeakValue);
+      std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;});
+      fSinglePEWave.resize(MaxBinsFFT, 0);
+      // If use channel dependent param filter
+      if(fUseParamFilter)
       {
-        if(fSinglePEChannels[i]==channelNumber) 
+        double GaussFilterPower = fPMTCalibrationDatabaseService->getGaussFilterPower(channelNumber);
+        double GaussFilterWC = fPMTCalibrationDatabaseService->getGaussFilterWC(channelNumber);
+        fFilterParamChannel = {GaussFilterWC, GaussFilterPower};
+        for(size_t k=0; k<fFilterParamChannel.size(); k++)
         {
-          fSinglePEWave = fSinglePEWaveVector[i];
-          double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);});
-          double SinglePENormalization = std::abs(fPeakAmplitude[i]/SPEPeakValue);
-          std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;});
-          fSinglePEWave.resize(MaxBinsFFT, 0);
-          AnalyseChannel = true;
-          // If use channel dependent param filter
-          if(fUseParamFilterInidividualChannel)
-          {
-            fFilterParamChannel = fFilterParamVector[i];
-            fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str());
-            for(size_t k=0; k<fFilterParamChannel.size(); k++)
-            {
-              fFilterTF1->SetParameter(k, fFilterParamChannel[k]);
-            }
-
-
-            mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<<fFilter << " for channel " << channelNumber <<std::endl;
-          }
-          break;
+          fFilterTF1->SetParameter(k, fFilterParamChannel[k]);
         }
+          mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<<fFilter << " for channel " << channelNumber <<std::endl;
       }
-      if(AnalyseChannel == false) mf::LogError("OpDeconvolutionAlg") << " SER for channel " << channelNumber <<" not found in the file \n";
     }
-    if(!AnalyseChannel) continue;
     //Read waveform
     size_t wfsize=wf.Waveform().size();
     if(wfsize>MaxBinsFFT){
@@ -537,15 +502,13 @@ std::vector<TComplex> opdet::OpDeconvolutionAlgWiener::DeconvolutionKernel(size_
     }
   }
 
-
   if(fDebug){
     std::string name="h_wienerfilter_"+std::to_string(NDecoWf);
     TH1F * hs_wiener = tfs->make< TH1F >
       (name.c_str(),"Wiener Filter;Frequency Bin;Magnitude",size/2, 0, size/2);
     for(size_t k=0; k<size/2; k++)
       hs_wiener->SetBinContent(k, TComplex::Abs( kernel[k]*serfft[k] ) );
   }
-  if(fUseParamFilter && fUseParamFilterInidividualChannel) delete fFilterTF1;
   return kernel;
 }
 

sbndcode/OpDetReco/OpDeconvolution/Alg/opdeconvolution_alg_data.fcl

@@ -6,7 +6,6 @@ OpDeconvolutionAlgData:
   tool_type: "OpDeconvolutionAlgWienerData"
   Debug: false
   MaxFFTSizePow: 16
-  OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root"
   BaseSampleBins: 30
   BaseVarCut: 50
   SkipChannelList: []