Draft: Feature/blip energy drift correction

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc

@@ -219,6 +219,7 @@ namespace blip {
   BlipRecoAlg::~BlipRecoAlg()
   {
     delete fCaloAlg;
+    delete ElifetimeTool;
   }
 
 
@@ -268,11 +269,12 @@ namespace blip {
     fCylinderRadius     = pset.get<float>         ("CylinderRadius",      15);
 
     fCaloAlg            = new calo::CalorimetryAlg( pset.get<fhicl::ParameterSet>("CaloAlg") );
+    ElifetimeTool       = new calo::NormalizeDriftSQLite( pset.get<fhicl::ParameterSet>("NormalizeDrift"));
     fCaloPlane          = pset.get<int>           ("CaloPlane",           2);
     fCalodEdx           = pset.get<float>         ("CalodEdx",            2.8);
     fESTAR_p0           = pset.get<float>         ("ESTAR_p0",            0.01730);
     fESTAR_p1           = pset.get<float>         ("ESTAR_p1",            0.00003479);
-    fLifetimeCorr       = pset.get<bool>          ("LifetimeCorrection",  false);
+    fLifetimeCorr       = pset.get<bool>          ("LifetimeCorrection",  true);
     fSCECorr            = pset.get<bool>          ("SCECorrection",       false);
     fYZUniformityCorr   = pset.get<bool>          ("YZUniformityCorrection",true);
     fModBoxA            = pset.get<float>         ("ModBoxA",             0.93);
@@ -311,7 +313,6 @@ namespace blip {
     std::cout<<"\n"
     <<"=========== BlipRecoAlg =========================\n"
     <<"Event "<<evt.id().event()<<" / run "<<evt.id().run()<<"\n";
-
     //=======================================
     // Reset things
     //=======================================
@@ -337,6 +338,8 @@ namespace blip {
     //auto const& detProp              = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
     //auto const& lifetime_provider   = art::ServiceHandle<lariov::UBElectronLifetimeService>()->GetProvider();
     //auto const& tpcCalib_provider   = art::ServiceHandle<lariov::TPCEnergyCalibService>()->GetProvider();
+    //Elifetime
+    ElifetimeTool->setup(evt);
 
     // -- geometry
     art::ServiceHandle<geo::Geometry> geom;
@@ -1141,12 +1144,26 @@ namespace blip {
       // ================================================================================
       float depEl   = std::max(0.0,(double)blip.Charge);
       float Efield  = kNominalEfield;
-
+      float recomb  = ModBoxRecomb(fCalodEdx,Efield);
+      blip.EnergyNoDriftCorrection   = depEl * (1./recomb) * kWion;
+      // METHOD 2: recombination factor using dE/dx from NIST tables (dE/dx = kinetic energy / CSDA)
+      float energy_estar = Q_to_E_ESTAR(depEl);
+      float energy_pstar = Q_to_E_PSTAR(depEl);
+      blip.EnergyESTARNoDriftCorrection  = energy_estar;
+      blip.EnergyPSTARNoDriftCorrection  = energy_pstar;
       // --- Lifetime correction ---
       // Ddisabled by default. Without knowing real T0 of a blip, attempting to 
       // apply this correction can do more harm than good! Note lifetime is in
       // units of 'ms', not microseconds, hence the 1E-3 conversion factor.
-      if( fLifetimeCorr && blip.Time>0 ) depEl *= exp( 1e-3*blip.Time/detProp.ElectronLifetime());
+      double tau = 0;
+      if(plist.size()==0) //Data and should use calibration db
+      {
+        EventTPCLifetimes = ElifetimeTool->GetRunInfo(evt.id().run());
+        if(blip.TPC==0) tau = EventTPCLifetimes.tau_E;
+        else tau = EventTPCLifetimes.tau_W; //west
+      }
+      else tau = detProp.ElectronLifetime();
+      if( fLifetimeCorr ) depEl *= exp( 1e-3*blip.Time/tau);
 
       // --- SCE corrections ---
       geo::Point_t point( blip.Position.X(),blip.Position.Y(),blip.Position.Z() );
@@ -1178,12 +1195,10 @@ namespace blip {
       }
 
       // METHOD 1: recombination factor from Mod Box model with a fixed dE/dx (fCalodEdx)
-      float recomb  = ModBoxRecomb(fCalodEdx,Efield);
       blip.Energy   = depEl * (1./recomb) * kWion;
-
       // METHOD 2: recombination factor using dE/dx from NIST tables (dE/dx = kinetic energy / CSDA)
-      float energy_estar = Q_to_E_ESTAR(depEl);
-      float energy_pstar = Q_to_E_PSTAR(depEl);
+      energy_estar = Q_to_E_ESTAR(depEl);
+      energy_pstar = Q_to_E_PSTAR(depEl); //reaculate after drift correction
       blip.EnergyESTAR = energy_estar;
       blip.EnergyPSTAR = energy_pstar;
       //std::cout<<"Calculating ESTAR energy dep...  "<<depEl<<", "<<Efield<<"\n";
@@ -1285,5 +1300,4 @@ namespace blip {
     printf("\n");
 
   }
-
 }

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h

@@ -57,6 +57,8 @@
 // Blip-specific utils
 #include "sbndcode/BlipRecoSBND/Utils/BlipUtils.h"
 #include "sbnobj/SBND/Blip/BlipDataTypes.h"
+//electron lifetime includes
+#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h" //Danger of double definition
 
 // ROOT stuff
 #include "TH1D.h"
@@ -187,6 +189,8 @@ namespace blip {
     bool                fYZUniformityCorr;
     float               fModBoxA;
     float               fModBoxB;
+    sbnd::calo::NormalizeDriftSQLite* ElifetimeTool;
+    sbnd::calo::NormalizeDriftSQLite::RunInfo EventTPCLifetimes;
 
     // --- Splines for calo
     TSpline3* spline_PSTAR;

sbndcode/BlipRecoSBND/Alg/CMakeLists.txt

@@ -36,6 +36,7 @@ art_make_library(LIBRARY_NAME  sbndcode_BlipRecoAlg
                            sbnobj::Common_CRT
                            sbndcode_BlipUtils
 			               sbnobj::SBND_Blip
+                        sbndcode_GIMME_LIFETIMES
 )
 
 install_headers()

sbndcode/BlipRecoSBND/BlipAna_module.cc

@@ -256,6 +256,9 @@ class BlipAnaTreeDataStruct
   float blip_energy[kMaxBlips];       // blip reco energy [MeV]
   float blip_energy_estar[kMaxBlips]; // blip reco energy [MeV], based on ESTAR
   float blip_energy_pstar[kMaxBlips]; // blip reco energy [MeV], based on PSTAR
+  float blip_energy_nodriftcorr[kMaxBlips];       // blip reco energy [MeV]
+  float blip_energy_estar_nodriftcorr[kMaxBlips]; // blip reco energy [MeV], based on ESTAR
+  float blip_energy_pstar_nodriftcorr[kMaxBlips]; // blip reco energy [MeV], based on PSTAR
   float blip_energyTrue[kMaxBlips];   // blip truth energy [MeV]
   float blip_yzcorr[kMaxBlips];       // YZ uniformity correction factor (already applied)
   bool  blip_isMC[kMaxBlips];         // blip is matched to MC particle
@@ -408,6 +411,9 @@ class BlipAnaTreeDataStruct
     FillWith(blip_energy,     -999);
     FillWith(blip_energy_estar, -999);
     FillWith(blip_energy_pstar, -999);
+    FillWith(blip_energy_nodriftcorr,     -999);
+    FillWith(blip_energy_estar_nodriftcorr, -999);
+    FillWith(blip_energy_pstar_nodriftcorr, -999);
     FillWith(blip_energyTrue, -999);
     FillWith(blip_yzcorr,     -9);
     FillWith(blip_proxtrkdist,-99);
@@ -522,6 +528,9 @@ class BlipAnaTreeDataStruct
     evtTree->Branch("blip_energy",blip_energy,"blip_energy[nblips]/F");
     evtTree->Branch("blip_energy_estar",blip_energy_estar,"blip_energy_estar[nblips]/F");
     evtTree->Branch("blip_energy_pstar",blip_energy_pstar,"blip_energy_pstar[nblips]/F");
+    evtTree->Branch("blip_energy_nodriftcorr",blip_energy_nodriftcorr,"blip_energy_nodriftcorr[nblips]/F");
+    evtTree->Branch("blip_energy_estar_nodriftcorr",blip_energy_estar_nodriftcorr,"blip_energy_estar_nodriftcorr[nblips]/F");
+    evtTree->Branch("blip_energy_pstar_nodriftcorr",blip_energy_pstar_nodriftcorr,"blip_energy_pstar_nodriftcorr[nblips]/F");
     //evtTree->Branch("blip_yzcorr",blip_yzcorr,"blip_yzcorr[nblips]/F");
     //evtTree->Branch("blip_energyTrue",blip_energyTrue,"blip_energyTrue[nblips]/F");
     evtTree->Branch("blip_incylinder",blip_incylinder,"blip_incylinder[nblips]/O");
@@ -1369,6 +1378,9 @@ void BlipAna::analyze(const art::Event& evt)
     fData->blip_energy[i]      = blp.Energy;
     fData->blip_energy_estar[i]= blp.EnergyESTAR;
     fData->blip_energy_pstar[i]= blp.EnergyPSTAR;
+    fData->blip_energy_nodriftcorr[i]      = blp.EnergyNoDriftCorrection;
+    fData->blip_energy_estar_nodriftcorr[i]= blp.EnergyESTARNoDriftCorrection;
+    fData->blip_energy_pstar_nodriftcorr[i]= blp.EnergyPSTARNoDriftCorrection;
 
     //fData->blip_yzcorr[i]     = tpcCalib.YZdqdxCorrection(fCaloPlane,blp.Position.Y(),blp.Position.Z());
 

sbndcode/BlipRecoSBND/blipreco_configs.fcl

@@ -29,7 +29,7 @@ sbnd_blipalg:
   CalodEdx:             2.8             #// assumed dE/dx for MeV-scale recombination corrections
   ESTAR_p0:             0.01730         #// p0 for electrons to energy conversion function using ESTAR table
   ESTAR_p1:		0.00003479	#// p1 for electrons to	energy conversion function using ESTAR table
-  LifetimeCorrection:   false           #// correct blips for lifetime (assumes T0 = 0)
+  LifetimeCorrection:   true           #// correct blips for lifetime (assumes T0 = 0)
   SCECorrection:        false           #// correct blips for space charge (requires known XYZ)
   CaloAlg:              @local::sbnd_calorimetryalgmc
 

sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt

@@ -17,9 +17,11 @@ set(TOOL_LIBRARIES
   wda::wda
 )
 
-cet_build_plugin(NormalizeDriftSQLite art::tool LIBRARIES ${TOOL_LIBRARIES})
+#cet_build_plugin(NormalizeDriftSQLite art::tool LIBRARIES ${TOOL_LIBRARIES})
 cet_build_plugin(NormalizeYZ art::tool LIBRARIES ${TOOL_LIBRARIES})
-
+cet_make_library(LIBRARY_NAME sbndcode_GIMME_LIFETIMES
+                 SOURCE    NormalizeDriftSQLite_tool.cc
+                 LIBRARIES  ${TOOL_LIBRARIES})
 install_headers()
 install_fhicl()
 install_source()

sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc

@@ -1,72 +1,4 @@
-// Framework Includes
-#include "art/Framework/Core/EDProducer.h"
-#include "art/Framework/Principal/Event.h"
-#include "art/Framework/Principal/Handle.h"
-#include "art/Framework/Services/Registry/ServiceHandle.h"
-#include "art/Persistency/Common/PtrMaker.h"
-#include "art/Utilities/ToolMacros.h"
-#include "cetlib_except/exception.h"
-#include "cetlib/cpu_timer.h"
-#include "fhiclcpp/ParameterSet.h"
-#include "messagefacility/MessageLogger/MessageLogger.h"
-
-#include "larevt/CalibrationDBI/Providers/DBFolder.h"
-
-// Tool include
-#include "larreco/Calorimetry/INormalizeCharge.h"
-
-// Services
-#include "lardata/DetectorInfoServices/DetectorClocksService.h"
-
-// Lab helpers
-#include "wda.h"
-
-// C++
-#include <string>
-#include <optional>
-#include <cassert>
-
-namespace sbnd {
-  namespace calo {
-
-class NormalizeDriftSQLite : public INormalizeCharge
-{
-public:
-  NormalizeDriftSQLite(fhicl::ParameterSet const &pset);
-
-  void configure(const fhicl::ParameterSet& pset) override;
-  void setup(const art::Event& e) override;
-  double Normalize(double dQdx, const art::Event &e, const recob::Hit &h, const geo::Point_t &location, const geo::Vector_t &direction, double t0) override;
-
-private:
-  // Configuration
-  std::string fDBFileName;
-  std::string fDBTag;
-  bool fVerbose;
-
-  lariov::DBFolder fDB;
-
-  std::optional<detinfo::DetectorClocksData> fClockData; // need delayed construction
-
-  // Class to hold data from DB
-  class RunInfo {
-  public:
-    double tau_E;
-    double tau_W;
-  };
-
-  // Helpers
-  RunInfo GetRunInfo(uint64_t run);
-
-  // Cache run requests
-  std::map<uint32_t, RunInfo> fRunInfos;
-};
-
-DEFINE_ART_CLASS_TOOL(NormalizeDriftSQLite)
-
-  } // end namespace calo
-} // end namespace sbnd
-
+#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h"
 
 sbnd::calo::NormalizeDriftSQLite::NormalizeDriftSQLite(fhicl::ParameterSet const &pset):
   fDBFileName(pset.get<std::string>("DBFileName")),

sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h

@@ -0,0 +1,70 @@
+#ifndef NORMALIZEDRIFTSQLITE_TOOL_H
+#define NORMALIZEDRIFTSQLITE_TOOL_H
+// Framework Includes
+//#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "art/Persistency/Common/PtrMaker.h"
+#include "art/Utilities/ToolMacros.h"
+#include "cetlib_except/exception.h"
+#include "cetlib/cpu_timer.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "larevt/CalibrationDBI/Providers/DBFolder.h"
+
+// Tool include
+#include "larreco/Calorimetry/INormalizeCharge.h"
+
+// Services
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+
+// Lab helpers
+//#include "wda.h"
+
+// C++
+#include <string>
+#include <optional>
+#include <cassert>
+
+namespace sbnd {
+  namespace calo {
+
+    class NormalizeDriftSQLite : public INormalizeCharge
+    {
+    public:
+      NormalizeDriftSQLite(fhicl::ParameterSet const &pset);
+
+      void configure(const fhicl::ParameterSet& pset) override;
+      void setup(const art::Event& e) override;
+      double Normalize(double dQdx, const art::Event &e, const recob::Hit &h, const geo::Point_t &location, const geo::Vector_t &direction, double t0) override;
+      // Class to hold data from DB
+      class RunInfo {
+      public:
+	double tau_E;
+	double tau_W;
+      };
+
+      // Helpers
+      RunInfo GetRunInfo(uint64_t run);
+
+    private:
+      // Configuration
+      std::string fDBFileName;
+      std::string fDBTag;
+      bool fVerbose;
+
+      lariov::DBFolder fDB;
+
+      std::optional<detinfo::DetectorClocksData> fClockData; // need delayed construction
+
+      // Cache run requests
+      std::map<uint32_t, RunInfo> fRunInfos;
+    };
+
+    DEFINE_ART_CLASS_TOOL(NormalizeDriftSQLite)
+
+      } // end namespace calo
+} // end namespace sbnd
+#endif

sbndcode/JobConfigurations/standard/reco/reco2_data.fcl

@@ -63,8 +63,9 @@ physics.producers.opt0finderSCE.CaloProducer:     "pandoraSCECaloData"
 physics.producers.pandoraSCEShower.PFParticleLabel:              "pandoraSCE"
 physics.producers.pandoraSCEShowerSBN.PFParticleLabel:           "pandoraSCE"
 
-
+physics.producers.blipreco.BlipAlg.CaloAlg: @local::sbnd_calorimetryalgdata
 
 physics.caloskimana: [ caloskim, crtana ]
 physics.ana:            @local::sbnd_reco2_analyzer_sequence
+physics.analyzers.blipana.BlipAlg.CaloAlg: @local::sbnd_calorimetryalgdata
 physics.end_paths: [stream1, caloskimana, ana ]

sbndcode/OpDetSim/CMakeLists.txt

@@ -9,7 +9,7 @@ cet_build_plugin( sbndPDMapAlg art::tool
                     lardataobj::RawData
                     lardata::DetectorInfoServices_DetectorClocksServiceStandard_service
                     larpandora::LArPandoraInterface
-                    sbncode::OpDet_PDMapAlgSimple_tool
+                    sbncode::OpDet_PDMapAlg
                     sbndcode_Utilities_SignalShapingServiceSBND_service
                     nurandom::RandomUtils_NuRandomService_service
                     art::Framework_Services_Optional_RandomNumberGenerator_service

ups/product_deps

@@ -253,7 +253,7 @@ wpdir   product_dir     wire-cell-cfg
 #
 ####################################
 product		version		qual	flags		<table_format=2>
-sbncode		v10_21_00       -
+sbncode		v10_21_02       -
 cetmodules	v3_24_01	-	only_for_build
 sbnd_data	v01_42_00	-	
 sbndutil        v10_06_01    -       optional