evgen::RadioGen Class Reference

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

Inheritance diagram for evgen::RadioGen:

Public Member Functions
	RadioGen (fhicl::ParameterSet const &pset)
Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDProducer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept
	Producer (fhicl::ParameterSet const &)
	Producer (Producer const &)=delete
	Producer (Producer &&)=delete
Producer &	operator= (Producer const &)=delete
Producer &	operator= (Producer &&)=delete
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)
Public Member Functions inherited from art::Modifier
	~Modifier () noexcept
	Modifier ()
	Modifier (Modifier const &)=delete
	Modifier (Modifier &&)=delete
Modifier &	operator= (Modifier const &)=delete
Modifier &	operator= (Modifier &&)=delete
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
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)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

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)
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::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

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >
Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)
Protected Member Functions inherited from art::ModuleBase
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 ()

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

Member Typedef Documentation

typedef double evgen::RadioGen::td_Mass

private

Definition at line 195 of file RadioGen_module.cc.

typedef int evgen::RadioGen::ti_PDGID

private

Definition at line 194 of file RadioGen_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 280 of file RadioGen_module.cc.

    : 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>{}->createEngine(*this, pset, "Seed"))
  {
    produces< std::vector<simb::MCTruth> >();
    produces< sumdata::RunData, art::InRun >();
    
    
    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 (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")    ;}
      else if(nuclideName=="60Co" ){readfile("60Co","Cobalt_60.root")   ;}
      else if(nuclideName=="85Kr" ){readfile("85Kr","Krypton_85.root")  ;}
      else if(nuclideName=="40K"  ){readfile("40K","Potassium_40.root") ;}
      else if(nuclideName=="232Th"){readfile("232Th","Thorium_232.root");}
      else if(nuclideName=="238U" ){readfile("238U","Uranium_238.root") ;}
      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"); //Each possible beta decay mode of Ar42 is given it's own .root file for now.
        readfile("42Ar_2", "Argon_42_2.root"); //This allows us to know which decay chain to follow for the dexcitation gammas.
        readfile("42Ar_3", "Argon_42_3.root"); //The dexcitation gammas are not included in the root files as we want to
        readfile("42Ar_4", "Argon_42_4.root"); //probabilistically simulate the correct coincident gammas, which we cannot guarantee
        readfile("42Ar_5", "Argon_42_5.root"); //by sampling a histogram.
        continue;
      } //Ar42  is handeled separately later
      else{
        std::string searchName = nuclideName;
        searchName+=".root";
        readfile(nuclideName,searchName);
      }
    }
  }

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

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

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

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

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

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

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

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

private

Returns the start and end of the readout window.

Definition at line 979 of file RadioGen_module.cc.

                                                            {
    // 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 606 of file RadioGen_module.cc.

  {
    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)};
  }

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

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

  {
    ///unique_ptr allows ownership to be transferred to the art::Event after the put statement
    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  )

private

Definition at line 623 of file RadioGen_module.cc.

  {
    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 = "Radionuclides/";
    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;
    }
  }

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

private

Definition at line 835 of file RadioGen_module.cc.

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

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

  {
    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; }
  }

Member Data Documentation

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

private

Definition at line 259 of file RadioGen_module.cc.

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

private

Definition at line 258 of file RadioGen_module.cc.

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

private

Definition at line 261 of file RadioGen_module.cc.

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

private

Definition at line 260 of file RadioGen_module.cc.

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

private

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

Definition at line 240 of file RadioGen_module.cc.

CLHEP::HepRandomEngine& evgen::RadioGen::fEngine

private

Definition at line 266 of file RadioGen_module.cc.

bool evgen::RadioGen::fIsFirstSignalSpecial

private

Definition at line 249 of file RadioGen_module.cc.

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

private

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

Definition at line 239 of file RadioGen_module.cc.

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

private

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

Definition at line 238 of file RadioGen_module.cc.

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

private

Beginning of time window to simulate in ns.

Definition at line 241 of file RadioGen_module.cc.

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

private

End of time window to simulate in ns.

Definition at line 242 of file RadioGen_module.cc.

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

private

Bottom corner x position (cm) in world coordinates.

Definition at line 243 of file RadioGen_module.cc.

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

private

Top corner x position (cm) in world coordinates.

Definition at line 246 of file RadioGen_module.cc.

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

private

Bottom corner y position (cm) in world coordinates.

Definition at line 244 of file RadioGen_module.cc.

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

private

Top corner y position (cm) in world coordinates.

Definition at line 247 of file RadioGen_module.cc.

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

private

Bottom corner z position (cm) in world coordinates.

Definition at line 245 of file RadioGen_module.cc.

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

private

Top corner z position (cm) in world coordinates.

Definition at line 248 of file RadioGen_module.cc.

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

private

Definition at line 263 of file RadioGen_module.cc.

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

private

Definition at line 262 of file RadioGen_module.cc.

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

private

Definition at line 265 of file RadioGen_module.cc.

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

private

Definition at line 264 of file RadioGen_module.cc.

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

private

Definition at line 257 of file RadioGen_module.cc.

int evgen::RadioGen::trackidcounter

private

Serial number for the MC track ID.

Definition at line 250 of file RadioGen_module.cc.

---

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

• larsim/larsim/EventGenerator/RadioGen_module.cc