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

Inheritance diagram for evgen::RadioGen:

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 193 of file RadioGen_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer

inherited

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 205 of file RadioGen_module.cc.

typedef int evgen::RadioGen::ti_PDGID

private

Definition at line 203 of file RadioGen_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 288 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.

     : EDProducer{pset}
     , fX0{pset.get<std::vector<double>>("X0", {})}
     , fY0{pset.get<std::vector<double>>("Y0", {})}
     , fZ0{pset.get<std::vector<double>>("Z0", {})}
     , fX1{pset.get<std::vector<double>>("X1", {})}
     , fY1{pset.get<std::vector<double>>("Y1", {})}
     , fZ1{pset.get<std::vector<double>>("Z1", {})}
     , fIsFirstSignalSpecial{pset.get<bool>("IsFirstSignalSpecial", false)}
     // create a default random engine; obtain the random seed from NuRandomService,
     // unless overridden in configuration with key "Seed"
     , fEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(createEngine(0),
                                                                                  pset,
                                                                                  "Seed"))
   {
     produces<std::vector<simb::MCTruth>>();
     produces<sumdata::RunData, art::InRun>();
 
     auto const fPathToFile = pset.get<std::string>("PathToFile", "Radionuclides/");
     auto const nuclide = pset.get<std::vector<std::string>>("Nuclide");
     auto const material = pset.get<std::vector<std::string>>("Material");
     auto const Bq = pset.get<std::vector<double>>("BqPercc");
     auto t0 = pset.get<std::vector<double>>("T0", {});
     auto t1 = pset.get<std::vector<double>>("T1", {});
 
     if (fPathToFile.find_last_of('/') !=
         fPathToFile.length() - 1) { //Checking that last character in path is a forward slash
       throw cet::exception("RadioGen")
         << "Last entry of the path to radiological files needs to be a forward slash ('/')\n";
     }
 
     if (fT0.empty() || fT1.empty()) { // better be both empty...
       if (!fT0.empty() || !fT1.empty()) {
         throw art::Exception(art::errors::Configuration)
           << "RadioGen T0 and T1 need to be both non-empty, or both empty"
              " (now T0 has "
           << fT0.size() << " entries and T1 has " << fT0.size() << ")\n";
       }
       auto const [defaultT0, defaultT1] = defaulttimewindow();
       t0.push_back(defaultT0);
       t1.push_back(defaultT1);
     }
 
     //
     // First, the materials are assigned to the coordinates explicitly
     // configured; then, each of the remaining materials is assigned to one
     // of the named volumes.
     // TODO: reaction to mismatches is "not so great".
     //
     mf::LogInfo("RadioGen") << "Configuring activity:";
     for (std::size_t iCoord : util::counter(fX0.size())) {
 
       fNuclide.push_back(nuclide.at(iCoord));
       fMaterial.push_back(material.at(iCoord));
       fBq.push_back(Bq.at(iCoord));
       // replicate the last timing if none is specified
       fT0.push_back(t0[std::min(iCoord, t0.size() - 1U)]);
       fT1.push_back(t1[std::min(iCoord, t1.size() - 1U)]);
 
       dumpNuclideSettings(mf::LogVerbatim("RadioGen"), iCoord);
 
     } // for
 
     std::size_t iNext = fX0.size();
     auto const volumes = pset.get<std::vector<std::string>>("Volumes", {});
     for (auto&& [iVolume, volName] : util::enumerate(volumes)) {
       // this expands the coordinate vectors
       auto const nVolumes = addvolume(volName);
       if (nVolumes == 0) {
         throw art::Exception(art::errors::Configuration)
           << "No volume named '" << volName << "' was found!\n";
       }
 
       std::size_t const iVolFirst = fNuclide.size();
       for (auto iVolInstance : util::counter(nVolumes)) {
         fNuclide.push_back(nuclide.at(iNext));
         fMaterial.push_back(material.at(iNext));
         fBq.push_back(Bq.at(iNext));
         // replicate the last timing if none is specified
         fT0.push_back(t0[std::min(iNext, t0.size() - 1U)]);
         fT1.push_back(t1[std::min(iNext, t1.size() - 1U)]);
 
         dumpNuclideSettings(mf::LogVerbatim("RadioGen"), iVolFirst + iVolInstance, volName);
       }
       ++iNext;
     } // for
 
     // check for consistency of vector sizes
     unsigned int nsize = fNuclide.size();
     if (nsize == 0) {
       throw art::Exception(art::errors::Configuration) << "No nuclide configured!\n";
     }
 
     if (fMaterial.size() != nsize)
       throw cet::exception("RadioGen") << "Different size Material vector and Nuclide vector\n";
     if (fBq.size() != nsize)
       throw cet::exception("RadioGen") << "Different size Bq vector and Nuclide vector\n";
     if (fT0.size() != nsize)
       throw cet::exception("RadioGen") << "Different size T0 vector and Nuclide vector\n";
     if (fT1.size() != nsize)
       throw cet::exception("RadioGen") << "Different size T1 vector and Nuclide vector\n";
     if (fX0.size() != nsize)
       throw cet::exception("RadioGen") << "Different size X0 vector and Nuclide vector\n";
     if (fY0.size() != nsize)
       throw cet::exception("RadioGen") << "Different size Y0 vector and Nuclide vector\n";
     if (fZ0.size() != nsize)
       throw cet::exception("RadioGen") << "Different size Z0 vector and Nuclide vector\n";
     if (fX1.size() != nsize)
       throw cet::exception("RadioGen") << "Different size X1 vector and Nuclide vector\n";
     if (fY1.size() != nsize)
       throw cet::exception("RadioGen") << "Different size Y1 vector and Nuclide vector\n";
     if (fZ1.size() != nsize)
       throw cet::exception("RadioGen") << "Different size Z1 vector and Nuclide vector\n";
 
     for (std::string& nuclideName : fNuclide) {
       if (nuclideName == "39Ar") { readfile("39Ar", "Argon_39.root", fPathToFile); }
       else if (nuclideName == "60Co") {
         readfile("60Co", "Cobalt_60.root", fPathToFile);
       }
       else if (nuclideName == "85Kr") {
         readfile("85Kr", "Krypton_85.root", fPathToFile);
       }
       else if (nuclideName == "40K") {
         readfile("40K", "Potassium_40.root", fPathToFile);
       }
       else if (nuclideName == "232Th") {
         readfile("232Th", "Thorium_232.root", fPathToFile);
       }
       else if (nuclideName == "238U") {
         readfile("238U", "Uranium_238.root", fPathToFile);
       }
       else if (nuclideName == "222Rn") {
         continue;
       } //Rn222 is handled separately later
       else if (nuclideName == "59Ni") {
         continue;
       } //Rn222 is handled separately later
       else if (nuclideName == "42Ar") {
         readfile(
           "42Ar_1",
           "Argon_42_1.root",
           fPathToFile); //Each possible beta decay mode of Ar42 is given it's own .root file for now.
         readfile(
           "42Ar_2",
           "Argon_42_2.root",
           fPathToFile); //This allows us to know which decay chain to follow for the dexcitation gammas.
         readfile(
           "42Ar_3",
           "Argon_42_3.root",
           fPathToFile); //The dexcitation gammas are not included in the root files as we want to
         readfile(
           "42Ar_4",
           "Argon_42_4.root",
           fPathToFile); //probabilistically simulate the correct coincident gammas, which we cannot guarantee
         readfile("42Ar_5", "Argon_42_5.root", fPathToFile); //by sampling a histogram.
         continue;
       } //Ar42  is handeled separately later
       else {
         std::string searchName = nuclideName;
         searchName += ".root";
         readfile(nuclideName, searchName, fPathToFile);
       }
     }
   }

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 480 of file RadioGen_module.cc.

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

