Light Calorimetry

CMakeLists.txt

@@ -15,7 +15,7 @@
 
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
-set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.14.02)
+set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.14.02.01)
 find_package(cetmodules REQUIRED)
 project(sbndcode LANGUAGES CXX)
 

sbndcode/Calorimetry/CMakeLists.txt

@@ -4,6 +4,8 @@ set (MODULE_LIBRARIES
                 larsim::LegacyLArG4
                 larcorealg::Geometry
                 larcore::Geometry_Geometry_service
+                larevt::CalibrationDBI_Providers
+                wda::wda
                 lardata::Utilities
                 lardataobj::RecoBase
                 art::Framework_Core

sbndcode/Calorimetry/LightCaloProducer_module.cc

@@ -44,6 +44,10 @@
 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
 
+// Calibration database includes for electron lifetime
+#include "larevt/CalibrationDBI/Providers/DBFolder.h"
+#include "wda.h"
+
 // SBND includes
 #include "sbndcode/OpDetSim/sbndPDMapAlg.hh"
 #include "sbnobj/Common/Reco/LightCalo.h"
@@ -121,6 +125,15 @@ class sbnd::LightCaloProducer : public art::EDProducer {
   // Returns the weighted mean of the light vector 
   double CalcMean(std::vector<double> total_light, std::vector<double> total_err);
 
+  // Struct to hold electron lifetime data from database
+  struct ELifetimeInfo {
+    double tau_tpc0;
+    double tau_tpc1;
+  };
+
+  // Helper to get electron lifetime from database
+  ELifetimeInfo GetELifetimeFromDB(uint64_t run);
+
   // fcl parameters 
   std::vector<std::string> fopflash_producer_v;
   std::vector<std::string> fopflash_ara_producer_v;
@@ -159,6 +172,12 @@ class sbnd::LightCaloProducer : public art::EDProducer {
   double fxarapucavuv_viseff;
   double fxarapucavis_eff; 
 
+  bool fuse_elifetime_db;
+  std::string felifetime_db_file;
+  std::string felifetime_db_tag;
+  std::unique_ptr<lariov::DBFolder> felifetime_db;
+  std::map<uint32_t, ELifetimeInfo> felifetime_cache;
+
   std::vector<float> fopdet_vuv_eff;
   std::vector<float> fopdet_vis_eff;
   std::vector<int>   fopdet_mask;
@@ -237,6 +256,14 @@ sbnd::LightCaloProducer::LightCaloProducer(fhicl::ParameterSet const& p)
   fxarapucavuv_viseff = p.get<double>("XArapucaVUVVISEff");
   fxarapucavis_eff    = p.get<double>("XArapucaVISEff");
 
+  // Electron lifetime database configuration
+  fuse_elifetime_db = p.get<bool>("UseELifetimeDB", false);
+  if (fuse_elifetime_db) {
+    felifetime_db_file = p.get<std::string>("ELifetimeDBFile");
+    felifetime_db_tag  = p.get<std::string>("ELifetimeDBTag");
+    felifetime_db = std::make_unique<lariov::DBFolder>(felifetime_db_file, "", "", felifetime_db_tag, true, false);
+  }
+
   // fill efficiency vectors 
   fopdet_vuv_eff.resize(nchan,0.);
   fopdet_vis_eff.resize(nchan,0.);
@@ -269,21 +296,22 @@ sbnd::LightCaloProducer::LightCaloProducer(fhicl::ParameterSet const& p)
   geom = lar::providerFrom<geo::Geometry>();
 
   art::ServiceHandle<art::TFileService> fs;
-  _tree = fs->make<TTree>("lightcalo","");
-  _tree->Branch("run",           &_run,          "run/I");
-  _tree->Branch("subrun",        &_subrun,       "subrun/I");
-  _tree->Branch("event",         &_event,        "event/I");
-  _tree->Branch("pfpid",         &_pfpid,        "pfpid/I");
-  _tree->Branch("opflash_time",  &_opflash_time, "opflash_time/D");
-
-  _tree->Branch("rec_gamma",    "std::vector<double>", &_rec_gamma);
-  _tree->Branch("dep_pe",       "std::vector<double>", &_dep_pe);
-  _tree->Branch("visibility",   "std::vector<double>", &_visibility);
-
-  _tree->Branch("slice_Q",       &_slice_Q,      "slice_Q/D");
-  _tree->Branch("slice_L",       &_slice_L,      "slice_L/D");
-  _tree->Branch("slice_E",       &_slice_E,      "slice_E/D");
-
+  if (ffill_tree){
+    _tree = fs->make<TTree>("lightcalo","");
+    _tree->Branch("run",           &_run,          "run/I");
+    _tree->Branch("subrun",        &_subrun,       "subrun/I");
+    _tree->Branch("event",         &_event,        "event/I");
+    _tree->Branch("pfpid",         &_pfpid,        "pfpid/I");
+    _tree->Branch("opflash_time",  &_opflash_time, "opflash_time/D");
+
+    _tree->Branch("rec_gamma",    "std::vector<double>", &_rec_gamma);
+    _tree->Branch("dep_pe",       "std::vector<double>", &_dep_pe);
+    _tree->Branch("visibility",   "std::vector<double>", &_visibility);
+
+    _tree->Branch("slice_Q",       &_slice_Q,      "slice_Q/D");
+    _tree->Branch("slice_L",       &_slice_L,      "slice_L/D");
+    _tree->Branch("slice_E",       &_slice_E,      "slice_E/D");
+  }
   // Call appropriate produces<>() functions here.
   produces<std::vector<sbn::LightCalo>>();
   produces<art::Assns<recob::Slice,   sbn::LightCalo>>();
@@ -325,21 +353,21 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
   ::art::Handle<std::vector<recob::Slice>> slice_h;
   e.getByLabel(fslice_producer, slice_h);
   if(!slice_h.isValid() || slice_h->empty()){
-    std::cout << "don't have good slices!" << std::endl;
+    std::cout << "[LightCaloProducer] : " << fslice_producer << " doesn't have good slices!" << std::endl;
     return;
   }
 
   ::art::Handle<std::vector<recob::PFParticle>> pfp_h;
   e.getByLabel(fslice_producer, pfp_h);
   if(!pfp_h.isValid() || pfp_h->empty()) {
-    std::cout << "don't have good PFParticle!" << std::endl;
+    std::cout << "[LightCaloProducer] : " << fslice_producer << " doesn't have good PFParticles!" << std::endl;
     return;
   }
 
   ::art::Handle<std::vector<recob::SpacePoint>> spacepoint_h;
   e.getByLabel(fslice_producer, spacepoint_h);
   if(!spacepoint_h.isValid() || spacepoint_h->empty()) {
-    std::cout << "don't have good SpacePoints!" << std::endl;
+    std::cout << "[LightCaloProducer] : " << fslice_producer << " don't have good SpacePoints!" << std::endl;
     return;
   }
 
@@ -350,7 +378,7 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
     ::art::Handle<std::vector<recob::OpFlash>> flash_h;
     e.getByLabel(fopflash_producer_v[i], flash_h);
     if (!flash_h.isValid() || flash_h->empty()) {
-      std::cout << "don't have good PMT flashes from producer " << fopflash_producer_v[i] << std::endl;
+      std::cout << "[LightCaloProducer] : " << "don't have good PMT flashes from producer " << fopflash_producer_v[i] << std::endl;
     }
     else{
       if (fopflash_producer_v[i].find("tpc0") != std::string::npos)
@@ -376,7 +404,7 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
     ::art::Handle<std::vector<sbn::TPCPMTBarycenterMatch>> bcfm_h;
     e.getByLabel(fbcfm_producer, bcfm_h);
     if(!bcfm_h.isValid() || bcfm_h->empty()) {
-      std::cout << "don't have good barycenter matches!" << std::endl;
+      std::cout << "[LightCaloProducer] : " << "don't have good barycenter matches!" << std::endl;
       return;
     }
     std::vector<art::Ptr<sbn::TPCPMTBarycenterMatch>> bcfm_v;
@@ -385,15 +413,15 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
     CollectMatches(bcfm_h, bcfm_v, fbcfm_producer, e, match_slices_v, match_op0, match_op1,
                     [this](art::Ptr<sbn::TPCPMTBarycenterMatch> bcfm) {
                       if (bcfm->flashTime > fopflash_max || bcfm->flashTime < fopflash_min) return false;
-                      if (bcfm->score < 0.02) return false;
+                      if (bcfm->score < fbcfmscore_cut) return false;
                       return true;
                     });
   }
   else if (fuse_opt0){
     ::art::Handle<std::vector<sbn::OpT0Finder>> opt0_h;
     e.getByLabel(fopt0_producer, opt0_h);
     if(!opt0_h.isValid() || opt0_h->empty()) {
-      std::cout << "don't have good OpT0Finder matches!" << std::endl;
+      std::cout << "[LightCaloProducer] : " << "don't have good OpT0Finder matches!" << std::endl;
       return;
     }
     std::vector<art::Ptr<sbn::OpT0Finder>> opt0_v;
@@ -430,7 +458,7 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
       ::art::Handle<std::vector<recob::OpFlash>> flash_ara_h;
       e.getByLabel(fopflash_ara_producer_v[i], flash_ara_h);
       if (!flash_ara_h.isValid() || flash_ara_h->empty()) {
-        std::cout << "don't have good X-ARAPUCA flashes from producer " << fopflash_ara_producer_v[i] << std::endl;
+        std::cout << "[LightCaloProducer] : " << "don't have good X-ARAPUCA flashes from producer " << fopflash_ara_producer_v[i] << std::endl;
       }
       else{
         if (fopflash_ara_producer_v[i].find("tpc0") != std::string::npos)
@@ -441,6 +469,17 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
     }
   }
 
+  // Get electron lifetime (from database if enabled, otherwise from detector properties)
+  double elifetime_tpc0, elifetime_tpc1;
+  if (fuse_elifetime_db) {
+    ELifetimeInfo elife_info = GetELifetimeFromDB(e.id().run());
+    elifetime_tpc0 = elife_info.tau_tpc0;
+    elifetime_tpc1 = elife_info.tau_tpc1;
+  } else {
+    elifetime_tpc0 = det_prop.ElectronLifetime();
+    elifetime_tpc1 = det_prop.ElectronLifetime();
+  }
+
   for (size_t n_slice=0; n_slice < match_slices_v.size(); n_slice++){
     // initialize tree variables 
     _pfpid = -1;
@@ -467,8 +506,9 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
 
     for (size_t i=0; i < slice_hits_v.size(); i++){
       auto hit = slice_hits_v[i];
-      auto drift_time = hit->PeakTime()*0.5 - clock_data.TriggerOffsetTPC(); 
-      double atten_correction = std::exp(drift_time/det_prop.ElectronLifetime()); // exp(us/us)
+      auto drift_time = clock_data.TPCTick2TrigTime(hit->PeakTime());
+      double elifetime = (hit->WireID().TPC == 0) ? elifetime_tpc0 : elifetime_tpc1;
+      double atten_correction = std::exp(drift_time/elifetime); // exp(us/us)
       auto hit_plane = hit->View();
       plane_charge.at(hit_plane) += hit->Integral()*atten_correction*(1/fcal_area_const.at(hit_plane));
       plane_hits.at(hit_plane)++; 
@@ -494,9 +534,9 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
           if (hit->View() !=bestHits) continue;
           const auto &position(sp->XYZ());
           geo::Point_t xyz(position[0],position[1],position[2]);
-          // correct for e- attenuation 
-          auto drift_time = hit->PeakTime()*0.5 - clock_data.TriggerOffsetTPC(); 
-          double atten_correction = std::exp(drift_time/det_prop.ElectronLifetime()); // exp(us/us)
+          auto drift_time = clock_data.TPCTick2TrigTime(hit->PeakTime());
+          double elifetime = (hit->WireID().TPC == 0) ? elifetime_tpc0 : elifetime_tpc1;
+          double atten_correction = std::exp(drift_time/elifetime); // exp(us/us)
           double charge = (1/fcal_area_const.at(hit->View()))*atten_correction*hit->Integral();
           if (xyz.X() < 0) has_sp0 = true;
           else             has_sp1 = true;
@@ -515,7 +555,7 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
     // find matching flashes if slice has reco spacepoints in both TPCs
     if ((has_sp0 && has_sp1) && (opflash0.isNull() || opflash1.isNull())){
       if (opflash0.isNull() && opflash1.isNull()){
-        if (fverbose) std::cout << "No opflashes found for slice with spacepoints in both TPCs." << std::endl;
+        if (fverbose) std::cout << "[LightCaloProducer] : " << "No opflashes found for slice with spacepoints in both TPCs." << std::endl;
         return;
       }
       else if (opflash0.isNull()) opflash0 = FindMatchingFlash(flash0_v, opflash1->Time());
@@ -531,7 +571,7 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
       flash_time = opflash1->Time(); 
     }
 
-    if (flash_in_0==false && flash_in_1==false && fverbose)
+    if (flash_in_0==false && flash_in_1==false)
       return;
 
     // get total L count
@@ -597,9 +637,9 @@ void sbnd::LightCaloProducer::CalculateCalorimetry(art::Event& e,
     _slice_E = (_slice_L + _slice_Q)*1e-9*g4param->Wph(); // GeV, Wph = 19.5 eV  
 
     if (fverbose){
-      std::cout << "charge: " << _slice_Q << std::endl;
-      std::cout << "light:  " << _slice_L << std::endl;
-      std::cout << "energy: " << _slice_E << std::endl;
+      std::cout << "[LightCaloProducer] : " << "charge: " << _slice_Q << std::endl;
+      std::cout << "[LightCaloProducer] : " << "light:  " << _slice_L << std::endl;
+      std::cout << "[LightCaloProducer] : " << "energy: " << _slice_E << std::endl;
     }
 
     sbn::LightCalo lightcalo{_slice_Q,_slice_L,_slice_E,bestHits};
@@ -657,7 +697,7 @@ void sbnd::LightCaloProducer::CollectMatches(const art::Handle<std::vector<Match
     auto slice = it->first;
     auto flash_v = it->second;
     if (flash_v.size() > 2){
-      std::cout << "more than two opflash matched to this slice!" << std::endl;
+      std::cout << "[LightCaloProducer] : " << "more than two opflash matched to this slice!" << std::endl;
       continue;
     }
     bool found_opflash0 = false;
@@ -692,7 +732,7 @@ void sbnd::LightCaloProducer::CollectMatches(const art::Handle<std::vector<Match
 art::Ptr<recob::OpFlash> sbnd::LightCaloProducer::FindMatchingFlash(const std::vector<art::Ptr<recob::OpFlash>> &flash_v,
                                                                     double ref_time)
 {
-  double best_dt = 1.e9;
+  double best_dt = std::numeric_limits<double>::max();  
   art::Ptr<recob::OpFlash> best;
   for (size_t i = 0; i < flash_v.size(); ++i) {
     art::Ptr<recob::OpFlash> cand = flash_v[i];
@@ -778,4 +818,33 @@ double sbnd::LightCaloProducer::CalcMean(std::vector<double> total_gamma, std::v
   return total_mean; 
 }
 
+sbnd::LightCaloProducer::ELifetimeInfo sbnd::LightCaloProducer::GetELifetimeFromDB(uint64_t run) {
+  // Check cache first
+  if (felifetime_cache.count(run)) {
+    return felifetime_cache.at(run);
+  }
+
+  // Query database - translate run into fake "timestamp" 
+  // (same convention as NormalizeDriftSQLite)
+  felifetime_db->UpdateData((run + 1000000000) * 1000000000);
+
+  ELifetimeInfo info;
+  double tau_E, tau_W;
+  felifetime_db->GetNamedChannelData(0, "etau_sce_spatial_east", tau_E);
+  felifetime_db->GetNamedChannelData(0, "etau_sce_spatial_west", tau_W);
+
+  info.tau_tpc0 = tau_E*1e3; // the db value is in ms, convert to us
+  info.tau_tpc1 = tau_W*1e3; // the db value is in ms, convert to us
+
+  if (fverbose) {
+    std::cout << "[LightCaloProducer] : Electron lifetime from DB for run " << run << std::endl;
+    std::cout << "[LightCaloProducer] : TPC0 (East): " << info.tau_tpc0 << " us" << std::endl;
+    std::cout << "[LightCaloProducer] : TPC1 (West): " << info.tau_tpc1 << " us" << std::endl;
+  }
+
+  // Cache the result
+  felifetime_cache[run] = info;
+  return info;
+}
+
 DEFINE_ART_MODULE(sbnd::LightCaloProducer)

sbndcode/Calorimetry/lightcalo_sbnd_data.fcl

@@ -1,11 +1,16 @@
 #include "lightcalo_sbnd.fcl"
+#include "calibration_database_GlobalTags_sbnd.fcl"
 
 BEGIN_PROLOG
 sbnd_lightcalo_data: @local::sbnd_lightcalo
 sbnd_lightcalo_data.CalAreaConstants: @local::sbnd_calorimetryalgdata.CalAreaConstants
 sbnd_lightcalo_data.OpFlashMin: -0.8 # assumption is light-triggered data
 sbnd_lightcalo_data.OpFlashMax: 2.0  # assumption is light-triggered data
 
+sbnd_lightcalo_data.UseELifetimeDB:  true
+sbnd_lightcalo_data.ELifetimeDBFile: "tpc_elifetime"
+sbnd_lightcalo_data.ELifetimeDBTag:  @local::SBND_Calibration_GlobalTags.tpc_elifetime_data
+
 sbnd_lightcalo_data_sce: @local::sbnd_lightcalo_data
 sbnd_lightcalo_data_sce.SliceProducer:       "pandoraSCE"
 sbnd_lightcalo_data_sce.OpT0FinderProducer:  "opt0finderSCE"

sbndcode/Calorimetry/run_lightcalo.fcl

@@ -28,7 +28,7 @@ outputs:
     out1:
     {
         module_type: RootOutput
-        fileName: "lightcalo-art.root"
+        fileName: "%ifb_lightcalo-art.root"
         dataTier: "reconstructed"
         compressionLevel:1
     }

sbndcode/JobConfigurations/standard/reco/reco2_data.fcl

@@ -2,8 +2,8 @@
 #include "sbnd_tpcpmt3dbarycentermatching_config.fcl"
 #include "sbnd_lightpropagationcorrection_config.fcl"
 #include "frameshift_sbnd_data.fcl"
-#include "standard_reco2_sbnd.fcl"
 #include "lightcalo_sbnd_data.fcl"
+#include "standard_reco2_sbnd.fcl"
 
 services:
 {

ups/product_deps

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