module to produce single or multiple specified particles in the detector More...

Inheritance diagram for evgen::GaisserParam:

Public Types
using	ModuleType = EDProducer

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

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

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

void	fillProductDescriptions ()

void	registerProducts (ProductDescriptions &productsToRegister)

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

std::unique_ptr< Worker >	makeWorker (WorkerParams const &wp)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Protected Member Functions
ConsumesCollector &	consumesCollector ()

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

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

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

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

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

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

Private Types
typedef std::map< double, TH1 * >	dhist_Map

typedef std::map< double, TH1 * >::iterator	dhist_Map_it

Private Member Functions
void	produce (art::Event &evt) override

void	beginJob () override

void	beginRun (art::Run &run) override

void	SampleOne (unsigned int i, simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)

void	Sample (simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)

std::pair< double, double >	GetThetaAndEnergy (double rand1, double rand2)

void	MakePDF ()

void	ResetMap ()

double	GaisserMuonFlux_Integrand (Double_t x, Double_t par)

double	GaisserFlux (double e, double theta)

std::vector< double >	GetBinning (const TAxis *axis, bool finalEdge=true)

Private Attributes
TFile *	m_File

dhist_Map *	m_PDFmap

TH1 *	m_thetaHist

CLHEP::HepRandomEngine &	fEngine
	art-managed random-number engine More...

int	fMode

bool	fPadOutVectors

std::vector< int >	fPDG
	PDG code of particles to generate. More...

int	fCharge
	Charge. More...

std::string	fInputDir
	Input Directory. More...

double	fEmin
	Minimum Kinetic Energy (GeV) More...

double	fEmax
	Maximum Kinetic Energy (GeV) More...

double	fEmid
	Energy to go from low to high (GeV) More...

int	fEBinsLow
	Number of low energy Bins. More...

int	fEBinsHigh
	Number of high energy Bins. More...

double	fThetamin
	Minimum theta. More...

double	fThetamax
	Maximum theta. More...

int	fThetaBins
	Number of theta Bins. More...

double	fXHalfRange
	Max X position. More...

double	fYInput
	Max Y position. More...

double	fZHalfRange
	Max Z position. More...

double	fT0
	Central t position (ns) in world coordinates. More...

double	fSigmaT
	Variation in t position (ns) More...

int	fTDist
	How to distribute t (gaus, or uniform) More...

bool	fSetParam
	Which version of Gaissers Param. More...

bool	fSetRead
	Whether to Read. More...

bool	fSetWrite
	Whether to Write. More...

bool	fSetReWrite
	Whether to ReWrite pdfs. More...

double	fEpsilon
	Minimum integration sum.... More...

double	fCryoBoundaries [6]

double	xNeg = 0

double	xPos = 0

double	zNeg = 0

double	zPos = 0

double	fCenterX = 0

double	fCenterZ = 0

TTree *	fTree

double	Time

double	Momentum

double	KinEnergy

double	Gamma

double	Energy

double	Theta

double	Phi

double	pnorm

double	XPosition

double	YPosition

double	ZPosition

double	DirCosineX

double	DirCosineY

double	DirCosineZ

Static Private Attributes
static constexpr int	kGAUS = 1

Detailed Description

module to produce single or multiple specified particles in the detector

Definition at line 61 of file GaisserParam_module.cc.

Member Typedef Documentation

typedef std::map<double, TH1*> evgen::GaisserParam::dhist_Map

private

Definition at line 73 of file GaisserParam_module.cc.

typedef std::map<double, TH1*>::iterator evgen::GaisserParam::dhist_Map_it

private

Definition at line 74 of file GaisserParam_module.cc.

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 26 of file Producer.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 147 of file GaisserParam_module.cc.

