LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
evgen::RadioGen Class Reference

Module to generate particles created by radiological decay, patterned off of SingleGen. More...

Inheritance diagram for evgen::RadioGen:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 RadioGen (fhicl::ParameterSet const &pset)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
void fillProductDescriptions ()
 
void registerProducts (ProductDescriptions &productsToRegister)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
std::unique_ptr< Worker > makeWorker (WorkerParams const &wp)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Protected Member Functions

ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Private Types

typedef int ti_PDGID
 
typedef double td_Mass
 

Private Member Functions

void produce (art::Event &evt)
 
void beginRun (art::Run &run)
 
std::size_t addvolume (std::string const &volumeName)
 
void SampleOne (unsigned int i, simb::MCTruth &mct)
 Checks that the node represents a box well aligned to world frame axes. More...
 
TLorentzVector dirCalc (double p, double m)
 
void readfile (std::string nuclide, std::string const &filename, std::string const &PathToFile)
 
void samplespectrum (std::string nuclide, int &itype, double &t, double &m, double &p)
 
void Ar42Gamma2 (std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
 
void Ar42Gamma3 (std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
 
void Ar42Gamma4 (std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
 
void Ar42Gamma5 (std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
 
double samplefromth1d (TH1D &hist)
 
template<typename Stream >
void dumpNuclideSettings (Stream &&out, std::size_t iNucl, std::string const &volumeName={}) const
 Prints the settings for the specified nuclide and volume. More...
 
std::pair< double, double > defaulttimewindow () const
 Returns the start and end of the readout window. More...
 

Private Attributes

std::vector< std::string > fNuclide
 List of nuclides to simulate. Example: "39Ar". More...
 
std::vector< std::string > fMaterial
 List of regexes of materials in which to generate the decays. Example: "LAr". More...
 
std::vector< double > fBq
 Radioactivity in Becquerels (decay per sec) per cubic cm. More...
 
std::vector< double > fT0
 Beginning of time window to simulate in ns. More...
 
std::vector< double > fT1
 End of time window to simulate in ns. More...
 
std::vector< double > fX0
 Bottom corner x position (cm) in world coordinates. More...
 
std::vector< double > fY0
 Bottom corner y position (cm) in world coordinates. More...
 
std::vector< double > fZ0
 Bottom corner z position (cm) in world coordinates. More...
 
std::vector< double > fX1
 Top corner x position (cm) in world coordinates. More...
 
std::vector< double > fY1
 Top corner y position (cm) in world coordinates. More...
 
std::vector< double > fZ1
 Top corner z position (cm) in world coordinates. More...
 
bool fIsFirstSignalSpecial
 
int trackidcounter
 Serial number for the MC track ID. More...
 
std::string fPathToFile
 
std::vector< std::string > spectrumname
 
std::vector< std::unique_ptr< TH1D > > alphaspectrum
 
std::vector< double > alphaintegral
 
std::vector< std::unique_ptr< TH1D > > betaspectrum
 
std::vector< double > betaintegral
 
std::vector< std::unique_ptr< TH1D > > gammaspectrum
 
std::vector< double > gammaintegral
 
std::vector< std::unique_ptr< TH1D > > neutronspectrum
 
std::vector< double > neutronintegral
 
CLHEP::HepRandomEngine & fEngine
 

Detailed Description

Module to generate particles created by radiological decay, patterned off of SingleGen.

The module generates the products of radioactive decay of some known nuclides. Each nuclide decay producing a single particle (with exceptions, as for example argon(A=42)), whose spectrum is specified in a ROOT file to be found in the FW_SEARCH_PATH path, and it can be a alpha, beta, gamma or neutron. In case multiple decay channels are possible, each decay is stochastically chosen weighting the channels according to the integral of their spectrum. The normalization of the spectrum is otherwise ignored.

Nuclides can be added by making the proper distributions available in a file called after the name used for it in the Nuclide configuration parameter (e.g 14C.root if the nuclide key is 14C): check the existing nuclide files for examples.

A special treatment is encoded for argon(A=42) and radon(A=222) (and, "temporary", for nichel(A=59) ).

Decays happen only in volumes specified in the configuration, and with a rate also specified in the configuration. Volumes are always box-shaped, and can be specified by coordinates or by name. In addition, within each volume decays will be generated only in the subvolumes matching the specified materials.

Note
All beta decays emit positrons.

Configuration parameters

  • X0, Y0, Z0, X1, Y1, Z1 (lists of real numbers, optional): if specified, they describe the box volumes where to generate decays; all lists must have the same size, and each entry i defines the box between coordinates (X0[i], Y0[i], Z0[i]) and (X1[i], Y1[i], Z1[i]) expressed in world coordinates and in centimeters;
  • Volumes (list of strings, optional): if specified, each entry represents all the volumes in the geometry whose name exactly matches the entry; the volumes are added after the ones explicitly listed by their coordinates (configuration parameters X0, Y0, Z0, X1, Y1 and Z1); each volume name counts as one entry, even if it expands to multiple volumes; Volumes can be omitted if volumes are specified with the coordinate parameters;
  • Nuclide (list of strings, mandatory): the list of decaying nuclides, one per volume entry; note that each of the elements in X0 (and the other 5 coordinate parameters) and each of the elements in Volumes counts as one entry, even for entries of the latter which match multiple volumes (in that case, all matching volumes are assigned the same nuclide parameters); since documentation is never maintained, refer to the code for a list of the supported materials; a subset of them is: 39Ar, 60Co, 85Kr, 40K, 232Th, 238U, 222Rn, 59Ni and 42Ar;
  • Material (list of regular expressions, mandatory): for each nuclide, the name of the materials allowed to decay in this mode; this name is a regular expression (C++ default std::regex), which needs to match the name of the material as specified in the detector geometry (usually in GDML format); a material is mandatory for each nuclide; if no material selection is desired, the all-matching pattern ".*" can be used;
  • BqPercc (list of real numbers, mandatory): activity of the nuclides, in the same order as they are in Nuclide parameter; each value is specified as becquerel per cubic centimeter;
  • T0, T1 (list of real values, optional): range of time when the decay may happen, in simulation time; each nuclide is assigned a time range; differently from the other parameters, the list of times can be shorter than the number of nuclides, in which case all the nuclides with no matching time range will be assigned the last range specified; as a specially special case, if no time is specified at all, the time window is assigned as from one readout window (detinfo::DetectorPropertiesData::ReadOutWindowSize()) before the nominal hardware trigger time (detinfo::DetectorClocksData::TriggerTime()) up to a number of ticks after the trigger time equivalent to the full simulated TPC waveform (detinfo::DetectorPropertiesData::NumberTimeSamples()); this makes it a quite poor default, so you may want to avoid it.

Definition at line 194 of file RadioGen_module.cc.

Member Typedef Documentation

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 26 of file Producer.h.

typedef double evgen::RadioGen::td_Mass
private

Definition at line 206 of file RadioGen_module.cc.

typedef int evgen::RadioGen::ti_PDGID
private

Definition at line 204 of file RadioGen_module.cc.

Constructor & Destructor Documentation

evgen::RadioGen::RadioGen ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 289 of file RadioGen_module.cc.

References addvolume(), art::errors::Configuration, util::counter(), art::detail::EngineCreator::createEngine(), defaulttimewindow(), dumpNuclideSettings(), util::enumerate(), fBq, fEngine, fIsFirstSignalSpecial, fMaterial, fNuclide, fPathToFile, fT0, fT1, fX0, fX1, fY0, fY1, fZ0, fZ1, readfile(), and t1.

290  : EDProducer{pset}
291  , fX0{pset.get<std::vector<double>>("X0", {})}
292  , fY0{pset.get<std::vector<double>>("Y0", {})}
293  , fZ0{pset.get<std::vector<double>>("Z0", {})}
294  , fX1{pset.get<std::vector<double>>("X1", {})}
295  , fY1{pset.get<std::vector<double>>("Y1", {})}
296  , fZ1{pset.get<std::vector<double>>("Z1", {})}
297  , fIsFirstSignalSpecial{pset.get<bool>("IsFirstSignalSpecial", false)}
298  // create a default random engine; obtain the random seed from NuRandomService,
299  // unless overridden in configuration with key "Seed"
300  , fEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(createEngine(0),
301  pset,
302  "Seed"))
303  {
304  produces<std::vector<simb::MCTruth>>();
305  produces<sumdata::RunData, art::InRun>();
306 
307  auto const fPathToFile = pset.get<std::string>("PathToFile", "Radionuclides/");
308  auto const nuclide = pset.get<std::vector<std::string>>("Nuclide");
309  auto const material = pset.get<std::vector<std::string>>("Material");
310  auto const Bq = pset.get<std::vector<double>>("BqPercc");
311  auto t0 = pset.get<std::vector<double>>("T0", {});
312  auto t1 = pset.get<std::vector<double>>("T1", {});
313 
314  if (fPathToFile.find_last_of('/') !=
315  fPathToFile.length() - 1) { //Checking that last character in path is a forward slash
316  throw cet::exception("RadioGen")
317  << "Last entry of the path to radiological files needs to be a forward slash ('/')\n";
318  }
319 
320  if (fT0.empty() || fT1.empty()) { // better be both empty...
321  if (!fT0.empty() || !fT1.empty()) {
323  << "RadioGen T0 and T1 need to be both non-empty, or both empty"
324  " (now T0 has "
325  << fT0.size() << " entries and T1 has " << fT0.size() << ")\n";
326  }
327  auto const [defaultT0, defaultT1] = defaulttimewindow();
328  t0.push_back(defaultT0);
329  t1.push_back(defaultT1);
330  }
331 
332  //
333  // First, the materials are assigned to the coordinates explicitly
334  // configured; then, each of the remaining materials is assigned to one
335  // of the named volumes.
336  // TODO: reaction to mismatches is "not so great".
337  //
338  mf::LogInfo("RadioGen") << "Configuring activity:";
339  for (std::size_t iCoord : util::counter(fX0.size())) {
340 
341  fNuclide.push_back(nuclide.at(iCoord));
342  fMaterial.push_back(material.at(iCoord));
343  fBq.push_back(Bq.at(iCoord));
344  // replicate the last timing if none is specified
345  fT0.push_back(t0[std::min(iCoord, t0.size() - 1U)]);
346  fT1.push_back(t1[std::min(iCoord, t1.size() - 1U)]);
347 
348  dumpNuclideSettings(mf::LogVerbatim("RadioGen"), iCoord);
349 
350  } // for
351 
352  std::size_t iNext = fX0.size();
353  auto const volumes = pset.get<std::vector<std::string>>("Volumes", {});
354  for (auto&& [iVolume, volName] : volumes | ranges::views::enumerate) {
355  // this expands the coordinate vectors
356  auto const nVolumes = addvolume(volName);
357  if (nVolumes == 0) {
359  << "No volume named '" << volName << "' was found!\n";
360  }
361 
362  std::size_t const iVolFirst = fNuclide.size();
363  for (auto iVolInstance : util::counter(nVolumes)) {
364  fNuclide.push_back(nuclide.at(iNext));
365  fMaterial.push_back(material.at(iNext));
366  fBq.push_back(Bq.at(iNext));
367  // replicate the last timing if none is specified
368  fT0.push_back(t0[std::min(iNext, t0.size() - 1U)]);
369  fT1.push_back(t1[std::min(iNext, t1.size() - 1U)]);
370 
371  dumpNuclideSettings(mf::LogVerbatim("RadioGen"), iVolFirst + iVolInstance, volName);
372  }
373  ++iNext;
374  } // for
375 
376  // check for consistency of vector sizes
377  unsigned int nsize = fNuclide.size();
378  if (nsize == 0) {
379  throw art::Exception(art::errors::Configuration) << "No nuclide configured!\n";
380  }
381 
382  if (fMaterial.size() != nsize)
383  throw cet::exception("RadioGen") << "Different size Material vector and Nuclide vector\n";
384  if (fBq.size() != nsize)
385  throw cet::exception("RadioGen") << "Different size Bq vector and Nuclide vector\n";
386  if (fT0.size() != nsize)
387  throw cet::exception("RadioGen") << "Different size T0 vector and Nuclide vector\n";
388  if (fT1.size() != nsize)
389  throw cet::exception("RadioGen") << "Different size T1 vector and Nuclide vector\n";
390  if (fX0.size() != nsize)
391  throw cet::exception("RadioGen") << "Different size X0 vector and Nuclide vector\n";
392  if (fY0.size() != nsize)
393  throw cet::exception("RadioGen") << "Different size Y0 vector and Nuclide vector\n";
394  if (fZ0.size() != nsize)
395  throw cet::exception("RadioGen") << "Different size Z0 vector and Nuclide vector\n";
396  if (fX1.size() != nsize)
397  throw cet::exception("RadioGen") << "Different size X1 vector and Nuclide vector\n";
398  if (fY1.size() != nsize)
399  throw cet::exception("RadioGen") << "Different size Y1 vector and Nuclide vector\n";
400  if (fZ1.size() != nsize)
401  throw cet::exception("RadioGen") << "Different size Z1 vector and Nuclide vector\n";
402 
403  for (std::string& nuclideName : fNuclide) {
404  if (nuclideName == "39Ar") { readfile("39Ar", "Argon_39.root", fPathToFile); }
405  else if (nuclideName == "60Co") {
406  readfile("60Co", "Cobalt_60.root", fPathToFile);
407  }
408  else if (nuclideName == "85Kr") {
409  readfile("85Kr", "Krypton_85.root", fPathToFile);
410  }
411  else if (nuclideName == "40K") {
412  readfile("40K", "Potassium_40.root", fPathToFile);
413  }
414  else if (nuclideName == "232Th") {
415  readfile("232Th", "Thorium_232.root", fPathToFile);
416  }
417  else if (nuclideName == "238U") {
418  readfile("238U", "Uranium_238.root", fPathToFile);
419  }
420  else if (nuclideName == "222Rn") {
421  continue;
422  } //Rn222 is handled separately later
423  else if (nuclideName == "59Ni") {
424  continue;
425  } //Rn222 is handled separately later
426  else if (nuclideName == "42Ar") {
427  readfile(
428  "42Ar_1",
429  "Argon_42_1.root",
430  fPathToFile); //Each possible beta decay mode of Ar42 is given it's own .root file for now.
431  readfile(
432  "42Ar_2",
433  "Argon_42_2.root",
434  fPathToFile); //This allows us to know which decay chain to follow for the dexcitation gammas.
435  readfile(
436  "42Ar_3",
437  "Argon_42_3.root",
438  fPathToFile); //The dexcitation gammas are not included in the root files as we want to
439  readfile(
440  "42Ar_4",
441  "Argon_42_4.root",
442  fPathToFile); //probabilistically simulate the correct coincident gammas, which we cannot guarantee
443  readfile("42Ar_5", "Argon_42_5.root", fPathToFile); //by sampling a histogram.
444  continue;
445  } //Ar42 is handeled separately later
446  else {
447  std::string searchName = nuclideName;
448  searchName += ".root";
449  readfile(nuclideName, searchName, fPathToFile);
450  }
451  }
452  }
code to link reconstructed objects back to the MC truth information
base_engine_t & createEngine(seed_t seed)
TTree * t1
Definition: plottest35.C:26
std::vector< double > fZ0
Bottom corner z position (cm) in world coordinates.
std::vector< std::string > fNuclide
List of nuclides to simulate. Example: "39Ar".
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::pair< double, double > defaulttimewindow() const
Returns the start and end of the readout window.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
void dumpNuclideSettings(Stream &&out, std::size_t iNucl, std::string const &volumeName={}) const
Prints the settings for the specified nuclide and volume.
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:49
std::vector< double > fT1
End of time window to simulate in ns.
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:292
std::vector< double > fY1
Top corner y position (cm) in world coordinates.
std::vector< double > fZ1
Top corner z position (cm) in world coordinates.
void readfile(std::string nuclide, std::string const &filename, std::string const &PathToFile)
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::size_t addvolume(std::string const &volumeName)
std::string fPathToFile
CLHEP::HepRandomEngine & fEngine
std::vector< double > fT0
Beginning of time window to simulate in ns.
std::vector< std::string > fMaterial
List of regexes of materials in which to generate the decays. Example: "LAr".
std::vector< double > fX0
Bottom corner x position (cm) in world coordinates.
std::vector< double > fBq
Radioactivity in Becquerels (decay per sec) per cubic cm.
std::vector< double > fX1
Top corner x position (cm) in world coordinates.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fY0
Bottom corner y position (cm) in world coordinates.

Member Function Documentation

std::size_t evgen::RadioGen::addvolume ( std::string const &  volumeName)
private

Adds all volumes with the specified name to the coordinates. Volumes must be boxes aligned to the world frame axes.

Definition at line 481 of file RadioGen_module.cc.

References geo::ROOTGeometryNavigator::apply(), fX0, fX1, fY0, fY1, fZ0, fZ1, and geo::origin().

Referenced by RadioGen().

482  {
483 
484  auto const& geom = *(lar::providerFrom<geo::Geometry>());
485 
486  std::vector<geo::GeoNodePath> volumePaths;
487  auto findVolume = [&volumePaths, volumeName](auto& path) {
488  if (path.current()->GetVolume()->GetName() == volumeName) volumePaths.push_back(path);
489  return true;
490  };
491 
492  geo::ROOTGeometryNavigator navigator{*(geom.ROOTGeoManager())};
493 
494  navigator.apply(findVolume);
495 
496  for (geo::GeoNodePath const& path : volumePaths) {
497 
498  //
499  // find the coordinates of the volume in local coordinates
500  //
501  TGeoShape const* pShape = path.current()->GetVolume()->GetShape();
502  auto pBox = dynamic_cast<TGeoBBox const*>(pShape);
503  if (!pBox) {
504  throw cet::exception("RadioGen") << "Volume '" << path.current()->GetName() << "' is a "
505  << pShape->IsA()->GetName() << ", not a TGeoBBox.\n";
506  }
507 
508  geo::Point_t const origin = geo::vect::makeFromCoords<geo::Point_t>(pBox->GetOrigin());
509  geo::Vector_t const diag = {pBox->GetDX(), pBox->GetDY(), pBox->GetDZ()};
510 
511  //
512  // convert to world coordinates
513  //
514 
515  auto const trans = path.currentTransformation<geo::TransformationMatrix>();
516 
517  geo::Point_t min, max;
518  trans.Transform(origin - diag, min);
519  trans.Transform(origin + diag, max);
520 
521  //
522  // add to the coordinates
523  //
524  fX0.push_back(std::min(min.X(), max.X()));
525  fY0.push_back(std::min(min.Y(), max.Y()));
526  fZ0.push_back(std::min(min.Z(), max.Z()));
527  fX1.push_back(std::max(min.X(), max.X()));
528  fY1.push_back(std::max(min.Y(), max.Y()));
529  fZ1.push_back(std::max(min.Z(), max.Z()));
530 
531  } // for
532 
533  return volumePaths.size();
534  } // RadioGen::addvolume()
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
std::vector< double > fZ0
Bottom corner z position (cm) in world coordinates.
Executes an operation on all the nodes of the ROOT geometry.
bool apply(GeoNodePath &path, Op &&op) const
Applies the specified operation to all nodes under the path.
std::vector< double > fY1
Top corner y position (cm) in world coordinates.
std::vector< double > fZ1
Top corner z position (cm) in world coordinates.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
Representation of a node and its ancestry.
Definition: GeoNodePath.h:34
std::vector< double > fX0
Bottom corner x position (cm) in world coordinates.
std::vector< double > fX1
Top corner x position (cm) in world coordinates.
ROOT::Math::Transform3D TransformationMatrix
Type of transformation matrix used in geometry.
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:229
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fY0
Bottom corner y position (cm) in world coordinates.
void evgen::RadioGen::Ar42Gamma2 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &  v_prods)
private

Definition at line 934 of file RadioGen_module.cc.

References dirCalc().

Referenced by Ar42Gamma3(), Ar42Gamma4(), Ar42Gamma5(), and SampleOne().

935  {
936  ti_PDGID pdgid = 22;
937  td_Mass m = 0.0; //we are writing gammas
938  std::vector<double> vd_p = {.0015246}; //Momentum in GeV
939  for (auto p : vd_p) {
940  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
941  }
942  }
TLorentzVector dirCalc(double p, double m)
void evgen::RadioGen::Ar42Gamma3 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &  v_prods)
private

Definition at line 944 of file RadioGen_module.cc.

References Ar42Gamma2(), and dirCalc().

Referenced by SampleOne().

945  {
946  ti_PDGID pdgid = 22;
947  td_Mass m = 0.0; //we are writing gammas
948  std::vector<double> vd_p = {.0003126};
949  for (auto p : vd_p) {
950  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
951  }
952  Ar42Gamma2(v_prods);
953  }
TLorentzVector dirCalc(double p, double m)
void Ar42Gamma2(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
void evgen::RadioGen::Ar42Gamma4 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &  v_prods)
private

Definition at line 955 of file RadioGen_module.cc.

References Ar42Gamma2(), dirCalc(), and fEngine.

Referenced by Ar42Gamma5(), and SampleOne().

956  {
957  CLHEP::RandFlat flat(fEngine);
958  ti_PDGID pdgid = 22;
959  td_Mass m = 0.0; //we are writing gammas
960  double chan1 = (0.052 / (0.052 + 0.020));
961  if (flat.fire() < chan1) {
962  std::vector<double> vd_p = {.0008997}; //Momentum in GeV
963  for (auto p : vd_p) {
964  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
965  }
966  Ar42Gamma2(v_prods);
967  }
968  else {
969  std::vector<double> vd_p = {.0024243}; //Momentum in GeV
970  for (auto p : vd_p) {
971  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
972  }
973  }
974  }
TLorentzVector dirCalc(double p, double m)
void Ar42Gamma2(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
CLHEP::HepRandomEngine & fEngine
void evgen::RadioGen::Ar42Gamma5 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &  v_prods)
private

Definition at line 976 of file RadioGen_module.cc.

References Ar42Gamma2(), Ar42Gamma4(), dirCalc(), and fEngine.

Referenced by SampleOne().

977  {
978  CLHEP::RandFlat flat(fEngine);
979  ti_PDGID pdgid = 22;
980  td_Mass m = 0.0; //we are writing gammas
981  double chan1 = (0.0033 / (0.0033 + 0.0201 + 0.041));
982  double chan2 = (0.0201 / (0.0033 + 0.0201 + 0.041));
983  double chanPick = flat.fire();
984  if (chanPick < chan1) {
985  std::vector<double> vd_p = {0.000692, 0.001228}; //Momentum in GeV
986  for (auto p : vd_p) {
987  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
988  }
989  Ar42Gamma2(v_prods);
990  }
991  else if (chanPick < (chan1 + chan2)) {
992  std::vector<double> vd_p = {0.0010212}; //Momentum in GeV
993  for (auto p : vd_p) {
994  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
995  }
996  Ar42Gamma4(v_prods);
997  }
998  else {
999  std::vector<double> vd_p = {0.0019208}; //Momentum in GeV
1000  for (auto p : vd_p) {
1001  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
1002  }
1003  Ar42Gamma2(v_prods);
1004  }
1005  }
TLorentzVector dirCalc(double p, double m)
void Ar42Gamma2(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
CLHEP::HepRandomEngine & fEngine
void Ar42Gamma4(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
void evgen::RadioGen::beginRun ( art::Run run)
privatevirtual

Reimplemented from art::EDProducer.

Definition at line 455 of file RadioGen_module.cc.

References geo::GeometryCore::DetectorName(), art::fullRun(), and art::Run::put().

456  {
458  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
459  }
constexpr auto fullRun()
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:140
ROOT libraries.
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

62  {
63  return collector_.consumes<T, BT>(tag);
64  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ProductToken< T > consumes(InputTag const &)
ConsumesCollector & art::ModuleBase::consumesCollector ( )
protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

58  {
59  return collector_;
60  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

76  {
77  collector_.consumesMany<T, BT>();
78  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

69  {
70  return collector_.consumesView<T, BT>(tag);
71  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ViewToken< Element > consumesView(InputTag const &)
std::pair< double, double > evgen::RadioGen::defaulttimewindow ( ) const
private

Returns the start and end of the readout window.

Definition at line 1044 of file RadioGen_module.cc.

References DEFINE_ART_MODULE, detinfo::makeDetectorTimings(), detinfo::DetectorPropertiesData::NumberTimeSamples(), and detinfo::DetectorPropertiesData::ReadOutWindowSize().

Referenced by RadioGen().

1045  {
1046  // values are in simulation time scale (and therefore in [ns])
1047  using namespace detinfo::timescales;
1048 
1049  auto const& timings = detinfo::makeDetectorTimings(
1051  detinfo::DetectorPropertiesData const& detInfo =
1053 
1054  //
1055  // we take a number of (TPC electronics) ticks before the trigger time,
1056  // and we go to a number of ticks after the trigger time;
1057  // that shift is one readout window size
1058  //
1059 
1060  auto const trigTimeTick = timings.toTick<electronics_tick>(timings.TriggerTime());
1061  electronics_time_ticks const beforeTicks{-static_cast<int>(detInfo.ReadOutWindowSize())};
1062  electronics_time_ticks const afterTicks{detInfo.NumberTimeSamples()};
1063 
1064  return {double(timings.toTimeScale<simulation_time>(trigTimeTick + beforeTicks)),
1065  double(timings.toTimeScale<simulation_time>(trigTimeTick + afterTicks))};
1066  } // RadioGen::defaulttimewindow()
timescale_traits< ElectronicsTimeCategory >::tick_t electronics_tick
A point on the electronics time scale expressed in its ticks.
timescale_traits< ElectronicsTimeCategory >::tick_interval_t electronics_time_ticks
An interval on the electronics time scale expressed in its ticks.
timescale_traits< SimulationTimeCategory >::time_point_t simulation_time
A point in time on the simulation time scale.
detinfo::DetectorTimings makeDetectorTimings(detinfo::DetectorClocksData const &detClocks)
Namespace including different time scales as defined in LArSoft.
TLorentzVector evgen::RadioGen::dirCalc ( double  p,
double  m 
)
private

Definition at line 679 of file RadioGen_module.cc.

References fEngine.

Referenced by Ar42Gamma2(), Ar42Gamma3(), Ar42Gamma4(), Ar42Gamma5(), and SampleOne().

680  {
681  CLHEP::RandFlat flat(fEngine);
682  // isotropic production angle for the decay product
683  double costheta = (2.0 * flat.fire() - 1.0);
684  if (costheta < -1.0) costheta = -1.0;
685  if (costheta > 1.0) costheta = 1.0;
686  double const sintheta = sqrt(1.0 - costheta * costheta);
687  double const phi = 2.0 * M_PI * flat.fire();
688  return TLorentzVector{p * sintheta * std::cos(phi),
689  p * sintheta * std::sin(phi),
690  p * costheta,
691  std::sqrt(p * p + m * m)};
692  }
CLHEP::HepRandomEngine & fEngine
void art::detail::Producer::doBeginJob ( SharedResources const &  resources)
inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

23  {
24  setupQueues(resources);
25  ProcessingFrame const frame{ScheduleID{}};
26  beginJobWithFrame(frame);
27  }
virtual void setupQueues(SharedResources const &)=0
virtual void beginJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doBeginRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

66  {
67  auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
68  ProcessingFrame const frame{mc.scheduleID()};
69  beginRunWithFrame(r, frame);
70  r.commitProducts();
71  return true;
72  }
TRandom r
Definition: spectrum.C:23
virtual void beginRunWithFrame(Run &, ProcessingFrame const &)=0
static RangeSet forRun(RunID)
Definition: RangeSet.cc:51
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

86  {
87  auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
88  ProcessingFrame const frame{mc.scheduleID()};
89  beginSubRunWithFrame(sr, frame);
90  sr.commitProducts();
91  return true;
92  }
virtual void beginSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
static RangeSet forSubRun(SubRunID)
Definition: RangeSet.cc:57
void art::detail::Producer::doEndJob ( )
inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

31  {
32  ProcessingFrame const frame{ScheduleID{}};
33  endJobWithFrame(frame);
34  }
virtual void endJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doEndRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

76  {
77  auto r = rp.makeRun(mc, rp.seenRanges());
78  ProcessingFrame const frame{mc.scheduleID()};
79  endRunWithFrame(r, frame);
80  r.commitProducts();
81  return true;
82  }
TRandom r
Definition: spectrum.C:23
virtual void endRunWithFrame(Run &, ProcessingFrame const &)=0
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

96  {
97  auto sr = srp.makeSubRun(mc, srp.seenRanges());
98  ProcessingFrame const frame{mc.scheduleID()};
99  endSubRunWithFrame(sr, frame);
100  sr.commitProducts();
101  return true;
102  }
virtual void endSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
bool art::detail::Producer::doEvent ( EventPrincipal ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

110  {
111  auto e = ep.makeEvent(mc);
112  ++counts_run;
113  ProcessingFrame const frame{mc.scheduleID()};
114  produceWithFrame(e, frame);
115  e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
116  ++counts_passed;
117  return true;
118  }
bool const checkPutProducts_
Definition: Producer.h:70
Float_t e
Definition: plot.C:35
virtual void produceWithFrame(Event &, ProcessingFrame const &)=0
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

45  {
46  ProcessingFrame const frame{ScheduleID{}};
48  }
virtual void respondToCloseInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

59  {
60  ProcessingFrame const frame{ScheduleID{}};
62  }
virtual void respondToCloseOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

38  {
39  ProcessingFrame const frame{ScheduleID{}};
41  }
virtual void respondToOpenInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

52  {
53  ProcessingFrame const frame{ScheduleID{}};
55  }
virtual void respondToOpenOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
template<typename Stream >
void evgen::RadioGen::dumpNuclideSettings ( Stream &&  out,
std::size_t  iNucl,
std::string const &  volumeName = {} 
) const
private

Prints the settings for the specified nuclide and volume.

Definition at line 1030 of file RadioGen_module.cc.

References fBq, fMaterial, fNuclide, fT0, fT1, fX0, fX1, fY0, fY1, fZ0, and fZ1.

Referenced by RadioGen().

1033  {
1034 
1035  out << "[#" << iNucl << "] " << fNuclide[iNucl] << " (" << fBq[iNucl] << " Bq/cm^3)"
1036  << " in " << fMaterial[iNucl] << " from " << fT0[iNucl] << " to " << fT1[iNucl]
1037  << " ns in ( " << fX0[iNucl] << ", " << fY0[iNucl] << ", " << fZ0[iNucl] << ") to ( "
1038  << fX1[iNucl] << ", " << fY1[iNucl] << ", " << fZ1[iNucl] << ")";
1039  if (!volumeName.empty()) out << " (\"" << volumeName << "\")";
1040 
1041  } // RadioGen::dumpNuclideSettings()
std::vector< double > fZ0
Bottom corner z position (cm) in world coordinates.
std::vector< std::string > fNuclide
List of nuclides to simulate. Example: "39Ar".
std::vector< double > fT1
End of time window to simulate in ns.
std::vector< double > fY1
Top corner y position (cm) in world coordinates.
std::vector< double > fZ1
Top corner z position (cm) in world coordinates.
std::vector< double > fT0
Beginning of time window to simulate in ns.
std::vector< std::string > fMaterial
List of regexes of materials in which to generate the decays. Example: "LAr".
std::vector< double > fX0
Bottom corner x position (cm) in world coordinates.
std::vector< double > fBq
Radioactivity in Becquerels (decay per sec) per cubic cm.
std::vector< double > fX1
Top corner x position (cm) in world coordinates.
std::vector< double > fY0
Bottom corner y position (cm) in world coordinates.
void art::Modifier::fillProductDescriptions ( )
inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

11  {
13  }
void fillDescriptions(ModuleDescription const &md)
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const
inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

44  {
45  return collector_.getConsumables();
46  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables() const
std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const &  wp)
inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

38  {
39  return doMakeWorker(wp);
40  }
virtual std::unique_ptr< Worker > doMakeWorker(WorkerParams const &wp)=0
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

83  {
84  return collector_.mayConsume<T, BT>(tag);
85  }
ProductToken< T > mayConsume(InputTag const &)
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

97  {
98  collector_.mayConsumeMany<T, BT>();
99  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

90  {
91  return collector_.mayConsumeView<T, BT>(tag);
92  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ViewToken< Element > mayConsumeView(InputTag const &)
ModuleDescription const & art::ModuleBase::moduleDescription ( ) const
inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

14  {
15  if (md_.has_value()) {
16  return *md_;
17  }
18 
20  "There was an error while calling moduleDescription().\n"}
21  << "The moduleDescription() base-class member function cannot be called\n"
22  "during module construction. To determine which module is "
23  "responsible\n"
24  "for calling it, find the '<module type>:<module "
25  "label>@Construction'\n"
26  "tag in the message prefix above. Please contact artists@fnal.gov\n"
27  "for guidance.\n";
28  }
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void evgen::RadioGen::produce ( art::Event evt)
privatevirtual

unique_ptr allows ownership to be transferred to the art::Event after the put statement

Implements art::EDProducer.

Definition at line 462 of file RadioGen_module.cc.

References fNuclide, simb::kSingleParticle, MF_LOG_DEBUG, art::Event::put(), SampleOne(), simb::MCTruth::SetOrigin(), and trackidcounter.

463  {
465  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
466 
467  simb::MCTruth truth;
469 
470  trackidcounter = -1;
471  for (unsigned int i = 0; i < fNuclide.size(); ++i) {
472  SampleOne(i, truth);
473  } //end loop over nuclides
474 
475  MF_LOG_DEBUG("RadioGen") << truth;
476  truthcol->push_back(truth);
477  evt.put(std::move(truthcol));
478  }
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
std::vector< std::string > fNuclide
List of nuclides to simulate. Example: "39Ar".
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
single particles thrown at the detector
Definition: MCTruth.h:26
int trackidcounter
Serial number for the MC track ID.
void SampleOne(unsigned int i, simb::MCTruth &mct)
Checks that the node represents a box well aligned to world frame axes.
#define MF_LOG_DEBUG(id)
Event generator information.
Definition: MCTruth.h:32
void evgen::RadioGen::readfile ( std::string  nuclide,
std::string const &  filename,
std::string const &  PathToFile 
)
private

Definition at line 696 of file RadioGen_module.cc.

References alphaintegral, alphaspectrum, betaintegral, betaspectrum, f, fNuclide, gammaintegral, gammaspectrum, neutronintegral, neutronspectrum, spectrumname, and y.

Referenced by RadioGen().

699  {
700  bool found{false};
701  std::regex const re_argon{"42Ar.*"};
702  for (size_t i = 0; i < fNuclide.size(); i++) {
703  if (
704  fNuclide[i] ==
705  nuclide) { //This check makes sure that the nuclide we are searching for is in fact in our fNuclide list. Ar42 handeled separately.
706  found = true;
707  break;
708  } //End If nuclide is in our list. Next is the special case of Ar42
709  else if (std::regex_match(nuclide, re_argon) && fNuclide[i] == "42Ar") {
710  found = true;
711  break;
712  }
713  }
714  if (!found) return;
715 
716  Bool_t addStatus = TH1::AddDirectoryStatus();
717  TH1::AddDirectory(
718  kFALSE); // cloned histograms go in memory, and aren't deleted when files are closed.
719  // be sure to restore this state when we're out of the routine.
720 
721  spectrumname.push_back(nuclide);
722  cet::search_path sp("FW_SEARCH_PATH");
723  std::string fn2 = PathToFile;
724  fn2 += filename;
725  std::string fullname;
726  sp.find_file(fn2, fullname);
727  struct stat sb;
728  if (fullname.empty() || stat(fullname.c_str(), &sb) != 0)
729  throw cet::exception("RadioGen")
730  << "Input spectrum file " << fn2 << " not found in FW_SEARCH_PATH!\n";
731 
732  TFile f(fullname.c_str(), "READ");
733  TGraph* alphagraph = (TGraph*)f.Get("Alphas");
734  TGraph* betagraph = (TGraph*)f.Get("Betas");
735  TGraph* gammagraph = (TGraph*)f.Get("Gammas");
736  TGraph* neutrongraph = (TGraph*)f.Get("Neutrons");
737 
738  if (alphagraph) {
739  int np = alphagraph->GetN();
740  double* y = alphagraph->GetY();
741  std::string name;
742  name = "RadioGen_";
743  name += nuclide;
744  name += "_Alpha";
745  auto alphahist = std::make_unique<TH1D>(name.c_str(), "Alpha Spectrum", np, 0, np);
746  for (int i = 0; i < np; i++) {
747  alphahist->SetBinContent(i + 1, y[i]);
748  alphahist->SetBinError(i + 1, 0);
749  }
750  alphaintegral.push_back(alphahist->Integral());
751  alphaspectrum.push_back(move(alphahist));
752  }
753  else {
754  alphaintegral.push_back(0);
755  alphaspectrum.push_back(nullptr);
756  }
757 
758  if (betagraph) {
759  int np = betagraph->GetN();
760 
761  double* y = betagraph->GetY();
762  std::string name;
763  name = "RadioGen_";
764  name += nuclide;
765  name += "_Beta";
766  auto betahist = std::make_unique<TH1D>(name.c_str(), "Beta Spectrum", np, 0, np);
767  for (int i = 0; i < np; i++) {
768  betahist->SetBinContent(i + 1, y[i]);
769  betahist->SetBinError(i + 1, 0);
770  }
771  betaintegral.push_back(betahist->Integral());
772  betaspectrum.push_back(move(betahist));
773  }
774  else {
775  betaintegral.push_back(0);
776  betaspectrum.push_back(nullptr);
777  }
778 
779  if (gammagraph) {
780  int np = gammagraph->GetN();
781  double* y = gammagraph->GetY();
782  std::string name;
783  name = "RadioGen_";
784  name += nuclide;
785  name += "_Gamma";
786  auto gammahist = std::make_unique<TH1D>(name.c_str(), "Gamma Spectrum", np, 0, np);
787  for (int i = 0; i < np; i++) {
788  gammahist->SetBinContent(i + 1, y[i]);
789  gammahist->SetBinError(i + 1, 0);
790  }
791  gammaintegral.push_back(gammahist->Integral());
792  gammaspectrum.push_back(move(gammahist));
793  }
794  else {
795  gammaintegral.push_back(0);
796  gammaspectrum.push_back(nullptr);
797  }
798 
799  if (neutrongraph) {
800  int np = neutrongraph->GetN();
801  double* y = neutrongraph->GetY();
802  std::string name;
803  name = "RadioGen_";
804  name += nuclide;
805  name += "_Neutron";
806  auto neutronhist = std::make_unique<TH1D>(name.c_str(), "Neutron Spectrum", np, 0, np);
807  for (int i = 0; i < np; i++) {
808  neutronhist->SetBinContent(i + 1, y[i]);
809  neutronhist->SetBinError(i + 1, 0);
810  }
811  neutronintegral.push_back(neutronhist->Integral());
812  neutronspectrum.push_back(move(neutronhist));
813  }
814  else {
815  neutronintegral.push_back(0);
816  neutronspectrum.push_back(nullptr);
817  }
818 
819  f.Close();
820  TH1::AddDirectory(addStatus);
821 
822  double total =
823  alphaintegral.back() + betaintegral.back() + gammaintegral.back() + neutronintegral.back();
824  if (total > 0) {
825  alphaintegral.back() /= total;
826  betaintegral.back() /= total;
827  gammaintegral.back() /= total;
828  neutronintegral.back() /= total;
829  }
830  }
std::vector< std::unique_ptr< TH1D > > gammaspectrum
std::vector< std::string > fNuclide
List of nuclides to simulate. Example: "39Ar".
Float_t y
Definition: compare.C:6
std::vector< double > neutronintegral
std::vector< std::string > spectrumname
TFile f
Definition: plotHisto.C:6
std::vector< std::unique_ptr< TH1D > > alphaspectrum
std::vector< double > alphaintegral
std::vector< double > betaintegral
std::vector< std::unique_ptr< TH1D > > neutronspectrum
std::vector< double > gammaintegral
std::vector< std::unique_ptr< TH1D > > betaspectrum
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void art::Modifier::registerProducts ( ProductDescriptions productsToRegister)
inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

17  {
18  ProductRegistryHelper::registerProducts(productsToRegister,
20  }
void registerProducts(ProductDescriptions &productsToRegister, ModuleDescription const &md)
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
double evgen::RadioGen::samplefromth1d ( TH1D &  hist)
private

Definition at line 896 of file RadioGen_module.cc.

References fEngine, nbinsx, and x.

Referenced by samplespectrum().

897  {
898  CLHEP::RandFlat flat(fEngine);
899 
900  int nbinsx = hist.GetNbinsX();
901  std::vector<double> partialsum;
902  partialsum.resize(nbinsx + 1);
903  partialsum[0] = 0;
904 
905  for (int i = 1; i <= nbinsx; i++) {
906  double hc = hist.GetBinContent(i);
907  if (hc < 0)
908  throw cet::exception("RadioGen") << "Negative bin: " << i << " " << hist.GetName() << "\n";
909  partialsum[i] = partialsum[i - 1] + hc;
910  }
911  double integral = partialsum[nbinsx];
912  if (integral == 0) return 0;
913  // normalize to unit sum
914  for (int i = 1; i <= nbinsx; i++)
915  partialsum[i] /= integral;
916 
917  double r1 = flat.fire();
918  int ibin = TMath::BinarySearch(nbinsx, &(partialsum[0]), r1);
919  Double_t x = hist.GetBinLowEdge(ibin + 1);
920  if (r1 > partialsum[ibin]) {
921  x += hist.GetBinWidth(ibin + 1) * (r1 - partialsum[ibin]) /
922  (partialsum[ibin + 1] - partialsum[ibin]);
923  }
924  return x;
925  }
Float_t x
Definition: compare.C:6
TH2F * hist
Definition: plot.C:134
Int_t nbinsx
Definition: plot.C:23
CLHEP::HepRandomEngine & fEngine
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgen::RadioGen::SampleOne ( unsigned int  i,
simb::MCTruth mct 
)
private

Checks that the node represents a box well aligned to world frame axes.

Definition at line 539 of file RadioGen_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), Ar42Gamma2(), Ar42Gamma3(), Ar42Gamma4(), Ar42Gamma5(), dirCalc(), energy, fBq, fEngine, fIsFirstSignalSpecial, fMaterial, fNuclide, fT0, fT1, fX0, fX1, fY0, fY1, fZ0, fZ1, part, geo::GeometryCore::ROOTGeoManager(), samplespectrum(), and trackidcounter.

Referenced by produce().

540  {
542  TGeoManager* geomanager = geo->ROOTGeoManager();
543 
544  CLHEP::RandFlat flat(fEngine);
545  CLHEP::RandPoisson poisson(fEngine);
546 
547  // figure out how many decays to generate, assuming that the entire prism consists of the radioactive material.
548  // we will skip over decays in other materials later.
549 
550  double rate =
551  fabs(fBq[i] * (fT1[i] - fT0[i]) * (fX1[i] - fX0[i]) * (fY1[i] - fY0[i]) * (fZ1[i] - fZ0[i])) /
552  1.0E9;
553  long ndecays = poisson.shoot(rate);
554 
555  std::regex const re_material{fMaterial[i]};
556  for (unsigned int idecay = 0; idecay < ndecays; idecay++) {
557  // generate just one particle at a time. Need a little recoding if a radioactive
558  // decay generates multiple daughters that need simulation
559  // uniformly distributed in position and time
560  //
561  // JStock: Leaving this as a single position for the decay products. For now I will assume they all come from the same spot.
562  TLorentzVector pos(
563  fX0[i] + flat.fire() * (fX1[i] - fX0[i]),
564  fY0[i] + flat.fire() * (fY1[i] - fY0[i]),
565  fZ0[i] + flat.fire() * (fZ1[i] - fZ0[i]),
566  (idecay == 0 && fIsFirstSignalSpecial) ? 0 : (fT0[i] + flat.fire() * (fT1[i] - fT0[i])));
567 
568  // discard decays that are not in the proper material
569  std::string volmaterial =
570  geomanager->FindNode(pos.X(), pos.Y(), pos.Z())->GetMedium()->GetMaterial()->GetName();
571  if (!std::regex_match(volmaterial, re_material)) continue;
572 
573  //Moved pdgid into the next statement, so that it is localized.
574  // electron=11, photon=22, alpha = 1000020040, neutron = 2112
575 
576  //JStock: Allow us to have different particles from the same decay. This requires multiple momenta.
577  std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>
578  v_prods; //(First is for PDGID, second is mass, third is Momentum)
579 
580  if (fNuclide[i] == "222Rn") // Treat 222Rn separately
581  {
582  double p = 0;
583  double t = 0.00548952;
584  td_Mass m = m_alpha;
585  ti_PDGID pdgid = 1000020040; //td_Mass = double. ti_PDGID = int;
586  double energy = t + m;
587  double p2 = energy * energy - m * m;
588  if (p2 > 0)
589  p = TMath::Sqrt(p2);
590  else
591  p = 0;
592  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
593  } //End special case RN222
594  else if (
595  fNuclide[i] ==
596  "59Ni") { //Treat 59Ni Calibration Source separately (as I haven't made a spectrum for it, and ultimately it should be handeled with multiple particle outputs.
597  double p = 0.008997;
598  td_Mass m = 0;
599  ti_PDGID pdgid =
600  22; // td_Mas=double. ti_PDFID=int. Assigning p directly, as t=p for gammas.
601  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
602  } //end special case Ni59 calibration source
603  else if (fNuclide[i] == "42Ar") { // Spot for special treatment of Ar42.
604  double p = 0;
605  double t = 0;
606  td_Mass m = 0;
607  ti_PDGID pdgid = 0; //td_Mass = double. ti_PDGID = int;
608  double bSelect = flat.fire(); //Make this a random number from 0 to 1.
609  if (bSelect < 0.819) { //beta channel 1. No Gamma. beta Q value 3525.22 keV
610  samplespectrum("42Ar_1", pdgid, t, m, p);
611  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
612  //No gamma here.
613  }
614  else if (bSelect < 0.9954) { //beta channel 2. 1 Gamma (1524.6 keV). beta Q value 2000.62
615  samplespectrum("42Ar_2", pdgid, t, m, p);
616  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
617  Ar42Gamma2(v_prods);
618  }
619  else if (
620  bSelect <
621  0.9988) { //beta channel 3. 1 Gamma Channel. 312.6 keV + gamma 2. beta Q value 1688.02 keV
622  samplespectrum("42Ar_3", pdgid, t, m, p);
623  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
624  Ar42Gamma3(v_prods);
625  }
626  else if (
627  bSelect <
628  0.9993) { //beta channel 4. 2 Gamma Channels. Either 899.7 keV (i 0.052) + gamma 2 or 2424.3 keV (i 0.020). beta Q value 1100.92 keV
629  samplespectrum("42Ar_4", pdgid, t, m, p);
630  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
631  Ar42Gamma4(v_prods);
632  }
633  else { //beta channel 5. 3 gamma channels. 692.0 keV + 1228.0 keV + Gamma 2 (i 0.0033) ||OR|| 1021.2 keV + gamma 4 (i 0.0201) ||OR|| 1920.8 keV + gamma 2 (i 0.041). beta Q value 79.82 keV
634  samplespectrum("42Ar_5", pdgid, t, m, p);
635  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
636  Ar42Gamma5(v_prods);
637  }
638  //Add beta.
639  //Call gamma function for beta mode.
640  }
641  else { //General Case.
642  double p = 0;
643  double t = 0;
644  td_Mass m = 0;
645  ti_PDGID pdgid = 0; //td_Mass = double. ti_PDGID = int;
646  samplespectrum(fNuclide[i], pdgid, t, m, p);
647  std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom =
648  std::make_tuple(pdgid, m, dirCalc(p, m));
649  v_prods.push_back(partMassMom);
650  } //end else (not RN or other special case
651 
652  //JStock: Modify this to now loop over the v_prods.
653  for (auto prodEntry : v_prods) {
654  // set track id to a negative serial number as these are all primary particles and have id <= 0
655  int trackid = trackidcounter;
656  ti_PDGID pdgid = std::get<0>(prodEntry);
657  td_Mass m = std::get<1>(prodEntry);
658  TLorentzVector pvec = std::get<2>(prodEntry);
659  trackidcounter--;
660  std::string primary("primary");
661 
662  // alpha particles need a little help since they're not in the TDatabasePDG table
663  // // so don't rely so heavily on default arguments to the MCParticle constructor
664  if (pdgid == 1000020040) {
665  simb::MCParticle part(trackid, pdgid, primary, -1, m, 1);
666  part.AddTrajectoryPoint(pos, pvec);
667  mct.Add(part);
668  } // end "If alpha"
669  else {
670  simb::MCParticle part(trackid, pdgid, primary);
671  part.AddTrajectoryPoint(pos, pvec);
672  mct.Add(part);
673  } // end All standard cases.
674  } //End Loop over all particles produces in this single decay.
675  }
676  }
TLorentzVector dirCalc(double p, double m)
std::vector< double > fZ0
Bottom corner z position (cm) in world coordinates.
std::vector< std::string > fNuclide
List of nuclides to simulate. Example: "39Ar".
void Ar42Gamma5(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
TGeoManager * ROOTGeoManager() const
Access to the ROOT geometry description manager.
TString part[npart]
Definition: Style.C:32
std::vector< double > fT1
End of time window to simulate in ns.
std::vector< double > fY1
Top corner y position (cm) in world coordinates.
std::vector< double > fZ1
Top corner z position (cm) in world coordinates.
int trackidcounter
Serial number for the MC track ID.
double energy
Definition: plottest35.C:25
void Ar42Gamma3(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
void samplespectrum(std::string nuclide, int &itype, double &t, double &m, double &p)
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
void Ar42Gamma2(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
CLHEP::HepRandomEngine & fEngine
std::vector< double > fT0
Beginning of time window to simulate in ns.
void Ar42Gamma4(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
std::vector< std::string > fMaterial
List of regexes of materials in which to generate the decays. Example: "LAr".
std::vector< double > fX0
Bottom corner x position (cm) in world coordinates.
std::vector< double > fBq
Radioactivity in Becquerels (decay per sec) per cubic cm.
ROOT libraries.
std::vector< double > fX1
Top corner x position (cm) in world coordinates.
std::vector< double > fY0
Bottom corner y position (cm) in world coordinates.
void evgen::RadioGen::samplespectrum ( std::string  nuclide,
int &  itype,
double &  t,
double &  m,
double &  p 
)
private

Definition at line 832 of file RadioGen_module.cc.

References alphaintegral, alphaspectrum, betaintegral, betaspectrum, e, fEngine, gammaintegral, gammaspectrum, neutronspectrum, samplefromth1d(), and spectrumname.

Referenced by SampleOne().

833  {
834  CLHEP::RandFlat flat(fEngine);
835 
836  int inuc = -1;
837  for (size_t i = 0; i < spectrumname.size(); i++) {
838  if (nuclide == spectrumname[i]) {
839  inuc = i;
840  break;
841  }
842  }
843  if (inuc == -1) {
844  t = 0; // throw an exception in the future
845  itype = 0;
846  throw cet::exception("RadioGen") << "Ununderstood nuclide: " << nuclide << "\n";
847  }
848 
849  double rtype = flat.fire();
850 
851  itype = -1;
852  m = 0;
853  p = 0;
854  for (
855  int itry = 0; itry < 10;
856  itry++) // maybe a tiny normalization issue with a sum of 0.99999999999 or something, so try a few times.
857  {
858  if (rtype <= alphaintegral[inuc] && alphaspectrum[inuc] != nullptr) {
859  itype = 1000020040; // alpha
860  m = m_alpha;
861  t = samplefromth1d(*alphaspectrum[inuc]) / 1000000.0;
862  }
863  else if (rtype <= alphaintegral[inuc] + betaintegral[inuc] && betaspectrum[inuc] != nullptr) {
864  itype = 11; // beta
865  m = m_e;
866  t = samplefromth1d(*betaspectrum[inuc]) / 1000000.0;
867  }
868  else if (rtype <= alphaintegral[inuc] + betaintegral[inuc] + gammaintegral[inuc] &&
869  gammaspectrum[inuc] != nullptr) {
870  itype = 22; // gamma
871  m = 0;
872  t = samplefromth1d(*gammaspectrum[inuc]) / 1000000.0;
873  }
874  else if (neutronspectrum[inuc] != nullptr) {
875  itype = 2112;
876  m = m_neutron;
877  t = samplefromth1d(*neutronspectrum[inuc]) / 1000000.0;
878  }
879  if (itype >= 0) break;
880  }
881  if (itype == -1) {
882  throw cet::exception("RadioGen")
883  << "Normalization problem with nuclide: " << nuclide << "\n";
884  }
885  double e = t + m;
886  p = e * e - m * m;
887  if (p >= 0) { p = TMath::Sqrt(p); }
888  else {
889  p = 0;
890  }
891  }
std::vector< std::unique_ptr< TH1D > > gammaspectrum
std::vector< std::string > spectrumname
std::vector< std::unique_ptr< TH1D > > alphaspectrum
double samplefromth1d(TH1D &hist)
std::vector< double > alphaintegral
std::vector< double > betaintegral
CLHEP::HepRandomEngine & fEngine
std::vector< std::unique_ptr< TH1D > > neutronspectrum
Float_t e
Definition: plot.C:35
std::vector< double > gammaintegral
std::vector< std::unique_ptr< TH1D > > betaspectrum
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  md)
inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

32  {
33  md_ = md;
34  }
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

50  {
51  // Now that we know we have seen all the consumes declarations,
52  // sort the results for fast lookup later.
53  collector_.sortConsumables(current_process_name);
54  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
void sortConsumables(std::string const &current_process_name)

Member Data Documentation

std::vector<double> evgen::RadioGen::alphaintegral
private

Definition at line 269 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<std::unique_ptr<TH1D> > evgen::RadioGen::alphaspectrum
private

Definition at line 268 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<double> evgen::RadioGen::betaintegral
private

Definition at line 271 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<std::unique_ptr<TH1D> > evgen::RadioGen::betaspectrum
private

Definition at line 270 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<double> evgen::RadioGen::fBq
private

Radioactivity in Becquerels (decay per sec) per cubic cm.

Definition at line 250 of file RadioGen_module.cc.

Referenced by dumpNuclideSettings(), RadioGen(), and SampleOne().

CLHEP::HepRandomEngine& evgen::RadioGen::fEngine
private
bool evgen::RadioGen::fIsFirstSignalSpecial
private

Definition at line 259 of file RadioGen_module.cc.

Referenced by RadioGen(), and SampleOne().

std::vector<std::string> evgen::RadioGen::fMaterial
private

List of regexes of materials in which to generate the decays. Example: "LAr".

Definition at line 249 of file RadioGen_module.cc.

Referenced by dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<std::string> evgen::RadioGen::fNuclide
private

List of nuclides to simulate. Example: "39Ar".

Definition at line 247 of file RadioGen_module.cc.

Referenced by dumpNuclideSettings(), produce(), RadioGen(), readfile(), and SampleOne().

std::string evgen::RadioGen::fPathToFile
private

Definition at line 261 of file RadioGen_module.cc.

Referenced by RadioGen().

std::vector<double> evgen::RadioGen::fT0
private

Beginning of time window to simulate in ns.

Definition at line 251 of file RadioGen_module.cc.

Referenced by dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fT1
private

End of time window to simulate in ns.

Definition at line 252 of file RadioGen_module.cc.

Referenced by dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fX0
private

Bottom corner x position (cm) in world coordinates.

Definition at line 253 of file RadioGen_module.cc.

Referenced by addvolume(), dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fX1
private

Top corner x position (cm) in world coordinates.

Definition at line 256 of file RadioGen_module.cc.

Referenced by addvolume(), dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fY0
private

Bottom corner y position (cm) in world coordinates.

Definition at line 254 of file RadioGen_module.cc.

Referenced by addvolume(), dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fY1
private

Top corner y position (cm) in world coordinates.

Definition at line 257 of file RadioGen_module.cc.

Referenced by addvolume(), dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fZ0
private

Bottom corner z position (cm) in world coordinates.

Definition at line 255 of file RadioGen_module.cc.

Referenced by addvolume(), dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::fZ1
private

Top corner z position (cm) in world coordinates.

Definition at line 258 of file RadioGen_module.cc.

Referenced by addvolume(), dumpNuclideSettings(), RadioGen(), and SampleOne().

std::vector<double> evgen::RadioGen::gammaintegral
private

Definition at line 273 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<std::unique_ptr<TH1D> > evgen::RadioGen::gammaspectrum
private

Definition at line 272 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<double> evgen::RadioGen::neutronintegral
private

Definition at line 275 of file RadioGen_module.cc.

Referenced by readfile().

std::vector<std::unique_ptr<TH1D> > evgen::RadioGen::neutronspectrum
private

Definition at line 274 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

std::vector<std::string> evgen::RadioGen::spectrumname
private

Definition at line 267 of file RadioGen_module.cc.

Referenced by readfile(), and samplespectrum().

int evgen::RadioGen::trackidcounter
private

Serial number for the MC track ID.

Definition at line 260 of file RadioGen_module.cc.

Referenced by produce(), and SampleOne().


The documentation for this class was generated from the following file: