Inheritance diagram for evgen::RadioGen:

Public Types
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

template<typename UserConfig , typename KeysToIgnore = void>
using	Table = ProducerBase::Table< UserConfig, KeysToIgnore >

Public Member Functions
	RadioGen (fhicl::ParameterSet const &pset)

virtual	~RadioGen ()

void	produce (art::Event &evt)

void	beginRun (art::Run &run)

void	reconfigure (fhicl::ParameterSet const &p)

template<typename PROD , BranchType B = InEvent>
ProductID	getProductID (std::string const &instanceName={}) const

template<typename PROD , BranchType B>
ProductID	getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const

bool	modifiesEvent () const

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename T , art::BranchType BT>
art::ProductToken< T >	consumes (InputTag const &it)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , art::BranchType BT>
art::ViewToken< T >	consumesView (InputTag const &it)

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename T , art::BranchType BT>
art::ProductToken< T >	mayConsume (InputTag const &it)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , art::BranchType BT>
art::ViewToken< T >	mayConsumeView (InputTag const &it)

base_engine_t &	createEngine (seed_t seed)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)

seed_t	get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)

Static Public Member Functions
static cet::exempt_ptr< Consumer >	non_module_context ()

Protected Member Functions
CurrentProcessingContext const *	currentContext () const

void	validateConsumedProduct (BranchType const bt, ProductInfo const &pi)

void	prepareForJob (fhicl::ParameterSet const &pset)

void	showMissingConsumes () const

Private Types
typedef int	ti_PDGID

typedef double	td_Mass

Private Member Functions
void	SampleOne (unsigned int i, simb::MCTruth &mct)

TLorentzVector	dirCalc (double p, double m)

void	readfile (std::string nuclide, std::string 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)

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

int	trackidcounter
	Serial number for the MC track ID. More...

const double	m_e = 0.000510998928

const double	m_alpha = 3.727379240

const double	m_neutron = 0.9395654133

std::vector< std::string >	spectrumname

std::vector< TH1D * >	alphaspectrum

std::vector< double >	alphaintegral

std::vector< TH1D * >	betaspectrum

std::vector< double >	betaintegral

std::vector< TH1D * >	gammaspectrum

std::vector< double >	gammaintegral

std::vector< TH1D * >	neutronspectrum

std::vector< double >	neutronintegral

Detailed Description

Module to generate particles created by radiological decay, patterend off of SingleGen Currently it generates only in rectangular prisms oriented along the x,y,z axes

Definition at line 97 of file RadioGen_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 43 of file EDProducer.h.

typedef double evgen::RadioGen::td_Mass

private

Definition at line 111 of file RadioGen_module.cc.

typedef int evgen::RadioGen::ti_PDGID

private

Definition at line 110 of file RadioGen_module.cc.

using art::EDProducer::WorkerType = WorkerT<EDProducer>

inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 177 of file RadioGen_module.cc.

   {
 
     this->reconfigure(pset);
 
     // create a default random engine; obtain the random seed from NuRandomService,
     // unless overridden in configuration with key "Seed"
     art::ServiceHandle<rndm::NuRandomService>()
       ->createEngine(*this, pset, "Seed");
 
     produces< std::vector<simb::MCTruth> >();
     produces< sumdata::RunData, art::InRun >();
 
   }

evgen::RadioGen::~RadioGen ( )

virtual

Definition at line 193 of file RadioGen_module.cc.

194 {

195 }

Member Function Documentation

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

private

Definition at line 692 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){
       std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
       v_prods.push_back(prods);
     }
   }

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

private

Definition at line 700 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){
       std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
       v_prods.push_back(prods);
     }
     Ar42Gamma2(v_prods);    
   }

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

private

Definition at line 709 of file RadioGen_module.cc.

References art::RandomNumberGenerator::getEngine().

                                                                                           {
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine();
     CLHEP::RandFlat     flat(engine);
     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){
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(prods);
       }
       Ar42Gamma2(v_prods);
     }else{
       std::vector<double> vd_p = {.0024243};//Momentum in GeV
       for(auto p : vd_p){
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(prods);
       }
     }
   }

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

private

Definition at line 730 of file RadioGen_module.cc.

References DEFINE_ART_MODULE, and art::RandomNumberGenerator::getEngine().

                                                                                           {
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine();
     CLHEP::RandFlat     flat(engine);
     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){
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(prods);
       }
       Ar42Gamma2(v_prods);
     }else if (chanPick<(chan1+chan2)){
       std::vector<double> vd_p = {0.0010212};//Momentum in GeV
       for(auto p : vd_p){
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(prods);
       }
       Ar42Gamma4(v_prods);
     }else{
       std::vector<double> vd_p = {0.0019208};//Momentum in GeV
       for(auto p : vd_p){
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> prods = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(prods);
       }
       Ar42Gamma2(v_prods);
     }
   }

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

virtual

Reimplemented from art::EDProducer.

Definition at line 257 of file RadioGen_module.cc.

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

   {
 
     // grab the geometry object to see what geometry we are using
 
     art::ServiceHandle<geo::Geometry> geo;
     std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetectorName()));
     run.put(std::move(runcol));
 
     return;
   }

template<typename T , BranchType = InEvent>

ProductToken<T> art::Consumer::consumes ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::consumes ( InputTag const & it )

inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::consumesMany ( )

inherited

Definition at line 162 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::Consumer::consumesView ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::consumesView ( InputTag const & it )

inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t seed )

inherited

Definition at line 26 of file EngineCreator.cc.

References art::EngineCreator::rng().

Referenced by evgen::CosmicsGen::CosmicsGen(), rndm::NuRandomService::createEngine(), cluster::fuzzyCluster::fuzzyCluster(), cluster::HoughLineFinder::HoughLineFinder(), art::MixFilter< T >::initEngine_(), larg4::LArG4::LArG4(), evgen::LightSource::LightSource(), evgen::NeutronOsc::NeutronOsc(), evgen::NucleonDecay::NucleonDecay(), opdet::OpMCDigi::OpMCDigi(), opdet::OptDetDigitizer::OptDetDigitizer(), phot::PhotonLibraryPropagation::PhotonLibraryPropagation(), detsim::SimDriftElectrons::SimDriftElectrons(), evgen::SingleGen::SingleGen(), evgen::SNNueAr40CCGen::SNNueAr40CCGen(), ToyOneShowerGen::ToyOneShowerGen(), and trkf::Track3DKalman::Track3DKalman().

 {
   return rng()->createEngine(placeholder_schedule_id(), seed);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make
	)

inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make,
		label_t const &	engine_label
	)

inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
 }

CurrentProcessingContext const * art::EDProducer::currentContext ( ) const

protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

   {
     return current_context_.get();
   }

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

private

Definition at line 416 of file RadioGen_module.cc.

References art::RandomNumberGenerator::getEngine().

                                                     {
       art::ServiceHandle<art::RandomNumberGenerator> rng;
       CLHEP::HepRandomEngine &engine = rng->getEngine();
       CLHEP::RandFlat  flat(engine);
       // 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 sintheta = sqrt(1.0-costheta*costheta);
       double phi = 2.0*M_PI*flat.fire();
         TLorentzVector pvec(p*sintheta*std::cos(phi),
           p*sintheta*std::sin(phi),
           p*costheta,
           std::sqrt(p*p+m*m));
       return pvec;
   }

EngineCreator::seed_t EngineCreator::get_seed_value	(	fhicl::ParameterSet const &	pset,
		char const	key[] = `"seed"`,
		seed_t const	implicit_seed = `-1`
	)

inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

 {
   auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
   return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
 }

template<typename PROD , BranchType B>

ProductID art::EDProducer::getProductID ( std::string const & instanceName = {} ) const

inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

   {
     return ProducerBase::getProductID<PROD, B>(moduleDescription_,
                                                instanceName);
   }

template<typename PROD , BranchType B>

ProductID art::ProducerBase::getProductID	(	ModuleDescription const &	moduleDescription,
		std::string const &	instanceName
	)		const

inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

   {
     auto const& pd =
       get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
     return pd.productID();
   }

template<typename T , BranchType = InEvent>

ProductToken<T> art::Consumer::mayConsume ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::mayConsume ( InputTag const & it )

inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::mayConsumeMany ( )

inherited

Definition at line 205 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const & it )

inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

bool art::ProducerBase::modifiesEvent ( ) const

inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

     {
       return true;
     }

cet::exempt_ptr< art::Consumer > art::Consumer::non_module_context ( )

staticinherited

Definition at line 76 of file Consumer.cc.

 {
   static Consumer invalid{InvalidTag{}};
   return &invalid;
 }

void art::Consumer::prepareForJob ( fhicl::ParameterSet const & pset )

protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

 {
   if (!moduleContext_)
     return;
 
   pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
   for (auto& consumablesPerBranch : consumables_) {
     cet::sort_all(consumablesPerBranch);
   }
 }

void evgen::RadioGen::produce ( art::Event & evt )

virtual

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

Implements art::EDProducer.

Definition at line 270 of file RadioGen_module.cc.

References simb::kSingleParticle, LOG_DEBUG, art::Event::put(), and simb::MCTruth::SetOrigin().

   {
 
     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
 
     LOG_DEBUG("RadioGen") << truth;
     truthcol->push_back(truth);
     evt.put(std::move(truthcol));
     return;
   }

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

private

Definition at line 436 of file RadioGen_module.cc.

References f, and y.

   {
     int ifound = 0;
     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.
         ifound = 1;
         break;
       } //End If nuclide is in our list. Next is the special case of Ar42
       else if ( (std::regex_match(nuclide, std::regex( "42Ar.*" )) ) && fNuclide[i]=="42Ar" ){ 
         ifound = 1;
         break;
       }
     }
     if (ifound == 0) 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";
       TH1D *alphahist = (TH1D*) new 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);
       }
       alphaspectrum.push_back(alphahist);
       alphaintegral.push_back(alphahist->Integral());
     }
     else
     {
       alphaspectrum.push_back(0);
       alphaintegral.push_back(0);
     }
     
     
     if (betagraph)
     {
       int np = betagraph->GetN();
       
       double *y = betagraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Beta";
       TH1D *betahist = (TH1D*) new 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);
       }
       betaspectrum.push_back(betahist);
       betaintegral.push_back(betahist->Integral());
     }
     else
     {
       betaspectrum.push_back(0);
       betaintegral.push_back(0);
     }
     
     if (gammagraph)
     {
       int np = gammagraph->GetN();
       double *y = gammagraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Gamma";
       TH1D *gammahist = (TH1D*) new 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);
       }
       gammaspectrum.push_back(gammahist);
       gammaintegral.push_back(gammahist->Integral());
     }
     else
     {
       gammaspectrum.push_back(0);
       gammaintegral.push_back(0);
     }
     
     if (neutrongraph)
     {
       int np = neutrongraph->GetN();
       double *y = neutrongraph->GetY();
       std::string name;
       name = "RadioGen_";
       name += nuclide;
       name += "_Neutron";
       TH1D *neutronhist = (TH1D*) new 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);
       }
       neutronspectrum.push_back(neutronhist);
       neutronintegral.push_back(neutronhist->Integral());
     }
     else
     {
       neutronspectrum.push_back(0);
       neutronintegral.push_back(0);
     }
     
     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 evgen::RadioGen::reconfigure ( fhicl::ParameterSet const & p )

Definition at line 198 of file RadioGen_module.cc.

