evgen::RadioGen Class Reference

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...
bool	fIsFirstSignalSpecial
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 178 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 194 of file RadioGen_module.cc.

  {
  }

Member Function Documentation

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

private

Definition at line 696 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 704 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 713 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 734 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 259 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 419 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.

Referenced by art::RootOutput::beginSubRun(), art::OutputModule::doBeginRun(), art::OutputModule::doBeginSubRun(), art::OutputModule::doEndRun(), art::OutputModule::doEndSubRun(), art::ProducingService::doPostReadEvent(), art::ProducingService::doPostReadRun(), art::ProducingService::doPostReadSubRun(), art::OutputModule::doWriteEvent(), art::ProcessPackage< L >::postScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::Run >::End::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::End::postScheduleSignal(), art::ProcessPackage< L >::preScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::preScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::preScheduleSignal(), art::EventProcessor::readEvent(), art::EventProcessor::readRun(), art::EmptyEvent::readRun_(), art::EventProcessor::readSubRun(), and art::EmptyEvent::readSubRun_().

{
  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 272 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 439 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 199 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");
    fIsFirstSignalSpecial      = p.get< bool >("IsFirstSignalSpecial", false);

    // 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 658 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 295 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]),
          (idecay==0 && fIsFirstSignalSpecial) ? 0 : ( fT0[i] + flat.fire()*(fT1[i] - fT0[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 587 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 164 of file RadioGen_module.cc.

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

private

Definition at line 163 of file RadioGen_module.cc.

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

private

Definition at line 166 of file RadioGen_module.cc.

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

private

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

bool evgen::RadioGen::fIsFirstSignalSpecial

private

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

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

private

Definition at line 167 of file RadioGen_module.cc.

const double evgen::RadioGen::m_alpha = 3.727379240

private

Definition at line 159 of file RadioGen_module.cc.

const double evgen::RadioGen::m_e = 0.000510998928

private

Definition at line 158 of file RadioGen_module.cc.

const double evgen::RadioGen::m_neutron = 0.9395654133

private

Definition at line 160 of file RadioGen_module.cc.

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

private

Definition at line 170 of file RadioGen_module.cc.

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

private

Definition at line 169 of file RadioGen_module.cc.

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

private

Definition at line 162 of file RadioGen_module.cc.

int evgen::RadioGen::trackidcounter

private

Serial number for the MC track ID.

Definition at line 151 of file RadioGen_module.cc.

---

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

• /cvmfs/larsoft.opensciencegrid.org/products/larsim/v07_09_00/source/larsim/EventGenerator/RadioGen_module.cc