From 9e2397d6af42139b9f4317712cf4ea71d78d9733 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Thu, 6 Nov 2025 01:35:54 -0600 Subject: [PATCH 01/28] Implement data-based gaussian gain fluctuations --- sbndcode/OpDetSim/PMTAlg/CMakeLists.txt | 17 ++++- .../OpDetSim/PMTAlg/PMTGainFluctuations.hh | 9 ++- .../PMTAlg/PMTGaussianGainFluctuation_tool.cc | 62 +++++++++++++++++++ .../PMTAlg/pmtgainfluctuations_config.fcl | 5 ++ 4 files changed, 91 insertions(+), 2 deletions(-) create mode 100644 sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc diff --git a/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt b/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt index 115aa504d..9a6003224 100644 --- a/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt +++ b/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt @@ -5,6 +5,13 @@ cet_build_plugin(PMTGainFluctuations1Dynode art::tool nurandom::RandomUtils_NuRandomService_service ) +cet_build_plugin(PMTGaussianGainFluctuation art::tool + SOURCE + PMTGaussianGainFluctuation_tool.cc + LIBRARIES + nurandom::RandomUtils_NuRandomService_service + larcore::Geometry_Geometry_service +) cet_build_plugin(PMTNonLinearityTF1 art::tool SOURCE @@ -13,8 +20,16 @@ cet_build_plugin(PMTNonLinearityTF1 art::tool ROOT::Hist ) +#cet_build_plugin(PMTNonLinearityTF1ChannelByChannel art::tool +# SOURCE +# PMTNonLinearityTF1ChannelByChannel_tool.cc +# LIBRARIES +# art::Framework_Services_Registry +# ROOT::Hist +# larcore::Geometry_Geometry_service +#) install_headers() install_fhicl() install_source() -FILE(GLOB fcl_files *.fcl) +FILE(GLOB fcl_files *.fcl) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh b/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh index 6d5a635f9..1c111e593 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh +++ b/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh @@ -20,7 +20,14 @@ public: virtual ~PMTGainFluctuations() noexcept = default; //Returns fluctuated factor for SPR - virtual double GainFluctuation(unsigned int npe, CLHEP::HepRandomEngine* eng) = 0; + virtual double GainFluctuation(unsigned int npe, CLHEP::HepRandomEngine* eng){ + // Default implementation, can be overridden + return npe; + } + virtual double GainFluctuation(int ch, unsigned int npe, CLHEP::HepRandomEngine* eng) { + // Default implementation, can be overridden + return GainFluctuation(npe, eng); + } }; #endif diff --git a/sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc new file mode 100644 index 000000000..c413cd109 --- /dev/null +++ b/sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc @@ -0,0 +1,62 @@ +//////////////////////////////////////////////////////////////////////// +// Specific class tool for PMTGainFluctuations +// File: PMTGaussianGainFluctuation.hh +// Base class: PMTGainFluctuations.hh +// Algorithm based on function +// 'multiplicationStageGain(unsigned int i /* = 1 */) const' +// in icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +//////////////////////////////////////////////////////////////////////// + +#include "fhiclcpp/ParameterSet.h" +#include "art/Utilities/ToolMacros.h" +#include "art/Utilities/make_tool.h" +#include "art/Utilities/ToolConfigTable.h" +#include "nurandom/RandomUtils/NuRandomService.h" +#include "CLHEP/Random/RandGaussQ.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" + +#include + +#include "sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" + +namespace opdet { + class PMTGaussianGainFluctuation; +} + + +class opdet::PMTGaussianGainFluctuation : opdet::PMTGainFluctuations { +public: + + //Configuration parameters + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + }; + + explicit PMTGaussianGainFluctuation(art::ToolConfigTable const& config); + + //Returns fluctuated factor for SPR + double GainFluctuation(int ch, unsigned int npe, CLHEP::HepRandomEngine* eng) override; + +private: + + //PMTCalibrationDatabase service + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; +}; + + +opdet::PMTGaussianGainFluctuation::PMTGaussianGainFluctuation(art::ToolConfigTable const& config) +{ + fPMTCalibrationDatabaseService = lar::providerFrom(); +} + +double opdet::PMTGaussianGainFluctuation::GainFluctuation(int ch, unsigned int npe, CLHEP::HepRandomEngine* eng){ + double spe_amp = fPMTCalibrationDatabaseService->getSPEAmplitude(ch); + double ChannelGainFluctuation = fPMTCalibrationDatabaseService->getSPEAmplitudeStd(ch); + return CLHEP::RandGaussQ::shoot(eng, npe, std::sqrt(npe)*(ChannelGainFluctuation/spe_amp)); +} + +DEFINE_ART_CLASS_TOOL(opdet::PMTGaussianGainFluctuation) diff --git a/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl b/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl index 5b611b525..59e59b535 100644 --- a/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl +++ b/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl @@ -18,4 +18,9 @@ FirstDynodeGainFluctuations: Gain: 1e7 # total typical PMT gain } +GaussianGainFluctuations: +{ + tool_type: "PMTGaussianGainFluctuation" +} + END_PROLOG From 1ae5c2a1aa872f71f756775581197deb8ba0c0ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 10 Nov 2025 04:01:30 -0600 Subject: [PATCH 02/28] Add channel-by-channel non linearity --- sbndcode/OpDetSim/PMTAlg/CMakeLists.txt | 16 +-- sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh | 9 +- ...PMTNonLinearityTF1ChannelByChannel_tool.cc | 123 ++++++++++++++++++ .../PMTAlg/pmtnonlinearity_config.fcl | 11 ++ 4 files changed, 150 insertions(+), 9 deletions(-) create mode 100644 sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc diff --git a/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt b/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt index 9a6003224..e97172251 100644 --- a/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt +++ b/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt @@ -20,14 +20,14 @@ cet_build_plugin(PMTNonLinearityTF1 art::tool ROOT::Hist ) -#cet_build_plugin(PMTNonLinearityTF1ChannelByChannel art::tool -# SOURCE -# PMTNonLinearityTF1ChannelByChannel_tool.cc -# LIBRARIES -# art::Framework_Services_Registry -# ROOT::Hist -# larcore::Geometry_Geometry_service -#) +cet_build_plugin(PMTNonLinearityTF1ChannelByChannel art::tool + SOURCE + PMTNonLinearityTF1ChannelByChannel_tool.cc + LIBRARIES + art::Framework_Services_Registry + ROOT::Hist + larcore::Geometry_Geometry_service +) install_headers() install_fhicl() diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh index 13b4bf04e..762c8c2a2 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh @@ -20,7 +20,14 @@ public: virtual ~PMTNonLinearity() noexcept = default; //Returns rescaled number of PE - virtual double NObservedPE(size_t bin, std::vector & pe_vector) = 0; + virtual double NObservedPE(size_t bin, std::vector & pe_vector){ + return pe_vector[bin]; + } + virtual double NObservedPE(int opch, size_t bin, std::vector & pe_vector){ + // Default implementation calls the non-channel-dependent version + return NObservedPE(bin, pe_vector); + } + }; #endif diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc new file mode 100644 index 000000000..02128fece --- /dev/null +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc @@ -0,0 +1,123 @@ +//////////////////////////////////////////////////////////////////////// +// Specific class tool for PMTNonLinearity +// File: PMTNonLinearityTF1_tool.hh +// Base class: PMTNonLinearity.hh +//////////////////////////////////////////////////////////////////////// + +#include "fhiclcpp/ParameterSet.h" +#include "art/Utilities/ToolMacros.h" +#include "art/Utilities/ToolConfigTable.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" + +#include "larcore/Geometry/WireReadout.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" + +#include +#include +#include +#include "TF1.h" + +#include "sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh" + + +namespace opdet { + class PMTNonLinearityTF1ChannelByChannel; +} + + +class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { +public: + + //Configuration parameters + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom attenuationForm { + Name("AttenuationForm"), + Comment("Non linearity functional form") + }; + + fhicl::Sequence> attenuationFormParams { + Name("AttenuationFormParams"), + Comment("Parameters for non linearity functional form") + }; + + fhicl::Atom attenuationPreTime { + Name("AttenuationPreTime"), + Comment("For the non-linear attenuation, consider photons arriving given by this time window. In nanoseconds.") + }; + + fhicl::Sequence nonLinearRange { + Name("NonLinearRange"), + Comment("Non linear range. Assume linear response/completely saturated response out of this range. In #PE.") + }; + + }; + + explicit PMTNonLinearityTF1ChannelByChannel(art::ToolConfigTable const& config); + + //Returns rescaled #pe after non linearity + double NObservedPE(int opch, size_t bin, std::vector & pe_vector) override; + +private: + //Configuration parameters + std::string fAttenuationForm; + std::vector> fAttenuationFormParams; + unsigned int fAttenuationPreTime; + std::vector fNonLinearRange; + + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); + + //TF1 for non linearity function + std::map fNonLinearTF1Map; + + // Vector to store non linearity attenuation values + std::vector> fPEAttenuation_V; + std::vector fPEAttenuation; + std::vector fPESaturationValue_V; + int fPESaturationValue; + +}; + + +opdet::PMTNonLinearityTF1ChannelByChannel::PMTNonLinearityTF1ChannelByChannel(art::ToolConfigTable const& config) + : fAttenuationForm { config().attenuationForm() } + , fAttenuationFormParams { config().attenuationFormParams() } + , fAttenuationPreTime { config().attenuationPreTime() } + , fNonLinearRange { config().nonLinearRange() } +{ + //Initialize TF1 for each channel with channel-dependent parameters + for(size_t opch=0; opchSetParameter(k, fAttenuationFormParams[opch][k]); + } + } + + // Initialize attenuation vector + fPEAttenuation_V.resize(fWireReadout.NOpChannels()); + for(size_t opch=0; opchEval(pe)/pe; + } + } + + fPESaturationValue_V.resize(fWireReadout.NOpChannels(), 1); + for(size_t opch=0; opchEval(fNonLinearRange[1])); + } +} + +double opdet::PMTNonLinearityTF1ChannelByChannel::NObservedPE(int opch, size_t bin, std::vector & pe_vector){ + size_t start_bin = bin-fAttenuationPreTime; + if(fAttenuationPreTime<0) start_bin=0; + unsigned int npe_acc = std::accumulate(pe_vector.begin()+start_bin,pe_vector.begin()+bin+1, 0); + + if(npe_acc Date: Mon, 10 Nov 2025 05:58:52 -0600 Subject: [PATCH 03/28] Add SER, gain fluct and nonlinearity per channel --- sbndcode/OpDetSim/DigiPMTSBNDAlg.cc | 209 ++++++++++++++++++++++------ sbndcode/OpDetSim/DigiPMTSBNDAlg.hh | 53 ++++++- 2 files changed, 216 insertions(+), 46 deletions(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index 53ec18074..e373fd79b 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -17,6 +17,8 @@ namespace opdet { , fPMTCoatedVISEff_tpc1(fParams.PMTCoatedVISEff_tpc1 / fParams.larProp->ScintPreScale()) , fPMTUncoatedEff_tpc1(fParams.PMTUncoatedEff_tpc1/ fParams.larProp->ScintPreScale()) // , fSinglePEmodel(fParams.SinglePEmodel) + , fUseDataNoise(fParams.UseDataNoise) + , fOpDetNoiseFile(fParams.OpDetNoiseFile) , fEngine(fParams.engine) , fFlatGen(*fEngine) , fPoissonQGen(*fEngine) @@ -51,59 +53,103 @@ namespace opdet { sp.find_file(fParams.PMTDataFile, fname); TFile* file = TFile::Open(fname.c_str(), "READ"); + std::vector>* SinglePEVec_p; + std::vector* fSinglePEChannels_p; + std::vector* fPeakAmplitude_p; + + file->GetObject("SERChannels", fSinglePEChannels_p); + file->GetObject("SinglePEVec", SinglePEVec_p); + file->GetObject("PeakAmplitude", fPeakAmplitude_p); + + fSinglePEWaveVector = *SinglePEVec_p; + fSinglePEChannels = *fSinglePEChannels_p; + fPeakAmplitude = *fPeakAmplitude_p; + // TPB emission time histogram for pmt_coated histogram std::vector* timeTPB_p; file->GetObject("timeTPB", timeTPB_p); fTimeTPB = std::make_unique (*fEngine, timeTPB_p->data(), timeTPB_p->size()); + // PMT calibration database service + fPMTCalibrationDatabaseService = lar::providerFrom(); + fPMTHDOpticalWaveformsPtr = art::make_tool(fParams.HDOpticalWaveformParams); + + //Resize the SER vector to the number of channels + fSinglePEWave_HD.resize(320); //shape of single pulse - if (fParams.PMTSinglePEmodel) { + if (fParams.PMTSinglePEmodel=="testbench") { mf::LogDebug("DigiPMTSBNDAlg") << " using testbench pe response"; std::vector* SinglePEVec_p; file->GetObject("SinglePEVec_HD", SinglePEVec_p); fSinglePEWave = *SinglePEVec_p; - // Prepare HD waveforms - fPMTHDOpticalWaveformsPtr = art::make_tool(fParams.HDOpticalWaveformParams); - fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD,fSinglePEWave); - - pulsesize = fSinglePEWave_HD[0].size(); + // Create the same SER for all channels + for(size_t i=0; iproduceSER_HD(fSinglePEWave_HD[i], fSinglePEWave); + } + pulsesize = fSinglePEWave_HD[0][0].size(); mf::LogDebug("DigiPMTSBNDAlg")<<"HD wvfs size: "<getReconstructChannel(i)) continue; + fSinglePEWave = fPMTCalibrationDatabaseService->getSER(i); + double SPEAmplitude = fPMTCalibrationDatabaseService->getSPEAmplitude(i); + double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);}); + double SinglePENormalization = std::abs(SPEAmplitude/SPEPeakValue); + std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;}); + if(fSinglePEWave.size()==0) continue; + if(fSinglePEWave.size()>0) pulsesize = fSinglePEWave.size(); + fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD[i], fSinglePEWave); + } + } + else{ + throw cet::exception("DigiPMTSBNDAlg") << "Wrong PMTSinglePEmodel configuration: " << fParams.PMTSinglePEmodel << std::endl; + } if(fParams.MakeGainFluctuations){ fPMTGainFluctuationsPtr = art::make_tool(fParams.GainFluctuationsParams); } if(fParams.SimulateNonLinearity){ - fPMTNonLinearityPtr = art::make_tool(fParams.NonLinearityParams); - } + fPMTNonLinearityPtr = art::make_tool(fParams.NonLinearityParams); + } // infer pulse polarity from SER peak sign double minADC_SinglePE = *min_element(fSinglePEWave.begin(), fSinglePEWave.end()); double maxADC_SinglePE = *max_element(fSinglePEWave.begin(), fSinglePEWave.end()); - fPositivePolarity = std::abs(maxADC_SinglePE) > std::abs(minADC_SinglePE); - + fPositivePolarity = std::abs(maxADC_SinglePE) > std::abs(minADC_SinglePE); + + // get ADC saturation value // currently assumes all dynamic range for PE (no overshoot) fADCSaturation = (fPositivePolarity ? fParams.PMTBaseline + fParams.PMTADCDynamicRange : fParams.PMTBaseline - fParams.PMTADCDynamicRange); + file->Close(); + // Initialize noise file + std::string fname_noise; + cet::search_path sp_noise("FW_SEARCH_PATH"); + if(fParams.UseDataNoise){ + std::cout << " Trying to open file " << fParams.OpDetNoiseFile << std::endl; + sp_noise.find_file(fParams.OpDetNoiseFile, fname_noise); + noise_file = TFile::Open(fname_noise.c_str(), "READ"); + } - file->Close(); } // end constructor - - DigiPMTSBNDAlg::~DigiPMTSBNDAlg(){} + DigiPMTSBNDAlg::~DigiPMTSBNDAlg(){ + if(fParams.UseDataNoise) noise_file->Close(); + } void DigiPMTSBNDAlg::ConstructWaveformUncoatedPMT( @@ -186,7 +232,7 @@ namespace opdet { //photon time in ns (w.r.t. the waveform start time a.k.a t_min) tphoton = ttsTime + simphotons[i].Time - t_min + ttpb + fParams.CableTime; - // store the pgoton time if it's within the readout window + // store the photon time if it's within the readout window if(tphoton > 0 && tphoton < nPE_v.size()) nPE_v[(size_t)tphoton]++; } } @@ -194,16 +240,22 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -265,17 +317,22 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } //Adding noise and saturation - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -320,16 +377,21 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -398,17 +460,22 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } //Adding noise and saturation - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -438,8 +505,9 @@ namespace opdet { } - void DigiPMTSBNDAlg::AddSPE(size_t time, std::vector& wave, double npe) + void DigiPMTSBNDAlg::AddSPE(size_t time, std::vector& wave, int ch, double npe) { + if(!fPMTCalibrationDatabaseService->getReconstructChannel(ch)) return; // time bin HD (double precision) // used to gert the time-shifted SER double time_bin_hd = fSampling*time; @@ -454,11 +522,11 @@ namespace opdet { // simulate gain fluctuations double npe_anode = npe; if(fParams.MakeGainFluctuations) - npe_anode=fPMTGainFluctuationsPtr->GainFluctuation(npe, fEngine); - + //npe_anode=fPMTGainFluctuationsPtr->GainFluctuation(npe, fEngine); + npe_anode=fPMTGainFluctuationsPtr->GainFluctuation(ch, npe, fEngine); // add SER to the waveform std::transform(min_it, max_it, - fSinglePEWave_HD[wvf_shift].begin(), min_it, + fSinglePEWave_HD[ch][wvf_shift].begin(), min_it, [npe_anode](auto a, auto b) { return a+npe_anode*b; }); } @@ -474,7 +542,7 @@ namespace opdet { } - void DigiPMTSBNDAlg::AddLineNoise(std::vector& wave) + void DigiPMTSBNDAlg::AddLineNoise(std::vector& wave, int ch) { // TODO: after running the profiler I can see that this is where // most cycles are being used. Potentially some improvement could @@ -495,7 +563,66 @@ namespace opdet { } - void DigiPMTSBNDAlg::AddDarkNoise(std::vector& wave) + void DigiPMTSBNDAlg::AddDataNoise(std::vector& wave, int ch) + { + // Get the noise waveform from data + // Get all the noise waveforms that correspond to the channel we are using + // Select a random noise waveform + // Compare the length of the noise and the simulation waveform: + // if the simulation length is larger then choose another random noise waveform to fill up the missing items + + if (!noise_file || noise_file->IsZombie()) { + throw cet::exception("DigiPMTSBNDAlg") << " PMT Noise file could not be opened " << std::endl; + return; + } + + TList *keys = noise_file->GetListOfKeys(); + if (!keys || keys->IsEmpty()) { + throw cet::exception("DigiPMTSBNDAlg") << " PMT Noise file is empty " << std::endl; + return; + } + + std::string opChannelName = "opchannel_" + std::to_string(ch) + "_pmt"; + std::vector keylist; + TIter nextkey(keys); + TKey *key; + while ((key = (TKey*)nextkey())) { + std::string channel_name = key->GetName(); + if (channel_name.find(opChannelName) != std::string::npos) { + keylist.push_back(key); + } + } + + if (keylist.empty()) { + throw cet::exception("DigiPMTSBNDAlg") << " PMT Noise file has no data for channel " << ch << std::endl; + return; + } + + std::vector noise_wform; + int waveBins = wave.size(); + int currentSize = 0; + + // Llenar el vector de ruido hasta igualar la longitud del waveform + while (currentSize < waveBins) { + int noiseWformIdx = static_cast(fEngine->flat() * keylist.size()); + TH1 *noiseHist = (TH1*)keylist[noiseWformIdx]->ReadObj(); + + for (int i = 1; i <= noiseHist->GetNbinsX(); i++) { + noise_wform.push_back(noiseHist->GetBinContent(i)); + currentSize += 1; + if (currentSize >= waveBins) break; + } + delete noiseHist; + } + + // Agregar ruido al waveform + for (size_t i = 0; i < wave.size(); i++) { + wave[i] += noise_wform[i] + fParams.PMTBaseline; + } + } + + + void DigiPMTSBNDAlg::AddDarkNoise(std::vector& wave, int ch) { double timeBin; // Multiply by 10^9 since fParams.DarkNoiseRate is in Hz (conversion from s to ns) @@ -503,7 +630,7 @@ namespace opdet { double darkNoiseTime = fExponentialGen.fire(mean); while(darkNoiseTime < wave.size()) { timeBin = std::round(darkNoiseTime); - if(timeBin < wave.size()) {AddSPE(fSamplingPeriod*timeBin, wave);} + if(timeBin < wave.size()) {AddSPE(fSamplingPeriod*timeBin, wave, ch);} // Find next time to add dark noise darkNoiseTime += fExponentialGen.fire(mean); } @@ -541,7 +668,6 @@ namespace opdet { return t_min; } - // TODO: this function is not being used anywhere! ~icaza double DigiPMTSBNDAlg::FindMinimumTimeLite( sim::SimPhotonsLite const& litesimphotons, @@ -596,6 +722,7 @@ namespace opdet { fBaseConfig.PMTCoatedVISEff_tpc1 = config.pmtcoatedVISEff_tpc1(); fBaseConfig.PMTUncoatedEff_tpc1 = config.pmtuncoatedEff_tpc1(); fBaseConfig.PMTSinglePEmodel = config.PMTsinglePEmodel(); + fBaseConfig.PositivePolarity = config.PositivePolarity(); fBaseConfig.PMTRiseTime = config.pmtriseTime(); fBaseConfig.PMTFallTime = config.pmtfallTime(); fBaseConfig.PMTMeanAmplitude = config.pmtmeanAmplitude(); @@ -605,6 +732,8 @@ namespace opdet { fBaseConfig.TTS = config.tts(); fBaseConfig.CableTime = config.cableTime(); fBaseConfig.PMTDataFile = config.pmtDataFile(); + fBaseConfig.UseDataNoise = config.UseDataNoise(); + fBaseConfig.OpDetNoiseFile = config.OpDetNoiseFile(); fBaseConfig.MakeGainFluctuations = config.gainFluctuationsParams.get_if_present(fBaseConfig.GainFluctuationsParams); fBaseConfig.SimulateNonLinearity = config.nonLinearityParams.get_if_present(fBaseConfig.NonLinearityParams); config.hdOpticalWaveformParams.get_if_present(fBaseConfig.HDOpticalWaveformParams); diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh index ed9190dd4..471d57fa0 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh @@ -41,8 +41,14 @@ #include "sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh" #include "sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh" #include "sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms.hh" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" + #include "TFile.h" +#include "TKey.h" +#include "TH1.h" +#include "TH1D.h" namespace opdet { @@ -69,10 +75,14 @@ namespace opdet { double PMTCoatedVISEff_tpc1; //PMT (coated) efficiency for reflected (VIS) light (TPC1) double PMTUncoatedEff_tpc1; //PMT (uncoated) efficiency (TPC1) std::string PMTDataFile; //File containing timing emission structure for TPB, and single PE profile from data - bool PMTSinglePEmodel; //Model for single pe response, false for ideal, true for test bench meas + std::string OpDetNoiseFile; + std::string PMTSinglePEmodel; //Model for single pe response, false for ideal, true for test bench meas bool MakeGainFluctuations; //Fluctuate PMT gain fhicl::ParameterSet GainFluctuationsParams; bool SimulateNonLinearity; //Fluctuate PMT gain + bool PositivePolarity; + bool OscillateAfterPulse; + bool UseDataNoise; fhicl::ParameterSet NonLinearityParams; fhicl::ParameterSet HDOpticalWaveformParams; @@ -124,6 +134,8 @@ namespace opdet { return fParams.PMTBaseline; } + void AddOscillationAfterPulse( std::vector& wave); + private: ConfigurationParameters_t fParams; @@ -138,10 +150,14 @@ namespace opdet { double fPMTUncoatedEff_tpc1; bool fPositivePolarity; int fADCSaturation; + bool fUseDataNoise; + std::string fOpDetNoiseFile; double sigma1; double sigma2; + TFile* noise_file; + const double transitTimeSpread_frac = 2.0 * std::sqrt(2.0 * std::log(2.0)); CLHEP::HepRandomEngine* fEngine; //!< Reference to art-managed random-number engine @@ -151,6 +167,9 @@ namespace opdet { CLHEP::RandExponential fExponentialGen; std::unique_ptr fTimeTPB; // histogram for getting the TPB emission time for coated PMTs + //PMTCalibrationDatabase service + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; + //PMTFluctuationsAlg std::unique_ptr fPMTGainFluctuationsPtr; //HDWaveforms @@ -159,14 +178,19 @@ namespace opdet { //PMTNonLinearity std::unique_ptr fPMTNonLinearityPtr; - void AddSPE(size_t time, std::vector& wave, double npe = 1); // add single pulse to auxiliary waveform + void AddSPE(size_t time, std::vector& wave, int ch ,double npe = 1); // add single pulse to auxiliary waveform void Pulse1PE(std::vector& wave); double Transittimespread(double fwhm); std::vector fSinglePEWave; // single photon pulse vector - std::vector> fSinglePEWave_HD; // single photon pulse vector + std::vector>> fSinglePEWave_HD; // single photon pulse vector int pulsesize; //size of 1PE waveform std::unordered_map< raw::Channel_t, std::vector > fFullWaveforms; + + + std::vector> fSinglePEWaveVector; + std::vector fSinglePEChannels; + std::vector fPeakAmplitude; void CreatePDWaveformUncoatedPMT( sim::SimPhotons const& SimPhotons, @@ -193,8 +217,9 @@ namespace opdet { std::unordered_map& DirectPhotonsMap, std::unordered_map& ReflectedPhotonsMap); void CreateSaturation(std::vector& wave);//Including saturation effects (dynamic range) - void AddLineNoise(std::vector& wave); //add noise to baseline - void AddDarkNoise(std::vector& wave); //add dark noise + void AddLineNoise(std::vector& wave, int ch); //add noise to baseline + void AddDataNoise(std::vector& wave, int ch); //add noise from data to baseline + void AddDarkNoise(std::vector& wave, int ch); //add dark noise double FindMinimumTime( sim::SimPhotons const&, int ch, @@ -299,7 +324,7 @@ namespace opdet { Comment("PMT (uncoated) detection efficiency (TPC1)") }; - fhicl::Atom PMTsinglePEmodel { + fhicl::Atom PMTsinglePEmodel { Name("PMTSinglePEmodel"), Comment("Model used for single PE response of PMT. =0 is ideal, =1 is testbench") }; @@ -309,6 +334,21 @@ namespace opdet { Comment("File containing timing emission distribution for TPB and single pe pulse from data") }; + fhicl::Atom PositivePolarity { + Name("PositivePolarity"), + Comment("Polarity of the waveforms") + }; + + fhicl::Atom UseDataNoise { + Name("UseDataNoise"), + Comment("Use data noise from file") + }; + + fhicl::Atom OpDetNoiseFile { + Name("OpDetNoiseFile"), + Comment("Use data noise from this file") + }; + fhicl::OptionalDelegatedParameter gainFluctuationsParams { Name("GainFluctuationsParams"), Comment("Parameters used for SinglePE response fluctuations") @@ -324,6 +364,7 @@ namespace opdet { Comment("Parameters used for high definition waveform") }; + }; //struct Config DigiPMTSBNDAlgMaker(Config const& config); //Constructor From f3d8beeaf4dc892f458eaabe90173981b339b7fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 11 Nov 2025 04:38:28 -0600 Subject: [PATCH 04/28] Fix TF1 memory leak --- .../PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc index 02128fece..6b7f8569f 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc @@ -67,11 +67,10 @@ class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { std::vector> fAttenuationFormParams; unsigned int fAttenuationPreTime; std::vector fNonLinearRange; - geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); //TF1 for non linearity function - std::map fNonLinearTF1Map; + std::map> fNonLinearTF1Map; // Vector to store non linearity attenuation values std::vector> fPEAttenuation_V; @@ -90,7 +89,7 @@ opdet::PMTNonLinearityTF1ChannelByChannel::PMTNonLinearityTF1ChannelByChannel(ar { //Initialize TF1 for each channel with channel-dependent parameters for(size_t opch=0; opch(("NonLinearTF1_"+std::to_string(opch)).c_str(), fAttenuationForm.c_str()); for(size_t k=0; kSetParameter(k, fAttenuationFormParams[opch][k]); } From e606f655cbd4b01dc52d9499f5b5c39e6f1eb5d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Thu, 13 Nov 2025 05:04:49 -0600 Subject: [PATCH 05/28] Add intialization method for using calibration database --- ...PMTNonLinearityTF1ChannelByChannel_tool.cc | 26 +++++++++++-------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc index 6b7f8569f..9c5d23039 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc @@ -19,7 +19,8 @@ #include "TF1.h" #include "sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh" - +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" namespace opdet { class PMTNonLinearityTF1ChannelByChannel; @@ -38,11 +39,6 @@ class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { Name("AttenuationForm"), Comment("Non linearity functional form") }; - - fhicl::Sequence> attenuationFormParams { - Name("AttenuationFormParams"), - Comment("Parameters for non linearity functional form") - }; fhicl::Atom attenuationPreTime { Name("AttenuationPreTime"), @@ -62,9 +58,10 @@ class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { double NObservedPE(int opch, size_t bin, std::vector & pe_vector) override; private: + + void ConfigureChannelByChannelTF1s(); //Configuration parameters std::string fAttenuationForm; - std::vector> fAttenuationFormParams; unsigned int fAttenuationPreTime; std::vector fNonLinearRange; geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); @@ -78,21 +75,27 @@ class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { std::vector fPESaturationValue_V; int fPESaturationValue; + //PMTCalibrationDatabase service + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; }; opdet::PMTNonLinearityTF1ChannelByChannel::PMTNonLinearityTF1ChannelByChannel(art::ToolConfigTable const& config) : fAttenuationForm { config().attenuationForm() } - , fAttenuationFormParams { config().attenuationFormParams() } , fAttenuationPreTime { config().attenuationPreTime() } , fNonLinearRange { config().nonLinearRange() } { + fPMTCalibrationDatabaseService = lar::providerFrom(); +} + + +void opdet::PMTNonLinearityTF1ChannelByChannel::ConfigureChannelByChannelTF1s(){ + //Initialize TF1 for each channel with channel-dependent parameters for(size_t opch=0; opch(("NonLinearTF1_"+std::to_string(opch)).c_str(), fAttenuationForm.c_str()); - for(size_t k=0; kSetParameter(k, fAttenuationFormParams[opch][k]); - } + fNonLinearTF1Map[opch]->SetParameter(0, fPMTCalibrationDatabaseService->getNonLineatiryPESat(opch)); + fNonLinearTF1Map[opch]->SetParameter(1, fPMTCalibrationDatabaseService->getNonLineatiryAlpha(opch)); } // Initialize attenuation vector @@ -110,6 +113,7 @@ opdet::PMTNonLinearityTF1ChannelByChannel::PMTNonLinearityTF1ChannelByChannel(ar } } + double opdet::PMTNonLinearityTF1ChannelByChannel::NObservedPE(int opch, size_t bin, std::vector & pe_vector){ size_t start_bin = bin-fAttenuationPreTime; if(fAttenuationPreTime<0) start_bin=0; From 211159344785aadbad0d89debe2f4c2bc6d4ec61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Thu, 13 Nov 2025 05:08:36 -0600 Subject: [PATCH 06/28] Add configuration function for old version --- sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1_tool.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1_tool.cc index c93b1762f..069e01df6 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1_tool.cc +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1_tool.cc @@ -79,7 +79,9 @@ opdet::PMTNonLinearityTF1::PMTNonLinearityTF1(art::ToolConfigTable const , fAttenuationFormParams { config().attenuationFormParams() } , fAttenuationPreTime { config().attenuationPreTime() } , fNonLinearRange { config().nonLinearRange() } -{ +{} + +void opdet::PMTNonLinearityTF1::ConfigureNonLinearity(){ fNonLinearTF1 = new TF1("NonLinearTF1", fAttenuationForm.c_str()); for(size_t k=0; kSetParameter(k, fAttenuationFormParams[k]); From 67cabffa6e8655f8bf342e1d2a9b64e238c07c71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Thu, 13 Nov 2025 05:09:05 -0600 Subject: [PATCH 07/28] Call configuration function before using the tool --- sbndcode/OpDetSim/DigiPMTSBNDAlg.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index e373fd79b..e90bb1dd4 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -123,6 +123,7 @@ namespace opdet { } if(fParams.SimulateNonLinearity){ fPMTNonLinearityPtr = art::make_tool(fParams.NonLinearityParams); + fPMTNonLinearityPtr->ConfigureNonLinearity(); } // infer pulse polarity from SER peak sign From b5ada101665d971646fd7212b9879b23ae4a235f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Thu, 13 Nov 2025 05:20:39 -0600 Subject: [PATCH 08/28] Declare configure method on base class --- sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh | 3 +++ .../PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc | 7 +++---- sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1_tool.cc | 5 +++-- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh index 762c8c2a2..5b00f3b56 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh @@ -18,6 +18,9 @@ class opdet::PMTNonLinearity { public: //Constructor virtual ~PMTNonLinearity() noexcept = default; + + //Configure non-linearity tool + virtual void ConfigureNonLinearity() = 0; //Returns rescaled number of PE virtual double NObservedPE(size_t bin, std::vector & pe_vector){ diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc index 9c5d23039..c7e217b91 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc @@ -56,10 +56,9 @@ class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { //Returns rescaled #pe after non linearity double NObservedPE(int opch, size_t bin, std::vector & pe_vector) override; - + void ConfigureNonLinearity() override; + private: - - void ConfigureChannelByChannelTF1s(); //Configuration parameters std::string fAttenuationForm; unsigned int fAttenuationPreTime; @@ -89,7 +88,7 @@ opdet::PMTNonLinearityTF1ChannelByChannel::PMTNonLinearityTF1ChannelByChannel(ar } -void opdet::PMTNonLinearityTF1ChannelByChannel::ConfigureChannelByChannelTF1s(){ +void opdet::PMTNonLinearityTF1ChannelByChannel::ConfigureNonLinearity(){ //Initialize TF1 for each channel with channel-dependent parameters for(size_t opch=0; opch & pe_vector) override; + void ConfigureNonLinearity() override; private: //Configuration parameters @@ -66,7 +67,7 @@ class opdet::PMTNonLinearityTF1 : opdet::PMTNonLinearity { std::vector fNonLinearRange; //TF1 for non linearity function - TF1 *fNonLinearTF1; + std::unique_ptr fNonLinearTF1; // Vector to store non linearity attenuation values std::vector fPEAttenuation_V; @@ -82,7 +83,7 @@ opdet::PMTNonLinearityTF1::PMTNonLinearityTF1(art::ToolConfigTable const {} void opdet::PMTNonLinearityTF1::ConfigureNonLinearity(){ - fNonLinearTF1 = new TF1("NonLinearTF1", fAttenuationForm.c_str()); + fNonLinearTF1 = std::make_unique("NonLinearTF1", fAttenuationForm.c_str()); for(size_t k=0; kSetParameter(k, fAttenuationFormParams[k]); } From 8e0354caac04ac7fed488b7f71aaac8657dcbeb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 18 Nov 2025 04:11:25 -0600 Subject: [PATCH 09/28] Remove non linearity parameters from fcl --- sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl | 1 - 1 file changed, 1 deletion(-) diff --git a/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl b/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl index ec9268bdc..344293a24 100644 --- a/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl +++ b/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl @@ -18,7 +18,6 @@ PMTNonLinearityTF1ChannelByChannel: #AttenuationForm: "sqrt( (sqrt(1+8*(x/[0])**[1])-1)/( 4*(x/[0])**[1] ) )" #AttenuationFormParams: [78.15, 1.96] AttenuationForm: "x/sqrt(1+(x/[0])**[1])" - AttenuationFormParams: [[269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.1], [280.0, 2.2], [160.0, 1.5], [280.0, 2.1], [220.0, 2.0], [250.0, 1.9], [280.0, 2.2], [160.0, 1.5], [190.0, 1.5], [190.0, 1.9], [250.0, 1.8], [280.0, 1.5], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [370.0, 2.2], [280.0, 2.2], [220.0, 1.8], [190.0, 1.8], [250.0, 1.5], [250.0, 2.1], [269, 1.84], [269, 1.84], [400.0, 2.2], [220.0, 1.7], [280.0, 1.9], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [280.0, 2.2], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.0], [250.0, 1.9], [280.0, 2.2], [250.0, 2.0], [269, 1.84], [280.0, 2.2], [190.0, 1.8], [250.0, 2.0], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.2], [220.0, 1.6], [269, 1.84], [130.0, 1.5], [190.0, 1.8], [250.0, 2.1], [280.0, 2.2], [220.0, 2.0], [160.0, 1.5], [220.0, 2.0], [310.0, 2.0], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.1], [190.0, 1.9], [190.0, 1.8], [160.0, 1.7], [190.0, 1.6], [250.0, 1.9], [280.0, 2.2], [310.0, 2.2], [269, 1.84], [280.0, 2.2], [250.0, 2.0], [160.0, 2.2], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.0], [269, 1.84], [269, 1.84], [250.0, 2.0], [190.0, 1.9], [269, 1.84], [269, 1.84], [269, 1.84], [280.0, 2.1], [310.0, 2.2], [269, 1.84], [370.0, 2.2], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.0], [220.0, 2.0], [190.0, 1.5], [310.0, 2.1], [220.0, 2.1], [269, 1.84], [220.0, 2.1], [310.0, 2.2], [269, 1.84], [269, 1.84], [160.0, 1.5], [280.0, 2.2], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [220.0, 2.0], [310.0, 2.2], [220.0, 1.8], [190.0, 2.0], [370.0, 2.2], [190.0, 1.9], [250.0, 2.1], [130.0, 1.5], [269, 1.84], [220.0, 1.8], [160.0, 1.6], [160.0, 1.7], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84], [269, 1.84]] AttenuationPreTime: 4 #ns NonLinearRange: [100, 701] } From 05846ef29050fb33752ddc95a9b183ea153cfdc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 18 Nov 2025 04:12:51 -0600 Subject: [PATCH 10/28] Update database configuration --- .../PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl b/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl index c2462179f..6b2a9d13a 100644 --- a/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl +++ b/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl @@ -3,8 +3,8 @@ BEGIN_PROLOG sbnd_pmtcalibrationdatabaseservice : { service_provider: "PMTCalibrationDatabaseService" CorrectionTags: { - PMTCalibrationDatabaseTag: "v1r1" - DatabaseTimeStamp: 1757601071000000000 + PMTCalibrationDatabaseTag: "v1r2" + DatabaseTimeStamp: 1763157679000000000 TableName: "pds_calibration" SERLength: 550 } From fd7926de2c019805fe2ced9b22042e43f3cafb32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 05:25:30 -0600 Subject: [PATCH 11/28] Add PMT pulse oscillation module --- .../PMTPulseOscillations/CMakeLists.txt | 28 +++ .../PMTPulseOscillation_module.cc | 161 ++++++++++++++++++ .../PMTPulseOscillations/job/CMakeLists.txt | 1 + .../job/run_pmtpulseoscillation.fcl | 60 +++++++ .../job/sbnd_pmtpulseoscillation_config.fcl | 19 +++ 5 files changed, 269 insertions(+) create mode 100644 sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt create mode 100644 sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc create mode 100644 sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt create mode 100644 sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl create mode 100644 sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt b/sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt new file mode 100644 index 000000000..2eb94cd3f --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt @@ -0,0 +1,28 @@ +set( + MODULE_LIBRARIES + sbndcode_OpDetSim + larcore::Geometry_Geometry_service + lardataobj::Simulation + lardata::Utilities + lardataobj::RawData + sbndcode_Utilities_SignalShapingServiceSBND_service + nurandom::RandomUtils_NuRandomService_service + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Optional_RandomNumberGenerator_service + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib + cetlib_except::cetlib_except + CLHEP::CLHEP + ${ROOT_BASIC_LIB_LIST} +) + +cet_build_plugin(PMTPulseOscillation art::module SOURCE PMTPulseOscillation_module.cc LIBRARIES ${MODULE_LIBRARIES}) + +install_headers() +install_fhicl() +install_source() +add_subdirectory(job) +##FILE(GLOB fcl_files *.fcl) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc new file mode 100644 index 000000000..614b709a9 --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc @@ -0,0 +1,161 @@ +//////////////////////////////////////////////////////////////////////// +// Class: SBNDOpDeconvolution +// Plugin Type: producer (art v3_06_03) +// File: SBNDOpDeconvolution_module.cc +// +// Generated at Tue Jul 13 06:29:02 2021 by Francisco Nicolas-Arnaldos using cetskelgen +// from cetlib version v3_11_01. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art/Utilities/make_tool.h" +#include "art_root_io/TFileService.h" + +#include "CLHEP/Random/RandFlat.h" + +#include + +#include "lardataobj/RawData/OpDetWaveform.h" + +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" +#include "TF1.h" +#include "TH1F.h" + +namespace opdet { + class PMTPulseOscillation; +} + + +class opdet::PMTPulseOscillation : public art::EDProducer { + public: + explicit PMTPulseOscillation(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + PMTPulseOscillation(PMTPulseOscillation const&) = delete; + PMTPulseOscillation(PMTPulseOscillation&&) = delete; + PMTPulseOscillation& operator=(PMTPulseOscillation const&) = delete; + PMTPulseOscillation& operator=(PMTPulseOscillation&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + void CreateOscillatedWaveform(const raw::OpDetWaveform& oldWaveform , raw::OpDetWaveform& newWaveform ); + + + private: + + // Declare member data here. + std::string fInputLabel; + std::vector fPDTypes; + std::vector fElectronics; + std::string fOscillationFunction; + double fOscillationFrequency; + double fOscillationDampingConstant; + double fSamplingPeriod; + double fPulseOscillationThreshold; + int fOscillationOffset; + //PDS map + opdet::sbndPDMapAlg pdsmap; + std::vector fOscillationFunctionParams; + TF1 *fOscillationTF1; + + + CLHEP::HepRandomEngine* engine = nullptr; //!< Reference to art-managed random-number engine + CLHEP::RandFlat fFlatGen; + + +}; + + +opdet::PMTPulseOscillation::PMTPulseOscillation(fhicl::ParameterSet const& p) + : EDProducer{p}, + fFlatGen(engine) // , + // More initializers here. +{ + std::cout << " DB 1 " << std::endl; + // Call appropriate produces<>() functions here. + // Call appropriate consumes<>() for any products to be retrieved by this module. + fInputLabel = p.get< std::string >("InputLabel"); + fPDTypes = p.get< std::vector >("PDTypes"); + fElectronics = p.get< std::vector >("Electronics"); + fOscillationFunction = p.get("OscillationFunction"); + fSamplingPeriod = p.get("SamplingPeriod"); + fOscillationFrequency = p.get("OscillationFrequency"); + fOscillationDampingConstant = p.get("OscillationDampingConstant"); + fOscillationFunctionParams = p.get>("OscillationFunctionParams"); + fPulseOscillationThreshold = p.get("PulseOscillationThreshold"); + fOscillationOffset = p.get("OscillationOffset"); + fOscillationTF1 = new TF1("OscillationTemplate", fOscillationFunction.c_str()); + for(size_t i=0; iSetParameter(i, fOscillationFunctionParams[i]); + } + + produces< std::vector< raw::OpDetWaveform > >(); +} + +void opdet::PMTPulseOscillation::produce(art::Event& e) +{ + //Load the waveforms + art::Handle< std::vector< raw::OpDetWaveform > > wfHandle; + e.getByLabel(fInputLabel, wfHandle); + if (!wfHandle.isValid()) { + std::cout << " not finding the handles " << std::endl; + mf::LogError("PMTPulseOscillation")<<"Input waveforms with input label "< > OscillatedWf_VecPtr(std::make_unique< std::vector< raw::OpDetWaveform > > ()); + for(auto const& wf : *wfHandle){ + if(std::find(fPDTypes.begin(), fPDTypes.end(), pdsmap.pdType(wf.ChannelNumber()) ) != fPDTypes.end() ){ + raw::OpDetWaveform new_wf = wf; + CreateOscillatedWaveform(wf, new_wf); + OscillatedWf_VecPtr->push_back(new_wf); + } + else OscillatedWf_VecPtr->push_back(wf); + } + + e.put( std::move(OscillatedWf_VecPtr) ); +} + + +void opdet::PMTPulseOscillation::CreateOscillatedWaveform(const raw::OpDetWaveform& oldWaveform, raw::OpDetWaveform& newWaveform) +{ + std::cout << " DB 2 " << std::endl; + //Find the idx of the maximum (or minimum) point of the waveform + int starting_ttick = std::distance(oldWaveform.begin(), std::min_element(oldWaveform.begin(), oldWaveform.end())); + double baseline = accumulate( oldWaveform.begin(), oldWaveform.end(), 0.0) / oldWaveform.size(); + double pulse_max = abs(oldWaveform[starting_ttick]-baseline); + if(pulse_maxEval(pulse_max); + std::cout << " DB 3 " << std::endl; + + for(size_t i=starting_ttick; i(i)-static_cast(starting_ttick))*fSamplingPeriod; + double envelope = std::exp(-delta_time / fOscillationDampingConstant); + std::cout << " DB 5 " << std::endl; + double phase = 2.0 * M_PI * fOscillationFrequency * delta_time; + std::cout << " DB 6 " << std::endl; + double new_value = static_cast(newWaveform.Waveform()[i+fOscillationOffset]) + oscillation_amplitude * envelope * std::cos(phase+M_PI); + std::cout << " DB 7 " << std::endl; + // Add some dither to avoid quantization effects + std::cout << " throwing random number " << fFlatGen.fire(1.0) << std::endl; + std::cout << " DB 8 " << std::endl; + double dither = fFlatGen.fire(1.0) - 0.5; + std::cout << " DB 9 " << std::endl; + newWaveform.Waveform()[i+fOscillationOffset] = std::round(new_value + dither); + std::cout << " DB 10 " << std::endl; + if(envelope < 0.05) return; + } + std::cout << " DB 4 " << std::endl; +} + +DEFINE_ART_MODULE(opdet::PMTPulseOscillation) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt b/sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt new file mode 100644 index 000000000..eca2a52d6 --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt @@ -0,0 +1 @@ +install_fhicl() \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl b/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl new file mode 100644 index 000000000..4ba9d6fde --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl @@ -0,0 +1,60 @@ +#include "services_sbnd.fcl" +#include "messages_sbnd.fcl" +#include "sam_sbnd.fcl" +#include "sbnd_pmtpulseoscillation_config.fcl" + + +process_name: pmtpulseoscillation + +services: +{ + # Load the service that manages root files for histograms. + TFileService: { fileName: "hist.root" } + RandomNumberGenerator: {} # required by fuzzyCluster + message: @local::sbnd_message_services_prod # from messages_sbnd.fcl +} + +#source is now a root file +source: +{ + module_type: RootInput + maxEvents: -1 # Number of events to create +} + +outputs: +{ + out1: + { + module_type: RootOutput + fileName: "%ifb_%p.root" + dataTier: "reconstructed" + compressionLevel: 1 + } +} + +physics: +{ + + producers: + { + # random number saver + rns: { module_type: RandomNumberSaver } + pmtpulseoscillation: @local::sbnd_pmtpulseoscillation + } + + + # define the producer and filter modules for this path, order matters, + # filters reject all following items. see lines starting physics.producers below + reco: [rns, pmtpulseoscillation ] + + stream1: [out1] + + # trigger_paths is a keyword and contains the paths that modify the art::event, + # ie filters and producers + trigger_paths: [reco] + + + # end_paths is a keyword and contains the paths that do not modify the art::Event, + # ie analyzers and output streams. these all run simultaneously + end_paths: [stream1] +} \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl b/sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl new file mode 100644 index 000000000..143a9762b --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl @@ -0,0 +1,19 @@ +BEGIN_PROLOG + +sbnd_pmtpulseoscillation: +{ + # Parameters for ideal SER simulation + module_type: "PMTPulseOscillation" + InputLabel: "opdaq" + PDTypes: ["pmt_coated", "pmt_uncoated"] + Electronics: [] + OscillationFrequency: 0.000178571 + OscillationDampingConstant: 50000 + SamplingPeriod: 2 + OscillationFunction: "[0]*exp([1] * x)" + OscillationFunctionParams: [0.20292146, 0.00079746] + PulseOscillationThreshold: 70 + OscillationOffset: 0 +} + +END_PROLOG \ No newline at end of file From bbadc030b0a549807103026c316a0dbdd0b9eb15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 05:25:57 -0600 Subject: [PATCH 12/28] Modify MC workflow for realistic PDS --- .../standard/reco/config/workflow_reco1.fcl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl index fbc2ea6d7..9a375942d 100644 --- a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl +++ b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl @@ -2,9 +2,9 @@ # #include "caldata_sbnd.fcl" #include "hitfindermodules_sbnd.fcl" -#include "opdeconvolution_sbnd.fcl" -#include "sbnd_ophitfinder_deco.fcl" -#include "sbnd_flashfinder_deco.fcl" +#include "opdeconvolution_sbnd_data.fcl" +#include "sbnd_ophitfinder_deco_data.fcl" +#include "sbnd_flashfinder_deco_data.fcl" #include "crtrecoproducers_sbnd.fcl" #include "mlreco_sbnd.fcl" #include "cluster_sbnd.fcl" From ef093364840fd5749e621569df48873f228a6ebc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 06:13:52 -0600 Subject: [PATCH 13/28] Fix random number generator --- .../PMTPulseOscillation_module.cc | 27 +++++-------------- .../job/run_pmtpulseoscillation.fcl | 4 ++- 2 files changed, 9 insertions(+), 22 deletions(-) diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc index 614b709a9..faaf754d5 100644 --- a/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc +++ b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc @@ -20,6 +20,7 @@ #include "art_root_io/TFileService.h" #include "CLHEP/Random/RandFlat.h" +#include "nurandom/RandomUtils/NuRandomService.h" #include @@ -67,21 +68,16 @@ class opdet::PMTPulseOscillation : public art::EDProducer { opdet::sbndPDMapAlg pdsmap; std::vector fOscillationFunctionParams; TF1 *fOscillationTF1; - - - CLHEP::HepRandomEngine* engine = nullptr; //!< Reference to art-managed random-number engine - CLHEP::RandFlat fFlatGen; - - + CLHEP::HepRandomEngine& fEngine; }; opdet::PMTPulseOscillation::PMTPulseOscillation(fhicl::ParameterSet const& p) : EDProducer{p}, - fFlatGen(engine) // , + fEngine(art::ServiceHandle{}->registerAndSeedEngine( + createEngine(0, "HepJamesRandom", "pulseosc"), "HepJamesRandom", "pulseosc", p, "Seed")) // More initializers here. { - std::cout << " DB 1 " << std::endl; // Call appropriate produces<>() functions here. // Call appropriate consumes<>() for any products to be retrieved by this module. fInputLabel = p.get< std::string >("InputLabel"); @@ -128,34 +124,23 @@ void opdet::PMTPulseOscillation::produce(art::Event& e) void opdet::PMTPulseOscillation::CreateOscillatedWaveform(const raw::OpDetWaveform& oldWaveform, raw::OpDetWaveform& newWaveform) { - std::cout << " DB 2 " << std::endl; //Find the idx of the maximum (or minimum) point of the waveform int starting_ttick = std::distance(oldWaveform.begin(), std::min_element(oldWaveform.begin(), oldWaveform.end())); double baseline = accumulate( oldWaveform.begin(), oldWaveform.end(), 0.0) / oldWaveform.size(); double pulse_max = abs(oldWaveform[starting_ttick]-baseline); if(pulse_maxEval(pulse_max); - std::cout << " DB 3 " << std::endl; - for(size_t i=starting_ttick; i(i)-static_cast(starting_ttick))*fSamplingPeriod; double envelope = std::exp(-delta_time / fOscillationDampingConstant); - std::cout << " DB 5 " << std::endl; double phase = 2.0 * M_PI * fOscillationFrequency * delta_time; - std::cout << " DB 6 " << std::endl; double new_value = static_cast(newWaveform.Waveform()[i+fOscillationOffset]) + oscillation_amplitude * envelope * std::cos(phase+M_PI); - std::cout << " DB 7 " << std::endl; // Add some dither to avoid quantization effects - std::cout << " throwing random number " << fFlatGen.fire(1.0) << std::endl; - std::cout << " DB 8 " << std::endl; - double dither = fFlatGen.fire(1.0) - 0.5; - std::cout << " DB 9 " << std::endl; + std::cout << CLHEP::RandFlat::shoot(&fEngine,1.0) << std::endl; + double dither = CLHEP::RandFlat::shoot(&fEngine,1.0) - 0.5; newWaveform.Waveform()[i+fOscillationOffset] = std::round(new_value + dither); - std::cout << " DB 10 " << std::endl; if(envelope < 0.05) return; } - std::cout << " DB 4 " << std::endl; } DEFINE_ART_MODULE(opdet::PMTPulseOscillation) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl b/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl index 4ba9d6fde..4e727f6dd 100644 --- a/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl @@ -1,4 +1,5 @@ #include "services_sbnd.fcl" +#include "simulationservices_sbnd.fcl" #include "messages_sbnd.fcl" #include "sam_sbnd.fcl" #include "sbnd_pmtpulseoscillation_config.fcl" @@ -12,6 +13,7 @@ services: TFileService: { fileName: "hist.root" } RandomNumberGenerator: {} # required by fuzzyCluster message: @local::sbnd_message_services_prod # from messages_sbnd.fcl + NuRandomService: @local::sbnd_seedservice } #source is now a root file @@ -38,7 +40,7 @@ physics: producers: { # random number saver - rns: { module_type: RandomNumberSaver } + rns: { module_type: "RandomNumberSaver" } pmtpulseoscillation: @local::sbnd_pmtpulseoscillation } From 3f7d0abf4b7560581fe957a04710ef5ac0a4ae05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 08:21:12 -0600 Subject: [PATCH 14/28] Revert "Update database configuration" This reverts commit 05846ef29050fb33752ddc95a9b183ea153cfdc3. --- .../PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl b/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl index 6b2a9d13a..c2462179f 100644 --- a/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl +++ b/sbndcode/Calibration/PDSDatabaseInterface/pmtcalibrationdatabase_sbnd.fcl @@ -3,8 +3,8 @@ BEGIN_PROLOG sbnd_pmtcalibrationdatabaseservice : { service_provider: "PMTCalibrationDatabaseService" CorrectionTags: { - PMTCalibrationDatabaseTag: "v1r2" - DatabaseTimeStamp: 1763157679000000000 + PMTCalibrationDatabaseTag: "v1r1" + DatabaseTimeStamp: 1757601071000000000 TableName: "pds_calibration" SERLength: 550 } From 943695fa104675cc16cbfdabb7f2d71d16c1136a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 09:01:37 -0600 Subject: [PATCH 15/28] Modify HDOpticalWaveform configuration for data time sampling --- sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl b/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl index 7ee6220b5..5a5e2a23f 100644 --- a/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl +++ b/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl @@ -11,7 +11,7 @@ IncludeHDOpticalWaveforms_PMT: { tool_type: "HDOpticalWaveforms" - ResolutionIncrease: 2. + ResolutionIncrease: 1. } END_PROLOG From 21d819146b81c359cea547985b3bb7660dc7991f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 09:14:37 -0600 Subject: [PATCH 16/28] Remove polarity as fcl-configurable parameter --- sbndcode/OpDetSim/DigiPMTSBNDAlg.hh | 5 ----- 1 file changed, 5 deletions(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh index 471d57fa0..bd1638914 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh @@ -334,11 +334,6 @@ namespace opdet { Comment("File containing timing emission distribution for TPB and single pe pulse from data") }; - fhicl::Atom PositivePolarity { - Name("PositivePolarity"), - Comment("Polarity of the waveforms") - }; - fhicl::Atom UseDataNoise { Name("UseDataNoise"), Comment("Use data noise from file") From d326923a2c546e581f3fcdf215b817deef3c1dc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 09:15:49 -0600 Subject: [PATCH 17/28] Remove polarity as fcl-configurable parameter --- sbndcode/OpDetSim/DigiPMTSBNDAlg.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index e90bb1dd4..ed93663ae 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -723,7 +723,6 @@ namespace opdet { fBaseConfig.PMTCoatedVISEff_tpc1 = config.pmtcoatedVISEff_tpc1(); fBaseConfig.PMTUncoatedEff_tpc1 = config.pmtuncoatedEff_tpc1(); fBaseConfig.PMTSinglePEmodel = config.PMTsinglePEmodel(); - fBaseConfig.PositivePolarity = config.PositivePolarity(); fBaseConfig.PMTRiseTime = config.pmtriseTime(); fBaseConfig.PMTFallTime = config.pmtfallTime(); fBaseConfig.PMTMeanAmplitude = config.pmtmeanAmplitude(); From 8c3b39184a7c333605d7f5aa2ffe64f08fa53f32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 09:45:41 -0600 Subject: [PATCH 18/28] Fix noise baseline and cleanup --- sbndcode/OpDetSim/DigiPMTSBNDAlg.cc | 26 +++++--------------------- sbndcode/OpDetSim/DigiPMTSBNDAlg.hh | 8 +++----- 2 files changed, 8 insertions(+), 26 deletions(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index ed93663ae..743172a20 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -47,36 +47,24 @@ namespace opdet { fSampling = fSampling / 1000.0; //in GHz, to cancel with ns fSamplingPeriod = 1./fSampling; - std::string fname; cet::search_path sp("FW_SEARCH_PATH"); sp.find_file(fParams.PMTDataFile, fname); TFile* file = TFile::Open(fname.c_str(), "READ"); - std::vector>* SinglePEVec_p; - std::vector* fSinglePEChannels_p; - std::vector* fPeakAmplitude_p; - - file->GetObject("SERChannels", fSinglePEChannels_p); - file->GetObject("SinglePEVec", SinglePEVec_p); - file->GetObject("PeakAmplitude", fPeakAmplitude_p); - - fSinglePEWaveVector = *SinglePEVec_p; - fSinglePEChannels = *fSinglePEChannels_p; - fPeakAmplitude = *fPeakAmplitude_p; - // TPB emission time histogram for pmt_coated histogram std::vector* timeTPB_p; file->GetObject("timeTPB", timeTPB_p); fTimeTPB = std::make_unique (*fEngine, timeTPB_p->data(), timeTPB_p->size()); - - // PMT calibration database service + + // PMT calibration database service fPMTCalibrationDatabaseService = lar::providerFrom(); fPMTHDOpticalWaveformsPtr = art::make_tool(fParams.HDOpticalWaveformParams); + int NOpChannels = fWireReadout.NOpChannels(); //Resize the SER vector to the number of channels - fSinglePEWave_HD.resize(320); + fSinglePEWave_HD.resize(NOpChannels); //shape of single pulse if (fParams.PMTSinglePEmodel=="testbench") { mf::LogDebug("DigiPMTSBNDAlg") << " using testbench pe response"; @@ -125,23 +113,19 @@ namespace opdet { fPMTNonLinearityPtr = art::make_tool(fParams.NonLinearityParams); fPMTNonLinearityPtr->ConfigureNonLinearity(); } - // infer pulse polarity from SER peak sign double minADC_SinglePE = *min_element(fSinglePEWave.begin(), fSinglePEWave.end()); double maxADC_SinglePE = *max_element(fSinglePEWave.begin(), fSinglePEWave.end()); fPositivePolarity = std::abs(maxADC_SinglePE) > std::abs(minADC_SinglePE); - // get ADC saturation value // currently assumes all dynamic range for PE (no overshoot) fADCSaturation = (fPositivePolarity ? fParams.PMTBaseline + fParams.PMTADCDynamicRange : fParams.PMTBaseline - fParams.PMTADCDynamicRange); file->Close(); - // Initialize noise file std::string fname_noise; cet::search_path sp_noise("FW_SEARCH_PATH"); if(fParams.UseDataNoise){ - std::cout << " Trying to open file " << fParams.OpDetNoiseFile << std::endl; sp_noise.find_file(fParams.OpDetNoiseFile, fname_noise); noise_file = TFile::Open(fname_noise.c_str(), "READ"); } @@ -618,7 +602,7 @@ namespace opdet { // Agregar ruido al waveform for (size_t i = 0; i < wave.size(); i++) { - wave[i] += noise_wform[i] + fParams.PMTBaseline; + wave[i] += noise_wform[i]; } } diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh index bd1638914..712d65eb2 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh @@ -37,6 +37,7 @@ #include "lardata/DetectorInfoServices/DetectorClocksServiceStandard.h" #include "lardataobj/Simulation/SimPhotons.h" #include "lardata/DetectorInfoServices/LArPropertiesService.h" +#include "larcore/Geometry/WireReadout.h" #include "sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh" #include "sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh" @@ -186,11 +187,8 @@ namespace opdet { std::vector>> fSinglePEWave_HD; // single photon pulse vector int pulsesize; //size of 1PE waveform std::unordered_map< raw::Channel_t, std::vector > fFullWaveforms; - - - std::vector> fSinglePEWaveVector; - std::vector fSinglePEChannels; - std::vector fPeakAmplitude; + + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); void CreatePDWaveformUncoatedPMT( sim::SimPhotons const& SimPhotons, From 72155095cefb787c226384cf7462880470ad76cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 09:46:23 -0600 Subject: [PATCH 19/28] Update default configuration --- sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl | 2 +- sbndcode/OpDetSim/digi_pmt_sbnd.fcl | 12 ++++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl b/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl index 344293a24..5fdeca0ef 100644 --- a/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl +++ b/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl @@ -19,7 +19,7 @@ PMTNonLinearityTF1ChannelByChannel: #AttenuationFormParams: [78.15, 1.96] AttenuationForm: "x/sqrt(1+(x/[0])**[1])" AttenuationPreTime: 4 #ns - NonLinearRange: [100, 701] + NonLinearRange: [100, 7000] } END_PROLOG diff --git a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl index a5d6ca7d2..10a49f837 100644 --- a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl +++ b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl @@ -15,8 +15,8 @@ sbnd_digipmt_alg: PMTChargeToADC: -25.97 #charge to adc factor # Parameters for test bench SER simulation - PMTSinglePEmodel: true #false for ideal PMT response, true for test bench measured response - PMTDataFile: "OpDetSim/digi_pmt_sbnd_v2int0.root" # located in sbnd_data + PMTSinglePEmodel: "data" + PMTDataFile: "OpDetSim/digi_pmt_sbnd_v2int0.root" # located in sbnd_data # Time delays TTS: 2.4 #Transit Time Spread in ns @@ -27,6 +27,10 @@ sbnd_digipmt_alg: PMTBaseline: 14745 #in ADC PMTBaselineRMS: 2.6 #in ADC + #Sample noise from data + UseDataNoise: true + OpDetNoiseFile: "PMTNoiseTemplate.root" + # Dark counts PMTDarkNoiseRate: 1000.0 #in Hz @@ -41,11 +45,11 @@ sbnd_digipmt_alg: # Simulate gain fluctuations # (comment-out to skip gain fluctuations simulation) - GainFluctuationsParams: @local::FirstDynodeGainFluctuations + GainFluctuationsParams: @local::GaussianGainFluctuations # Simulate PMT non linear response # (comment-out to skip non linearities simulation) - NonLinearityParams: @local::PMTNonLinearityTF1 + NonLinearityParams: @local::PMTNonLinearityTF1ChannelByChannel HDOpticalWaveformParamsPMT: @local::IncludeHDOpticalWaveforms_PMT } From caf4bbdcf00823fe6b08da83c79c81b0d976d929 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 12:22:35 -0600 Subject: [PATCH 20/28] Remove cout --- .../OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc index faaf754d5..3f514fb17 100644 --- a/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc +++ b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc @@ -136,7 +136,6 @@ void opdet::PMTPulseOscillation::CreateOscillatedWaveform(const raw::OpDetWavefo double phase = 2.0 * M_PI * fOscillationFrequency * delta_time; double new_value = static_cast(newWaveform.Waveform()[i+fOscillationOffset]) + oscillation_amplitude * envelope * std::cos(phase+M_PI); // Add some dither to avoid quantization effects - std::cout << CLHEP::RandFlat::shoot(&fEngine,1.0) << std::endl; double dither = CLHEP::RandFlat::shoot(&fEngine,1.0) - 0.5; newWaveform.Waveform()[i+fOscillationOffset] = std::round(new_value + dither); if(envelope < 0.05) return; From f5b2fc970efca30622fa4c850cee82bcd11aa832 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Wed, 19 Nov 2025 12:31:27 -0600 Subject: [PATCH 21/28] Modify standard config to include realistic PMT MC and major fcl refactoring --- .../standard/reco/config/workflow_reco1.fcl | 18 +++---- .../standard/reco/reco1_data.fcl | 6 +-- .../standard/standard_detsim_sbnd.fcl | 4 +- .../job/opdeconvolution_sbnd.fcl | 31 +++++++++++ .../job/opdeconvolution_sbnd_data.fcl | 24 --------- .../job/run_decohitfinder_data.fcl | 5 +- .../job/sbnd_flashfinder_deco.fcl | 53 +++++++++++++++++++ .../job/sbnd_flashfinder_deco_data.fcl | 27 ---------- .../job/sbnd_ophitfinder_deco.fcl | 43 +++++++++++++++ .../job/sbnd_ophitfinder_deco_data.fcl | 36 ------------- .../OpFlash/job/sbnd_flashfinder_data.fcl | 36 ------------- 11 files changed, 144 insertions(+), 139 deletions(-) delete mode 100644 sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl delete mode 100644 sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl delete mode 100644 sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl delete mode 100644 sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl diff --git a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl index 9a375942d..5441cc73f 100644 --- a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl +++ b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl @@ -2,9 +2,9 @@ # #include "caldata_sbnd.fcl" #include "hitfindermodules_sbnd.fcl" -#include "opdeconvolution_sbnd_data.fcl" -#include "sbnd_ophitfinder_deco_data.fcl" -#include "sbnd_flashfinder_deco_data.fcl" +#include "opdeconvolution_sbnd.fcl" +#include "sbnd_ophitfinder_deco.fcl" +#include "sbnd_flashfinder_deco.fcl" #include "crtrecoproducers_sbnd.fcl" #include "mlreco_sbnd.fcl" #include "cluster_sbnd.fcl" @@ -18,18 +18,18 @@ sbnd_reco1_producers:{ rns: { module_type: RandomNumberSaver } ### optical deconvolution - opdecopmt: @local::SBNDOpDeconvolutionPMT - opdecoxarapuca: @local::SBNDOpDeconvolutionXARAPUCA + opdecopmt: @local::SBNDOpDeconvolutionPMT_realisticMC + opdecoxarapuca: @local::SBNDOpDeconvolutionXARAPUCA ### optical hit finders # ophit: @local::sbnd_hit_finder - ophitpmt: @local::SBNDDecoOpHitFinderPMT - ophitxarapuca: @local::SBNDDecoOpHitFinderXArapuca + ophitpmt: @local::SBNDDecoOpHitFinderPMT_realisticMC + ophitxarapuca: @local::SBNDDecoOpHitFinderXArapuca ### flash finders # opflash: @local::sbnd_opflash - opflashtpc0: @local::SBNDDecoSimpleFlashTPC0 - opflashtpc1: @local::SBNDDecoSimpleFlashTPC1 + opflashtpc0: @local::SBNDDecoSimpleFlashTPC0_realisticMC + opflashtpc1: @local::SBNDDecoSimpleFlashTPC1_realisticMC ## opflash(arapucas): @local::sbnd_opflash opflashtpc0xarapuca: @local::SBNDDecoSimpleFlashTPC0Arapuca diff --git a/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl index 2c4a44285..b759287e1 100644 --- a/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl @@ -1,7 +1,7 @@ #include "wcsp_data_sbnd.fcl" -#include "opdeconvolution_sbnd_data.fcl" -#include "sbnd_ophitfinder_deco_data.fcl" -#include "sbnd_flashfinder_deco_data.fcl" +#include "opdeconvolution_sbnd.fcl" +#include "sbnd_ophitfinder_deco.fcl" +#include "sbnd_flashfinder_deco.fcl" #include "crt_channel_map_service.fcl" #include "crt_calib_service.fcl" #include "wfalign_sbnd_data.fcl" diff --git a/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl b/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl index 7c0a8e0a4..dea189cad 100644 --- a/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl +++ b/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl @@ -34,6 +34,7 @@ #include "wcsimsp_sbnd.fcl" #include "crtsimmodules_sbnd.fcl" #include "opdetdigitizer_sbnd.fcl" +#include "sbnd_pmtpulseoscillation_config.fcl" #include "rootoutput_sbnd.fcl" process_name: DetSim @@ -66,10 +67,11 @@ physics: simtpc2d: @local::sbnd_wcls_simsp crtsim: @local::sbnd_crtsim opdaq: @local::sbnd_opdetdigitizer + pmtpulseoscillation: @local::sbnd_pmtpulseoscillation } #define the producer and filter modules for this path, order matters, - simulate: [rns, simtpc2d, crtsim, opdaq] + simulate: [rns, simtpc2d, crtsim, opdaq, pmtpulseoscillation] #define the output stream, there could be more than one if using filters stream1: [ out1 ] diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl index 9599a707c..0d7790bf5 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl @@ -1,4 +1,5 @@ #include "opdeconvolution_alg.fcl" +#include "opdeconvolution_alg_data.fcl" BEGIN_PROLOG @@ -11,6 +12,8 @@ SBNDOpDeconvolution: OpDecoAlg: @local::OpDeconvolutionAlg } +###### PMT IDEAL MC ###### + SBNDOpDeconvolutionPMT: @local::SBNDOpDeconvolution SBNDOpDeconvolutionPMT.PDTypes: ["pmt_coated", "pmt_uncoated"] SBNDOpDeconvolutionPMT.Electronics: [""] @@ -20,6 +23,34 @@ SBNDOpDeconvolutionPMT.OpDecoAlg.FilterParams: [0.049, 2] #Freq in GHz SBNDOpDeconvolutionPMT.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 SBNDOpDeconvolutionPMT.OpDecoAlg.DecoWaveformPrecision: 0.005 + +###### PMT DATA ###### + +SBNDOpDeconvolutionPMT_data: @local::SBNDOpDeconvolution +SBNDOpDeconvolutionPMT_data.OpDecoAlg: @local::OpDeconvolutionAlgData +SBNDOpDeconvolutionPMT_data.InputLabel: "wfalign:PMTChannels" +SBNDOpDeconvolutionPMT_data.PDTypes: ["pmt_coated", "pmt_uncoated"] +SBNDOpDeconvolutionPMT_data.Electronics: [""] +SBNDOpDeconvolutionPMT_data.OpDecoAlg.OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root" +SBNDOpDeconvolutionPMT_data.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 +SBNDOpDeconvolutionPMT_data.OpDecoAlg.DecoWaveformPrecision: 0.005 +SBNDOpDeconvolutionPMT_data.OpDecoAlg.SkipChannelList: [39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302] + +###### PMT REALISTIC MC ###### + +SBNDOpDeconvolutionPMT_realisticMC: @local::SBNDOpDeconvolution +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg: @local::OpDeconvolutionAlgData +SBNDOpDeconvolutionPMT_realisticMC.InputLabel: "pmtpulseoscillation" +SBNDOpDeconvolutionPMT_realisticMC.PDTypes: ["pmt_coated", "pmt_uncoated"] +SBNDOpDeconvolutionPMT_realisticMC.Electronics: [""] +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root" +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.DecoWaveformPrecision: 0.005 +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.SkipChannelList: [39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302] + + +###### XA IDEAL MC ###### + SBNDOpDeconvolutionXARAPUCA: @local::SBNDOpDeconvolution SBNDOpDeconvolutionXARAPUCA.PDTypes: ["xarapuca_vuv", "xarapuca_vis"] SBNDOpDeconvolutionXARAPUCA.Electronics: ["daphne"] diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl deleted file mode 100644 index 66ad23bc2..000000000 --- a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl +++ /dev/null @@ -1,24 +0,0 @@ -#include "opdeconvolution_alg_data.fcl" - -BEGIN_PROLOG - -SBNDOpDeconvolution: -{ - module_type: "SBNDOpDeconvolution" - InputLabel: "wfalign:PMTChannels" - PDTypes: [] - Electronics: [] - OpDecoAlg: @local::OpDeconvolutionAlgData -} - -SBNDOpDeconvolutionPMT_data: @local::SBNDOpDeconvolution -SBNDOpDeconvolutionPMT_data.PDTypes: ["pmt_coated", "pmt_uncoated"] -SBNDOpDeconvolutionPMT_data.Electronics: [""] -SBNDOpDeconvolutionPMT_data.OpDecoAlg.OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root" -#SBNDOpDeconvolutionPMT_data.OpDecoAlg.UseParamFilter: true -#SBNDOpDeconvolutionPMT_data.OpDecoAlg.FilterParams: [0.049, 2] #Freq in GHz -SBNDOpDeconvolutionPMT_data.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 -SBNDOpDeconvolutionPMT_data.OpDecoAlg.DecoWaveformPrecision: 0.005 -SBNDOpDeconvolutionPMT_data.OpDecoAlg.SkipChannelList: [39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302] - -END_PROLOG diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl index f3d64fb84..c5fae55b8 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl @@ -3,9 +3,9 @@ #include "messages_sbnd.fcl" #include "sam_sbnd.fcl" -#include "opdeconvolution_sbnd_data.fcl" +#include "opdeconvolution_sbnd.fcl" #include "sbnd_flashfinder_deco_data.fcl" -#include "sbnd_ophitfinder_deco_data.fcl" +#include "sbnd_ophitfinder_deco.fcl" process_name: DecoRecoData @@ -54,7 +54,6 @@ physics: opflashtpc1 ] - stream1: [out1] trigger_paths: [reco] diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl index fd231e3f8..41d453e49 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl @@ -2,7 +2,11 @@ BEGIN_PROLOG + ####OpFlash finder for PMT deconvolved waveforms##### + + +###### PMT IDEAL MC ###### ###TPC0 SBNDDecoSimpleFlashTPC0: @local::SBNDSimpleFlashTPC0 SBNDDecoSimpleFlashTPC0.PECalib.SPEAreaGain: 200 @@ -22,8 +26,57 @@ SBNDDecoSimpleFlashTPC1.ReadoutDelay: 0.135 //cable time delay in us SBNDDecoSimpleFlashTPC1.CorrectLightPropagation: true +###### PMT DATA ###### +###TPC0 +SBNDDecoSimpleFlashTPC0_data: @local::SBNDSimpleFlashTPC0 +SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC0_data.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC0_data.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC0_data.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC0_data.UseT0Tool: true +SBNDDecoSimpleFlashTPC0_data.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC0_data.CorrectLightPropagation: true +SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data_v2.root" + +#TPC1 +SBNDDecoSimpleFlashTPC1_data: @local::SBNDSimpleFlashTPC1 +SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC1_data.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC1_data.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC1_data.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC1_data.UseT0Tool: true +SBNDDecoSimpleFlashTPC1_data.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC1_data.CorrectLightPropagation: true + + +###### PMT REALISTIC MC ###### +###TPC0 +SBNDDecoSimpleFlashTPC0_realisticMC: @local::SBNDSimpleFlashTPC0 +SBNDDecoSimpleFlashTPC0_realisticMC.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC0_realisticMC.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC0_realisticMC.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC0_realisticMC.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC0_realisticMC.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC0_realisticMC.UseT0Tool: true +SBNDDecoSimpleFlashTPC0_realisticMC.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC0_realisticMC.CorrectLightPropagation: true + +#TPC1 +SBNDDecoSimpleFlashTPC1_realisticMC: @local::SBNDSimpleFlashTPC1 +SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC1_realisticMC.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC1_realisticMC.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC1_realisticMC.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC1_realisticMC.UseT0Tool: true +SBNDDecoSimpleFlashTPC1_realisticMC.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC1_realisticMC.CorrectLightPropagation: true +SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data_v2.root" +###### XA IDEAL MC ###### ####OpFlash finder for XArapucas deconvolved waveforms##### ###TPC0 diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl deleted file mode 100644 index 45b3aa338..000000000 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl +++ /dev/null @@ -1,27 +0,0 @@ -#include "sbnd_flashfinder_data.fcl" - -BEGIN_PROLOG - -####OpFlash finder for PMT deconvolved waveforms##### -###TPC0 -SBNDDecoSimpleFlashTPC0_data: @local::SBNDSimpleFlashTPC0_data -SBNDDecoSimpleFlashTPC0_data.PECalib.SPEAreaGain: 200 -SBNDDecoSimpleFlashTPC0_data.OpHitProducers: ["ophitpmt"] -SBNDDecoSimpleFlashTPC0_data.OpHitInputTime: "RiseTime" -SBNDDecoSimpleFlashTPC0_data.UseT0Tool: true -SBNDDecoSimpleFlashTPC0_data.ReadoutDelay: 0 //cable time delay in us -SBNDDecoSimpleFlashTPC0_data.CorrectLightPropagation: true -SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data.root" - -#TPC1 -SBNDDecoSimpleFlashTPC1_data: @local::SBNDSimpleFlashTPC1_data -SBNDDecoSimpleFlashTPC1_data.PECalib.SPEAreaGain: 200 -SBNDDecoSimpleFlashTPC1_data.OpHitProducers: ["ophitpmt"] -SBNDDecoSimpleFlashTPC1_data.OpHitInputTime: "RiseTime" -SBNDDecoSimpleFlashTPC1_data.UseT0Tool: true -SBNDDecoSimpleFlashTPC1_data.ReadoutDelay: 0 //cable time delay in us -SBNDDecoSimpleFlashTPC1_data.CorrectLightPropagation: true -SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data.root" - - -END_PROLOG diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl index 2911860a9..989fe0193 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl @@ -2,6 +2,9 @@ BEGIN_PROLOG +###### MC CONFIGURATION ###### + + #####OpHit finder for PMT deconvolved waveforms##### SBNDDecoOpHitFinderPMT: @local::sbnd_ophit_finder_pmt ####SPE area must be 1./DecoWaveformPrecision @@ -63,4 +66,44 @@ SBNDDecoOpHitFinderXArapuca.PedAlgoPset.NumSampleTail: 50 SBNDDecoOpHitFinderXArapuca.PedAlgoPset.Method:2 +###### DATA CONFIGURATION ###### + +#####OpHit finder for PMT deconvolved waveforms##### +SBNDDecoOpHitFinderPMT_data: @local::sbnd_ophit_finder_pmt +####SPE area must be 1./DecoWaveformPrecision +SBNDDecoOpHitFinderPMT_data.SPEArea: 200 +SBNDDecoOpHitFinderPMT_data.InputModule: "opdecopmt" +SBNDDecoOpHitFinderPMT_data.HitThreshold: 1 +SBNDDecoOpHitFinderPMT_data.RiseTimeCalculator: @local::sbnd_opreco_risetimecalculator + +#HitAlgoPset +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThreshold: 25 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdByChannel: true +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdVector: [50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 18.0, 35.0, 19.0, 21.0, 19.0, 20.0, 21.0, 27.0, 19.0, 27.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 19.0, 37.0, 50.0, 20.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 24.0, 19.0, 26.0, 19.0, 12.0, 30.0, 50.0, 25.0, 14.0, 12.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 50.0, 50.0, 50.0, 19.0, 19.0, 24.0, 21.0, 28.0, 19.0, 22.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 50.0, 19.0, 22.0, 35.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 30.0, 50.0, 20.0, 25.0, 23.0, 32.0, 20.0, 22.0, 20.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 23.0, 50.0, 19.0, 25.0, 19.0, 18.0, 19.0, 15.0, 21.0, 24.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 20.0, 30.0, 26.0, 18.0, 34.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 27.0, 50.0, 50.0, 21.0, 28.0, 50.0, 50.0, 50.0, 19.0, 19.0, 50.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 19.0, 30.0, 20.0, 18.0, 50.0, 18.0, 19.0, 13.0, 50.0, 25.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 25.0, 22.0, 19.0, 29.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 23.0, 29.0, 19.0, 27.0, 20.0, 25.0, 22.0, 50.0, 19.0, 25.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0] +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NSigmaThreshold: 3.4 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndADCThreshold: 8 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndNSigmaThreshold: 0.2 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.MinPulseWidth: 4 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Name: "SlidingWindow" +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPostSample: 6 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPreSample: 3 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.PositivePolarity: true +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailADCThreshold: 2 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailNSigmaThreshold: 2 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Verbosity: false + +#BaselinePset +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Name:"Edges" +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleFront:200 +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleTail:200 +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Method:2 + + +###### REALISTICMC CONFIGURATION ###### + +#####OpHit finder for PMT deconvolved waveforms##### +SBNDDecoOpHitFinderPMT_realisticMC: @local::SBNDDecoOpHitFinderPMT_data +####SPE area must be 1./DecoWaveformPrecision +SBNDDecoOpHitFinderPMT_realisticMC.InputModule: "pmtpulseoscillation" + END_PROLOG diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl deleted file mode 100644 index 02eadf426..000000000 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl +++ /dev/null @@ -1,36 +0,0 @@ -#include "ophit_finder_sbnd.fcl" - -BEGIN_PROLOG - -#####OpHit finder for PMT deconvolved waveforms##### -SBNDDecoOpHitFinderPMT_data: @local::sbnd_ophit_finder_pmt -####SPE area must be 1./DecoWaveformPrecision -SBNDDecoOpHitFinderPMT_data.SPEArea: 200 -SBNDDecoOpHitFinderPMT_data.InputModule: "opdecopmt" -SBNDDecoOpHitFinderPMT_data.HitThreshold: 1 -SBNDDecoOpHitFinderPMT_data.RiseTimeCalculator: @local::sbnd_opreco_risetimecalculator - -#HitAlgoPset -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThreshold: 25 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdByChannel: true -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdVector: [50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 18.0, 35.0, 19.0, 21.0, 19.0, 20.0, 21.0, 27.0, 19.0, 27.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 19.0, 37.0, 50.0, 20.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 24.0, 19.0, 26.0, 19.0, 12.0, 30.0, 50.0, 25.0, 14.0, 12.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 50.0, 50.0, 50.0, 19.0, 19.0, 24.0, 21.0, 28.0, 19.0, 22.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 50.0, 19.0, 22.0, 35.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 30.0, 50.0, 20.0, 25.0, 23.0, 32.0, 20.0, 22.0, 20.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 23.0, 50.0, 19.0, 25.0, 19.0, 18.0, 19.0, 15.0, 21.0, 24.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 20.0, 30.0, 26.0, 18.0, 34.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 27.0, 50.0, 50.0, 21.0, 28.0, 50.0, 50.0, 50.0, 19.0, 19.0, 50.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 19.0, 30.0, 20.0, 18.0, 50.0, 18.0, 19.0, 13.0, 50.0, 25.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 25.0, 22.0, 19.0, 29.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 23.0, 29.0, 19.0, 27.0, 20.0, 25.0, 22.0, 50.0, 19.0, 25.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0] -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NSigmaThreshold: 3.4 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndADCThreshold: 8 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndNSigmaThreshold: 0.2 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.MinPulseWidth: 4 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Name: "SlidingWindow" -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPostSample: 6 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPreSample: 3 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.PositivePolarity: true -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailADCThreshold: 2 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailNSigmaThreshold: 2 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Verbosity: false - -#BaselinePset -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Name:"Edges" -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleFront:200 -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleTail:200 -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Method:2 - - -END_PROLOG diff --git a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl b/sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl deleted file mode 100644 index 19c4ccfa1..000000000 --- a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl +++ /dev/null @@ -1,36 +0,0 @@ -#include "sbnd_flashalgo.fcl" -#include "sbnd_flashcalib.fcl" -#include "sbnd_flashgeoalgo.fcl" -#include "sbnd_flasht0algo.fcl" -#include "sbnd_driftestimatoralgo.fcl" - -BEGIN_PROLOG - -SBNDSimpleFlash: -{ - module_type : "SBNDFlashFinder" - FlashFinderAlgo : "SimpleFlashAlgo" - AlgoConfig : @local::SimpleFlashStandard - OpHitProducers : ["ophitpmt","ophitarapuca"] - OpFlashProducer : "opflash" - PECalib : @local::NoCalib - OpHitInputTime : "PeakTime" - FlashGeoConfig : @local::FlashGeoThreshold - UseT0Tool : false - FlashT0Config : @local::FlashT0SelectedChannels - ReadoutDelay : 0. //in us - CorrectLightPropagation : false - DriftEstimatorConfig : @local::DriftEstimatorPMTRatio -} - -SBNDSimpleFlashTPC0_data: @local::SBNDSimpleFlash -SBNDSimpleFlashTPC0_data.AlgoConfig: @local::SimpleFlashTPC0 -SBNDSimpleFlashTPC0_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" -SBNDSimpleFlashTPC0_data.DriftEstimatorConfig.DataCalibration: true - -SBNDSimpleFlashTPC1_data: @local::SBNDSimpleFlash -SBNDSimpleFlashTPC1_data.AlgoConfig: @local::SimpleFlashTPC1 -SBNDSimpleFlashTPC1_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" -SBNDSimpleFlashTPC1_data.DriftEstimatorConfig.DataCalibration: true - -END_PROLOG From 1e662117519eac92030a07ce9b41bb082510946c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 05:04:22 -0600 Subject: [PATCH 22/28] Fix comment on PMT response --- sbndcode/OpDetSim/DigiPMTSBNDAlg.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh index 712d65eb2..ee2c68852 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh @@ -324,7 +324,7 @@ namespace opdet { fhicl::Atom PMTsinglePEmodel { Name("PMTSinglePEmodel"), - Comment("Model used for single PE response of PMT. =0 is ideal, =1 is testbench") + Comment("Model used for single PE response of PMT. Accepted strings are \"testbench\", \"data\", and \"ideal\"") }; fhicl::Atom pmtDataFile { From d05f478e33c43490558855392eb068f6d9ea2638 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 05:10:54 -0600 Subject: [PATCH 23/28] Remove CRT includes --- sbndcode/JobConfigurations/standard/reco/reco1_data.fcl | 2 -- 1 file changed, 2 deletions(-) diff --git a/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl index eeeafe464..4b6ae0e4d 100755 --- a/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl @@ -3,8 +3,6 @@ #include "opdeconvolution_sbnd.fcl" #include "sbnd_ophitfinder_deco.fcl" #include "sbnd_flashfinder_deco.fcl" -#include "crt_channel_map_service.fcl" -#include "crt_calib_service.fcl" #include "wfalign_sbnd_data.fcl" #include "standard_reco1_sbnd.fcl" From 6db4328c51d3204c9eb3d41685f8c9cef631b150 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 06:00:44 -0600 Subject: [PATCH 24/28] Treat correctly the on/off calibrated/uncalibrated PMTs --- sbndcode/OpDetSim/DigiPMTSBNDAlg.cc | 59 ++++++++++++++++++++++------- sbndcode/OpDetSim/DigiPMTSBNDAlg.hh | 2 + 2 files changed, 48 insertions(+), 13 deletions(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index 743172a20..6fc26b7ee 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -61,10 +61,10 @@ namespace opdet { // PMT calibration database service fPMTCalibrationDatabaseService = lar::providerFrom(); fPMTHDOpticalWaveformsPtr = art::make_tool(fParams.HDOpticalWaveformParams); - int NOpChannels = fWireReadout.NOpChannels(); //Resize the SER vector to the number of channels fSinglePEWave_HD.resize(NOpChannels); + //shape of single pulse if (fParams.PMTSinglePEmodel=="testbench") { mf::LogDebug("DigiPMTSBNDAlg") << " using testbench pe response"; @@ -89,17 +89,51 @@ namespace opdet { Pulse1PE(fSinglePEWave); } else if (fParams.PMTSinglePEmodel=="data"){ - // Read the SER from the calibration database (channel independent) + //Build the average data PMT waveform for channels with uncalibrated SER + int dataSERSize; + // Retrieve the data SER vector size from the first calibrated channel in the database + for(size_t i=0; igetReconstructChannel(i)){ + dataSERSize = fPMTCalibrationDatabaseService->getSER(i).size(); + break; + } + } + fAverageDataSER.resize(dataSERSize); + fAverageDataSPEAmplitude = 0; + int nChannelsWithDataSER = 0; for(size_t i=0; igetReconstructChannel(i)) continue; - fSinglePEWave = fPMTCalibrationDatabaseService->getSER(i); - double SPEAmplitude = fPMTCalibrationDatabaseService->getSPEAmplitude(i); - double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);}); - double SinglePENormalization = std::abs(SPEAmplitude/SPEPeakValue); - std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;}); - if(fSinglePEWave.size()==0) continue; - if(fSinglePEWave.size()>0) pulsesize = fSinglePEWave.size(); - fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD[i], fSinglePEWave); + for (size_t j = 0; j < fAverageDataSER.size() && j < fPMTCalibrationDatabaseService->getSER(i).size(); ++j) { fAverageDataSER[j] += fPMTCalibrationDatabaseService->getSER(i)[j]; } + fAverageDataSPEAmplitude += fPMTCalibrationDatabaseService->getSPEAmplitude(i); + nChannelsWithDataSER++; + } + std::transform(fAverageDataSER.begin(), fAverageDataSER.end(), fAverageDataSER.begin(), [nChannelsWithDataSER](double val) {return val / nChannelsWithDataSER;}); + double PeakValue = *std::max_element(fAverageDataSER.begin(), fAverageDataSER.end(), [](double a, double b) {return std::abs(a) < std::abs(b);}); + fAverageDataSPEAmplitude = fAverageDataSPEAmplitude / nChannelsWithDataSER; + // Normalise the average data SER to the average SPE amplitude + double AverageDataPENormalization = std::abs(fAverageDataSPEAmplitude/PeakValue); + std::transform(fAverageDataSER.begin(), fAverageDataSER.end(), fAverageDataSER.begin(), [AverageDataPENormalization](double val) {return val * AverageDataPENormalization;}); + // Read the SER from the calibration database (channel independent) + for(size_t i=0; igetOnPMT(i)){ + // PMT is OFF + continue; + } + else if( fPMTCalibrationDatabaseService->getOnPMT(i) && !fPMTCalibrationDatabaseService->getReconstructChannel(i)){ + // PMT is ON but does not have a calibrated SER. These are only relevant for trigger emulation, not used in downstream reconstruciton + fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD[i], fAverageDataSER); + } + else{ + // PMT is ON and has a calibrated SER + fSinglePEWave = fPMTCalibrationDatabaseService->getSER(i); + double SPEAmplitude = fPMTCalibrationDatabaseService->getSPEAmplitude(i); + double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);}); + double SinglePENormalization = std::abs(SPEAmplitude/SPEPeakValue); + std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;}); + if(fSinglePEWave.size()==0) continue; + if(fSinglePEWave.size()>0) pulsesize = fSinglePEWave.size(); + fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD[i], fSinglePEWave); + } } } else{ @@ -129,7 +163,6 @@ namespace opdet { sp_noise.find_file(fParams.OpDetNoiseFile, fname_noise); noise_file = TFile::Open(fname_noise.c_str(), "READ"); } - } // end constructor DigiPMTSBNDAlg::~DigiPMTSBNDAlg(){ @@ -491,8 +524,8 @@ namespace opdet { void DigiPMTSBNDAlg::AddSPE(size_t time, std::vector& wave, int ch, double npe) - { - if(!fPMTCalibrationDatabaseService->getReconstructChannel(ch)) return; + { + if(!fPMTCalibrationDatabaseService->getOnPMT(ch)) return; // time bin HD (double precision) // used to gert the time-shifted SER double time_bin_hd = fSampling*time; diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh index ee2c68852..7a7d7bcd0 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh @@ -184,6 +184,8 @@ namespace opdet { double Transittimespread(double fwhm); std::vector fSinglePEWave; // single photon pulse vector + std::vector fAverageDataSER; // single photon pulse vector + double fAverageDataSPEAmplitude; std::vector>> fSinglePEWave_HD; // single photon pulse vector int pulsesize; //size of 1PE waveform std::unordered_map< raw::Channel_t, std::vector > fFullWaveforms; From 54714692b00656c1f8f606c45253d74c3e6f7862 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 06:07:28 -0600 Subject: [PATCH 25/28] Change PMT detection efficiencies --- sbndcode/OpDetSim/digi_pmt_sbnd.fcl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl index 10a49f837..a9ad05338 100644 --- a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl +++ b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl @@ -35,13 +35,13 @@ sbnd_digipmt_alg: PMTDarkNoiseRate: 1000.0 #in Hz # Detection efficiencies - PMTCoatedVUVEff_tpc0: 0.035 #PMT coated detection efficiency for direct (VUV) light - PMTCoatedVISEff_tpc0: 0.03882 #PMT coated detection efficiency for reflected (VIS) light - PMTUncoatedEff_tpc0: 0.03758 #PMT uncoated detection efficiency + PMTCoatedVUVEff_tpc0: 0.0315 #PMT coated detection efficiency for direct (VUV) light + PMTCoatedVISEff_tpc0: 0.03493 #PMT coated detection efficiency for reflected (VIS) light + PMTUncoatedEff_tpc0: 0.03382 #PMT uncoated detection efficiency - PMTCoatedVUVEff_tpc1: 0.03 #PMT coated detection efficiency for direct (VUV) light - PMTCoatedVISEff_tpc1: 0.03882 #PMT coated detection efficiency for reflected (VIS) light - PMTUncoatedEff_tpc1: 0.03758 #PMT uncoated detection efficiency + PMTCoatedVUVEff_tpc1: 0.027 #PMT coated detection efficiency for direct (VUV) light + PMTCoatedVISEff_tpc1: 0.03493 #PMT coated detection efficiency for reflected (VIS) light + PMTUncoatedEff_tpc1: 0.03382 #PMT uncoated detection efficiency # Simulate gain fluctuations # (comment-out to skip gain fluctuations simulation) From 4a6aeb5779787a190aff41bb8337bcb0ca5f60ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 06:09:33 -0600 Subject: [PATCH 26/28] Translate comments on the code --- sbndcode/OpDetSim/DigiPMTSBNDAlg.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index 6fc26b7ee..345a80bb2 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -620,7 +620,7 @@ namespace opdet { int waveBins = wave.size(); int currentSize = 0; - // Llenar el vector de ruido hasta igualar la longitud del waveform + // Fill noise vector until it matches the waveform length while (currentSize < waveBins) { int noiseWformIdx = static_cast(fEngine->flat() * keylist.size()); TH1 *noiseHist = (TH1*)keylist[noiseWformIdx]->ReadObj(); @@ -633,7 +633,7 @@ namespace opdet { delete noiseHist; } - // Agregar ruido al waveform + // Add noise to the waveform for (size_t i = 0; i < wave.size(); i++) { wave[i] += noise_wform[i]; } From ed1d23deacd07f9f716dd96cec6a25ac87e71e53 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 12:30:34 -0600 Subject: [PATCH 27/28] Add PMTPulseOscillation subdir --- sbndcode/OpDetSim/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/sbndcode/OpDetSim/CMakeLists.txt b/sbndcode/OpDetSim/CMakeLists.txt index 796328b68..35c22cf77 100755 --- a/sbndcode/OpDetSim/CMakeLists.txt +++ b/sbndcode/OpDetSim/CMakeLists.txt @@ -177,6 +177,7 @@ install_source() cet_enable_asserts() add_subdirectory(PMTAlg) add_subdirectory(HDWvf) +add_subdirectory(PMTPulseOscillations) # install sbnd_pds_mapping.json with mapping of the photon detectors install_fw(LIST sbnd_pds_mapping.json) From dc4ef3ed3be97a24c13491de9306d8b4ef731c6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Fri, 21 Nov 2025 16:09:07 -0600 Subject: [PATCH 28/28] Fix path within sbnd_data --- sbndcode/OpDetSim/digi_pmt_sbnd.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl index a9ad05338..801220e7a 100644 --- a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl +++ b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl @@ -29,7 +29,7 @@ sbnd_digipmt_alg: #Sample noise from data UseDataNoise: true - OpDetNoiseFile: "PMTNoiseTemplate.root" + OpDetNoiseFile: "OpDetReco/PMTNoiseTemplate.root" # Dark counts PMTDarkNoiseRate: 1000.0 #in Hz