UPDATED: Adding capability to deconvolve using channel-by-channel templates

duneopdet/OpticalDetector/Deconvolution/DecoAnalysis_module.cc

@@ -3,12 +3,17 @@
 // This analyzer module for Deconvolution creates histograms
 // with information from OpDetWaveforms and OpWaveforms.
 // @authors     : Daniele Guffanti, Maritza Delgado, Sergio Manthey Corchado
-// @created     : 2022 
-//=========================================================================== 
+// @created     : 2022
+//
+// Modified: Oct 7, 2024, Viktor Pec
+//     Added raw waveforms to the output. Redid directory structure
+//     inside the root file. Set the maximum number of histograms to
+//     save to 400.
+//     ===========================================================================
 
 #ifndef DecoAnalysis_h
 #define DecoAnalysis_h
- 
+
 // LArSoft includes
 #include "larcore/Geometry/Geometry.h"
 #include "lardataobj/RawData/OpDetWaveform.h"
@@ -24,19 +29,20 @@
 #include "art/Framework/Principal/Handle.h"
 #include "canvas/Persistency/Common/Ptr.h"
 #include "canvas/Persistency/Common/PtrVector.h"
+#include "canvas/Persistency/Common/FindOne.h"
 #include "art/Framework/Services/Registry/ServiceHandle.h"
 #include "art_root_io/TFileService.h"
 #include "art_root_io/TFileDirectory.h"
 #include "messagefacility/MessageLogger/MessageLogger.h"
-  
+
 // ROOT includes
 #include "TH1.h"
 #include "THStack.h"
 #include "TF1.h"
 #include "TLorentzVector.h"
 #include "TVector3.h"
 #include "TTree.h"
-  
+
 // C++ Includes
 #include <map>
 #include <vector>
@@ -48,102 +54,136 @@
 
 namespace opdet {
   class DecoAnalysis : public art::EDAnalyzer{
-    public:
-      // Standard constructor and destructor for an ART module.
-      DecoAnalysis(const fhicl::ParameterSet&);
-      virtual ~DecoAnalysis();
-      // The analyzer routine, called once per event. 
-      void analyze (const art::Event&); 
-
-    private:
-      // The parameters we'll read from the .fcl file.
-      std::string fInputModuleDeco;          // Input tag for OpDet collection
-      std::string fInputModuleDigi;          // Input tag for OpDet collection
-      std::string fInstanceName;             // Input tag for OpDet collection
-      double fSampleFreq;                    // in MHz
+  public:
+    // Standard constructor and destructor for an ART module.
+    DecoAnalysis(const fhicl::ParameterSet&);
+    virtual ~DecoAnalysis();
+    // The analyzer routine, called once per event.
+    void analyze (const art::Event&);
+
+  private:
+    // The parameters we'll read from the .fcl file.
+    std::string fInputModuleDeco;          // Input tag for OpDet collection
+    std::string fInputModuleDigi;          // Input tag for OpDet collection
+    std::string fInstanceNameDeco;             // Input tag for OpDet collection
+    std::string fInstanceNameDigi;             // Input tag for OpDet collection
+    double fSampleFreq;                    // in MHz
     //float fTimeBegin;                      // in us
     // float fTimeEnd;                        // in us
-  };  
-} 
 
-#endif 
+    static const size_t MAX_WFMS = 400;
+  };
+}
+
+#endif
 namespace opdet {
   DEFINE_ART_MODULE(DecoAnalysis)
 }
