Merge pull request #681 from SBNSoftware/feature/acastill_flash_drift_distance

Feature/acastill flash drift distance

Author: nathanielerowe

sbndcode/JobConfigurations/standard/reco/reco1_data.fcl

@@ -1,6 +1,7 @@
 #include "wcsp_data_sbnd.fcl"
 #include "opdeconvolution_sbnd_data.fcl"
 #include "sbnd_ophitfinder_deco_data.fcl"
+#include "sbnd_flashfinder_deco_data.fcl"
 #include "crt_channel_map_service.fcl"
 #include "crt_calib_service.fcl"
 #include "standard_reco1_sbnd.fcl"
@@ -23,8 +24,8 @@ physics.producers:
     crttracks:      @local::crttrackproducer_data_sbnd
     opdecopmt: @local::SBNDOpDeconvolutionPMT_data
     ophitpmt:  @local::SBNDDecoOpHitFinderPMT_data
-    opflashtpc0:   @local::SBNDDecoSimpleFlashTPC0
-    opflashtpc1:   @local::SBNDDecoSimpleFlashTPC1
+    opflashtpc0:   @local::SBNDDecoSimpleFlashTPC0_data
+    opflashtpc1:   @local::SBNDDecoSimpleFlashTPC1_data
 }
 
 physics.reco1: [sptpc2d, gaushit, numberofhitsfilter, cluster3d, crtstrips, 

sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl

@@ -4,7 +4,7 @@
 #include "sam_sbnd.fcl"
 
 #include "opdeconvolution_sbnd_data.fcl"
-#include "sbnd_flashfinder_deco.fcl"
+#include "sbnd_flashfinder_deco_data.fcl"
 #include "sbnd_ophitfinder_deco_data.fcl"
 
 process_name: DecoRecoData
@@ -43,8 +43,8 @@ physics:
   {
     opdecopmt:     @local::SBNDOpDeconvolutionPMT_data
     ophitpmt:  @local::SBNDDecoOpHitFinderPMT_data
-    opflashtpc0:   @local::SBNDDecoSimpleFlashTPC0
-    opflashtpc1:   @local::SBNDDecoSimpleFlashTPC1
+    opflashtpc0:   @local::SBNDDecoSimpleFlashTPC0_data
+    opflashtpc1:   @local::SBNDDecoSimpleFlashTPC1_data
   }
 
   reco: [

sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl

@@ -0,0 +1,27 @@
+#include "sbnd_flashfinder_data.fcl"
+
+BEGIN_PROLOG
+
+####OpFlash finder for PMT deconvolved waveforms#####
+###TPC0
+SBNDDecoSimpleFlashTPC0_data: @local::SBNDSimpleFlashTPC0_data
+SBNDDecoSimpleFlashTPC0_data.PECalib.SPEAreaGain: 200
+SBNDDecoSimpleFlashTPC0_data.OpHitProducers: ["ophitpmt"]
+SBNDDecoSimpleFlashTPC0_data.OpHitInputTime: "RiseTime"
+SBNDDecoSimpleFlashTPC0_data.UseT0Tool: true
+SBNDDecoSimpleFlashTPC0_data.ReadoutDelay: 0.135 //cable time delay in us
+SBNDDecoSimpleFlashTPC0_data.CorrectLightPropagation: true
+SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data.root"
+
+#TPC1
+SBNDDecoSimpleFlashTPC1_data: @local::SBNDSimpleFlashTPC1_data
+SBNDDecoSimpleFlashTPC1_data.PECalib.SPEAreaGain: 200
+SBNDDecoSimpleFlashTPC1_data.OpHitProducers: ["ophitpmt"]
+SBNDDecoSimpleFlashTPC1_data.OpHitInputTime: "RiseTime"
+SBNDDecoSimpleFlashTPC1_data.UseT0Tool: true
+SBNDDecoSimpleFlashTPC1_data.ReadoutDelay: 0.135 //cable time delay in us
+SBNDDecoSimpleFlashTPC1_data.CorrectLightPropagation: true
+SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data.root"
+
+
+END_PROLOG

sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc

@@ -54,6 +54,10 @@ namespace lightana
         fhicl::Comment("Group velocity for VIS photons")
       };
 
+      fhicl::Atom<bool> fDataCalibration {
+        fhicl::Name("DataCalibration"),
+        fhicl::Comment("Running on data")
+      };
     };
 
     // Default constructor
@@ -79,6 +83,9 @@ namespace lightana
     double fVGroupVIS;
     double fVGroupVUV_I;
 
+    // Whether to use data calibration for drift distance estimation
+    bool fDataCalibration;
+
     //Geo properties
     double fDriftDistance;
     double fVISLightPropTime;
@@ -168,28 +175,47 @@ namespace lightana
       }
     }
 
+    double pmtratio=-1.;
     // compute PMTRatio metric
