diff --git a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc b/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc index 9fb94f9d3..787dd4b10 100644 --- a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc +++ b/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 ("CylinderRadius", 15); fCaloAlg = new calo::CalorimetryAlg( pset.get("CaloAlg") ); + ElifetimeTool = new calo::NormalizeDriftSQLite( pset.get("NormalizeDrift")); fCaloPlane = pset.get ("CaloPlane", 2); fCalodEdx = pset.get ("CalodEdx", 2.8); fESTAR_p0 = pset.get ("ESTAR_p0", 0.01730); fESTAR_p1 = pset.get ("ESTAR_p1", 0.00003479); - fLifetimeCorr = pset.get ("LifetimeCorrection", false); + fLifetimeCorr = pset.get ("LifetimeCorrection", true); fSCECorr = pset.get ("SCECorrection", false); fYZUniformityCorr = pset.get ("YZUniformityCorrection",true); fModBoxA = pset.get ("ModBoxA", 0.93); @@ -311,7 +313,6 @@ namespace blip { std::cout<<"\n" <<"=========== BlipRecoAlg =========================\n" <<"Event "<()->DataFor(evt); //auto const& lifetime_provider = art::ServiceHandle()->GetProvider(); //auto const& tpcCalib_provider = art::ServiceHandle()->GetProvider(); + //Elifetime + ElifetimeTool->setup(evt); // -- geometry art::ServiceHandle 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... "<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()); diff --git a/sbndcode/BlipRecoSBND/blipreco_configs.fcl b/sbndcode/BlipRecoSBND/blipreco_configs.fcl index 01329758a..abf03428f 100644 --- a/sbndcode/BlipRecoSBND/blipreco_configs.fcl +++ b/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 diff --git a/sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt b/sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt index efa89c74a..297c1dfd3 100644 --- a/sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt +++ b/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() diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc index 1a3529c01..638ccc0f7 100644 --- a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc +++ b/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 -#include -#include - -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 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 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("DBFileName")), diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h new file mode 100644 index 000000000..aca48044c --- /dev/null +++ b/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 +#include +#include + +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 fClockData; // need delayed construction + + // Cache run requests + std::map fRunInfos; + }; + + DEFINE_ART_CLASS_TOOL(NormalizeDriftSQLite) + + } // end namespace calo +} // end namespace sbnd +#endif diff --git a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl index 9d7a95dba..45b5bbb0f 100755 --- a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl +++ b/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 ] diff --git a/sbndcode/OpDetSim/CMakeLists.txt b/sbndcode/OpDetSim/CMakeLists.txt index 35c22cf77..78740226d 100755 --- a/sbndcode/OpDetSim/CMakeLists.txt +++ b/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 diff --git a/ups/product_deps b/ups/product_deps index 4efe01d36..1b16c7d38 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -253,7 +253,7 @@ wpdir product_dir wire-cell-cfg # #################################### product version qual flags -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