Merge branch 'develop' into bugfix/wirecellgrouped

Author: nathanielerowe

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: []

sbndcode/WireCell/cfg/pgrapher/experiment/sbnd/dnnroi.jsonnet

@@ -0,0 +1,112 @@
+// This produces a function to configure DNN-ROI for one APA given
+// anode and torch service (ts) objects.
+// 
+// The prefix is prepended to all internal node names if uniqueness
+// beyond anode ID is needed.  The output_scale allows for an ad-hoc
+// scaling of dnnroi output.  The U and W planes will go through
+// dnnroi while hte W plane will be shunted.  What comes out will be a
+// unified frame with frame tag "dnnspN" where "N" is the anode ID.
+
+local wc = import "wirecell.jsonnet";
+local pg = import "pgraph.jsonnet";
+
+
+function (anode, ts_p0, ts_p1, prefix="dnnroi", output_scale=1.0, nticks=3428, tick_per_slice=4, nchunks=1)
+    local apaid = anode.data.ident;
+    local prename = prefix + std.toString(apaid);
+    local intags = ['loose_lf%d'%apaid, 'mp2_roi%d'%apaid,
+                     'mp3_roi%d'%apaid];
+
+    local dnnroi_u = pg.pnode({
+        type: "DNNROIFinding",
+        name: prename+"u",
+        data: {
+            anode: wc.tn(anode),
+            plane: 0,
+            intags: intags,
+            decon_charge_tag: "decon_charge%d" %apaid,
+            outtag: "dnnsp%du"%apaid,
+            input_scale: 0.00025, // 1/4000
+            output_scale: output_scale,
+            forward: wc.tn(ts_p0),
+            tick_per_slice: tick_per_slice,
+            nticks: nticks,
+            nchunks: nchunks
+        }
+    }, nin=1, nout=1, uses=[ts_p0, anode]);
+    local dnnroi_v = pg.pnode({
+        type: "DNNROIFinding",
+        name: prename+"v",
+        data: {
+            anode: wc.tn(anode),
+            plane: 1,
+            intags: intags,
+            decon_charge_tag: "decon_charge%d" %apaid,
+            outtag: "dnnsp%dv"%apaid,
+            input_scale: 0.00025, // 1/4000
+            output_scale: output_scale,
+            forward: wc.tn(ts_p1),
+            tick_per_slice: tick_per_slice,
+            nticks: nticks,
+            nchunks: nchunks
+        }
+    }, nin=1, nout=1, uses=[ts_p1, anode]);
+    local dnnroi_w = pg.pnode({
+        type: "PlaneSelector",
+        name: prename+"w",
+        data: {
+            anode: wc.tn(anode),
+            plane: 2,
+            tags: ["gauss%d"%apaid],
+            tag_rules: [{
+                frame: {".*":"DNNROIFinding"},
+                trace: {["gauss%d"%apaid]:"dnnsp%dw"%apaid},
+            }],
+        }
+    }, nin=1, nout=1, uses=[anode]);
+
+    local dnnpipes = [dnnroi_u, dnnroi_v, dnnroi_w];
+    local dnnfanout = pg.pnode({
+        type: "FrameFanout",
+        name: prename,
+        data: {
+            multiplicity: 3
+        }
+    }, nin=1, nout=3);
+
+    local dnnfanin = pg.pnode({
+        type: "FrameFanin",
+        name: prename,
+        data: {
+            multiplicity: 3,
+            tag_rules: [{
+                frame: {".*": "dnnsp%d%s" % [apaid,plane]},
+                trace: {".*": "dnnsp%d%s" % [apaid,plane]},
+            } for plane in ["u", "v", "w"]]
+        },
+    }, nin=3, nout=1);
+
+    local retagger = pg.pnode({
+      type: "Retagger",
+      name: 'dnnroi%d' % apaid,
+      data: {
+        // Note: retagger keeps tag_rules an array to be like frame fanin/fanout.
+        tag_rules: [{
+          // Retagger also handles "frame" and "trace" like fanin/fanout
+          // merge separately all traces like gaussN to gauss.
+          frame: {
+            ".*": "dnnsp%d" % apaid
+          },
+          merge: {
+            ".*": "dnnsp%d" % apaid
+          },
+        }],
+      },
+    }, nin=1, nout=1);
+    
+    pg.intern(innodes=[dnnfanout],
+              outnodes=[retagger],
+              centernodes=dnnpipes+[dnnfanin],
+              edges=[pg.edge(dnnfanout, dnnpipes[ind], ind, 0) for ind in [0,1,2]] +
+              [pg.edge(dnnpipes[ind], dnnfanin, 0, ind) for ind in [0,1,2]] +
+              [pg.edge(dnnfanin, retagger, 0, 0)])