References fhicl::ParameterSet::get().

   {
     // do not put seed in reconfigure because we don't want to reset 
     // the seed midstream -- same as SingleGen
 
     fNuclide       = p.get< std::vector<std::string>>("Nuclide");
     fMaterial      = p.get< std::vector<std::string>>("Material");
     fBq            = p.get< std::vector<double> >("BqPercc");
     fT0            = p.get< std::vector<double> >("T0");
     fT1            = p.get< std::vector<double> >("T1");
     fX0            = p.get< std::vector<double> >("X0");
     fY0            = p.get< std::vector<double> >("Y0");
     fZ0            = p.get< std::vector<double> >("Z0");
     fX1            = p.get< std::vector<double> >("X1");
     fY1            = p.get< std::vector<double> >("Y1");
     fZ1            = p.get< std::vector<double> >("Z1");
 
     // check for consistency of vector sizes
     
     unsigned int nsize = fNuclide.size();
     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 handeled separately later
       else if(nuclideName=="59Ni"){continue;} //Rn222 is handeled 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);
       }
     }
     
     return;
   }

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

private

Definition at line 655 of file RadioGen_module.cc.

References art::RandomNumberGenerator::getEngine(), and x.

   {
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine();
     CLHEP::RandFlat  flat(engine);
     
     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

Definition at line 293 of file RadioGen_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), energy, art::RandomNumberGenerator::getEngine(), part, and geo::GeometryCore::ROOTGeoManager().

                                                           {
 
     art::ServiceHandle<geo::Geometry> geo;
     TGeoManager *geomanager = geo->ROOTGeoManager(); 
 
     // get the random number generator service and make some CLHEP generators
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine();
     CLHEP::RandFlat     flat(engine);
     CLHEP::RandPoisson  poisson(engine);
 
     // 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);
 
     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]),
           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();
       //mf::LogDebug("RadioGen") << "Position: " << pos.X() << " " << pos.Y() << " " << pos.Z() << " " << pos.T() << " " << volmaterial << " " << fMaterial[i] << std::endl;
       if ( ! std::regex_match(volmaterial, std::regex(fMaterial[i])) ) continue;
       //mf::LogDebug("RadioGen") << "Decay accepted" << std::endl;
       
       //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;
         //Make TLorentzVector and push it to the back of v_prods.
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(partMassMom);
       }//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.
         std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
         v_prods.push_back(partMassMom);
       }//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);
           std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
           v_prods.push_back(partMassMom);
           //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);
           std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
           v_prods.push_back(partMassMom);
           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);
           std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
           v_prods.push_back(partMassMom);
           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);
           std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
           v_prods.push_back(partMassMom);
           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);
           std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom = std::make_tuple(pdgid, m, dirCalc(p,m));
           v_prods.push_back(partMassMom);
           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 584 of file RadioGen_module.cc.

References e, and art::RandomNumberGenerator::getEngine().

   {
     
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine();
     CLHEP::RandFlat  flat(engine);
     
     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] != 0)
       {
         itype = 1000020040; // alpha
         m = m_alpha;
         t = samplefromth1d(alphaspectrum[inuc])/1000000.0;
       }
       else if (rtype <= alphaintegral[inuc]+betaintegral[inuc] && betaspectrum[inuc] != 0)
       {
         itype = 11; // beta
         m = m_e;
         t = samplefromth1d(betaspectrum[inuc])/1000000.0;
       }
       else if ( rtype <= alphaintegral[inuc] + betaintegral[inuc] + gammaintegral[inuc] && gammaspectrum[inuc] != 0)
       {
         itype = 22; // gamma
         m = 0;
         t = samplefromth1d(gammaspectrum[inuc])/1000000.0;
       }
       else if( neutronspectrum[inuc] != 0)
       {
         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::Consumer::showMissingConsumes ( ) const

protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

 {
   if (!moduleContext_)
     return;
 
   // If none of the branches have missing consumes statements, exit early.
   if (std::all_of(cbegin(missingConsumes_),
                   cend(missingConsumes_),
                   [](auto const& perBranch) { return perBranch.empty(); }))
     return;
 
   constexpr cet::HorizontalRule rule{60};
   mf::LogPrint log{"MTdiagnostics"};
   log << '\n'
       << rule('=') << '\n'
       << "The following consumes (or mayConsume) statements are missing from\n"
       << module_context(moduleDescription_) << '\n'
       << rule('-') << '\n';
 
   cet::for_all_with_index(
     missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
       for (auto const& pi : perBranch) {
         log << "  "
             << assemble_consumes_statement(static_cast<BranchType>(i), pi)
             << '\n';
       }
     });
   log << rule('=');
 }