-    double PECoated=0, PEUncoated=0;
-    for(size_t boxID=0; boxID<fPDSBoxIDs.size(); boxID++){
-      //we need the uncoated PMT in each window and at least one coated
-      if( BoxMap_NUncoatedCh[boxID]==1 && BoxMap_NCoatedCh[boxID]>=1){
-        double CoWeight = 1./BoxMap_NCoatedCh[boxID];
-        PECoated+=CoWeight * BoxMap_PECoated[boxID];
-        PEUncoated+=BoxMap_PEUncoated[boxID];
+    if(fDataCalibration)
+    {
+      double TotalPE=0;
+      double RatioPerBoxWeight=0;
+      for(size_t boxID=0; boxID<fPDSBoxIDs.size(); boxID++){
+        double RatioPerBox;
+        double PECoated=0, PEUncoated=0;
+        //we need the uncoated PMT in each window and at least one coated
+        if( BoxMap_NUncoatedCh[boxID]==1 && BoxMap_NCoatedCh[boxID]>=1){
+          PECoated+= BoxMap_PECoated[boxID];
+          PEUncoated+=BoxMap_PEUncoated[boxID];
+          TotalPE += PECoated + PEUncoated;
+          RatioPerBox = (PEUncoated/PECoated)*BoxMap_NCoatedCh[boxID];
+          RatioPerBoxWeight += RatioPerBox * (PECoated + PEUncoated);
+        }
       }
+      if(TotalPE!=0) pmtratio = RatioPerBoxWeight/TotalPE;
+    }
+    else{
+      double PECoated=0, PEUncoated=0;
+      for(size_t boxID=0; boxID<fPDSBoxIDs.size(); boxID++){
+        //we need the uncoated PMT in each window and at least one coated
+        if( BoxMap_NUncoatedCh[boxID]==1 && BoxMap_NCoatedCh[boxID]>=1){
+          double CoWeight = 1./BoxMap_NCoatedCh[boxID];
+          PECoated+=CoWeight * BoxMap_PECoated[boxID];
+          PEUncoated+=BoxMap_PEUncoated[boxID];
+        }
+      }
+      if(PECoated!=0) pmtratio = PEUncoated/PECoated;
     }
 
-    if(PECoated!=0){
-      double pmtratio = PEUncoated/PECoated;
-
+    if(pmtratio!=-1){
       double drift_distance;
       if(pmtratio<=fPMTRatioCal[0])
         drift_distance=fDriftCal[0];
       else if(pmtratio>=fPMTRatioCal[fNCalBins-1])
         drift_distance=fDriftCal[fNCalBins-1];
       else
         drift_distance=Interpolate(pmtratio);
-
       return drift_distance;
     }
     else return fDriftCal[fNCalBins-1]; 
@@ -216,11 +242,11 @@ namespace lightana
 
     size_t upix = std::upper_bound(fPMTRatioCal.begin(), fPMTRatioCal.end(), val)-fPMTRatioCal.begin();
 
-    double slope = ( fDriftCal[upix]-fDriftCal[upix] ) / ( fPMTRatioCal[upix]-fPMTRatioCal[upix-1] );
+    double slope = ( fDriftCal[upix]-fDriftCal[upix-1] ) / ( fPMTRatioCal[upix]-fPMTRatioCal[upix-1] );
 
     return fDriftCal[upix-1] + slope * ( val - fPMTRatioCal[upix-1] );
   }
 
 }
 
-DEFINE_ART_CLASS_TOOL(lightana::DriftEstimatorPMTRatio)
+DEFINE_ART_CLASS_TOOL(lightana::DriftEstimatorPMTRatio)

sbndcode/OpDetReco/OpFlash/job/sbnd_driftestimatoralgo.fcl

@@ -8,6 +8,7 @@ DriftEstimatorPMTRatio:
   CalibrationFile: "OpDetReco/PMTRatioCalibration.root"
   VGroupVUV: @local::sbnd_vuv_timing_parameterization.vuv_vgroup_mean
   VGroupVIS: @local::sbnd_vis_timing_parameterization.vis_vmean
+  DataCalibration: false
 }
 
 END_PROLOG

sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl

@@ -0,0 +1,36 @@
+#include "sbnd_flashalgo.fcl"
+#include "sbnd_flashcalib.fcl"
+#include "sbnd_flashgeoalgo.fcl"
+#include "sbnd_flasht0algo.fcl"
+#include "sbnd_driftestimatoralgo.fcl"
+
+BEGIN_PROLOG
+
+SBNDSimpleFlash:
+{
+  module_type     : "SBNDFlashFinder"
+  FlashFinderAlgo : "SimpleFlashAlgo"
+  AlgoConfig      : @local::SimpleFlashStandard
+  OpHitProducers  : ["ophitpmt","ophitarapuca"]
+  OpFlashProducer : "opflash"
+  PECalib         : @local::NoCalib
+  OpHitInputTime  : "PeakTime"
+  FlashGeoConfig  : @local::FlashGeoThreshold
+  UseT0Tool       : false
+  FlashT0Config   : @local::FlashT0SelectedChannels
+  ReadoutDelay    : 0. //in us
+  CorrectLightPropagation : false
+  DriftEstimatorConfig    : @local::DriftEstimatorPMTRatio
+}
+
+SBNDSimpleFlashTPC0_data: @local::SBNDSimpleFlash
+SBNDSimpleFlashTPC0_data.AlgoConfig: @local::SimpleFlashTPC0
+SBNDSimpleFlashTPC0_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio"
+SBNDSimpleFlashTPC0_data.DriftEstimatorConfig.DataCalibration: true
+
+SBNDSimpleFlashTPC1_data: @local::SBNDSimpleFlash
+SBNDSimpleFlashTPC1_data.AlgoConfig: @local::SimpleFlashTPC1
+SBNDSimpleFlashTPC1_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio"
+SBNDSimpleFlashTPC1_data.DriftEstimatorConfig.DataCalibration: true
+
+END_PROLOG