References art::detail::EngineCreator::createEngine(), fCharge, fEBinsHigh, fEBinsLow, fEmax, fEmid, fEmin, fEngine, fEpsilon, fInputDir, fMode, fPadOutVectors, fPDG, fSetParam, fSetRead, fSetReWrite, fSetWrite, fSigmaT, fT0, fTDist, fThetaBins, fThetamax, fThetamin, fXHalfRange, fYInput, and fZHalfRange.

     : art::EDProducer{pset}
     // create a default random engine; obtain the random seed from NuRandomService,
     // unless overridden in configuration with key "Seed"
     , fEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(createEngine(0),
                                                                                  pset,
                                                                                  "Seed"))
     , fMode{pset.get<int>("ParticleSelectionMode")}
     , fPadOutVectors{pset.get<bool>("PadOutVectors")}
     , fPDG{pset.get<std::vector<int>>("PDG")}
     , fCharge{pset.get<int>("Charge")}
     , fInputDir{pset.get<std::string>("InputDir")}
     , fEmin{pset.get<double>("Emin")}
     , fEmax{pset.get<double>("Emax")}
     , fEmid{pset.get<double>("Emid")}
     , fEBinsLow{pset.get<int>("EBinsLow")}
     , fEBinsHigh{pset.get<int>("EBinsHigh")}
     , fThetamin{pset.get<double>("Thetamin")}
     , fThetamax{pset.get<double>("Thetamax")}
     , fThetaBins{pset.get<int>("ThetaBins")}
     , fXHalfRange{pset.get<double>("XHalfRange")}
     , fYInput{pset.get<double>("YInput")}
     , fZHalfRange{pset.get<double>("ZHalfRange")}
     , fT0{pset.get<double>("T0")}
     , fSigmaT{pset.get<double>("SigmaT")}
     , fTDist{pset.get<int>("TDist")}
     , fSetParam{pset.get<bool>("SetParam")}
     , fSetRead{pset.get<bool>("SetRead")}
     , fSetWrite{pset.get<bool>("SetWrite")}
     , fSetReWrite{pset.get<bool>("SetReWrite")}
     , fEpsilon{pset.get<double>("Epsilon")}
   {
     produces<std::vector<simb::MCTruth>>();
     produces<sumdata::RunData, art::InRun>();
   }

Member Function Documentation

void evgen::GaisserParam::beginJob ( )

overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 184 of file GaisserParam_module.cc.

References DirCosineX, DirCosineY, DirCosineZ, Energy, fCenterX, fCenterZ, fCryoBoundaries, fTree, geo::GeometryCore::Iterate(), KinEnergy, Momentum, Phi, Theta, Time, xNeg, xPos, XPosition, YPosition, zNeg, zPos, and ZPosition.

   {
     //Work out center of cryostat(s)
     art::ServiceHandle<geo::Geometry const> geom;
     for (auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
       cryostat.Boundaries(fCryoBoundaries);
       if (xNeg > fCryoBoundaries[0]) xNeg = fCryoBoundaries[0];
       if (xPos < fCryoBoundaries[1]) xPos = fCryoBoundaries[1];
       if (zNeg > fCryoBoundaries[4]) zNeg = fCryoBoundaries[4];
       if (zPos < fCryoBoundaries[5]) zPos = fCryoBoundaries[5];
     }
     fCenterX = xNeg + (xPos - xNeg) / 2;
     fCenterZ = zNeg + (zPos - zNeg) / 2;
 
     // Make the Histograms....
     art::ServiceHandle<art::TFileService const> tfs;
     fTree = tfs->make<TTree>("Generator", "analysis tree");
     fTree->Branch("XPosition", &XPosition, "XPosition/D");
     fTree->Branch("YPosition", &YPosition, "YPosition/D");
     fTree->Branch("ZPosition", &ZPosition, "ZPosition/D");
     fTree->Branch("Time", &Time, "Time/D");
     fTree->Branch("Momentum", &Momentum, "Momentum/D");
     fTree->Branch("KinEnergy", &KinEnergy, "KinEnergy/D");
     fTree->Branch("Energy", &Energy, "Energy/D");
     fTree->Branch("DirCosineX", &DirCosineX, "DirCosineX/D");
     fTree->Branch("DirCosineY", &DirCosineY, "DirCosineY/D");
     fTree->Branch("DirCosineZ", &DirCosineZ, "DirCosineZ/D");
     fTree->Branch("Theta", &Theta, "Theta/D");
     fTree->Branch("Phi", &Phi, "Phi/D");
   }

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

overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 216 of file GaisserParam_module.cc.

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

   {
     // Check fcl parameters were set correctly
     if (fThetamax > M_PI / 2 + 0.01)
       throw cet::exception("GaisserParam")
         << "\nThetamax has to be less than " << M_PI / 2 << ", but was entered as " << fThetamax
         << ", this cause an error so leaving program now...\n\n";
     if (fThetamin < 0)
       throw cet::exception("GaisserParam")
         << "\nThetamin has to be more than 0, but was entered as " << fThetamin
         << ", this cause an error so leaving program now...\n\n"
         << std::endl;
     if (fThetamax < fThetamin)
       throw cet::exception("GaisserParam")
         << "\nMinimum angle is bigger than maximum angle....causes an error so leaving program "
            "now....\n\n"
         << std::endl;
 
     // grab the geometry object to see what geometry we are using
     art::ServiceHandle<geo::Geometry const> geo;
     run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
     MakePDF();
   }

template<typename T , BranchType BT>

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

protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

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

ConsumesCollector & art::ModuleBase::consumesCollector ( )

protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

   {
     return collector_;
   }

template<typename T , BranchType BT>

void art::ModuleBase::consumesMany ( )

protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

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

template<typename Element , BranchType = InEvent>

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

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::consumesView ( InputTag const & tag )

inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

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

void art::detail::Producer::doBeginJob ( SharedResources const & resources )

inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

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

bool art::detail::Producer::doBeginRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

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

bool art::detail::Producer::doBeginSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

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

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

inherited

Definition at line 30 of file Producer.cc.

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

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

bool art::detail::Producer::doEndRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

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

bool art::detail::Producer::doEndSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

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

bool art::detail::Producer::doEvent	(	EventPrincipal &	ep,
		ModuleContext const &	mc,
		std::atomic< std::size_t > &	counts_run,
		std::atomic< std::size_t > &	counts_passed,
		std::atomic< std::size_t > &	counts_failed
	)

inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

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

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

inherited

Definition at line 44 of file Producer.cc.

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

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

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

inherited

Definition at line 58 of file Producer.cc.

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

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

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

inherited

Definition at line 37 of file Producer.cc.

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

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

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

inherited

Definition at line 51 of file Producer.cc.

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

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

void art::Modifier::fillProductDescriptions ( )

inherited

Definition at line 10 of file Modifier.cc.

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

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

double evgen::GaisserParam::GaisserFlux	(	double	e,
		double	theta
	)

private

Definition at line 652 of file GaisserParam_module.cc.

References c1, and fSetParam.

Referenced by GaisserMuonFlux_Integrand().

   {
 
     double ct = cos(theta);
     double di;
     if (fSetParam) {
       double gamma = 2.7;
       double A = 0.14;
       double rc = 1.e-4; // fraction of prompt muons
       double c1 = sqrt(1. - (1. - pow(ct, 2.0)) / pow((1. + 32. / 6370.), 2.0)); // Earth curvature
       double deltae = 2.06e-3 * (1030. / c1 - 120.); // muon energy loss in atmosphere
       double em = e + deltae / 2.;
       double e1 = e + deltae;
       double pdec = pow((120. / (1030. / c1)), (1.04 / c1 / em)); // muon decay
       di = A * pow(e1, -gamma) *
            (1. / (1. + 1.1 * e1 * c1 / 115.) + 0.054 / (1. + 1.1 * e1 * c1 / 850.) + rc) * pdec;
     }
     else {
       double gamma = 2.7; // spectral index
       double A = 0.14;    // normalisation
       double C = 3.64;
       double gamma2 = 1.29;
       double ct_star = sqrt((pow(ct, 2) + 0.102573 * 0.102573 - 0.068287 * pow(ct, 0.958633) +
                              0.0407253 * pow(ct, 0.817285)) /
                             (1.0 + 0.102573 * 0.102573 - 0.068287 + 0.0407253));
       double eMod = e * (1.0 + C / (e * pow(ct_star, gamma2)));
       di = A * pow(eMod, -gamma) *
            (1. / (1. + 1.1 * e * ct_star / 115.) + 0.054 / (1. + 1.1 * e * ct_star / 850.));
     }
 
     return di;
   } // GaisserFlux

double evgen::GaisserParam::GaisserMuonFlux_Integrand	(	Double_t *	x,
		Double_t *	par
	)

private

Definition at line 686 of file GaisserParam_module.cc.

References GaisserFlux().

Referenced by MakePDF().

   {
 
     //---- calculate the flux
     double flux = 2.0 * M_PI * sin(x[1]) * GaisserFlux(x[0], x[1]);
 
     return flux;
   } // MuonFluxIntegrand

std::vector< double > evgen::GaisserParam::GetBinning	(	const TAxis *	axis,
		bool	finalEdge = `true`
	)

private

Definition at line 696 of file GaisserParam_module.cc.

Referenced by MakePDF().

   {
     std::vector<double> vbins;
     for (int ibin = 1; ibin <= axis->GetNbins() + 1; ibin++) {
       if (ibin <= axis->GetNbins())
         vbins.push_back(axis->GetBinLowEdge(ibin));
       else if (finalEdge)
         vbins.push_back(axis->GetBinLowEdge(ibin) + axis->GetBinWidth(ibin));
     }
     return vbins;
   } // Get Binning

std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const

inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

   {
     return collector_.getConsumables();
   }

std::pair< double, double > evgen::GaisserParam::GetThetaAndEnergy	(	double	rand1,
		double	rand2
	)

private

Definition at line 384 of file GaisserParam_module.cc.

References energy, m_PDFmap, and m_thetaHist.