void art::Consumer::validateConsumedProduct	(	BranchType const	bt,
		ProductInfo const &	pi
	)

protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

 {
   // Early exits if consumes tracking has been disabled or if the
   // consumed product is an allowed consumable.
   if (!moduleContext_)
     return;
 
   if (cet::binary_search_all(consumables_[bt], pi))
     return;
 
   if (requireConsumes_) {
     throw Exception(errors::ProductRegistrationFailure,
                     "Consumer: an error occurred during validation of a "
                     "retrieved product\n\n")
       << "The following consumes (or mayConsume) statement is missing from\n"
       << module_context(moduleDescription_) << ":\n\n"
       << "  " << assemble_consumes_statement(bt, pi) << "\n\n";
   }
 
   missingConsumes_[bt].insert(pi);
 }

Member Data Documentation

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

private

Definition at line 163 of file RadioGen_module.cc.

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

private

Definition at line 162 of file RadioGen_module.cc.

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

private

Definition at line 165 of file RadioGen_module.cc.

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

private

Definition at line 164 of file RadioGen_module.cc.

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

private

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

Definition at line 141 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 140 of file RadioGen_module.cc.

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

private

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

Definition at line 139 of file RadioGen_module.cc.

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

private

Beginning of time window to simulate in ns.

Definition at line 142 of file RadioGen_module.cc.

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

private

End of time window to simulate in ns.

Definition at line 143 of file RadioGen_module.cc.

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

private

Bottom corner x position (cm) in world coordinates.

Definition at line 144 of file RadioGen_module.cc.

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

private

Top corner x position (cm) in world coordinates.

Definition at line 147 of file RadioGen_module.cc.

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

private

Bottom corner y position (cm) in world coordinates.

Definition at line 145 of file RadioGen_module.cc.

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

private

Top corner y position (cm) in world coordinates.

Definition at line 148 of file RadioGen_module.cc.

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

private

Bottom corner z position (cm) in world coordinates.

Definition at line 146 of file RadioGen_module.cc.

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

private

Top corner z position (cm) in world coordinates.

Definition at line 149 of file RadioGen_module.cc.

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

private

Definition at line 167 of file RadioGen_module.cc.

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

private

Definition at line 166 of file RadioGen_module.cc.

const double evgen::RadioGen::m_alpha = 3.727379240

private

Definition at line 158 of file RadioGen_module.cc.

const double evgen::RadioGen::m_e = 0.000510998928

private

Definition at line 157 of file RadioGen_module.cc.

const double evgen::RadioGen::m_neutron = 0.9395654133

private

Definition at line 159 of file RadioGen_module.cc.

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

private

Definition at line 169 of file RadioGen_module.cc.

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

private

Definition at line 168 of file RadioGen_module.cc.

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

private

Definition at line 161 of file RadioGen_module.cc.

int evgen::RadioGen::trackidcounter

private

Serial number for the MC track ID.

Definition at line 150 of file RadioGen_module.cc.

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

larsim/v06_53_00/source/larsim/EventGenerator/RadioGen_module.cc

Public Types

Public Member Functions

Static Public Member Functions

Protected Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation