diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index e30bf2cc7..8a1097e2f 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -16,6 +16,9 @@ // LArSoft libraries #include "lardataalg/DetectorInfo/DetectorTimings.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" +#include "larcorealg/Geometry/OpDetGeo.h" +#include "larcorealg/Geometry/geo_vectors_utils.h" // geo::vect::toPoint() #include "larcorealg/CoreUtils/counter.h" #include "larcorealg/CoreUtils/StdUtils.h" // util::begin(), util::end() @@ -221,7 +224,7 @@ auto icarus::opdet::PMTsimulationAlg::makeGainFluctuator(int channel) const if (fParams.useGainCalibDB && fParams.gainCalibProvider) { // DB Gaussian: per-channel SPE area and width from database // fNominalSPEArea is the integral of the SPR template, i.e. the mean area per PE - // fBiasConstant covers the bias in the integral definitions btw SPR template and official reco + // fBiasConstant covers the bias in the integral definitions btw SPR template and official reco (SBN-docdb-46261) double const speArea = fParams.gainCalibProvider->getSPEArea(channel); double const speFitWidth = fParams.gainCalibProvider->getSPEFitWidth(channel); // gainRatio = speArea / fNominalSPEArea: mean effective PEs per true PE @@ -300,8 +303,19 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform photons_used->clear(); photons_used->SetChannel(channel); } + // distance-dependent survival S(d): PMT center cached per channel + bool const applySd = fParams.distanceSurvival.apply; + geo::Point_t const PMTcenter = fParams.wireReadoutGeom->OpDetGeoFromOpChannel(channel).GetCenter(); + for(auto const& ph : photons) { - if (!KicksPhotoelectron()) continue; + if (applySd) { + // Bernoulli survival p = QE * S(d), same stream as KicksPhotoelectron + double const d_cm + = (PMTcenter - geo::vect::toPoint(ph.InitialPosition)).R(); + if (CLHEP::RandFlat::shoot(fParams.randomEngine) + >= fQE * distanceSurvival(d_cm)) continue; + } + else if (!KicksPhotoelectron()) continue; if (photons_used) photons_used->push_back(ph); // copy @@ -429,6 +443,7 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform // std::cout << "\tadded pes... " << channel << " " << diff.count() << std::endl; // start=std::chrono::high_resolution_clock::now(); + ApplyTailSuppression(waveform); AddPedestal(channel, waveformStartTS, waveform); if(fParams.darkNoiseRate > 0.0_Hz) AddDarkNoise(waveform, channel); @@ -713,6 +728,16 @@ bool icarus::opdet::PMTsimulationAlg::KicksPhotoelectron() const { return CLHEP::RandFlat::shoot(fParams.randomEngine) < fQE; } +//------------------------------------------------------------------------------ +double icarus::opdet::PMTsimulationAlg::distanceSurvival(double d_cm) const +{ + auto const& edges = fParams.distanceSurvival.binEdges; + auto const it = std::upper_bound(edges.begin(), edges.end(), d_cm); + if ((it == edges.begin()) || (it == edges.end())) return 1.0; + return fParams.distanceSurvival.factors[std::distance(edges.begin(), it) - 1]; +} // icarus::opdet::PMTsimulationAlg::distanceSurvival() + + // ----------------------------------------------------------------------------- void icarus::opdet::PMTsimulationAlg::AddPhotoelectrons( PulseSampling_t const& pulse, Waveform_t& wave, tick const time_bin, @@ -867,6 +892,26 @@ void icarus::opdet::PMTsimulationAlg::ApplySaturation( } // icarus::opdet::PMTsimulationAlg::ApplySaturation(range) +// ----------------------------------------------------------------------------- +void icarus::opdet::PMTsimulationAlg::ApplyTailSuppression(Waveform_t& waveform) const +{ + if (!fParams.tailSuppression.apply) return; + + // tau [ns] * fSampling [MHz] / 1000.0f = tau_ticks + float const beta = 1.0f / (fParams.tailSuppression.tau * fSampling.value() / 1000.0f); + float const epsilon = fParams.tailSuppression.epsilon; + int const pol = fParams.pulsePolarity; + + float Q = 0.0f; + for (auto& sample : waveform) { + float const s = pol * sample.value(); + Q = (1.0f - beta) * Q + beta * std::max(0.0f, s); + sample = ADCcount{ pol * std::max(0.0f, s - epsilon * Q) }; + } + +} // icarus::opdet::PMTsimulationAlg::ApplyTailSuppression() + + // ----------------------------------------------------------------------------- /// Forces `waveform` ADC within the `min` to `max` range (`max` included). void icarus::opdet::PMTsimulationAlg::ClipWaveform @@ -1012,9 +1057,48 @@ icarus::opdet::PMTsimulationAlgMaker::PMTsimulationAlgMaker fBaseConfig.beamGateTriggerNReps = config.BeamGateTriggerNReps(); fBaseConfig.discrimAlgo = config.getDiscriminationAlgo(); + // + // tail suppression + // + if (config.TailSuppression()) { + fBaseConfig.tailSuppression.apply = config.TailSuppression()->Apply(); + fBaseConfig.tailSuppression.epsilon = config.TailSuppression()->Epsilon(); + fBaseConfig.tailSuppression.tau = config.TailSuppression()->Tau(); + } + + // + // distance-dependent photon survival + // + if (config.DistanceSurvival()) { + fBaseConfig.distanceSurvival.apply = config.DistanceSurvival()->Apply(); + fBaseConfig.distanceSurvival.binEdges = config.DistanceSurvival()->BinEdges(); + fBaseConfig.distanceSurvival.factors = config.DistanceSurvival()->Factors(); + } + // // parameter checks // + if (fBaseConfig.distanceSurvival.apply) { + auto const& ds = fBaseConfig.distanceSurvival; + if (ds.binEdges.size() != ds.factors.size() + 1) { + throw cet::exception("PMTsimulationAlg") + << "DistanceSurvival: BinEdges must have exactly one more entry than" + " Factors (got " << ds.binEdges.size() << " edges, " + << ds.factors.size() << " factors)\n"; + } + if (!std::is_sorted(ds.binEdges.begin(), ds.binEdges.end())) { + throw cet::exception("PMTsimulationAlg") + << "DistanceSurvival: BinEdges must be sorted in increasing order\n"; + } + for (double const f: ds.factors) { + if ((f < 0.0) || (f > 1.0)) { + throw cet::exception("PMTsimulationAlg") + << "DistanceSurvival: survival factors must be in [0, 1] (got " + << f << ")\n"; + } + } + } // distance survival checks + if (std::abs(fBaseConfig.pulsePolarity) != 1.0) { throw cet::exception("PMTsimulationAlg") << "Pulse polarity settings can be only +1.0 or -1.0 (got: " @@ -1045,6 +1129,7 @@ std::unique_ptr icarus::opdet::PMTsimulationAlgMaker::operator()( std::uint64_t beamGateTimestamp, detinfo::LArProperties const& larProp, + geo::WireReadoutGeom const& wireReadoutGeom, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, @@ -1058,7 +1143,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( { return std::make_unique(makeParams( beamGateTimestamp, - larProp, clockData, timingDelays, gainCalibProvider, + larProp, wireReadoutGeom, clockData, timingDelays, gainCalibProvider, SPRfunction, pedestalGenerator, mainRandomEngine, darkNoiseRandomEngine, elecNoiseRandomEngine, trackSelectedPhotons @@ -1071,6 +1156,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( auto icarus::opdet::PMTsimulationAlgMaker::makeParams( std::uint64_t beamGateTimestamp, detinfo::LArProperties const& larProp, + geo::WireReadoutGeom const& wireReadoutGeom, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, @@ -1094,6 +1180,7 @@ auto icarus::opdet::PMTsimulationAlgMaker::makeParams( params.beamGateTimestamp = beamGateTimestamp; params.larProp = &larProp; + params.wireReadoutGeom = &wireReadoutGeom; params.clockData = &clockData; params.detTimings = detinfo::makeDetectorTimings(params.clockData); params.timingDelays = timingDelays; diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index 71b66c72c..df4786736 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -28,6 +28,7 @@ // LArSoft libraries #include "lardataobj/RawData/OpDetWaveform.h" #include "lardataobj/Simulation/SimPhotons.h" +#include "larcorealg/Geometry/fwd.h" // geo::WireReadoutGeom #include "lardataalg/DetectorInfo/LArProperties.h" #include "lardataalg/DetectorInfo/DetectorClocksData.h" #include "lardataalg/DetectorInfo/DetectorTimings.h" @@ -43,6 +44,7 @@ #include "messagefacility/MessageLogger/MessageLogger.h" #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/types/OptionalAtom.h" +#include "fhiclcpp/types/OptionalTable.h" #include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/Table.h" @@ -163,10 +165,15 @@ class icarus::opdet::OpDetWaveformMakerClass { * The algorithm creates simulated PMT waveforms as read out by ICARUS, * including the generation of trigger primitives. * Contributions to the waveforms include: - * * physical photons + * * physical photons (optionally gated by a distance-dependent survival + * probability, see below) * * dark noise * * electronics noise * + * In addition, an optional post-processing correction ("tail suppression", + * see below) may be applied to the signal waveform to reproduce the + * behaviour observed in data in the tail of large pulses. + * * The algorithm processes an optical channel at a time, independently * and uncorrelated from the other channels. * For each channel, multiple waveforms may be generated according to the @@ -213,10 +220,22 @@ class icarus::opdet::OpDetWaveformMakerClass { * of the photon converting to a photoelectron, the quantum efficiency * check here should be skipped by setting the efficiency to 1. * + * When enabled by `DistanceSurvival.Apply`, the photoelectron conversion is + * additionally gated by a distance-dependent survival probability @f$ S(d) @f$: + * each photon is kept with total probability @f$ QE \times S(d) @f$, where + * @f$ d @f$ is the straight-line distance from the photon emission point to + * the PMT center. This is drawn from the same random stream as the plain + * quantum-efficiency acceptance, so the two are mutually exclusive at + * configuration time (either `KicksPhotoelectron()` or the @f$ QE \times S(d) @f$ + * gate is used, never both). @f$ S(d) @f$ is a piecewise-constant function of + * the distance (see the "Distance-dependent photon survival" section below). + * Since it requires the photon emission position, this feature needs the full + * `sim::SimPhotons` input (positions), not the position-less `sim::SimPhotonsLite`. + * * For each converting photon, a photoelectron is added to the channel by * placing a template waveform shape into the channel waveform. * - * If enabled by `ApplyTimingDelays`, the timing correction service + * If enabled by `ApplyTimingDelays`, the timing correction service * `icarusDB::IPMTTimingCorrectionService` is used to simulate the chain of * delays between the photon hitting the photocathode and the time its signal * is digitized. The delay is dominated by the signal cable length (~200ns). @@ -320,6 +339,63 @@ class icarus::opdet::OpDetWaveformMakerClass { * Noise can be disabled by using the `PMTnoNoiseGeneratorTool`. * * + * Distance-dependent photon survival (SBN-docdb-48511) + * ---------------------------------------------------- + * + * When enabled via `DistanceSurvival.Apply`, each scintillation photon is + * accepted as a photoelectron with probability @f$ QE \times S(d) @f$ instead + * of the plain quantum-efficiency acceptance, where @f$ d @f$ is the + * straight-line distance [cm] from the photon emission point to the PMT + * center and @f$ S(d) \in [0, 1] @f$ is a survival factor. This provides an + * effective, distance-dependent correction to the light yield on top of the + * upstream photon propagation. + * + * @f$ S(d) @f$ is piecewise constant: it is defined by a set of bin edges + * (`DistanceSurvival.BinEdges`, sorted, in cm) and one survival factor per + * bin (`DistanceSurvival.Factors`, so `BinEdges` has exactly one more entry + * than `Factors`). `Factors[i]` applies to @f$ d \in + * [ \mathrm{BinEdges}[i], \mathrm{BinEdges}[i+1] ) @f$; for distances outside + * the covered range the survival is 1 (photon subject to plain `QE` only). + * The acceptance is drawn from the same random stream as the quantum + * efficiency, and the correction is uniform across all channels. + * + * Because the emission position is required, this feature is only available + * on the full `sim::SimPhotons` input path; if only the position-less + * `sim::SimPhotonsLite` is available the module raises a configuration error + * rather than silently dropping the correction (see `SimPMTIcarus`). + * + * + * Tail suppression (SBN-docdb-48157, SBN-docdb-48511) + * --------------------------------------------------- + * + * In the tail of large pulses, data waveforms sit systematically below the + * level predicted by linearly scaling the SPE template. This correction is an + * empirical, tunable model to reproduce that behaviour in simulation. + * + * When enabled via `TailSuppression.Apply`, a low-pass filter is run over the + * signal waveform (pre-pedestal) to track the recent integrated charge. The + * filter is a recursive (exponential moving) average of the polarity-corrected + * samples: + * @f[ + * Q(n) = (1 - \beta)\,Q(n-1) + \beta\,\max(0,\,s(n)) + * @f] + * where @f$ s(n) @f$ is the polarity-corrected sample and @f$ \beta @f$ is the + * filter smoothing factor. It is the inverse of the decay time constant + * `TailSuppression.Tau` expressed in ticks, i.e. + * @f$ \beta = 1 / (\tau \cdot f_{\mathrm{sampling}}) @f$ with @f$ \tau @f$ the + * `Tau` value [ns] and @f$ f_{\mathrm{sampling}} @f$ the sampling frequency + * [GHz, ticks/ns]. The charge state @f$ Q @f$ is initialized to zero at the start of + * each channel's full waveform and does not persist across channels. + * + * The tracked charge is then subtracted from each sample, scaled by a tunable + * strength parameter `TailSuppression.Epsilon`, with the result clamped so it + * does not cross the baseline: + * @f[ + * s'(n) = \max(0,\, s(n) - \varepsilon\,Q(n)) + * @f] + * (the polarity is restored before the sample is written back). The + * correction parameters are the same across all channels. + * * Configuration * ============== * @@ -347,6 +423,7 @@ class icarus::opdet::OpDetWaveformMakerClass { * @f$ \mu_{i} \propto (\Delta V_{i})^{k} @f$ (with @f$ \mu_{i} @f$ the * gain for stage @f$ i @f$, @f$ \Delta V_{i} @f$ the drop of potential * of that stage and @f$ k @f$ the parameter set by `dynodeK`. + * * * `DiscrimAlgo` (choice, default: `"CrossingThreshold"`): selects one of the * hard-coded discrimination algorithms used for zero suppression. * The suppression algorithm identifies some interesting time points and @@ -366,6 +443,29 @@ class icarus::opdet::OpDetWaveformMakerClass { * always contain just noise. * * + * Signal corrections + * ------------------ + * + * These optional top-level configuration tables control the corrections + * described in the "Distance-dependent photon survival" and "Tail suppression" + * sections above. Both are disabled by default. + * + * * `DistanceSurvival` (table, optional): distance-dependent photon survival + * @f$ S(d) @f$. Sub-parameters: + * * `Apply` (default: `false`): enable the correction. Requires full + * `sim::SimPhotons` input (photon positions). + * * `BinEdges` (default: empty): distance bin edges [cm], sorted; must have + * exactly one more entry than `Factors`. + * * `Factors` (default: empty): survival probability per bin, each in + * `[0, 1]`; survival is 1 for distances outside the covered range. + * * `TailSuppression` (table, optional): empirical tail suppression correction. + * Sub-parameters: + * * `Apply` (default: `false`): enable the correction. + * * `Epsilon` (default: `0.0`): correction strength @f$ \varepsilon @f$ + * (dimensionless); typical values 0.27 (Run 1/2). + * * `Tau` (default: `250.0`): charge-state decay time constant [ns]. + * + * * Random number generators * ------------------------- * @@ -505,6 +605,23 @@ class icarus::opdet::PMTsimulationAlg { }; // struct PMTspecs_t + struct TailSuppressionParams_t { + bool apply = false; + float epsilon = 0.0f; ///< correction strength [dimensionless] + float tau = 250.0f; ///< charge-state decay time constant [ns] + }; + + /// Distance-dependent photon survival: each scintillation + /// photon is kept with probability `QE x S(d)`, where `d` is the + /// straight-line distance from the photon emission point to the PMT + /// center. `factors[i]` applies to `d` in `[binEdges[i], binEdges[i+1])` + /// [cm]; survival is 1 outside the covered range. + struct DistanceSurvivalParams_t { + bool apply = false; + std::vector binEdges; ///< bin edges [cm], size = factors + 1 + std::vector factors; ///< survival probabilities, in [0, 1] + }; + /// @{ /// @name High level configuration parameters. @@ -531,6 +648,8 @@ class icarus::opdet::PMTsimulationAlg { hertz darkNoiseRate; float saturation; //equivalent to the number of p.e. that saturates the electronic signal PMTspecs_t PMTspecs; ///< PMT specifications. + TailSuppressionParams_t tailSuppression; ///< Causal subtractive droop correction parameters. + DistanceSurvivalParams_t distanceSurvival; ///< Distance-dependent photon survival S(d). bool doGainFluctuations; ///< Whether to simulate gain fluctuations. bool useGainCalibDB; ///< Whether to use per-channel DB gain for fluctuations. bool doTimingDelays; ///< Whether to simulate timing delays. @@ -544,6 +663,9 @@ class icarus::opdet::PMTsimulationAlg { detinfo::LArProperties const* larProp = nullptr; ///< LarProperties service provider. + /// OpDet geometry mapping (required if `distanceSurvival.apply`). + geo::WireReadoutGeom const* wireReadoutGeom = nullptr; + detinfo::DetectorClocksData const* clockData = nullptr; /// Timing delays service interfacing with database @@ -637,6 +759,10 @@ class icarus::opdet::PMTsimulationAlg { template void printConfiguration(Stream&& out, std::string indent = "") const; + /// Whether the distance-dependent photon survival S(d) is enabled + /// (requires full `sim::SimPhotons` input: photon positions). + bool appliesDistanceSurvival() const + { return fParams.distanceSurvival.apply; } /// Manages the conversion between names and values of `DiscriminationAlgo`. static util::MultipleChoiceSelection const @@ -971,6 +1097,10 @@ class icarus::opdet::PMTsimulationAlg { /// Returns a random response whether a photon generates a photoelectron. bool KicksPhotoelectron() const; + + /// Returns the distance survival probability S(d) for a photon travelling + /// a straight-line distance `d_cm` [cm] (1 outside the configured range). + double distanceSurvival(double d_cm) const; /// Returns the ADC range allowed for photoelectron saturation. std::pair saturationRange(ADCcount baseline) const; @@ -984,7 +1114,10 @@ class icarus::opdet::PMTsimulationAlg { /// Applies the configured photoelectron saturation on the `waveform`. void ApplySaturation(Waveform_t& waveform, ADCcount baseline) const; - + + /// Applies the causal subtractive tail suppression correction to the signal waveform. + void ApplyTailSuppression(Waveform_t& waveform) const; + /// Forces `waveform` ADC within the `min` to `max` range (`max` included). static void ClipWaveform(Waveform_t& waveform, ADCcount min, ADCcount max); @@ -1034,7 +1167,51 @@ class icarus::opdet::PMTsimulationAlgMaker { }; // struct PMTspecConfig - + struct TailSuppressionConfig { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom Apply { + Name("Apply"), + Comment("Enable causal subtractive tail suppression correction"), + false + }; + fhicl::Atom Epsilon { + Name("Epsilon"), + Comment("Correction strength [dimensionless]"), + 0.0f + }; + fhicl::Atom Tau { + Name("Tau"), + Comment("Charge-state decay time constant [ns]"), + 250.0f + }; + + }; // struct TailSuppressionConfig + + struct DistanceSurvivalConfig { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom Apply { + Name("Apply"), + Comment("Enable distance-dependent photon survival S(d)"), + false + }; + fhicl::Sequence BinEdges { + Name("BinEdges"), + Comment("Photon-to-PMT distance bin edges [cm]; one more than Factors"), + std::vector{} + }; + fhicl::Sequence Factors { + Name("Factors"), + Comment("Survival probability per distance bin (S = 1 outside range)"), + std::vector{} + }; + + }; // struct DistanceSurvivalConfig + + /// Main algorithm FHiCL configuration. struct Config { using Name = fhicl::Name; @@ -1112,6 +1289,22 @@ class icarus::opdet::PMTsimulationAlgMaker { 1U }; + // + // tail suppression + // + fhicl::OptionalTable TailSuppression { + Name("TailSuppression"), + Comment("PMT anode baseline-droop correction (causal subtractive model)") + }; + + // + // distance-dependent photon survival + // + fhicl::OptionalTable DistanceSurvival { + Name("DistanceSurvival"), + Comment("Distance-dependent scintillation photon survival S(d)") + }; + // // dark noise // @@ -1186,6 +1379,7 @@ class icarus::opdet::PMTsimulationAlgMaker { std::unique_ptr operator()( std::uint64_t beamGateTimestamp, detinfo::LArProperties const& larProp, + geo::WireReadoutGeom const& wireReadoutGeom, detinfo::DetectorClocksData const& detClocks, icarusDB::PMTTimingCorrections const* timingDelays, icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, @@ -1219,6 +1413,7 @@ class icarus::opdet::PMTsimulationAlgMaker { PMTsimulationAlg::ConfigurationParameters_t makeParams( std::uint64_t beamGateTimestamp, detinfo::LArProperties const& larProp, + geo::WireReadoutGeom const& wireReadoutGeom, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, @@ -1320,7 +1515,21 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration << '\n' << indent << "Gain at first stage: " << fParams.PMTspecs.firstStageGain() << '\n' << indent << "SPR nominal area: " << fNominalSPEArea << " ADC x tick" << '\n' << indent << "Bias ratio: " << fBiasConstant + << '\n' << indent << "Tail suppression: " + << std::boolalpha << fParams.tailSuppression.apply ; + if (fParams.tailSuppression.apply) { + out << " (epsilon=" << fParams.tailSuppression.epsilon + << ", tau=" << fParams.tailSuppression.tau << " ns)"; + } + out + << '\n' << indent << "Distance survival: " + << std::boolalpha << fParams.distanceSurvival.apply + ; + if (fParams.distanceSurvival.apply) { + out << " (" << fParams.distanceSurvival.factors.size() << " bins over " + << fParams.distanceSurvival.binEdges.size() << " edges)"; + } out << '\n' << indent << "Pedestal: " << fPedestalGen->toString(indent + " ", ""); @@ -1335,7 +1544,6 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration out << '\n' << indent << "Template photoelectron waveform settings:" << '\n'; wsp.dump(std::forward(out), indent + " "); - out << '\n' << indent << "Track used photons: " << std::boolalpha << fParams.trackSelectedPhotons << '\n'; diff --git a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl index dfd6948bc..aab91269f 100644 --- a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl +++ b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl @@ -44,6 +44,10 @@ # 20251111 (mvicenzi@bnl.gov) # introduced configuration for new Run1/2 SPR (SPRisolHitRun9271) # introduced configuration for new Run3/4 SPR (SPRisolHitRun12715) +# 20260618 (mvicenzi@bnl.gov) +# introduced tail suppression tune for run2 +# 20260722 (mvicenzi@bnl.gov) +# introduced distance-based survival probability for run2 #include "opticalproperties_icarus.fcl" @@ -402,6 +406,13 @@ icarus_pmtsimulationalg_202202_noise: { Gain: @local::icarus_settings_opdet.nominalPMTgain # total typical PMT gain } # PMTspecs + # Subtractive tail suppression correction + TailSuppression: { + Apply: false + Epsilon: 0.0 + Tau: 250.0 # [ns] + } + } # icarus_pmtsimulationalg_202202_noise @@ -469,6 +480,23 @@ icarus_pmtsimulationalg_run2_2025: { Gain: @local::icarus_settings_opdet_run2.nominalPMTgain # total PMT gain } + # Tail suppression effective correction + # Run2 tune from SBN-doc-48511 + TailSuppression: { + Apply: true + Epsilon: 0.27 # correction strength + Tau: 500.0 # [ns] + } + + # Distance-dependent scintillation photon survival S(d): each photon kept + # with probability QE x S(d), S = 1 outside range + # Run2 tune from SBN-doc-48511 + DistanceSurvival: { + Apply: true + BinEdges: [ 150, 200, 250, 300, 350, 400, 450 ] # [cm] + Factors: [ 0.97, 0.87, 0.82, 0.82, 0.82, 0.85 ] + } + } # icarus_pmtsimulationalg_run2_2025 @@ -504,6 +532,20 @@ icarus_pmtsimulationalg_run4_2025: { Gain: @local::icarus_settings_opdet_run4.nominalPMTgain # total PMT gain } + # TODO: requires tuning on Run-4 data + TailSuppression: { + Apply: false + Epsilon: 0.27 + Tau: 500.0 # [ns] + } + + # TODO: requires tuning on Run-4 data + DistanceSurvival: { + Apply: false + BinEdges: [ 150, 200, 250, 300, 350, 400, 450 ] # [cm] + Factors: [ 0.97, 0.87, 0.82, 0.82, 0.82, 0.85 ] + } + } # icarus_pmtsimulationalg_run4_2025 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index dce586625..633c1d4d4 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -23,6 +23,7 @@ // LArSoft libraries #include "larcore/CoreUtils/ServiceUtil.h" +#include "larcore/Geometry/WireReadout.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "lardata/DetectorInfoServices/LArPropertiesService.h" #include "lardataalg/DetectorInfo/DetectorTimings.h" @@ -490,6 +491,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) auto PMTsimulator = makePMTsimulator( e.time().value(), // using the event generation time as beam time stamp *(lar::providerFrom()), + art::ServiceHandle()->Get(), clockData, fDoTimingDelays ? lar::providerFrom() : nullptr, fUseGainCalibDB ? art::ServiceHandle()->provider() : nullptr, @@ -515,6 +517,16 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) // Prefer SimPhotons if available. // Make sure that there are parallel inputs for both formats; bool const useLitePhotons = !pmtVector.isValid(); + + // Distance-dependent survival needs photon positions: never silently + // skip the correction on the position-less SimPhotonsLite path. + if (useLitePhotons && PMTsimulator->appliesDistanceSurvival()) { + throw art::Exception(art::errors::Configuration) + << "DistanceSurvival is enabled but only sim::SimPhotonsLite input" + " is available ('" << fInputModuleName.encode() << "'):" + " photon positions are required. Run the photon propagation with" + " full sim::SimPhotons.\n"; + } // storage for the photons that are not used (but still required) std::vector fakePhotons;