Referenced by SampleOne().

   {
     if (rand1 < 0 || rand1 > 1)
       std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand1 << std::endl;
     if (rand2 < 0 || rand2 > 1)
       std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand2 << std::endl;
 
     int thetaBin = 0;
     m_thetaHist->GetBinWithContent(double(rand1), thetaBin, 1, m_thetaHist->GetNbinsX(), 1.0);
     if (m_thetaHist->GetBinContent(thetaBin) < rand1) thetaBin += 1;
     double drand1 =
       (m_thetaHist->GetBinContent(thetaBin) - rand1) /
       (m_thetaHist->GetBinContent(thetaBin) - m_thetaHist->GetBinContent(thetaBin - 1));
     double thetaLow = m_thetaHist->GetXaxis()->GetBinLowEdge(thetaBin);
     double thetaUp = m_thetaHist->GetXaxis()->GetBinUpEdge(thetaBin);
     double theta = drand1 * (thetaUp - thetaLow) + thetaLow;
 
     int energyBin = 0;
     TH1* energyHist = 0;
     bool notFound = true;
     for (dhist_Map_it mapit = m_PDFmap->begin(); mapit != m_PDFmap->end(); mapit++) {
       if (fabs(mapit->first + thetaLow) < 0.000001) {
         energyHist = mapit->second;
         notFound = false;
         break;
       }
     }
     if (notFound) std::cout << "GetThetaAndEnergy: ERROR:\tInvalid theta!" << std::endl;
     //  MSG("string = " << To_TString(thetaLow) << ", double = " << thetaLow );
 
     energyHist->GetBinWithContent(double(rand2), energyBin, 1, energyHist->GetNbinsX(), 1.0);
     if (energyHist->GetBinContent(energyBin) < rand2) energyBin += 1;
     double drand2 =
       (energyHist->GetBinContent(energyBin) - rand2) /
       (energyHist->GetBinContent(energyBin) - energyHist->GetBinContent(energyBin - 1));
     double energyLow = energyHist->GetXaxis()->GetBinLowEdge(energyBin);
     double energyUp = energyHist->GetXaxis()->GetBinUpEdge(energyBin);
     double energy = drand2 * (energyUp - energyLow) + energyLow;
     //  MSG("MuFlux::GetThetaEnergy()\te = " << energy*1000. );
 
     return std::make_pair(theta, energy);
   }

void evgen::GaisserParam::MakePDF ( )

private

Definition at line 428 of file GaisserParam_module.cc.

References fEBinsHigh, fEBinsLow, fEmax, fEmid, fEmin, fEpsilon, fInputDir, fSetRead, fSetReWrite, fSetWrite, fThetaBins, fThetamax, fThetamin, GaisserMuonFlux_Integrand(), GetBinning(), m_File, m_PDFmap, m_thetaHist, and ResetMap().