Referenced by RadioGen().

   {
 
     auto const& geom = *(lar::providerFrom<geo::Geometry>());
 
     std::vector<geo::GeoNodePath> volumePaths;
     auto findVolume = [&volumePaths, volumeName](auto& path) {
       if (path.current().GetVolume()->GetName() == volumeName) volumePaths.push_back(path);
       return true;
     };
 
     geo::ROOTGeometryNavigator navigator{*(geom.ROOTGeoManager())};
 
     navigator.apply(findVolume);
 
     for (geo::GeoNodePath const& path : volumePaths) {
 
       //
       // find the coordinates of the volume in local coordinates
       //
       TGeoShape const* pShape = path.current().GetVolume()->GetShape();
       auto pBox = dynamic_cast<TGeoBBox const*>(pShape);
       if (!pBox) {
         throw cet::exception("RadioGen") << "Volume '" << path.current().GetName() << "' is a "
                                          << pShape->IsA()->GetName() << ", not a TGeoBBox.\n";
       }
 
       geo::Point_t const origin = geo::vect::makeFromCoords<geo::Point_t>(pBox->GetOrigin());
       geo::Vector_t const diag = {pBox->GetDX(), pBox->GetDY(), pBox->GetDZ()};
 
       //
       // convert to world coordinates
       //
 
       auto const trans = path.currentTransformation<geo::TransformationMatrix>();
 
       geo::Point_t min, max;
       trans.Transform(origin - diag, min);
       trans.Transform(origin + diag, max);
 
       //
       // add to the coordinates
       //
       fX0.push_back(std::min(min.X(), max.X()));
       fY0.push_back(std::min(min.Y(), max.Y()));
       fZ0.push_back(std::min(min.Z(), max.Z()));
       fX1.push_back(std::max(min.X(), max.X()));
       fY1.push_back(std::max(min.Y(), max.Y()));
       fZ1.push_back(std::max(min.Z(), max.Z()));
 
     } // for
 
     return volumePaths.size();
   } // RadioGen::addvolume()