-  
+
 namespace opdet {
   //-----------------------------------------------------------------------
   // Constructor
   DecoAnalysis::DecoAnalysis(fhicl::ParameterSet const& pset)
-  : EDAnalyzer(pset){
+    : EDAnalyzer(pset){
     // Read the fcl-file
     fInputModuleDeco  =   pset.get< std::string >("InputModuleDeco");
     fInputModuleDigi  =   pset.get< std::string >("InputModuleDigi");
-    fInstanceName     =   pset.get< std::string >("InstanceName");
-
+    fInstanceNameDeco =   pset.get< std::string >("InstanceNameDeco");
+    fInstanceNameDigi =   pset.get< std::string >("InstanceNameDigi");
+
     // Obtain parameters from DetectorClocksService
-    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();     
+    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
     fSampleFreq = clockData.OpticalClock().Frequency();
 
     art::ServiceHandle< art::TFileService > tfs;
   }
-  
+
   //-----------------------------------------------------------------------
   // Destructor
-  DecoAnalysis::~DecoAnalysis() {  
+  DecoAnalysis::~DecoAnalysis() {
   }
-    
+
   //-----------------------------------------------------------------------
-  void DecoAnalysis::analyze(const art::Event& evt){ 
+  void DecoAnalysis::analyze(const art::Event& evt){
+    auto mfi = mf::LogInfo("DecoAnalysis::analyze()");
+
+
     // Map to store how many waveforms are on one optical channel
     std::map< int, int > mapChannelWF;
-  
+
     // Get deconvolved "OpWaveforms" from the event
     art::Handle< std::vector< recob::OpWaveform >> deconv;
-    evt.getByLabel(fInputModuleDeco, fInstanceName, deconv);
+    evt.getByLabel(fInputModuleDeco, fInstanceNameDeco, deconv);
     std::vector< recob::OpWaveform > OpWaveform;
     for (auto const& wf : *deconv) {OpWaveform.push_back(wf);}
-    
+
     // Get Waveforms "OpDetWaveforms" from the event
-    art::Handle< std::vector< raw::OpDetWaveform >> digi;
-    evt.getByLabel(fInputModuleDigi, fInstanceName, digi);
-    std::vector< raw::OpDetWaveform > OpDetWaveform;
-    for (auto const& wf : *digi) {OpDetWaveform.push_back(wf);}
+    //art::Handle< std::vector< raw::OpDetWaveform >> digi;
+    art::FindOne<raw::OpDetWaveform> fopDigi(deconv, evt, fInputModuleDeco);
 
     // Access ART's TFileService, which will handle creating and writing histograms for us
     art::ServiceHandle< art::TFileService > tfs;
-
+    auto evt_dir = tfs->mkdir(TString::Format("run_%d_evt_%d", evt.id().run(), evt.id().event()).Data());
+
     double firstWaveformTime = -9999;
-    for (auto const& waveform : *digi)
-    {
-      if (firstWaveformTime < 0) firstWaveformTime = waveform.TimeStamp();
-    }
+    for (long unsigned int i = 0; i < fopDigi.size(); ++i)
+      {
+	auto const &waveform = fopDigi.at(i).ref();
+	if (firstWaveformTime == -9999 || firstWaveformTime > waveform.TimeStamp())
+	  firstWaveformTime = waveform.TimeStamp();
+      }
+
+    std::map<int,std::pair<art::TFileDirectory, art::TFileDirectory> > chan_dir_map;
+
+
+    mfi<<"\n"
+       <<"Number of deconvolved waveforms in this event: "<<OpWaveform.size()
+       <<"\n";
 
-    for (int i = 0; i < int(OpWaveform.size()); i++){
-      raw::OpDetWaveform waveform = OpDetWaveform.at(i);
-      recob::OpWaveform decowaveform = OpWaveform.at(i);
+    size_t max_wfms = MAX_WFMS;
+    if (MAX_WFMS > OpWaveform.size())
+      max_wfms = OpWaveform.size();
+
+    mfi<<"Number of deconvolved waveforms to be saved in this event: "<<max_wfms<<"\n";
+
+    for (size_t i = 0; i < max_wfms; i++){
+      raw::OpDetWaveform const& waveform = fopDigi.at(i).ref();
+      recob::OpWaveform& decowaveform = OpWaveform.at(i);
       int channel = decowaveform.Channel();
 
       // Implement different end time for waveforms of variable length
       double startTime = waveform.TimeStamp() - firstWaveformTime;
       double endTime = double(waveform.size())/fSampleFreq + startTime;
 
-      // Make a name for the histogram
-      std::stringstream histName;
-      histName << "event_"      << evt.id().event() 
-              << "_opchannel_" << channel
-              << "_decowaveform_"  << mapChannelWF[channel];
-
+
+      // create channel subdir
+      if (mapChannelWF[channel] == 0) {
+	chan_dir_map.emplace( channel,
+			      std::pair(evt_dir.mkdir(TString::Format("ch%d/raw", channel).Data()),
+					evt_dir.mkdir(TString::Format("ch%d/deconv", channel).Data())));
+      }
+
+      auto &chan_dir = chan_dir_map.at(channel);
+
       // Increase counter for number of waveforms on this optical channel
-      mapChannelWF[channel]++;
+      auto nth_waveform = mapChannelWF[channel]++;
+
 
       // Create a new histogram
-      TH1D *decowaveformHist = tfs->make< TH1D >(histName.str().c_str(),  TString::Format(";t - %f (#mus);",firstWaveformTime), decowaveform.NSignal(), startTime, endTime);
+      TH1D *decowaveformHist = chan_dir.second.make< TH1D >(TString::Format("waveform_%d", nth_waveform),  TString::Format(";t - %f (#mus);",firstWaveformTime), decowaveform.NSignal(), startTime, endTime);
       // Copy values from the waveform into the histogram
       for(size_t tick = 0; tick < decowaveform.NSignal(); tick++){
-        // Fill histogram with waveform                        
-        decowaveformHist->SetBinContent(tick, decowaveform.Signal()[tick]);
-      }                
-    }                                      
-  } 
+        // Fill histogram with waveform
+        decowaveformHist->SetBinContent(tick+1, decowaveform.Signal()[tick]);
+      }
+
+      // Create a new histogram with raw waveform
+      TH1D *rawwaveformHist = chan_dir.first.make< TH1D >(TString::Format("waveform_%d", nth_waveform),  TString::Format(";t - %f (#mus);",firstWaveformTime), waveform.Waveform().size(), startTime, endTime);
+      // Copy values from the waveform into the histogram
+      for(size_t tick = 0; tick < waveform.Waveform().size(); tick++){
+        // Fill histogram with waveform
+        rawwaveformHist->SetBinContent(tick+1, waveform.Waveform()[tick]);
+      }
+    }
+  }
 } // namespace opdet

duneopdet/OpticalDetector/Deconvolution/Deconvolution.fcl

@@ -1,3 +1,5 @@
+#include "dune_opdet_channels.fcl"
+
 BEGIN_PROLOG
 
 ###################################################################
@@ -8,48 +10,118 @@ WfmWienerfilter: {
     Name: "Wiener"
     Cutoff: 1.      # Cutoff is not used in this filter
 }
-  
+
 WfmGaussfilter: {
     Name: "Gauss"
     Cutoff: 0.13     # In MHz.The cutoff frequency is defined by the standard deviation in the frequency domain.
 }                    # The cutoff value should be changed if signal smoothing is not observed.
- 
+
 ###################################################################
 
-dune_deconvolution:
+generic_dune_deconvolution:
 {
   module_type:       "Deconvolution"
-  InputModule:       "opdigi"              
-  InstanceName:      ""    
- 
+  InputModule:       "opdigi"
+  InstanceName:      ""
+
   #The LineNoiseRMS,PreTrigger,Pedestal, Samples and DigiDataFile below have the same values as the Digitizer.
 
   LineNoiseRMS:      2.6      # Pedestal RMS in [ADC] counts, likely an underestimate
-  PreTrigger:        100      # In [ticks] 
+  PreTrigger:        100      # In [ticks]
   Pedestal:          1500     # In [ADC]
   Samples:           1000     # Timewindow (ReadoutWindow) in [ticks]
   Scale:             1        # Scaling of resulting deconvolution signal.
-  DigiDataColumn:    0        # SPE template source file column.                     
-  DigiDataFile:      "SPE_DAPHNE2_FBK_2022.dat"    
-                     # The SPE template with undershoot and without pretrigger (in ADC*us), 
-
+  SPETemplateFileDataColumn: 0  # SPE template source file column.
+  SPETemplateFiles:          [] # The SPE template with undershoot and without pretrigger (in ADC*us),
+  SPETemplateMapChannels:    []
+  SPETemplateMapTemplates:   []
+  IgnoreChannels:       []
+
+  NoiseTemplateFiles:        []
+  NoiseTemplateMapChannels:  []
+  NoiseTemplateMapTemplates: []
+
+  InputPolarity: 1 # polarity of the input raw waveform
+
   AutoScale:             true    # Scaling based on SPE amplitude from template (Use "true" for Wiener Filter and
                                  # "false" for Gauss Filter). If set to false the value of Scale is used.
   ApplyPostBLCorrection: false   # Correct baseline after the deconvolution process(Use "false" for Wiener).
   PedestalBuffer:        20      # In [ticks], should always be smaller than PreTrigger.
   ApplyPostfilter:       true    # Filter the waveforms after "Wiener" deconvolution.
-  
+
   WfmFilter: @local::WfmWienerfilter     # Write the filter: "WfmWienerfilter" or "WfmGaussfilter"
   # If gauss filter is used, the following parameters should be changed:
   #WfmGaussfilter.Cutoff: 0.13; AutoScale: false; Scale: 1.; ApplyPostBLCorrection: true; ApplyPostfilter: false;
-  WfmPostfilter: @local::WfmGaussfilter  # Only available "Gauss" postfilter. 
+  WfmPostfilter: @local::WfmGaussfilter  # Only available "Gauss" postfilter.
 }
- 
-#Postfilter Cutoff 
+
+#Postfilter Cutoff
 dune_deconvolution.WfmPostfilter.Cutoff: 1.8     # Use this value only for postfilter.
 
+
+# For backwards compatibility. This table is currently used in many workflows for different detectors, data and MC.
+dune_deconvolution: @local::generic_dune_deconvolution
+dune_deconvolution.SPETemplateFiles: ["SPE_DAPHNE2_FBK_2022.dat"]
+# Set SPE template channel map: channel = -1 for all channels
+dune_deconvolution.SPETemplateMapChannels:  [-1]
+dune_deconvolution.SPETemplateMapTemplates: [ 0]
+
+
+##### Below, there are detector and data/MC spicific configurations. #####
+
+
+### FD ###
+dune_fd_deconvolution: @local::dune_deconvolution
+
+
 #By debbuging and review the values (SNR,H,S,N,G0,G1,G,V,v) of the Wiener filter:
-deconvolution_snr: @local::dune_deconvolution
+deconvolution_snr: @local::dune_fd_deconvolution
 deconvolution_snr.OutputProduct: "SNR"
 
-END_PROLOG
+
+### ProtoDUNE HD ###
+protodunehd_deconvolution: @local::dune_deconvolution
+protodunehd_deconvolution.Samples:        1024
+protodunehd_deconvolution.Pedestal:       8180
+protodunehd_deconvolution.PreTrigger:     50
+protodunehd_deconvolution.PedestalBuffer: 30
+protodunehd_deconvolution.InputModule: "pdhddaphne"
+protodunehd_deconvolution.InstanceName: "daq"
+protodunehd_deconvolution.InputPolarity: -1
+protodunehd_deconvolution.ApplyPostfilter: true
+protodunehd_deconvolution.WfmPostfilter.Cutoff: 1.5
+protodunehd_deconvolution.ApplyPostBLCorrection: true
+
+protodunehd_deconvolution.IgnoreChannels:      @local::protodunehd_pds_channels.IgnoreChannels
+
+# SPE Template
+protodunehd_deconvolution.SPETemplateFiles:        @local::protodunehd_pds_channels.SPETemplateFiles
+protodunehd_deconvolution.SPETemplateMapChannels:  @local::protodunehd_pds_channels.SPETemplateMapChannels
+protodunehd_deconvolution.SPETemplateMapTemplates: @local::protodunehd_pds_channels.SPETemplateMapTemplates
+# Noise template
+protodunehd_deconvolution.NoiseTemplateFiles:        @local::protodunehd_pds_channels.NoiseTemplateFiles
+protodunehd_deconvolution.NoiseTemplateMapChannels:  @local::protodunehd_pds_channels.NoiseTemplateMapChannels
+protodunehd_deconvolution.NoiseTemplateMapTemplates: @local::protodunehd_pds_channels.NoiseTemplateMapTemplates
+
+
+### ProtoDUNE HD MC ###
+protodunehd_mc_deconvolution: @local::protodunehd_deconvolution
+# values below are as per settings for the digitizer,table protodunehd_opdigi from opticaldetectormodules_dune.fcl
+protodunehd_mc_deconvolution.LineNoiseRMS: 4.5
+protodunehd_mc_deconvolution.Pedestal: 8200
+protodunehd_mc_deconvolution.PreTrigger: 128
+protodunehd_mc_deconvolution.InputModule: "opdigi"
+protodunehd_mc_deconvolution.InstanceName: ""
+protodunehd_mc_deconvolution.InputPolarity: 1
+
+protodunehd_mc_deconvolution.IgnoreChannels: []
+
+protodunehd_mc_deconvolution.SPETemplateFiles:        @local::protodunehd_pds_channels_mc.SPETemplateFiles
+protodunehd_mc_deconvolution.SPETemplateMapChannels:  @local::protodunehd_pds_channels_mc.SPETemplateMapChannels
+protodunehd_mc_deconvolution.SPETemplateMapTemplates: @local::protodunehd_pds_channels_mc.SPETemplateMapTemplates
+# no noise templates used in PDHD MC
+protodunehd_deconvolution.NoiseTemplateFiles:        []
+protodunehd_deconvolution.NoiseTemplateMapChannels:  []
+protodunehd_deconvolution.NoiseTemplateMapTemplates: []
+
+END_PROLOG