Referenced by beginRun().

   {
     std::cout << "In my function MakePDF" << std::endl;
     m_File = 0;
     m_PDFmap = 0;
     m_thetaHist = 0;
     double TotalMuonFlux = 0;
 
     if (m_PDFmap) {
       std::cout << "PDFMAP" << std::endl;
       if (m_PDFmap->size() > 0 && !fSetReWrite) {
         std::cout << "MakePDF: Map has already been initialised. " << std::endl;
         std::cout << "Do fSetReWrite - true if you really want to override this map." << std::endl;
         return;
       }
       std::cout << fSetReWrite << std::endl;
       if (fSetReWrite) ResetMap();
     }
     else {
       m_PDFmap = new dhist_Map;
       std::cout << "Making new dhist_Map called m_PDFmap....." << std::endl;
     }
 
     TF2* muonSpec = new TF2("muonSpec",
                             this,
                             &evgen::GaisserParam::GaisserMuonFlux_Integrand,
                             fEmin,
                             fEmax,
                             fThetamin,
                             fThetamax,
                             0,
                             "GaisserParam",
                             "GaisserMuonFlux_Integrand");
     //--------------------------------------------
     //------------ Compute the pdfs
 
     //---- compute pdf for the theta
     TotalMuonFlux = muonSpec->Integral(
       fEmin, fEmax, fThetamin, fThetamax, fEpsilon); // Work out the muon flux at the surface
     std::cout << "Surface flux of muons = " << TotalMuonFlux << " cm-2 s-1" << std::endl;
 
     //---- work out if we're reading a file, writing to file, or neither
     std::ostringstream pdfFile;
     pdfFile << "GaisserPDF_" << fEmin << "-" << fEmid << "-" << fEmax << "-" << fEBinsLow << "-"
             << fEBinsHigh << "-" << fThetamin << "-" << fThetamax << "-" << fThetaBins << ".root";
     std::string tmpfileName = pdfFile.str();
     std::replace(tmpfileName.begin(), tmpfileName.end(), '+', '0');
     if (tmpfileName == "GaisserPDF_0-100-100100-1000-10000-0-1.5708-100.root")
       tmpfileName = "GaisserPDF_DefaultBins.root";
     else if (tmpfileName == "GaisserPDF_0-100-4000-1000-1000-0-1.5708-100.root")
       tmpfileName = "GaisserPDF_LowEnergy.root";
     else if (tmpfileName == "GaisserPDF_4000-10000-100000-1000-10000-0-1.5708-100.root")
       tmpfileName = "GaisserPDF_MidEnergy.root";
     else if (tmpfileName == "GaisserPDF_100000-500000-1e007-10000-100000-0-1.5708-100.root")
       tmpfileName = "GaisserPDF_HighEnergy.root";
 
     std::ostringstream pdfFilePath;
     pdfFilePath << fInputDir << tmpfileName;
     std::string fileName = pdfFilePath.str();
 
     cet::search_path sp("FW_SEARCH_PATH");
     std::string fROOTfile; //return /lbne/data/0-100-1.57.root
     if (sp.find_file(tmpfileName, fROOTfile)) fileName = fROOTfile;
 
     std::cout << "File path; " << fileName << std::endl;
 
     if (fSetRead) {
       struct stat buffer;
       fSetRead = stat(fileName.c_str(), &buffer) == 0; // Check if file exists already
       if (!fSetRead)
         std::cout << "WARNING- " + fileName + " does not exist." << std::endl;
       else {
         std::cout << "Reading PDF from file " + fileName << std::endl;
         m_File = new TFile(fileName.c_str()); // Open file
         if (m_File->IsZombie() ||
             m_File->TestBit(TFile::kRecovered)) { // Check that file is not corrupted
           std::cout << "WARNING- " + fileName + " is corrupted or cannot be read." << std::endl;
           fSetRead = false;
         }
       }
     }
 
     if (fSetRead) { // If the file exists then read it....
       std::cout << "Now going to read file...." << std::endl;
       fSetWrite = false; // Don't want to write as already exists
       double thetalow = fThetamin;
       m_thetaHist = (TH1D*)m_File->Get("pdf_theta");
       for (int i = 1; i <= fThetaBins; i++) {
         std::ostringstream pdfEnergyHist;
         pdfEnergyHist << "pdf_energy_" << i;
         std::string pdfEnergyHiststr = pdfEnergyHist.str();
 
         TH1* pdf_hist = (TH1D*)m_File->Get(pdfEnergyHiststr.c_str()); // Get the bin
         m_PDFmap->insert(std::make_pair(-thetalow, pdf_hist)); //---- -ve theta for quicker sorting
         thetalow = double(i) * (fThetamax) / double(fThetaBins); //---- increment the value of theta
       }
     }      // ------------------------------------------
     else { // File doesn't exist so want to make it.....
            // ------------------------------------------
       std::cout << "Generating a new muon flux PDF" << std::endl;
       if (fSetWrite) {
         std::cout << "Writing to PDF to file " + fileName << std::endl;
         m_File = new TFile(fileName.c_str(), "RECREATE");
       }
 
       double dnbins_theta = double(fThetaBins);
       m_thetaHist = new TH1D("pdf_theta", "pdf_theta", fThetaBins, fThetamin, fThetamax);
       for (int i = 1; i <= fThetaBins; i++) {
         double di = i == 0 ? 0.1 : double(i);
         double theta = di * (fThetamax) / dnbins_theta;
         double int_i = muonSpec->Integral(fEmin, fEmax, fThetamin, theta, fEpsilon);
         m_thetaHist->SetBinContent(i, int_i / TotalMuonFlux);
       }
       if (fSetWrite) m_thetaHist->Write();
       std::cout
         << "theta PDF complete... now making the energy PDFs (this will take a little longer)... "
         << std::endl;
 
       //---- now compute the energy pdf
       double thetalow = fThetamin;
       for (int i = 1; i <= fThetaBins; i++) {
 
         double di = double(i);
         double theta = di * (fThetamax) / fThetaBins;
 
         //---- compute the total integral
         double int_tot = muonSpec->Integral(fEmin, fEmax, thetalow, theta, fEpsilon);
 
         //---- compute pdf for the low energy
         int nbins = fEBinsLow;
         TH1* pdf_lowenergy = new TH1D("pdf_lowenergy", "pdf_lowenergy", nbins, fEmin, fEmid);
         double dnbins = double(nbins);
         for (int j = 1; j <= nbins; j++) {
           double dj = double(j);
           double int_j = muonSpec->Integral(
             fEmin, fEmin + dj * (fEmid - fEmin) / dnbins, thetalow, theta, fEpsilon);
           pdf_lowenergy->SetBinContent(j, int_j / int_tot);
         }
 
         //---- compute pdf for the high energy
         nbins = fEBinsHigh;
         dnbins = double(nbins);
         TH1* pdf_highenergy = new TH1D("pdf_highenergy", "pdf_highenergy", nbins, fEmid, fEmax);
         for (int j = 1; j <= nbins; j++) {
           double dj = double(j);
           double int_j = muonSpec->Integral(
             fEmin, fEmid + dj * (fEmax - fEmid) / dnbins, thetalow, theta, fEpsilon);
           pdf_highenergy->SetBinContent(j, int_j / int_tot);
         }
 
         //---- now combine the two energy hists
         std::vector<double> vxbins = GetBinning(pdf_lowenergy->GetXaxis(), false);
         std::vector<double> vxbins2 = GetBinning(pdf_highenergy->GetXaxis());
         vxbins.insert(vxbins.end(), vxbins2.begin(), vxbins2.end());
 
         int ibin = 0;
         double* xbins = new double[vxbins.size()];
         for (std::vector<double>::const_iterator binit = vxbins.begin(); binit != vxbins.end();
              binit++, ibin++)
           xbins[ibin] = (*binit);
         TH1* pdf_energy = new TH1D("pdf_energy", "pdf_energy", vxbins.size() - 1, xbins);
         int ibin2 = 1;
         for (ibin = 1; ibin <= pdf_lowenergy->GetNbinsX(); ibin++, ibin2++) {
           double content = pdf_lowenergy->GetBinContent(ibin);
           if (ibin == 1) content = content - 0.00001;
           pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), content);
         }
         for (ibin = 1; ibin <= pdf_highenergy->GetNbinsX(); ibin++, ibin2++) {
           pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), pdf_highenergy->GetBinContent(ibin));
         }
 
         //---- and remove any negative bins
         std::ostringstream Clonestr;
         Clonestr << "pdf_energy_" << i;
         TH1* pdf_energy_noneg = (TH1D*)pdf_energy->Clone(Clonestr.str().c_str());
         pdf_energy_noneg->Reset();
 
         double PDF = 0.0;
         double lastPD = 0.0;
         int nSkip = 0;
         for (ibin = 1; ibin <= pdf_energy->GetNbinsX(); ibin++) {
           double newPD = pdf_energy->GetBinContent(ibin);
           double probDiff = newPD - lastPD;
           if (probDiff < 0) {
             if (ibin != pdf_energy->GetNbinsX()) {
               nSkip++;
               continue;
             }
             else
               probDiff = 0;
           }
           else
             PDF += probDiff;
           if (PDF > 1) PDF = 1;
           for (int iskip = 0; iskip <= nSkip; iskip++)
             pdf_energy_noneg->Fill(pdf_energy_noneg->GetBinCenter(ibin - iskip), PDF);
           nSkip = 0;
           lastPD = newPD;
         }
 
         //---- add this hist, increment thetalow and delete unwanted TH1s
         if (fSetWrite) pdf_energy_noneg->Write();
         m_PDFmap->insert(
           std::make_pair(-thetalow, pdf_energy_noneg)); //---- -ve theta for quicker sorting
 
         //---- increment the value of theta
         thetalow = theta;
 
         //---- free up memory from unwanted hists
         delete pdf_lowenergy;
         delete pdf_highenergy;
         delete pdf_energy;
 
         std::cout << "\r===> " << 100.0 * double(i) / double(fThetaBins) << "% complete... "
                   << std::endl;
       } // ThetaBins
       std::cout << "finished the energy pdfs." << std::endl;
     } //---- if(!m_doRead)
 
     delete muonSpec;
     return;
   } // Make PDF