void evgen::RadioGen::Ar42Gamma2 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> & v_prods )

private

Definition at line 933 of file RadioGen_module.cc.

References dirCalc().

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

   {
     ti_PDGID pdgid = 22;
     td_Mass m = 0.0;                       //we are writing gammas
     std::vector<double> vd_p = {.0015246}; //Momentum in GeV
     for (auto p : vd_p) {
       v_prods.emplace_back(pdgid, m, dirCalc(p, m));
     }
   }

void evgen::RadioGen::Ar42Gamma3 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> & v_prods )

private

Definition at line 943 of file RadioGen_module.cc.

References Ar42Gamma2(), and dirCalc().

Referenced by SampleOne().

   {
     ti_PDGID pdgid = 22;
     td_Mass m = 0.0; //we are writing gammas
     std::vector<double> vd_p = {.0003126};
     for (auto p : vd_p) {
       v_prods.emplace_back(pdgid, m, dirCalc(p, m));
     }
     Ar42Gamma2(v_prods);
   }

void evgen::RadioGen::Ar42Gamma4 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> & v_prods )

private

Definition at line 954 of file RadioGen_module.cc.

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

Referenced by Ar42Gamma5(), and SampleOne().

   {
     CLHEP::RandFlat flat(fEngine);
     ti_PDGID pdgid = 22;
     td_Mass m = 0.0; //we are writing gammas
     double chan1 = (0.052 / (0.052 + 0.020));
     if (flat.fire() < chan1) {
       std::vector<double> vd_p = {.0008997}; //Momentum in GeV
       for (auto p : vd_p) {
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       }
       Ar42Gamma2(v_prods);
     }
     else {
       std::vector<double> vd_p = {.0024243}; //Momentum in GeV
       for (auto p : vd_p) {
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       }
     }
   }

void evgen::RadioGen::Ar42Gamma5 ( std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> & v_prods )

private

Definition at line 975 of file RadioGen_module.cc.

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

Referenced by SampleOne().

   {
     CLHEP::RandFlat flat(fEngine);
     ti_PDGID pdgid = 22;
     td_Mass m = 0.0; //we are writing gammas
     double chan1 = (0.0033 / (0.0033 + 0.0201 + 0.041));
     double chan2 = (0.0201 / (0.0033 + 0.0201 + 0.041));
     double chanPick = flat.fire();
     if (chanPick < chan1) {
       std::vector<double> vd_p = {0.000692, 0.001228}; //Momentum in GeV
       for (auto p : vd_p) {
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       }
       Ar42Gamma2(v_prods);
     }
     else if (chanPick < (chan1 + chan2)) {
       std::vector<double> vd_p = {0.0010212}; //Momentum in GeV
       for (auto p : vd_p) {
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       }
       Ar42Gamma4(v_prods);
     }
     else {
       std::vector<double> vd_p = {0.0019208}; //Momentum in GeV
       for (auto p : vd_p) {
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       }
       Ar42Gamma2(v_prods);
     }
   }

