diff --git a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc b/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc index 9fb94f9d3..d19210e88 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 sbnd::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,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... "<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..6d5093cec 100644 --- a/sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt +++ b/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() diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Definitions.h b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Definitions.h new file mode 100644 index 000000000..4fae7e23a --- /dev/null +++ b/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 +#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 \ No newline at end of file diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc new file mode 100644 index 000000000..27151fcf8 --- /dev/null +++ b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc @@ -0,0 +1,80 @@ +sbnd::calo::NormalizeDriftSQLite::NormalizeDriftSQLite(fhicl::ParameterSet const &pset): + fDBFileName(pset.get("DBFileName")), + fDBTag(pset.get("DBTag")), + fVerbose(pset.get("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()->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; +} diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.cc b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.cc new file mode 100644 index 000000000..d20e8318a --- /dev/null +++ b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.cc @@ -0,0 +1,2 @@ +#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.h" +#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc" \ No newline at end of file diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.h b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_class.h new file mode 100644 index 000000000..dde2cfcc4 --- /dev/null +++ b/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 diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc index 1a3529c01..2a6bfa65e 100644 --- a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc +++ b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.cc @@ -1,151 +1,2 @@ -// 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 - - -sbnd::calo::NormalizeDriftSQLite::NormalizeDriftSQLite(fhicl::ParameterSet const &pset): - fDBFileName(pset.get("DBFileName")), - fDBTag(pset.get("DBTag")), - fVerbose(pset.get("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()->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; -} - +#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h" +#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Implementation.cc" \ No newline at end of file diff --git a/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h new file mode 100644 index 000000000..25b67a3cc --- /dev/null +++ b/sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_tool.h @@ -0,0 +1,5 @@ +#ifndef NORMALIZEDRIFTSQLITE_TOOL_H +#define NORMALIZEDRIFTSQLITE_TOOL_H +#include "sbndcode/Calibration/TPCCalorimetry/NormalizeDriftSQLite_Definitions.h" +DEFINE_ART_CLASS_TOOL(sbnd::calo::NormalizeDriftSQLite) +#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