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 sbnd::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,13 +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() );
       if( fSCECorr ) {
@@ -1178,12 +1194,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 +1299,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_class.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,10 +29,14 @@ 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
-
+  NormalizeDrift:{
+  DBFileName: tpc_elifetime
+  DBTag: @local::SBND_Calibration_GlobalTags.tpc_elifetime_data
+  Verbose: false
+  }
   VetoBadChannels:      false            #// don't cluster hits from "bad" channels
   BadChanProducer:      "nfspl1:badchannels" #// label for vector of bad channel IDs from signal processing
   #BadChanFile:          "blipreco_badchannels.txt"

sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt

@@ -19,7 +19,9 @@ set(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_class.cc
+                 LIBRARIES  ${TOOL_LIBRARIES})
 install_headers()
 install_fhicl()
 install_source()

sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Definitions.h

@@ -0,0 +1,67 @@
+// 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

sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc

@@ -0,0 +1,80 @@
+sbnd::calo::NormalizeDriftSQLite::NormalizeDriftSQLite(fhicl::ParameterSet const &pset):
+  fDBFileName(pset.get<std::string>("DBFileName")),
+  fDBTag(pset.get<std::string>("DBTag")),
+  fVerbose(pset.get<bool>("Verbose", false)),
+  fDB(fDBFileName, "", "", fDBTag, true, false)
+{}
+
+void sbnd::calo::NormalizeDriftSQLite::configure(const fhicl::ParameterSet& pset) {}
+
+void sbnd::calo::NormalizeDriftSQLite::setup(const art::Event& e) {
+  fClockData.emplace(art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e));
+}
+
+sbnd::calo::NormalizeDriftSQLite::RunInfo sbnd::calo::NormalizeDriftSQLite::GetRunInfo(uint64_t run) {
+  // check the cache
+  if (fRunInfos.count(run)) {
+    return fRunInfos.at(run);
+  }
+
+  // Look up the run
+  //
+  // Translate the run into a fake "timestamp"
+  fDB.UpdateData((run+1000000000)*1000000000);
+
+  RunInfo thisrun;
+
+  double this_tau_E, this_tau_W;
+  fDB.GetNamedChannelData(0, "etau_sce_spatial_east", this_tau_E);
+  fDB.GetNamedChannelData(0, "etau_sce_spatial_west", this_tau_W);
+  thisrun.tau_E = this_tau_E;
+  thisrun.tau_W = this_tau_W;
+
+  if (fVerbose) std::cout << "NormalizeDriftSQLite Tool -- Lifetime Data:" << "\nTPC East: " << thisrun.tau_E << "\nTPC West: " << thisrun.tau_W << std::endl;
+
+  // Set the cache
+  fRunInfos[run] = thisrun;
+
+  return thisrun;
+}
+
+double sbnd::calo::NormalizeDriftSQLite::Normalize(double dQdx, const art::Event &e, 
+    const recob::Hit &hit, const geo::Point_t &location, const geo::Vector_t &direction, double t0) {
+
+  if (!fClockData) {
+    std::cout << "Error: fClockData is not valid" << std::endl;
+    throw cet::exception("fClockData is not valid");
+  }
+
+  // Get the info
+  RunInfo runelifetime = GetRunInfo(e.id().runID().run());
+
+  // lookup the TPC
+  double thiselifetime = -1;
+  unsigned tpc = hit.WireID().TPC;
+  unsigned cryo = hit.WireID().Cryostat;
+
+  // East
+  if (cryo == 0 && tpc == 0) thiselifetime = runelifetime.tau_E;
+
+  // West
+  if (cryo == 0 && tpc == 1) thiselifetime = runelifetime.tau_W;
+
+  // Get the hit time
+  double thit = fClockData->TPCTick2TrigTime(hit.PeakTime()) - t0;
+  thit = thit * 1.e-3;
+
+  if (fVerbose) std::cout << "NormalizeDriftSQLite Tool -- Norm factor: " << exp(thit / thiselifetime) << " at TPC: " << tpc << " Cryo: " << cryo << " Time: " << thit << " Track T0: " << t0 << ", x: " << location.X() << std::endl;
+
+  // Scale
+  if (thiselifetime > 0) {
+    dQdx = dQdx*exp(thit / thiselifetime);
+  }
+  // Throw exception if thiselifetime is not updated to non-zero value
+  else {
+    std::cout << "sbnd::calo::NormalizeDriftSQLite::Normalize electron lifetime is not found for run " << e.id().runID().run() << std::endl;
+    throw cet::exception("Electron lifetime is not found");
+  }
+
+  return dQdx;
+}

sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.cc

@@ -0,0 +1,2 @@
+#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.h"
+#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc"

sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.h

@@ -0,0 +1,4 @@
+#ifndef NORMALIZEDRIFTSQLITE_CLASS_H
+#define NORMALIZEDRIFTSQLITE_CLASS_H
+#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Definitions.h"
+#endif