void evgen::RadioGen::beginRun ( art::Run & run )

privatevirtual

Reimplemented from art::EDProducer.

Definition at line 454 of file RadioGen_module.cc.

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

   {
     art::ServiceHandle<geo::Geometry const> geo;
     run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
   }

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().

   {
     return collector_.consumes<T, BT>(tag);
   }

ConsumesCollector & art::ModuleBase::consumesCollector ( )

protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

   {
     return collector_;
   }

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().

   {
     collector_.consumesMany<T, BT>();
   }

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().

   {
     return collector_.consumesView<T, BT>(tag);
   }

std::pair< double, double > evgen::RadioGen::defaulttimewindow ( ) const

private

Returns the start and end of the readout window.

Definition at line 1043 of file RadioGen_module.cc.

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

Referenced by RadioGen().

   {
     // values are in simulation time scale (and therefore in [ns])
     using namespace detinfo::timescales;
 
     auto const& timings = detinfo::makeDetectorTimings(
       art::ServiceHandle<detinfo::DetectorClocksService>()->DataForJob());
     detinfo::DetectorPropertiesData const& detInfo =
       art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataForJob();
 
     //
     // we take a number of (TPC electronics) ticks before the trigger time,
     // and we go to a number of ticks after the trigger time;
     // that shift is one readout window size
     //
 
     auto const trigTimeTick = timings.toTick<electronics_tick>(timings.TriggerTime());
     electronics_time_ticks const beforeTicks{-static_cast<int>(detInfo.ReadOutWindowSize())};
     electronics_time_ticks const afterTicks{detInfo.NumberTimeSamples()};
 
     return {double(timings.toTimeScale<simulation_time>(trigTimeTick + beforeTicks)),
             double(timings.toTimeScale<simulation_time>(trigTimeTick + afterTicks))};
   } // RadioGen::defaulttimewindow()

TLorentzVector evgen::RadioGen::dirCalc	(	double	p,
		double	m
	)

private

Definition at line 678 of file RadioGen_module.cc.