std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const & wp )

inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

   {
     return doMakeWorker(wp);
   }

template<typename T , BranchType BT>

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

protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

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

template<typename T , BranchType BT>

void art::ModuleBase::mayConsumeMany ( )

protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

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

template<typename Element , BranchType = InEvent>

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

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const & tag )

inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

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

ModuleDescription const & art::ModuleBase::moduleDescription ( ) const

inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

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

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

overrideprivatevirtual

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

Implements art::EDProducer.

Definition at line 241 of file GaisserParam_module.cc.

References fEngine, simb::kSingleParticle, MF_LOG_DEBUG, art::Event::put(), Sample(), and simb::MCTruth::SetOrigin().

   {
     std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
 
     simb::MCTruth truth;
     truth.SetOrigin(simb::kSingleParticle);
 
     Sample(truth, fEngine);
 
     MF_LOG_DEBUG("GaisserParam") << truth;
 
     truthcol->push_back(truth);
 
     evt.put(std::move(truthcol));
   }

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

inherited

Definition at line 16 of file Modifier.cc.

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

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

void evgen::GaisserParam::ResetMap ( )

private

Definition at line 709 of file GaisserParam_module.cc.

References DEFINE_ART_MODULE, m_PDFmap, and m_thetaHist.

Referenced by MakePDF().

   {
     if (m_thetaHist) delete m_thetaHist;
     if (m_PDFmap) {
       for (dhist_Map_it mapit = m_PDFmap->begin(); mapit != m_PDFmap->end(); mapit++) {
         if (mapit->second) delete mapit->second;
       }
       m_PDFmap->clear();
       delete m_PDFmap;
       std::cout << "Reset PDFmap and thetaHist..." << std::endl;
     }
   } // ResetMap

