56 #include "cetlib_except/exception.h" 66 #include "CLHEP/Random/RandPoissonQ.h" 68 #include "range/v3/view/enumerate.hpp" 87 Comment(
"SimEnergyDeposit label.")};
89 Comment(
"Simulate slow scintillation light, default true"),
92 Comment(
"Simulate slow scintillation light")};
94 Comment(
"Simulate reflected visible light")};
96 Name(
"IncludeAnodeReflections"),
97 Comment(
"Simulate anode reflections, default false"),
100 Comment(
"Simulate light propagation time")};
102 Name(
"GeoPropTimeOnly"),
103 Comment(
"Simulate light propagation time geometric approximation, default false"),
106 Name(
"UseLitePhotons"),
107 Comment(
"Store SimPhotonsLite/OpDetBTRs instead of SimPhotons")};
109 Comment(
"Photons cannot cross the cathode")};
111 Name(
"OnlyActiveVolume"),
112 Comment(
"PAR fast sim usually only for active volume, default true"),
115 Comment(
"Set to true if light is only supported in C:1")};
117 Comment(
"Tool describing scintillation time structure")};
119 Name(
"UseXeAbsorption"),
120 Comment(
"Use Xe absorption length instead of Ar, default false"),
124 Comment(
"Configuration for visible timing parameterization")};
127 Comment(
"Configuration for visibile visibility parameterization")};
136 const std::vector<double>& OpDetVisibilities,
137 const int NumPhotons)
const;
139 void AddOpDetBTR(std::vector<sim::OpDetBacktrackerRecord>& opbtr,
140 std::vector<int>& ChannelMap,
199 ,
fGeom(*(lar::providerFrom<geo::Geometry>()))
205 ,
fSimTag(config().SimulationLabel())
218 mf::LogInfo(
"PDFastSimPAR") <<
"Initializing PDFastSimPAR." << std::endl;
230 <<
"Propagation time simulation requested, but VUVTiming not specified." 236 <<
"Reflected light or anode reflections simulation requested, but VisHits not specified." 242 <<
"Reflected light propagation time simulation requested, but VISTiming not specified." 249 <<
" cryostats is configured" 250 <<
" , and semi-analytic model is requested for scintillation photon propagation." 251 <<
" THIS CONFIGURATION IS NOT SUPPORTED and it is open to bugs" 252 <<
" (e.g. scintillation may be detected only in cryostat #0)." 253 <<
"\nThis would be normally a fatal error, but it has been forcibly overridden." 255 << std::string(80,
'=');
259 <<
"Photon propagation via semi-analytic model is not supported yet" 260 <<
" on detectors with more than one cryostat.";
277 auto log =
mf::LogTrace(
"PDFastSimPAR") <<
"PDFastSimPAR: active volume boundaries from " 280 log <<
"\n - C:" << iCryo <<
": " << box.Min() <<
" -- " << box.Max() <<
" cm";
285 mf::LogInfo(
"PDFastSimPAR") <<
"Using Lite Photons";
286 produces<std::vector<sim::SimPhotonsLite>>();
287 produces<std::vector<sim::OpDetBacktrackerRecord>>();
288 if (fDoReflectedLight) {
289 mf::LogInfo(
"PDFastSimPAR") <<
"Storing Reflected Photons";
290 produces<std::vector<sim::SimPhotonsLite>>(
"Reflected");
291 produces<std::vector<sim::OpDetBacktrackerRecord>>(
"Reflected");
295 mf::LogInfo(
"PDFastSimPAR") <<
"Using Sim Photons";
296 produces<std::vector<sim::SimPhotons>>();
297 if (fDoReflectedLight) {
298 mf::LogInfo(
"PDFastSimPAR") <<
"Storing Reflected Photons";
299 produces<std::vector<sim::SimPhotons>>(
"Reflected");
306 if (
fNTPC > 1 && fDriftDistance < 50)
309 mf::LogInfo(
"PDFastSimPAR") <<
"PDFastSimPAR Initialization finish.\n" 310 <<
"Simulate using semi-analytic model for number of hits." 317 mf::LogTrace(
"PDFastSimPAR") <<
"PDFastSimPAR Module Producer" 318 <<
"EventID: " <<
event.event();
320 std::vector<int> PDChannelToSOCMapDirect(
fNOpChannels, -1);
321 std::vector<int> PDChannelToSOCMapReflect(
fNOpChannels, -1);
324 auto phlit = std::make_unique<std::vector<sim::SimPhotonsLite>>();
325 auto opbtr = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
326 auto phlit_ref = std::make_unique<std::vector<sim::SimPhotonsLite>>();
327 auto opbtr_ref = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
328 auto& dir_phlitcol(*phlit);
329 auto& ref_phlitcol(*phlit_ref);
331 auto phot = std::make_unique<std::vector<sim::SimPhotons>>();
332 auto phot_ref = std::make_unique<std::vector<sim::SimPhotons>>();
333 auto& dir_photcol(*
phot);
334 auto& ref_photcol(*phot_ref);
339 dir_phlitcol[i].OpChannel = i;
340 ref_phlitcol[i].OpChannel = i;
347 dir_photcol[i].fOpChannel = i;
348 ref_photcol[i].fOpChannel = i;
358 auto const& edeps = edepHandle;
366 for (
auto const& edepi : *edeps) {
369 int nphot_fast = edepi.NumFPhotons();
370 int nphot_slow = edepi.NumSPhotons();
372 num_fastph += nphot_fast;
373 num_slowph += nphot_slow;
377 int trackID = edepi.TrackID();
378 int nphot = edepi.NumPhotons();
379 double edeposit = edepi.Energy() / nphot;
380 double pos[3] = {edepi.MidPointX(), edepi.MidPointY(), edepi.MidPointZ()};
381 geo::Point_t const ScintPoint = {pos[0], pos[1], pos[2]};
389 std::vector<double> OpDetVisibilities;
390 fVisibilityModel->detectedDirectVisibilities(OpDetVisibilities, ScintPoint);
398 std::vector<double> OpDetVisibilitiesAnode;
399 fVisibilityModel->detectedReflectedVisibilities(OpDetVisibilitiesAnode, ScintPoint,
true);
404 for (
size_t i = 0; i < AnodeDetectedNumFast.size(); ++i) {
405 DetectedNumFast[i] += AnodeDetectedNumFast[i];
407 for (
size_t i = 0; i < AnodeDetectedNumSlow.size(); ++i) {
408 DetectedNumSlow[i] += AnodeDetectedNumSlow[i];
416 std::vector<double> OpDetVisibilitiesRefl;
417 fVisibilityModel->detectedReflectedVisibilities(OpDetVisibilitiesRefl, ScintPoint,
false);
424 for (
size_t Reflected = 0; Reflected <= DoReflected; ++Reflected) {
425 for (
size_t channel = 0; channel <
fNOpChannels; channel++) {
429 int ndetected_fast = DetectedNumFast[channel];
430 int ndetected_slow = DetectedNumSlow[channel];
432 ndetected_fast = ReflDetectedNumFast[channel];
433 ndetected_slow = ReflDetectedNumSlow[channel];
440 std::vector<double> transport_time;
442 transport_time.resize(ndetected_fast + ndetected_slow);
443 fPropTimeModel->propagationTime(transport_time, ScintPoint, channel, Reflected);
450 int n = ndetected_fast;
452 for (
int i = 0; i <
n; ++i) {
454 double dtime = edepi.StartT() +
fScintTime->fastScintTime();
456 int time =
static_cast<int>(std::round(dtime));
458 ++ref_phlitcol[channel].DetectedPhotons[time];
460 ++dir_phlitcol[channel].DetectedPhotons[time];
465 int n = ndetected_slow;
467 for (
int i = 0; i <
n; ++i) {
468 double dtime = edepi.StartT() +
fScintTime->slowScintTime();
470 int time =
static_cast<int>(std::round(dtime));
472 ++ref_phlitcol[channel].DetectedPhotons[time];
474 ++dir_phlitcol[channel].DetectedPhotons[time];
479 AddOpDetBTR(*opbtr_ref, PDChannelToSOCMapReflect, tmpbtr);
481 AddOpDetBTR(*opbtr, PDChannelToSOCMapDirect, tmpbtr);
490 photon.
Energy = 2.9 * CLHEP::eV;
492 photon.
Energy = 9.7 * CLHEP::eV;
495 int n = ndetected_fast;
497 for (
int i = 0; i <
n; ++i) {
499 double dtime = edepi.StartT() +
fScintTime->fastScintTime();
501 int time =
static_cast<int>(std::round(dtime));
504 ref_photcol[channel].insert(ref_photcol[channel].
end(), 1, photon);
506 dir_photcol[channel].insert(dir_photcol[channel].
end(), 1, photon);
510 int n = ndetected_slow;
512 for (
int i = 0; i <
n; ++i) {
513 double dtime = edepi.StartT() +
fScintTime->slowScintTime();
515 int time =
static_cast<int>(std::round(dtime));
518 ref_photcol[channel].insert(ref_photcol[channel].
end(), 1, photon);
520 dir_photcol[channel].insert(dir_photcol[channel].
end(), 1, photon);
528 mf::LogTrace(
"PDFastSimPAR") <<
"Total points: " << num_points
529 <<
", total fast photons: " << num_fastph
530 <<
", total slow photons: " << num_slowph
531 <<
"\ndetected fast photons: " << num_fastdp
532 <<
", detected slow photons: " << num_slowdp;
535 event.put(move(phlit));
536 event.put(move(opbtr));
538 event.put(move(phlit_ref),
"Reflected");
539 event.put(move(opbtr_ref),
"Reflected");
543 event.put(move(
phot));
552 std::vector<int>& ChannelMap,
556 if (ChannelMap[iChan] < 0) {
557 ChannelMap[iChan] = opbtr.size();
558 opbtr.emplace_back(std::move(btr));
561 size_t idtest = ChannelMap[iChan];
563 for (
auto const& timePDclockSDP : timePDclockSDPsMap) {
564 for (
auto const& sdp : timePDclockSDP.second) {
565 double xyz[3] = {sdp.x, sdp.y, sdp.z};
566 opbtr.at(idtest).AddScintillationPhotons(
567 sdp.trackID, timePDclockSDP.first, sdp.numPhotons, xyz, sdp.energy);
576 const std::vector<double>& OpDetVisibilities,
577 const int NumPhotons)
const 579 for (
size_t i = 0; i < OpDetVisibilities.size(); ++i) {
580 DetectedNumPhotons[i] =
fRandPoissPhot->fire(OpDetVisibilities[i] * NumPhotons);
595 if (
fNTPC == 2 && ((ScintPoint.X() < 0.) != (OpDetPoint.X() < 0.)) &&
612 std::vector<geo::Point_t> opDetCenter;
613 for (
size_t const i : ::ranges::views::ints(
size_t(0),
fNOpChannels)) {
615 opDetCenter.push_back(opDet.
GetCenter());
const bool fGeoPropTimeOnly
const std::vector< geo::BoxBoundedGeo > fActiveVolumes
base_engine_t & createEngine(seed_t seed)
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
const size_t fNOpChannels
void AddOpDetBTR(std::vector< sim::OpDetBacktrackerRecord > &opbtr, std::vector< int > &ChannelMap, const sim::OpDetBacktrackerRecord &btr) const
static std::vector< geo::BoxBoundedGeo > extractActiveLArVolume(geo::GeometryCore const &geom)
Utilities related to art service access.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const bool fUseLitePhotons
std::unique_ptr< PropagationTimeModel > fPropTimeModel
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
fhicl::Atom< bool > DoSlowComponent
Point_t const & GetCenter() const
EDProducer(fhicl::ParameterSet const &pset)
PDFastSimPAR(Parameters const &config)
All information of a photon entering the sensitive optical detector volume.
fhicl::Atom< bool > DoReflectedLight
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::unique_ptr< ScintTime > fScintTime
const bool fOnlyActiveVolume
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
const bool fIncludePropTime
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
const bool fUseXeAbsorption
const bool fOpaqueCathode
Energy deposited on a readout Optical Detector by simulated tracks.
geo::Point_t InitialPosition
Scintillation position in world coordinates [cm].
const bool fDoReflectedLight
fhicl::Atom< bool > UseLitePhotons
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
std::unique_ptr< CLHEP::RandPoissonQ > fRandPoissPhot
CLHEP::HepRandomEngine & fScintTimeEngine
Simulation objects for optical detectors.
int OpDetNum() const
Returns the readout Optical Detector this object describes.
const bool fDoFastComponent
#define DEFINE_ART_MODULE(klass)
const std::vector< geo::Point_t > fOpDetCenter
Definitions of geometry vector data types.
std::unique_ptr< SemiAnalyticalModel > fVisibilityModel
fhicl::Atom< bool > OpaqueCathode
T get(std::string const &key) const
fhicl::Atom< bool > IncludeAnodeReflections
void produce(art::Event &) override
The data type to uniquely identify a TPC.
Description of the physical geometry of one entire detector.
const art::InputTag fSimTag
const bool fDoSlowComponent
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
unsigned int NOpDets() const
Number of OpDets in the whole detector.
Provides a base class aware of world box coordinates.
void AddScintillationPhotons(TrackID_t trackID, timePDclock_t timePDclock, double numberPhotons, double const *xyz, double energy)
Add scintillation photons and energy to this OpticalDetector.
Encapsulate the geometry of an optical detector.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
fhicl::Atom< bool > UseXeAbsorption
General LArSoft Utilities.
fhicl::Atom< bool > GeoPropTimeOnly
CLHEP::HepRandomEngine & fPhotonEngine
double DriftDistance() const
Drift distance is defined as the distance between the anode and the cathode, in centimeters.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
std::vector< geo::Point_t > opDetCenters() const
bool SetInSD
Whether the photon reaches the sensitive detector.
contains information for a single step in the detector simulation
void detectedNumPhotons(std::vector< int > &DetectedNumPhotons, const std::vector< double > &OpDetVisibilities, const int NumPhotons) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
fhicl::Atom< art::InputTag > SimulationLabel
MaybeLogger_< ELseverityLevel::ELsev_success, true > LogTrace
fhicl::Atom< bool > IncludePropTime
float Energy
Scintillation photon energy [GeV].
fhicl::Atom< bool > OnlyActiveVolume
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
bool isScintInActiveVolume(geo::Point_t const &ScintPoint)
const bool fIncludeAnodeReflections
timePDclockSDPs_t const & timePDclockSDPsMap() const
Returns all the deposited energy information as stored.
fhicl::Atom< bool > OnlyOneCryostat
art framework interface to geometry description
geo::GeometryCore const & fGeom
fhicl::Atom< bool > DoFastComponent
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
Event finding and building.
int MotherTrackID
ID of the GEANT4 track causing the scintillation.
const bool fOnlyOneCryostat