References fEngine.

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

   {
     CLHEP::RandFlat flat(fEngine);
     // isotropic production angle for the decay product
     double costheta = (2.0 * flat.fire() - 1.0);
     if (costheta < -1.0) costheta = -1.0;
     if (costheta > 1.0) costheta = 1.0;
     double const sintheta = sqrt(1.0 - costheta * costheta);
     double const phi = 2.0 * M_PI * flat.fire();
     return TLorentzVector{p * sintheta * std::cos(phi),
                           p * sintheta * std::sin(phi),
                           p * costheta,
                           std::sqrt(p * p + m * m)};
   }

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().

   {
     setupQueues(resources);
     ProcessingFrame const frame{ScheduleID{}};
     beginJobWithFrame(frame);
   }

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().

   {
     auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

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().

   {
     auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

void art::detail::Producer::doEndJob ( )

inherited

Definition at line 30 of file Producer.cc.

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

   {
     ProcessingFrame const frame{ScheduleID{}};
     endJobWithFrame(frame);
   }

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().

   {
     auto r = rp.makeRun(mc, rp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

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().

   {
     auto sr = srp.makeSubRun(mc, srp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

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().

   {
     auto e = ep.makeEvent(mc);
     ++counts_run;
     ProcessingFrame const frame{mc.scheduleID()};
     produceWithFrame(e, frame);
     e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
     ++counts_passed;
     return true;
   }

void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const & fb )

inherited

Definition at line 44 of file Producer.cc.

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

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const & fb )

inherited

Definition at line 58 of file Producer.cc.

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

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseOutputFilesWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const & fb )

inherited

Definition at line 37 of file Producer.cc.

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

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const & fb )

inherited

Definition at line 51 of file Producer.cc.

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

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenOutputFilesWithFrame(fb, frame);
   }

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 1029 of file RadioGen_module.cc.

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

Referenced by RadioGen().

   {
 
     out << "[#" << iNucl << "]  " << fNuclide[iNucl] << " (" << fBq[iNucl] << " Bq/cm^3)"
         << " in " << fMaterial[iNucl] << " from " << fT0[iNucl] << " to " << fT1[iNucl]
         << " ns in ( " << fX0[iNucl] << ", " << fY0[iNucl] << ", " << fZ0[iNucl] << ") to ( "
         << fX1[iNucl] << ", " << fY1[iNucl] << ", " << fZ1[iNucl] << ")";
     if (!volumeName.empty()) out << " (\"" << volumeName << "\")";
 
   } // RadioGen::dumpNuclideSettings()

void art::Modifier::fillProductDescriptions ( )

inherited

Definition at line 10 of file Modifier.cc.

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

   {
     ProductRegistryHelper::fillDescriptions(moduleDescription());
   }

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().

   {
     return collector_.getConsumables();
   }

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.

   {
     return doMakeWorker(wp);
   }

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().

   {
     return collector_.mayConsume<T, BT>(tag);
   }

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().

   {
     collector_.mayConsumeMany<T, BT>();
   }

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().

   {
     return collector_.mayConsumeView<T, BT>(tag);
   }

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().

   {
     if (md_.has_value()) {
       return *md_;
     }
 
     throw Exception{errors::LogicError,
                     "There was an error while calling moduleDescription().\n"}
       << "The moduleDescription() base-class member function cannot be called\n"
          "during module construction.  To determine which module is "
          "responsible\n"
          "for calling it, find the '<module type>:<module "
          "label>@Construction'\n"
          "tag in the message prefix above.  Please contact artists@fnal.gov\n"
          "for guidance.\n";
   }

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 461 of file RadioGen_module.cc.

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

   {
     std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
 
     simb::MCTruth truth;
     truth.SetOrigin(simb::kSingleParticle);
 
     trackidcounter = -1;
     for (unsigned int i = 0; i < fNuclide.size(); ++i) {
       SampleOne(i, truth);
     } //end loop over nuclides
 
     MF_LOG_DEBUG("RadioGen") << truth;
     truthcol->push_back(truth);
     evt.put(std::move(truthcol));
   }

void evgen::RadioGen::readfile	(	std::string	nuclide,
		std::string const &	filename,
		std::string const &	PathToFile
	)

private

Definition at line 695 of file RadioGen_module.cc.

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

Referenced by RadioGen().

   {
     bool found{false};
     std::regex const re_argon{"42Ar.*"};
     for (size_t i = 0; i < fNuclide.size(); i++) {
       if (
         fNuclide[i] ==
         nuclide) { //This check makes sure that the nuclide we are searching for is in fact in our fNuclide list. Ar42 handeled separately.
         found = true;
         break;
       } //End If nuclide is in our list. Next is the special case of Ar42
       else if (std::regex_match(nuclide, re_argon) && fNuclide[i] == "42Ar") {
         found = true;
         break;
       }
     }
     if (!found) return;
 
     Bool_t addStatus = TH1::AddDirectoryStatus();
     TH1::AddDirectory(
       kFALSE); // cloned histograms go in memory, and aren't deleted when files are closed.
     // be sure to restore this state when we're out of the routine.
 
     spectrumname.push_back(nuclide);
     cet::search_path sp("FW_SEARCH_PATH");
     std::string fn2 = PathToFile;
     fn2 += filename;
     std::string fullname;
     sp.find_file(fn2, fullname);
     struct stat sb;
     if (fullname.empty() || stat(fullname.c_str(), &sb) != 0)
       throw cet::exception("RadioGen")
         << "Input spectrum file " << fn2 << " not found in FW_SEARCH_PATH!\n";
 
     TFile f(fullname.c_str(), "READ");
     TGraph* alphagraph = (TGraph*)f.Get("Alphas");
     TGraph* betagraph = (TGraph*)f.Get("Betas");
     TGraph* gammagraph = (TGraph*)f.Get("Gammas");
     TGraph* neutrongraph = (TGraph*)f.Get("Neutrons");
 
     if (alphagraph) {
       int np = alphagraph->GetN();
       double* y = alphagraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Alpha";
       auto alphahist = std::make_unique<TH1D>(name.c_str(), "Alpha Spectrum", np, 0, np);
       for (int i = 0; i < np; i++) {
         alphahist->SetBinContent(i + 1, y[i]);
         alphahist->SetBinError(i + 1, 0);
       }
       alphaintegral.push_back(alphahist->Integral());
       alphaspectrum.push_back(move(alphahist));
     }
     else {
       alphaintegral.push_back(0);
       alphaspectrum.push_back(nullptr);
     }
 
     if (betagraph) {
       int np = betagraph->GetN();
 
       double* y = betagraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Beta";
       auto betahist = std::make_unique<TH1D>(name.c_str(), "Beta Spectrum", np, 0, np);
       for (int i = 0; i < np; i++) {
         betahist->SetBinContent(i + 1, y[i]);
         betahist->SetBinError(i + 1, 0);
       }
       betaintegral.push_back(betahist->Integral());
       betaspectrum.push_back(move(betahist));
     }
     else {
       betaintegral.push_back(0);
       betaspectrum.push_back(nullptr);
     }
 
     if (gammagraph) {
       int np = gammagraph->GetN();
       double* y = gammagraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Gamma";
       auto gammahist = std::make_unique<TH1D>(name.c_str(), "Gamma Spectrum", np, 0, np);
       for (int i = 0; i < np; i++) {
         gammahist->SetBinContent(i + 1, y[i]);
         gammahist->SetBinError(i + 1, 0);
       }
       gammaintegral.push_back(gammahist->Integral());
       gammaspectrum.push_back(move(gammahist));
     }
     else {
       gammaintegral.push_back(0);
       gammaspectrum.push_back(nullptr);
     }
 
     if (neutrongraph) {
       int np = neutrongraph->GetN();
       double* y = neutrongraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Neutron";
       auto neutronhist = std::make_unique<TH1D>(name.c_str(), "Neutron Spectrum", np, 0, np);
       for (int i = 0; i < np; i++) {
         neutronhist->SetBinContent(i + 1, y[i]);
         neutronhist->SetBinError(i + 1, 0);
       }
       neutronintegral.push_back(neutronhist->Integral());
       neutronspectrum.push_back(move(neutronhist));
     }
     else {
       neutronintegral.push_back(0);
       neutronspectrum.push_back(nullptr);
     }
 
     f.Close();
     TH1::AddDirectory(addStatus);
 
     double total =
       alphaintegral.back() + betaintegral.back() + gammaintegral.back() + neutronintegral.back();
     if (total > 0) {
       alphaintegral.back() /= total;
       betaintegral.back() /= total;
       gammaintegral.back() /= total;
       neutronintegral.back() /= total;
     }
   }

void art::Modifier::registerProducts ( ProductDescriptions & productsToRegister )

inherited

Definition at line 16 of file Modifier.cc.

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

   {
     ProductRegistryHelper::registerProducts(productsToRegister,
                                             moduleDescription());
   }

double evgen::RadioGen::samplefromth1d ( TH1D & hist )

private

Definition at line 895 of file RadioGen_module.cc.

References fEngine, nbinsx, and x.

Referenced by samplespectrum().

   {
     CLHEP::RandFlat flat(fEngine);
 
     int nbinsx = hist.GetNbinsX();
     std::vector<double> partialsum;
     partialsum.resize(nbinsx + 1);
     partialsum[0] = 0;
 
     for (int i = 1; i <= nbinsx; i++) {
       double hc = hist.GetBinContent(i);
       if (hc < 0)
         throw cet::exception("RadioGen") << "Negative bin:  " << i << " " << hist.GetName() << "\n";
       partialsum[i] = partialsum[i - 1] + hc;
     }
     double integral = partialsum[nbinsx];
     if (integral == 0) return 0;
     // normalize to unit sum
     for (int i = 1; i <= nbinsx; i++)
       partialsum[i] /= integral;
 
     double r1 = flat.fire();
     int ibin = TMath::BinarySearch(nbinsx, &(partialsum[0]), r1);
     Double_t x = hist.GetBinLowEdge(ibin + 1);
     if (r1 > partialsum[ibin]) {
       x += hist.GetBinWidth(ibin + 1) * (r1 - partialsum[ibin]) /
            (partialsum[ibin + 1] - partialsum[ibin]);
     }
     return x;
   }

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 538 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().

   {
     art::ServiceHandle<geo::Geometry const> geo;
     TGeoManager* geomanager = geo->ROOTGeoManager();
 
     CLHEP::RandFlat flat(fEngine);
     CLHEP::RandPoisson poisson(fEngine);
 
     // figure out how many decays to generate, assuming that the entire prism consists of the radioactive material.
     // we will skip over decays in other materials later.
 
     double rate =
       fabs(fBq[i] * (fT1[i] - fT0[i]) * (fX1[i] - fX0[i]) * (fY1[i] - fY0[i]) * (fZ1[i] - fZ0[i])) /
       1.0E9;
     long ndecays = poisson.shoot(rate);
 
     std::regex const re_material{fMaterial[i]};
     for (unsigned int idecay = 0; idecay < ndecays; idecay++) {
       // generate just one particle at a time.  Need a little recoding if a radioactive
       // decay generates multiple daughters that need simulation
       // uniformly distributed in position and time
       //
       // JStock: Leaving this as a single position for the decay products. For now I will assume they all come from the same spot.
       TLorentzVector pos(
         fX0[i] + flat.fire() * (fX1[i] - fX0[i]),
         fY0[i] + flat.fire() * (fY1[i] - fY0[i]),
         fZ0[i] + flat.fire() * (fZ1[i] - fZ0[i]),
         (idecay == 0 && fIsFirstSignalSpecial) ? 0 : (fT0[i] + flat.fire() * (fT1[i] - fT0[i])));
 
       // discard decays that are not in the proper material
       std::string volmaterial =
         geomanager->FindNode(pos.X(), pos.Y(), pos.Z())->GetMedium()->GetMaterial()->GetName();
       if (!std::regex_match(volmaterial, re_material)) continue;
 
       //Moved pdgid into the next statement, so that it is localized.
       // electron=11, photon=22, alpha = 1000020040, neutron = 2112
 
       //JStock: Allow us to have different particles from the same decay. This requires multiple momenta.
       std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>
         v_prods; //(First is for PDGID, second is mass, third is Momentum)
 
       if (fNuclide[i] == "222Rn") // Treat 222Rn separately
       {
         double p = 0;
         double t = 0.00548952;
         td_Mass m = m_alpha;
         ti_PDGID pdgid = 1000020040; //td_Mass = double. ti_PDGID = int;
         double energy = t + m;
         double p2 = energy * energy - m * m;
         if (p2 > 0)
           p = TMath::Sqrt(p2);
         else
           p = 0;
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       } //End special case RN222
       else if (
         fNuclide[i] ==
         "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.
         double p = 0.008997;
         td_Mass m = 0;
         ti_PDGID pdgid =
           22; // td_Mas=double. ti_PDFID=int. Assigning p directly, as t=p for gammas.
         v_prods.emplace_back(pdgid, m, dirCalc(p, m));
       }                                 //end special case Ni59 calibration source
       else if (fNuclide[i] == "42Ar") { // Spot for special treatment of Ar42.
         double p = 0;
         double t = 0;
         td_Mass m = 0;
         ti_PDGID pdgid = 0;           //td_Mass = double. ti_PDGID = int;
         double bSelect = flat.fire(); //Make this a random number from 0 to 1.
         if (bSelect < 0.819) {        //beta channel 1. No Gamma. beta Q value 3525.22 keV
           samplespectrum("42Ar_1", pdgid, t, m, p);
           v_prods.emplace_back(pdgid, m, dirCalc(p, m));
           //No gamma here.
         }
         else if (bSelect < 0.9954) { //beta channel 2. 1 Gamma (1524.6 keV). beta Q value 2000.62
           samplespectrum("42Ar_2", pdgid, t, m, p);
           v_prods.emplace_back(pdgid, m, dirCalc(p, m));
           Ar42Gamma2(v_prods);
         }
         else if (
           bSelect <
           0.9988) { //beta channel 3. 1 Gamma Channel. 312.6 keV + gamma 2. beta Q value 1688.02 keV
           samplespectrum("42Ar_3", pdgid, t, m, p);
           v_prods.emplace_back(pdgid, m, dirCalc(p, m));
           Ar42Gamma3(v_prods);
         }
         else if (
           bSelect <
           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
           samplespectrum("42Ar_4", pdgid, t, m, p);
           v_prods.emplace_back(pdgid, m, dirCalc(p, m));
           Ar42Gamma4(v_prods);
         }
         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
           samplespectrum("42Ar_5", pdgid, t, m, p);
           v_prods.emplace_back(pdgid, m, dirCalc(p, m));
           Ar42Gamma5(v_prods);
         }
         //Add beta.
         //Call gamma function for beta mode.
       }
       else { //General Case.
         double p = 0;
         double t = 0;
         td_Mass m = 0;
         ti_PDGID pdgid = 0; //td_Mass = double. ti_PDGID = int;
         samplespectrum(fNuclide[i], pdgid, t, m, p);
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom =
           std::make_tuple(pdgid, m, dirCalc(p, m));
         v_prods.push_back(partMassMom);
       } //end else (not RN or other special case
 
       //JStock: Modify this to now loop over the v_prods.
       for (auto prodEntry : v_prods) {
         // set track id to a negative serial number as these are all primary particles and have id <= 0
         int trackid = trackidcounter;
         ti_PDGID pdgid = std::get<0>(prodEntry);
         td_Mass m = std::get<1>(prodEntry);
         TLorentzVector pvec = std::get<2>(prodEntry);
         trackidcounter--;
         std::string primary("primary");
 
         // alpha particles need a little help since they're not in the TDatabasePDG table
         // // so don't rely so heavily on default arguments to the MCParticle constructor
         if (pdgid == 1000020040) {
           simb::MCParticle part(trackid, pdgid, primary, -1, m, 1);
           part.AddTrajectoryPoint(pos, pvec);
           mct.Add(part);
         } // end "If alpha"
         else {
           simb::MCParticle part(trackid, pdgid, primary);
           part.AddTrajectoryPoint(pos, pvec);
           mct.Add(part);
         } // end All standard cases.
       }   //End Loop over all particles produces in this single decay.
     }
   }