void evgen::GaisserParam::Sample	(	simb::MCTruth &	mct,
		CLHEP::HepRandomEngine &	engine
	)

private

Definition at line 357 of file GaisserParam_module.cc.

References fMode, fPDG, and SampleOne().

Referenced by produce().

   {
     switch (fMode) {
     case 0: // List generation mode: every event will have one of each
             // particle species in the fPDG array
       for (unsigned int i = 0; i < fPDG.size(); ++i) {
         SampleOne(i, mct, engine);
       } //end loop over particles
       break;
     case 1: // Random selection mode: every event will exactly one particle
             // selected randomly from the fPDG array
     {
       CLHEP::RandFlat flat(engine);
 
       unsigned int i = flat.fireInt(fPDG.size());
       SampleOne(i, mct, engine);
     } break;
     default:
       mf::LogWarning("UnrecognizeOption")
         << "GaisserParam does not recognize ParticleSelectionMode " << fMode;
       break;
     } // switch on fMode
 
     return;
   }

void evgen::GaisserParam::SampleOne	(	unsigned int	i,
		simb::MCTruth &	mct,
		CLHEP::HepRandomEngine &	engine
	)

private

Definition at line 261 of file GaisserParam_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), DirCosineX, DirCosineY, DirCosineZ, Energy, fCenterX, fCenterZ, fCharge, fPDG, fSigmaT, fT0, fTDist, fTree, fXHalfRange, fYInput, fZHalfRange, Gamma, GetThetaAndEnergy(), kGAUS, KinEnergy, m_PDFmap, Momentum, part, Phi, pnorm, Theta, Time, x, XPosition, YPosition, and ZPosition.

Referenced by Sample().

   {
 
     CLHEP::RandFlat flat(engine);
     CLHEP::RandGaussQ gauss(engine);
 
     Time = Momentum = KinEnergy = Gamma = Energy = Theta = Phi = pnorm = 0.0;
     XPosition = YPosition = ZPosition = DirCosineX = DirCosineY = DirCosineZ = 0.0;
 
     TVector3 x;
     TVector3 pmom;
 
     // set track id to -i as these are all primary particles and have id <= 0
     int trackid = -1 * (i + 1);
     std::string primary("primary");
 
     // Work out whether particle/antiparticle, and mass...
     double m = 0.0;
     fPDG[i] = 13;
     if (fCharge == 0) {
       if (1.0 - 2.0 * flat.fire() > 0) fPDG[i] = -fPDG[i];
     }
     else if (fCharge == 1)
       fPDG[i] = 13;
     else if (fCharge == 2)
       fPDG[i] = -13;
     static TDatabasePDG pdgt;
     TParticlePDG* pdgp = pdgt.GetParticle(fPDG[i]);
     if (pdgp) m = pdgp->Mass();
 
     // Work out T0...
     if (fTDist == kGAUS) { Time = gauss.fire(fT0, fSigmaT); }
     else {
       Time = fT0 + fSigmaT * (2.0 * flat.fire() - 1.0);
     }
 
     // Work out Positioning.....
     x[0] = (1.0 - 2.0 * flat.fire()) * fXHalfRange + fCenterX;
     x[1] = fYInput;
     x[2] = (1.0 - 2.0 * flat.fire()) * fZHalfRange + fCenterZ;
 
     // Make Lorentz vector for x and t....
     TLorentzVector pos(x[0], x[1], x[2], Time);
 
     // Access the pdf map which has been loaded.....
     if (m_PDFmap) {
 
       //---- get the muon theta and energy from histograms using 2 random numbers
       std::pair<double, double> theta_energy; //---- muon theta and energy
       theta_energy = GetThetaAndEnergy(flat.fire(), flat.fire());
 
       //---- Set theta, phi
       Theta = theta_energy.first; // Angle returned by GetThetaAndEnergy between 0 and pi/2
       Phi = M_PI * (1.0 - 2.0 * flat.fire()); // Randomly generated angle between -pi and pi
 
       //---- Set KE, E, p
       KinEnergy = theta_energy.second; // Energy returned by GetThetaAndEnergy
       Gamma = 1 + (KinEnergy / m);
       Energy = Gamma * m;
       Momentum = std::sqrt(Energy * Energy - m * m); // Get momentum
 
       pmom[0] = Momentum * std::sin(Theta) * std::cos(Phi);
       pmom[1] = -Momentum * std::cos(Theta);
       pmom[2] = Momentum * std::sin(Theta) * std::sin(Phi);
 
       pnorm = std::sqrt(pmom[0] * pmom[0] + pmom[1] * pmom[1] + pmom[2] * pmom[2]);
       DirCosineX = pmom[0] / pnorm;
       DirCosineY = pmom[1] / pnorm;
       DirCosineZ = pmom[2] / pnorm;
     }
     else {
       std::cout << "MuFlux map hasn't been initialised, aborting...." << std::endl;
       return;
     }
 
     TLorentzVector pvec(pmom[0], pmom[1], pmom[2], Energy);
 
     simb::MCParticle part(trackid, fPDG[i], primary);
     part.AddTrajectoryPoint(pos, pvec);
     XPosition = x[0];
     YPosition = x[1];
     ZPosition = x[2];
 
     std::cout << "Theta: " << Theta << " Phi: " << Phi << " KineticEnergy: " << KinEnergy
               << " Energy: " << Energy << " Momentum: " << Momentum << std::endl;
     std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T()
               << std::endl;
     std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
     std::cout << "Normalised..." << DirCosineX << " " << DirCosineY << " " << DirCosineZ
               << std::endl;
 
     fTree->Fill();
     mct.Add(part);
   }

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

inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

   {
     md_ = md;
   }

void art::ModuleBase::sortConsumables ( std::string const & current_process_name )

inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

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

Member Data Documentation

double evgen::GaisserParam::DirCosineX

private

Definition at line 140 of file GaisserParam_module.cc.

Referenced by beginJob(), and SampleOne().

double evgen::GaisserParam::DirCosineY

private

Definition at line 140 of file GaisserParam_module.cc.

Referenced by beginJob(), and SampleOne().

double evgen::GaisserParam::DirCosineZ

private

Definition at line 140 of file GaisserParam_module.cc.

Referenced by beginJob(), and SampleOne().

double evgen::GaisserParam::Energy

private

Definition at line 139 of file GaisserParam_module.cc.

Referenced by beginJob(), and SampleOne().

double evgen::GaisserParam::fCenterX = 0

private

Definition at line 134 of file GaisserParam_module.cc.

Referenced by beginJob(), and SampleOne().

double evgen::GaisserParam::fCenterZ = 0

private

Definition at line 135 of file GaisserParam_module.cc.

Referenced by beginJob(), and SampleOne().

int evgen::GaisserParam::fCharge

private

Charge.

Definition at line 101 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and SampleOne().

double evgen::GaisserParam::fCryoBoundaries[6]

private

Definition at line 129 of file GaisserParam_module.cc.

Referenced by beginJob().

int evgen::GaisserParam::fEBinsHigh

private

Number of high energy Bins.

Definition at line 108 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

int evgen::GaisserParam::fEBinsLow

private

Number of low energy Bins.

Definition at line 107 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

double evgen::GaisserParam::fEmax

private

Maximum Kinetic Energy (GeV)

Definition at line 105 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

double evgen::GaisserParam::fEmid

private

Energy to go from low to high (GeV)

Definition at line 106 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

double evgen::GaisserParam::fEmin

private

Minimum Kinetic Energy (GeV)

Definition at line 104 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

CLHEP::HepRandomEngine& evgen::GaisserParam::fEngine

private

art-managed random-number engine

Definition at line 89 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and produce().

double evgen::GaisserParam::fEpsilon

private

Minimum integration sum....

Definition at line 126 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

std::string evgen::GaisserParam::fInputDir

private

Input Directory.

Definition at line 102 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and MakePDF().

int evgen::GaisserParam::fMode

private

Particle Selection Mode 0–generate a list of all particles, 1–generate a single particle selected randomly from the list

Definition at line 93 of file GaisserParam_module.cc.

Referenced by GaisserParam(), and Sample().

bool evgen::GaisserParam::fPadOutVectors

private

Select to pad out configuration vectors if they are not of of the same length as PDG false: don't pad out - all values need to specified true: pad out - default values assumed and printed out

Definition at line 96 of file GaisserParam_module.cc.

Referenced by GaisserParam().

std::vector<int> evgen::GaisserParam::fPDG

private

PDG code of particles to generate.

Definition at line 100 of file GaisserParam_module.cc.

Referenced by GaisserParam(), Sample(), and SampleOne().

bool evgen::GaisserParam::fSetParam

private

Which version of Gaissers Param.

Definition at line 122 of file GaisserParam_module.cc.

Referenced by GaisserFlux(), and GaisserParam().