gar::evgen::RadioGen Class Reference

Inheritance diagram for gar::evgen::RadioGen:

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)
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 Member Functions
void	SampleOne (unsigned int i, simb::MCTruth &mct)
void	readfile (std::string nuclide, std::string filename)
void	samplespectrum (std::string nuclide, int &itype, double &t, double &m, double &p)
double	samplefromth1d (TH1D *hist)

Private Attributes
std::vector< std::string >	fNuclide
	List of nuclides to simulate. Example: "39Ar". 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
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
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, patterend off of SingleGen Currently it generates only in rectangular prisms oriented along the x,y,z axes

Definition at line 74 of file RadioGen_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 138 of file RadioGen_module.cc.

                                                    : EDProducer{pset},
      fEngine(art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this,pset,"Seed"))
    {
      
      this->reconfigure(pset);
      
      produces< std::vector<simb::MCTruth> >();
      produces< sumdata::RunData, ::art::InRun >();
      
    }

gar::evgen::RadioGen::~RadioGen ( )

virtual

Definition at line 150 of file RadioGen_module.cc.

    {
    }

Member Function Documentation

void gar::evgen::RadioGen::beginRun

(

::art::Run &

run

)

Definition at line 194 of file RadioGen_module.cc.

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

void gar::evgen::RadioGen::produce

(

::art::Event &

evt

)

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

Definition at line 207 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));
      return;
    }

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

private

Definition at line 309 of file RadioGen_module.cc.

    {
      int ifound = 0;
      for (size_t i=0; i<fNuclide.size(); i++)
      {
        if (fNuclide[i] == nuclide)
        {
          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");
      
      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-1);
        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-1);
        
        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-1);
        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);
      }
      
      f.Close();
      TH1::AddDirectory(addStatus);
      
      double total = alphaintegral.back() + betaintegral.back() + gammaintegral.back();
      if (total>0)
      {
        alphaintegral.back() /= total;
        betaintegral.back() /= total;
        gammaintegral.back() /= total;
      }
    }

void gar::evgen::RadioGen::reconfigure

(

fhicl::ParameterSet const &

p

)

Definition at line 155 of file RadioGen_module.cc.

    {
        // 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");
      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 (  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";
      
      readfile("39Ar","Argon_39.root");
      readfile("60Co","Cobalt_60.root");
      readfile("85Kr","Krypton_85.root");
      readfile("40K","Potassium_40.root");
      readfile("232Th","Thorium_232.root");
      readfile("238U","Uranium_238.root");

      return;
    }

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

private

Definition at line 493 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 gar::evgen::RadioGen::SampleOne ( unsigned int  i, simb::MCTruth &  mct  )

private

Definition at line 229 of file RadioGen_module.cc.

                                                            {
      
     // get the random number generator service and make some CLHEP generators
      CLHEP::RandFlat   flat(fEngine);
      //CLHEP::RandGauss  gauss(fEngine);
      CLHEP::RandPoisson  poisson(fEngine);
      
      // figure out how many decays to generate
      
      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
        
        int pdgid=0;  // electron=11, photon=22, alpha = 1000020040
        double t = 0; // kinetic energy of particle GeV
        double m = 0; // mass of daughter particle GeV
        double p = 0; // generated momentum (GeV)
        
          // Treat 222Rn separately
        if (fNuclide[i] == "222Rn")
        {
          pdgid = 1000020040;
          t     = 0.00548952;
          m     = m_alpha;
          double energy = t + m;
          double p2     = energy*energy - m*m;
          if (p2 > 0) p = TMath::Sqrt(p2);
          else        p = 0;
        }
        else samplespectrum(fNuclide[i],pdgid,t,m,p);
        
        
          // uniformly distributed in position and time
        
        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]) );
        
          // 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));
        
          // set track id to a negative serial number as these are all primary particles and have id <= 0
        int trackid = trackidcounter;
        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);
        }
        else
        {
          simb::MCParticle part(trackid, pdgid, primary);
          part.AddTrajectoryPoint(pos, pvec);
          mct.Add(part);
        }
      }
    }

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

private

Definition at line 430 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] != 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 ( gammaspectrum[inuc] != 0)
        {
          itype = 22; // gamma
          m = 0;
          t = samplefromth1d(gammaspectrum[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> gar::evgen::RadioGen::alphaintegral

private

Definition at line 124 of file RadioGen_module.cc.

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

private

Definition at line 123 of file RadioGen_module.cc.

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

private

Definition at line 126 of file RadioGen_module.cc.

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

private

Definition at line 125 of file RadioGen_module.cc.

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

private

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

Definition at line 104 of file RadioGen_module.cc.

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

private

Definition at line 130 of file RadioGen_module.cc.

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

private

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

Definition at line 103 of file RadioGen_module.cc.

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

private

Beginning of time window to simulate in ns.

Definition at line 105 of file RadioGen_module.cc.

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

private

End of time window to simulate in ns.

Definition at line 106 of file RadioGen_module.cc.

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

private

Bottom corner x position (cm) in world coordinates.

Definition at line 107 of file RadioGen_module.cc.

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

private

Top corner x position (cm) in world coordinates.

Definition at line 110 of file RadioGen_module.cc.

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

private

Bottom corner y position (cm) in world coordinates.

Definition at line 108 of file RadioGen_module.cc.

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

private

Top corner y position (cm) in world coordinates.

Definition at line 111 of file RadioGen_module.cc.

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

private

Bottom corner z position (cm) in world coordinates.

Definition at line 109 of file RadioGen_module.cc.

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

private

Top corner z position (cm) in world coordinates.

Definition at line 112 of file RadioGen_module.cc.

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

private

Definition at line 128 of file RadioGen_module.cc.

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

private

Definition at line 127 of file RadioGen_module.cc.

const double gar::evgen::RadioGen::m_alpha = 3.727379240

private

Definition at line 120 of file RadioGen_module.cc.

const double gar::evgen::RadioGen::m_e = 0.000510998928

private

Definition at line 119 of file RadioGen_module.cc.

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

private

Definition at line 122 of file RadioGen_module.cc.

int gar::evgen::RadioGen::trackidcounter

private

Serial number for the MC track ID.

Definition at line 113 of file RadioGen_module.cc.

---

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

• garsoft/EventGenerator/RadioGen_module.cc