void evgen::RadioGen::samplespectrum	(	std::string	nuclide,
		int &	itype,
		double &	t,
		double &	m,
		double &	p
	)

private

Definition at line 831 of file RadioGen_module.cc.

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

Referenced by SampleOne().

   {
     CLHEP::RandFlat flat(fEngine);
 
     int inuc = -1;
     for (size_t i = 0; i < spectrumname.size(); i++) {
       if (nuclide == spectrumname[i]) {
         inuc = i;
         break;
       }
     }
     if (inuc == -1) {
       t = 0; // throw an exception in the future
       itype = 0;
       throw cet::exception("RadioGen") << "Ununderstood nuclide:  " << nuclide << "\n";
     }
 
     double rtype = flat.fire();
 
     itype = -1;
     m = 0;
     p = 0;
     for (
       int itry = 0; itry < 10;
       itry++) // maybe a tiny normalization issue with a sum of 0.99999999999 or something, so try a few times.
     {
       if (rtype <= alphaintegral[inuc] && alphaspectrum[inuc] != nullptr) {
         itype = 1000020040; // alpha
         m = m_alpha;
         t = samplefromth1d(*alphaspectrum[inuc]) / 1000000.0;
       }
       else if (rtype <= alphaintegral[inuc] + betaintegral[inuc] && betaspectrum[inuc] != nullptr) {
         itype = 11; // beta
         m = m_e;
         t = samplefromth1d(*betaspectrum[inuc]) / 1000000.0;
       }
       else if (rtype <= alphaintegral[inuc] + betaintegral[inuc] + gammaintegral[inuc] &&
                gammaspectrum[inuc] != nullptr) {
         itype = 22; // gamma
         m = 0;
         t = samplefromth1d(*gammaspectrum[inuc]) / 1000000.0;
       }
       else if (neutronspectrum[inuc] != nullptr) {
         itype = 2112;
         m = m_neutron;
         t = samplefromth1d(*neutronspectrum[inuc]) / 1000000.0;
       }
       if (itype >= 0) break;
     }
     if (itype == -1) {
       throw cet::exception("RadioGen")
         << "Normalization problem with nuclide:  " << nuclide << "\n";
     }
     double e = t + m;
     p = e * e - m * m;
     if (p >= 0) { p = TMath::Sqrt(p); }
     else {
       p = 0;
     }
   }

void art::ModuleBase::setModuleDescription ( ModuleDescription const & md )

inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

   {
     md_ = md;
   }

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().

   {
     // Now that we know we have seen all the consumes declarations,
     // sort the results for fast lookup later.
     collector_.sortConsumables(current_process_name);
   }