LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
evgen::GaisserParam Class Reference

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

Inheritance diagram for evgen::GaisserParam:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

typedef std::map< double, TH1 * > dhist_Map
 
typedef std::map< double, TH1 * >::iterator dhist_Map_it
 
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 GaisserParam (fhicl::ParameterSet const &pset)
 
virtual ~GaisserParam ()
 
void produce (art::Event &evt)
 
void beginJob ()
 
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_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (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< Consumernon_module_context ()
 

Public Attributes

TFile * m_File
 
dhist_Mapm_PDFmap
 
TH1 * m_thetaHist
 

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 Member Functions

void SampleOne (unsigned int i, simb::MCTruth &mct)
 
void Sample (simb::MCTruth &mct)
 
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

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 const int kGAUS = 1
 

Detailed Description

module to produce single or multiple specified particles in the detector

Definition at line 74 of file GaisserParam_module.cc.

Member Typedef Documentation

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

Definition at line 87 of file GaisserParam_module.cc.

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

Definition at line 88 of file GaisserParam_module.cc.

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.

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

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

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

Definition at line 177 of file GaisserParam_module.cc.

178  {
179  this->reconfigure(pset);
180 
181  // create a default random engine; obtain the random seed from NuRandomService,
182  // unless overridden in configuration with key "Seed"
184  ->createEngine(*this, pset, "Seed");
185 
186  produces< std::vector<simb::MCTruth> >();
187  produces< sumdata::RunData, art::InRun >();
188  }
base_engine_t & createEngine(seed_t seed)
void reconfigure(fhicl::ParameterSet const &p)
evgen::GaisserParam::~GaisserParam ( )
virtual

Definition at line 191 of file GaisserParam_module.cc.

192  {
193  }

Member Function Documentation

void evgen::GaisserParam::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 236 of file GaisserParam_module.cc.

References geo::GeometryCore::CryostatBoundaries(), art::TFileDirectory::make(), and geo::GeometryCore::Ncryostats().

237  {
238  //Work out center of cryostat(s)
240  for (unsigned int i=0; i < geom->Ncryostats() ; i++ ) {
242  if ( xNeg > fCryoBoundaries[0] ) xNeg = fCryoBoundaries[0];
243  if ( xPos < fCryoBoundaries[1] ) xPos = fCryoBoundaries[1];
244  if ( zNeg > fCryoBoundaries[4] ) zNeg = fCryoBoundaries[4];
245  if ( zPos < fCryoBoundaries[5] ) zPos = fCryoBoundaries[5];
246  }
247  fCenterX = xNeg + (xPos-xNeg)/2;
248  fCenterZ = zNeg + (zPos-zNeg)/2;
249 
250  // Make the Histograms....
252  /*
253  fPositionX = tfs->make<TH1D>("fPositionX" ,"Position (cm)" ,500,fCenterX-(fXHalfRange+10) ,fCenterX+(fXHalfRange+10));
254  fPositionY = tfs->make<TH1D>("fPositionY" ,"Position (cm)" ,500,-(fYInput+10),(fYInput+10));
255  fPositionZ = tfs->make<TH1D>("fPositionZ" ,"Position (cm)" ,500,fCenterZ-(fZHalfRange+10) ,fCenterZ+(fZHalfRange+10));
256  fTime = tfs->make<TH1D>("fTime" ,"Time (s)" ,500,0,1e6);
257  fMomentumHigh = tfs->make<TH1D>("fMomentumHigh","Momentum (GeV)",500,0,fEmax);
258  fMomentum = tfs->make<TH1D>("fMomentum" ,"Momentum (GeV)",500,0,100);
259  fEnergy = tfs->make<TH1D>("fEnergy" ,"Energy (GeV)" ,500,0,fEmax);
260 
261  fDirCosineX = tfs->make<TH1D>("fDirCosineX","Normalised Direction cosine",500,-1,1);
262  fDirCosineY = tfs->make<TH1D>("fDirCosineY","Normalised Direction cosine",500,-1,1);
263  fDirCosineZ = tfs->make<TH1D>("fDirCosineZ","Normalised Direction cosine",500,-1,1);
264 
265  fTheta = tfs->make<TH1D>("fTheta" ,"Angle (radians)",500,-365,365);
266  fPhi = tfs->make<TH1D>("fPhi" ,"Angle (radians)",500,-365,365);
267  */
268  fTree = tfs->make<TTree>("Generator","analysis tree");
269  fTree->Branch("XPosition" ,&XPosition ,"XPosition/D" );
270  fTree->Branch("YPosition" ,&YPosition ,"YPosition/D" );
271  fTree->Branch("ZPosition" ,&ZPosition ,"ZPosition/D" );
272  fTree->Branch("Time" ,&Time ,"Time/D" );
273  fTree->Branch("Momentum" ,&Momentum ,"Momentum/D" );
274  fTree->Branch("KinEnergy" ,&KinEnergy ,"KinEnergy/D" );
275  fTree->Branch("Energy" ,&Energy ,"Energy/D" );
276  fTree->Branch("DirCosineX",&DirCosineX,"DirCosineX/D");
277  fTree->Branch("DirCosineY",&DirCosineY,"DirCosineY/D");
278  fTree->Branch("DirCosineZ",&DirCosineZ,"DirCosineZ/D");
279  fTree->Branch("Theta" ,&Theta ,"Theta/D" );
280  fTree->Branch("Phi" ,&Phi ,"Phi/D" );
281  }
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
void CryostatBoundaries(double *boundaries, geo::CryostatID const &cid) const
Returns the boundaries of the specified cryostat.
void evgen::GaisserParam::beginRun ( art::Run run)
virtual

Reimplemented from art::EDProducer.

Definition at line 284 of file GaisserParam_module.cc.

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

285  {
286  // grab the geometry object to see what geometry we are using
288  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetectorName()));
289 
290  // Check fcl parameters were set correctly
291  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";
292  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;
293  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;
294 
295  run.put(std::move(runcol));
296  MakePDF ();
297 
298  return;
299  }
double fThetamax
Maximum theta.
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:148
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
double fThetamin
Minimum theta.
Namespace collecting geometry-related classes utilities.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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().

148 {
149  if (!moduleContext_)
150  return ProductToken<T>::invalid();
151 
152  consumables_[BT].emplace_back(ConsumableType::Product,
153  TypeID{typeid(T)},
154  it.label(),
155  it.instance(),
156  it.process());
157  return ProductToken<T>{it};
158 }
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 162 of file Consumer.h.

163 {
164  if (!moduleContext_)
165  return;
166 
167  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
168 }
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
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().

173 {
174  if (!moduleContext_)
175  return ViewToken<T>::invalid();
176 
177  consumables_[BT].emplace_back(ConsumableType::ViewElement,
178  TypeID{typeid(T)},
179  it.label(),
180  it.instance(),
181  it.process());
182  return ViewToken<T>{it};
183 }
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
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().

34 {
35  return rng()->createEngine(
36  placeholder_schedule_id(), seed, kind_of_engine_to_make);
37 }
long seed
Definition: chem4.cc:68
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
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().

43 {
44  return rng()->createEngine(
45  placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
46 }
long seed
Definition: chem4.cc:68
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
CurrentProcessingContext const * art::EDProducer::currentContext ( ) const
protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

121  {
122  return current_context_.get();
123  }
CPC_exempt_ptr current_context_
Definition: EDProducer.h:116
double evgen::GaisserParam::GaisserFlux ( double  e,
double  theta 
)
private

Definition at line 711 of file GaisserParam_module.cc.

References c1.

711  {
712 
713  double ct = cos(theta);
714  double di;
715  if(fSetParam){
716  // double gamma=2.77; // LVD spectrum: spectral index
717  // double A=1.84*0.14; // normalisation
718  double gamma = 2.7;
719  double A = 0.14;
720  double rc = 1.e-4; // fraction of prompt muons
721  double c1 = sqrt(1.-(1.-pow(ct,2.0))/pow( (1.+32./6370.) ,2.0)); // Earth curvature
722  double deltae = 2.06e-3 * (1030. / c1 - 120.); // muon energy loss in atmosphere
723  double em = e + deltae/2.;
724  double e1 = e + deltae;
725  double pdec = pow( (120. / (1030. / c1)), (1.04 / c1 / em)); // muon decay
726  di=A*pow(e1,-gamma)*(1./(1.+1.1*e1*c1/115.) + 0.054/(1.+1.1*e1*c1/850.) + rc)*pdec;
727  }
728  else{
729  double gamma=2.7; // spectral index
730  double A=0.14; // normalisation
731  double C = 3.64;
732  double gamma2 = 1.29;
733  double ct_star = sqrt( (pow(ct,2) + 0.102573*0.102573 - 0.068287*pow(ct,0.958633)
734  + 0.0407253*pow(ct,0.817285) )/(1.0 + 0.102573*0.102573 - 0.068287 + 0.0407253) );
735  double eMod = e*(1.0+C/(e*pow(ct_star,gamma2)));
736  di=A*pow(eMod,-gamma)*(1./(1.+1.1*e*ct_star/115.) + 0.054/(1.+1.1*e*ct_star/850.));
737  }
738 
739  return di;
740  } // GaisserFlux
TCanvas * c1
Definition: plotHisto.C:7
bool fSetParam
Which version of Gaissers Param.
Float_t e
Definition: plot.C:34
double evgen::GaisserParam::GaisserMuonFlux_Integrand ( Double_t *  x,
Double_t *  par 
)
private

Definition at line 743 of file GaisserParam_module.cc.

Referenced by MakePDF().

743  {
744 
745  //---- calculate the flux
746  double flux = 2.0*M_PI*sin(x[1])*GaisserFlux(x[0],x[1]);
747 
748  return flux;
749  } // MuonFluxIntegrand
Float_t x
Definition: compare.C:6
double GaisserFlux(double e, double theta)
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_().

52 {
53  auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
54  return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
55 }
std::vector< double > evgen::GaisserParam::GetBinning ( const TAxis *  axis,
bool  finalEdge = true 
)
private

Definition at line 752 of file GaisserParam_module.cc.

752  {
753  std::vector<double> vbins;
754  for(int ibin=1; ibin<=axis->GetNbins()+1; ibin++){
755  if(ibin<=axis->GetNbins()) vbins.push_back(axis->GetBinLowEdge(ibin));
756  else if(finalEdge) vbins.push_back(axis->GetBinLowEdge(ibin) + axis->GetBinWidth(ibin));
757  }
758  return vbins;
759  } // Get Binning
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_.

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
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().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
Int_t B
Definition: plot.C:25
std::pair< double, double > evgen::GaisserParam::GetThetaAndEnergy ( double  rand1,
double  rand2 
)
private

Definition at line 466 of file GaisserParam_module.cc.

References energy.

467  {
468  if(rand1 < 0 || rand1 > 1) std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand1 << std::endl;
469  if(rand2 < 0 || rand2 > 1) std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand2 << std::endl;
470 
471  int thetaBin = 0;
472  m_thetaHist->GetBinWithContent(double(rand1),thetaBin,1,m_thetaHist->GetNbinsX(),1.0);
473  if(m_thetaHist->GetBinContent(thetaBin) < rand1) thetaBin += 1;
474  double drand1 = (m_thetaHist->GetBinContent(thetaBin) - rand1)/(m_thetaHist->GetBinContent(thetaBin) - m_thetaHist->GetBinContent(thetaBin-1));
475  double thetaLow = m_thetaHist->GetXaxis()->GetBinLowEdge(thetaBin);
476  double thetaUp = m_thetaHist->GetXaxis()->GetBinUpEdge(thetaBin);
477  double theta = drand1*(thetaUp-thetaLow) + thetaLow;
478 
479  int energyBin = 0;
480  TH1* energyHist = 0;
481  bool notFound = true;
482  for(dhist_Map_it mapit=m_PDFmap->begin(); mapit!=m_PDFmap->end(); mapit++){
483  if( fabs(mapit->first+thetaLow)<0.000001 ) {
484  energyHist = mapit->second;
485  notFound = false;
486  break;
487  }
488  }
489  if(notFound) std::cout << "GetThetaAndEnergy: ERROR:\tInvalid theta!" << std::endl;
490  // MSG("string = " << To_TString(thetaLow) << ", double = " << thetaLow );
491 
492  energyHist->GetBinWithContent(double(rand2),energyBin,1,energyHist->GetNbinsX(),1.0);
493  if(energyHist->GetBinContent(energyBin) < rand2) energyBin += 1;
494  double drand2 = (energyHist->GetBinContent(energyBin) - rand2)/(energyHist->GetBinContent(energyBin) - energyHist->GetBinContent(energyBin-1));
495  double energyLow = energyHist->GetXaxis()->GetBinLowEdge(energyBin);
496  double energyUp = energyHist->GetXaxis()->GetBinUpEdge(energyBin);
497  double energy = drand2*(energyUp-energyLow) + energyLow;
498  // MSG("MuFlux::GetThetaEnergy()\te = " << energy*1000. );
499 
500  return std::make_pair(theta,energy);
501  }
double energy
Definition: plottest35.C:25
std::map< double, TH1 * >::iterator dhist_Map_it
void evgen::GaisserParam::MakePDF ( )
private

Definition at line 504 of file GaisserParam_module.cc.

References GaisserMuonFlux_Integrand().

505  {
506  std::cout << "In my function MakePDF" << std::endl;
507  m_File=0;
508  m_PDFmap=0;
509  m_thetaHist=0;
510  double TotalMuonFlux=0;
511 
512  if(m_PDFmap){
513  std::cout << "PDFMAP" << std::endl;
514  if(m_PDFmap->size()>0 && !fSetReWrite){
515  std::cout << "MakePDF: Map has already been initialised. " << std::endl;
516  std::cout << "Do fSetReWrite - true if you really want to override this map." << std::endl;
517  return;
518  }
519  std::cout << fSetReWrite << std::endl;
520  if(fSetReWrite) ResetMap();
521  }
522  else{
523  m_PDFmap = new dhist_Map;
524  std::cout << "Making new dhist_Map called m_PDFmap....." << std::endl;
525  }
526 
527  TF2* muonSpec = new TF2("muonSpec", this,
530  "GaisserParam", "GaisserMuonFlux_Integrand"
531  );
532  //--------------------------------------------
533  //------------ Compute the pdfs
534 
535  //---- compute pdf for the theta
536  TotalMuonFlux = muonSpec->Integral(fEmin, fEmax, fThetamin, fThetamax, fEpsilon ); // Work out the muon flux at the surface
537  std::cout << "Surface flux of muons = " << TotalMuonFlux << " cm-2 s-1" << std::endl;
538 
539  //---- work out if we're reading a file, writing to file, or neither
540  std::ostringstream pdfFile;
541  pdfFile << "GaisserPDF_"<<fEmin<<"-"<<fEmid<<"-"<<fEmax<<"-"<< fEBinsLow<<"-"<<fEBinsHigh<<"-"<<fThetamin<<"-"<<fThetamax<<"-"<<fThetaBins<<".root";
542  std::string tmpfileName = pdfFile.str();
543  std::replace(tmpfileName.begin(),tmpfileName.end(),'+','0');
544  if (tmpfileName == "GaisserPDF_0-100-100100-1000-10000-0-1.5708-100.root") tmpfileName = "GaisserPDF_DefaultBins.root";
545  else if (tmpfileName == "GaisserPDF_0-100-4000-1000-1000-0-1.5708-100.root") tmpfileName = "GaisserPDF_LowEnergy.root";
546  else if (tmpfileName == "GaisserPDF_4000-10000-100000-1000-10000-0-1.5708-100.root") tmpfileName = "GaisserPDF_MidEnergy.root";
547  else if (tmpfileName == "GaisserPDF_100000-500000-1e007-10000-100000-0-1.5708-100.root") tmpfileName = "GaisserPDF_HighEnergy.root";
548 
549  std::ostringstream pdfFilePath;
550  pdfFilePath << fInputDir << tmpfileName;
551  std::string fileName = pdfFilePath.str();
552 
553  // fTemplateFile = pset.get< std::string >("TemplateFile");
554  // //fCalorimetryModuleLabel = pset.get< std::string >("CalorimetryModuleLabel");
555 
556  cet::search_path sp("FW_SEARCH_PATH");
557  std::string fROOTfile; //return /lbne/data/0-100-1.57.root
558  if( sp.find_file(tmpfileName, fROOTfile) ) fileName = fROOTfile;
559 
560  std::cout << "File path; " << fileName << std::endl;
561 
562  if(fSetRead){
563  struct stat buffer;
564  fSetRead = stat(fileName.c_str(), &buffer) == 0; // Check if file exists already
565  if(!fSetRead) std::cout << "WARNING- "+fileName+" does not exist." << std::endl;
566  else{
567  std::cout << "Reading PDF from file "+fileName << std::endl;
568  m_File = new TFile(fileName.c_str()); // Open file
569  if(m_File->IsZombie() || m_File->TestBit(TFile::kRecovered)){ // Check that file is not corrupted
570  std::cout << "WARNING- "+fileName+" is corrupted or cannot be read." << std::endl;
571  fSetRead = false;
572  }
573  }
574  }
575 
576  if(fSetRead){ // If the file exists then read it....
577  std::cout << "Now going to read file...." << std::endl;
578  fSetWrite = false; // Don't want to write as already exists
579  double thetalow = fThetamin;
580  m_thetaHist = (TH1D*) m_File->Get("pdf_theta");
581  for(int i=1; i<=fThetaBins; i++){
582  std::ostringstream pdfEnergyHist;
583  pdfEnergyHist << "pdf_energy_"<<i;
584  std::string pdfEnergyHiststr = pdfEnergyHist.str();
585 
586  TH1* pdf_hist = (TH1D*) m_File->Get( pdfEnergyHiststr.c_str() ); // Get the bin
587  m_PDFmap->insert(std::make_pair(-thetalow,pdf_hist)); //---- -ve theta for quicker sorting
588  thetalow = double(i)*(fThetamax)/double(fThetaBins); //---- increment the value of theta
589  }
590  } // ------------------------------------------
591  else { // File doesn't exist so want to make it.....
592  // ------------------------------------------
593  std::cout << "Generating a new muon flux PDF" << std::endl;
594  if(fSetWrite){
595  std::cout << "Writing to PDF to file "+fileName << std::endl;
596  m_File = new TFile(fileName.c_str(),"RECREATE");
597  }
598 
599  double dnbins_theta = double(fThetaBins);
600  m_thetaHist = new TH1D("pdf_theta", "pdf_theta", fThetaBins, fThetamin, fThetamax);
601  for(int i=1; i<=fThetaBins; i++){
602  double di = i==0 ? 0.1 : double(i);
603  double theta = di*(fThetamax)/dnbins_theta;
604  double int_i = muonSpec->Integral(fEmin, fEmax, fThetamin, theta, fEpsilon);
605  m_thetaHist->SetBinContent(i, int_i/TotalMuonFlux);
606  }
607  if(fSetWrite) m_thetaHist->Write();
608  std::cout << "theta PDF complete... now making the energy PDFs (this will take a little longer)... " << std::endl;
609 
610  //---- now compute the energy pdf
611  double thetalow = fThetamin;
612  for(int i=1; i<=fThetaBins; i++){
613 
614  double di = double(i);
615  double theta = di*(fThetamax)/fThetaBins;
616 
617  //---- compute the total integral
618  double int_tot = muonSpec->Integral(fEmin, fEmax, thetalow, theta, fEpsilon);
619 
620  //---- compute pdf for the low energy
621  int nbins = fEBinsLow;
622  TH1* pdf_lowenergy = new TH1D("pdf_lowenergy", "pdf_lowenergy", nbins, fEmin, fEmid);
623  double dnbins = double(nbins);
624  for(int j=1; j<=nbins; j++){
625  double dj = double(j);
626  double int_j = muonSpec->Integral(fEmin, fEmin + dj*(fEmid-fEmin)/dnbins, thetalow, theta, fEpsilon);
627  // std::std::cout << j << "(" << m_emin << " --> " << m_emin + dj*m_emid/dnbins << ") = " << int_j/int_tot << std::std::endl;
628  // std::std::cout << j << "(" << m_emin << " --> " << m_emin + dj*(m_emid-m_emin)/dnbins << ") = " << int_j/int_tot << std::std::endl;
629  pdf_lowenergy->SetBinContent(j, int_j/int_tot);
630  }
631 
632  //---- compute pdf for the high energy
633  nbins = fEBinsHigh;
634  dnbins=double(nbins);
635  TH1* pdf_highenergy = new TH1D("pdf_highenergy", "pdf_highenergy", nbins, fEmid, fEmax);
636  for(int j=1; j<=nbins; j++){
637  double dj = double(j);
638  double int_j = muonSpec->Integral(fEmin, fEmid + dj*(fEmax-fEmid)/dnbins, thetalow, theta, fEpsilon);
639  // std::cout << j << "(" << m_emin << " --> " << m_emid + dj*(m_emax-m_emid)/dnbins << ") = " << int_j/int_tot << std::endl;
640  pdf_highenergy->SetBinContent(j, int_j/int_tot);
641  }
642 
643  //---- now combine the two energy hists
644  std::vector<double> vxbins = GetBinning(pdf_lowenergy->GetXaxis(),false);
645  std::vector<double> vxbins2 = GetBinning(pdf_highenergy->GetXaxis());
646  vxbins.insert(vxbins.end(), vxbins2.begin(), vxbins2.end());
647 
648  int ibin = 0;
649  double* xbins = new double[vxbins.size()];
650  for(std::vector<double>::const_iterator binit=vxbins.begin(); binit!=vxbins.end(); binit++, ibin++) xbins[ibin]=(*binit);
651  TH1* pdf_energy = new TH1D("pdf_energy", "pdf_energy", vxbins.size()-1, xbins);
652  int ibin2 = 1;
653  for(ibin = 1; ibin<=pdf_lowenergy->GetNbinsX(); ibin++, ibin2++){
654  double content = pdf_lowenergy->GetBinContent(ibin);
655  if(ibin == 1) content = content - 0.00001;
656  pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), content);
657  }
658  for(ibin = 1; ibin<=pdf_highenergy->GetNbinsX(); ibin++, ibin2++){
659  pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), pdf_highenergy->GetBinContent(ibin));
660  }
661 
662  //---- and remove any negative bins
663  std::ostringstream Clonestr;
664  Clonestr << "pdf_energy_"<<i;
665  TH1* pdf_energy_noneg = (TH1D*) pdf_energy->Clone( Clonestr.str().c_str() );
666  pdf_energy_noneg->Reset();
667 
668  double PDF = 0.0;
669  double lastPD = 0.0;
670  int nSkip=0;
671  for(ibin = 1; ibin<=pdf_energy->GetNbinsX(); ibin++){
672  double newPD = pdf_energy->GetBinContent(ibin);
673  double probDiff = newPD - lastPD;
674  if(probDiff<0){
675  if(ibin!=pdf_energy->GetNbinsX()){
676  nSkip++;
677  continue;
678  }
679  else probDiff = 0;
680  }
681  else PDF += probDiff;
682  if(PDF > 1) PDF = 1;
683  for(int iskip=0; iskip <= nSkip; iskip++) pdf_energy_noneg->Fill(pdf_energy_noneg->GetBinCenter(ibin-iskip), PDF);
684  nSkip=0;
685  lastPD = newPD;
686  }
687 
688 
689  //---- add this hist, increment thetalow and delete unwanted TH1s
690  if(fSetWrite) pdf_energy_noneg->Write();
691  m_PDFmap->insert(std::make_pair(-thetalow,pdf_energy_noneg)); //---- -ve theta for quicker sorting
692 
693  //---- increment the value of theta
694  thetalow = theta;
695 
696  //---- free up memory from unwanted hists
697  delete pdf_lowenergy;
698  delete pdf_highenergy;
699  delete pdf_energy;
700 
701  std::cout << "\r===> " << 100.0*double(i)/double(fThetaBins) << "% complete... " << std::endl;
702  } // ThetaBins
703  std::cout << "finished the energy pdfs." << std::endl;
704  }//---- if(!m_doRead)
705 
706  delete muonSpec;
707  return;
708  } // Make PDF
double fEpsilon
Minimum integration sum....
double fThetamax
Maximum theta.
int fThetaBins
Number of theta Bins.
double fEmid
Energy to go from low to high (GeV)
bool fSetReWrite
Whether to ReWrite pdfs.
std::string fInputDir
Input Directory.
int fEBinsLow
Number of low energy Bins.
intermediate_table::const_iterator const_iterator
std::map< double, TH1 * > dhist_Map
double fEmin
Minimum Kinetic Energy (GeV)
std::vector< double > GetBinning(const TAxis *axis, bool finalEdge=true)
double fEmax
Maximum Kinetic Energy (GeV)
int fEBinsHigh
Number of high energy Bins.
bool fSetRead
Whether to Read.
double fThetamin
Minimum theta.
double GaisserMuonFlux_Integrand(Double_t *x, Double_t *par)
bool fSetWrite
Whether to Write.
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().

191 {
192  if (!moduleContext_)
193  return ProductToken<T>::invalid();
194 
195  consumables_[BT].emplace_back(ConsumableType::Product,
196  TypeID{typeid(T)},
197  it.label(),
198  it.instance(),
199  it.process());
200  return ProductToken<T>{it};
201 }
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 205 of file Consumer.h.

206 {
207  if (!moduleContext_)
208  return;
209 
210  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
211 }
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
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().

216 {
217  if (!moduleContext_)
218  return ViewToken<T>::invalid();
219 
220  consumables_[BT].emplace_back(ConsumableType::ViewElement,
221  TypeID{typeid(T)},
222  it.label(),
223  it.instance(),
224  it.process());
225  return ViewToken<T>{it};
226 }
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

41  {
42  return true;
43  }
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().

90 {
91  if (!moduleContext_)
92  return;
93 
94  pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
95  for (auto& consumablesPerBranch : consumables_) {
96  cet::sort_all(consumablesPerBranch);
97  }
98 }
bool requireConsumes_
Definition: Consumer.h:137
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
void evgen::GaisserParam::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 302 of file GaisserParam_module.cc.

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

303  {
305  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
306 
307  simb::MCTruth truth;
309 
310  Sample(truth);
311 
312  LOG_DEBUG("GaisserParam") << truth;
313 
314  truthcol->push_back(truth);
315 
316  evt.put(std::move(truthcol));
317 
318  return;
319  }
void Sample(simb::MCTruth &mct)
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:78
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
single particles thrown at the detector
Definition: MCTruth.h:24
#define LOG_DEBUG(id)
Event generator information.
Definition: MCTruth.h:30
void evgen::GaisserParam::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 196 of file GaisserParam_module.cc.

References fhicl::ParameterSet::get().

197  {
198  // do not put seed in reconfigure because we don't want to reset
199  // the seed midstream
200 
201  fPadOutVectors = p.get< bool >("PadOutVectors");
202  fMode = p.get< int >("ParticleSelectionMode");
203  fPDG = p.get< std::vector<int> >("PDG");
204 
205  fCharge = p.get< int >("Charge");
206  fInputDir = p.get< std::string >("InputDir");
207 
208  fEmin = p.get< double >("Emin");
209  fEmax = p.get< double >("Emax");
210  fEmid = p.get< double >("Emid");
211  fEBinsLow = p.get< int >("EBinsLow");
212  fEBinsHigh = p.get< int >("EBinsHigh");
213 
214  fThetamin = p.get< double >("Thetamin");
215  fThetamax = p.get< double >("Thetamax");
216  fThetaBins = p.get< int >("ThetaBins");
217 
218  fXHalfRange = p.get< double >("XHalfRange");
219  fYInput = p.get< double >("YInput");
220  fZHalfRange = p.get< double >("ZHalfRange");
221 
222  fT0 = p.get< double >("T0");
223  fSigmaT = p.get< double >("SigmaT");
224  fTDist = p.get< int >("TDist");
225 
226  fSetParam = p.get< bool >("SetParam");
227  fSetRead = p.get< bool >("SetRead");
228  fSetWrite = p.get< bool >("SetWrite");
229  fSetReWrite = p.get< bool >("SetReWrite");
230  fEpsilon = p.get< double >("Epsilon");
231 
232  return;
233  }
double fEpsilon
Minimum integration sum....
double fThetamax
Maximum theta.
double fYInput
Max Y position.
int fThetaBins
Number of theta Bins.
double fEmid
Energy to go from low to high (GeV)
bool fSetReWrite
Whether to ReWrite pdfs.
std::string fInputDir
Input Directory.
int fEBinsLow
Number of low energy Bins.
double fEmin
Minimum Kinetic Energy (GeV)
double fSigmaT
Variation in t position (ns)
bool fSetParam
Which version of Gaissers Param.
double fEmax
Maximum Kinetic Energy (GeV)
double fXHalfRange
Max X position.
int fEBinsHigh
Number of high energy Bins.
std::vector< int > fPDG
PDG code of particles to generate.
bool fSetRead
Whether to Read.
double fThetamin
Minimum theta.
int fTDist
How to distribute t (gaus, or uniform)
double fT0
Central t position (ns) in world coordinates.
bool fSetWrite
Whether to Write.
double fZHalfRange
Max Z position.
void evgen::GaisserParam::ResetMap ( )
private

Definition at line 762 of file GaisserParam_module.cc.

References DEFINE_ART_MODULE.

762  {
763  if(m_thetaHist) delete m_thetaHist;
764  if(m_PDFmap){
765  for(dhist_Map_it mapit=m_PDFmap->begin(); mapit!=m_PDFmap->end(); mapit++){
766  if(mapit->second) delete mapit->second;
767  }
768  m_PDFmap->clear();
769  delete m_PDFmap;
770  std::cout << "Reset PDFmap and thetaHist..." << std::endl;
771  }
772  } // ResetMap
std::map< double, TH1 * >::iterator dhist_Map_it
void evgen::GaisserParam::Sample ( simb::MCTruth mct)
private

Definition at line 436 of file GaisserParam_module.cc.

References art::RandomNumberGenerator::getEngine().

437  {
438  switch (fMode) {
439  case 0: // List generation mode: every event will have one of each
440  // particle species in the fPDG array
441  for (unsigned int i=0; i<fPDG.size(); ++i) {
442  SampleOne(i,mct);
443  }//end loop over particles
444  break;
445  case 1: // Random selection mode: every event will exactly one particle
446  // selected randomly from the fPDG array
447  {
449  CLHEP::HepRandomEngine &engine = rng->getEngine();
450  CLHEP::RandFlat flat(engine);
451 
452  unsigned int i=flat.fireInt(fPDG.size());
453  SampleOne(i,mct);
454  }
455  break;
456  default:
457  mf::LogWarning("UnrecognizeOption") << "GaisserParam does not recognize ParticleSelectionMode "
458  << fMode;
459  break;
460  } // switch on fMode
461 
462  return;
463  }
base_engine_t & getEngine() const
void SampleOne(unsigned int i, simb::MCTruth &mct)
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< int > fPDG
PDG code of particles to generate.
void evgen::GaisserParam::SampleOne ( unsigned int  i,
simb::MCTruth mct 
)
private

Definition at line 324 of file GaisserParam_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), art::RandomNumberGenerator::getEngine(), part, and x.

324  {
325 
326  // get the random number generator service and make some CLHEP generators
328  CLHEP::HepRandomEngine &engine = rng->getEngine();
329  CLHEP::RandFlat flat(engine);
330  CLHEP::RandGaussQ gauss(engine);
331 
332  Time = Momentum = KinEnergy = Gamma = Energy = Theta = Phi = pnorm = 0.0;
334 
335  TVector3 x;
336  TVector3 pmom;
337 
338  // set track id to -i as these are all primary particles and have id <= 0
339  int trackid = -1*(i+1);
340  std::string primary("primary");
341 
342  // Work out whether particle/antiparticle, and mass...
343  double m =0.0;
344  fPDG[i] = 13;
345  if (fCharge == 0 ) {
346  if(1.0-2.0*flat.fire() > 0) fPDG[i]=-fPDG[i];
347  }
348  else if ( fCharge == 1 ) fPDG[i] = 13;
349  else if ( fCharge == 2 ) fPDG[i] = -13;
350  static TDatabasePDG pdgt;
351  TParticlePDG* pdgp = pdgt.GetParticle(fPDG[i]);
352  if (pdgp) m = pdgp->Mass();
353 
354  // Work out T0...
355  if(fTDist==kGAUS){
356  Time = gauss.fire(fT0, fSigmaT);
357  }
358  else {
359  Time = fT0 + fSigmaT*(2.0*flat.fire()-1.0);
360  }
361 
362  // Work out Positioning.....
363  x[0] = (1.0 - 2.0*flat.fire())*fXHalfRange + fCenterX;
364  x[1] = fYInput;
365  x[2] = (1.0 - 2.0*flat.fire())*fZHalfRange + fCenterZ;
366 
367  // Make Lorentz vector for x and t....
368  TLorentzVector pos(x[0], x[1], x[2], Time);
369 
370  // Access the pdf map which has been loaded.....
371  if(m_PDFmap) {
372 
373  //---- get the muon theta and energy from histograms using 2 random numbers
374  std::pair<double,double> theta_energy; //---- muon theta and energy
375  theta_energy = GetThetaAndEnergy(flat.fire(),flat.fire());
376 
377  //---- Set theta, phi
378  Theta = theta_energy.first; // Angle returned by GetThetaAndEnergy between 0 and pi/2
379  Phi = M_PI*( 1.0-2.0*flat.fire() ); // Randomly generated angle between -pi and pi
380 
381  //---- Set KE, E, p
382  KinEnergy = theta_energy.second; // Energy returned by GetThetaAndEnergy
383  Gamma = 1 + (KinEnergy/m);
384  Energy = Gamma * m;
385  Momentum = std::sqrt(Energy*Energy-m*m); // Get momentum
386 
387  pmom[0] = Momentum*std::sin(Theta)*std::cos(Phi);
388  pmom[1] = -Momentum*std::cos(Theta);
389  pmom[2] = Momentum*std::sin(Theta)*std::sin(Phi);
390 
391  pnorm = std::sqrt( pmom[0]*pmom[0] + pmom[1]*pmom[1] + pmom[2]*pmom[2] );
392  DirCosineX = pmom[0] / pnorm;
393  DirCosineY = pmom[1] / pnorm;
394  DirCosineZ = pmom[2] / pnorm;
395  }
396  else {
397  std::cout << "MuFlux map hasn't been initialised, aborting...." << std::endl;
398  return;
399  }
400 
401  TLorentzVector pvec(pmom[0], pmom[1], pmom[2], Energy );
402 
403  simb::MCParticle part(trackid, fPDG[i], primary);
404  part.AddTrajectoryPoint(pos, pvec);
405  /*
406  fTree->Branch();
407  fTree->Branch();
408  fTree->Branch();
409  fPositionX ->Fill (x[0]);
410  fPositionY ->Fill (x[1]);
411  fPositionZ ->Fill (x[2]);
412  fTime ->Fill (Time);
413  fMomentumHigh ->Fill (Momentum);
414  fMomentum ->Fill (Momentum);
415  fEnergy ->Fill (Energy);
416  fDirCosineX ->Fill (DirCosineX);
417  fDirCosineY ->Fill (DirCosineY);
418  fDirCosineZ ->Fill (DirCosineZ);
419  fTheta ->Fill (Theta*180/M_PI);
420  fPhi ->Fill (Phi *180/M_PI);
421  */
422  XPosition = x[0];
423  YPosition = x[1];
424  ZPosition = x[2];
425 
426  std::cout << "Theta: " << Theta << " Phi: " << Phi << " KineticEnergy: " << KinEnergy << " Energy: " << Energy << " Momentum: " << Momentum << std::endl;
427  std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
428  std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
429  std::cout << "Normalised..." << DirCosineX << " " << DirCosineY << " " << DirCosineZ << std::endl;
430 
431  fTree->Fill();
432  mct.Add(part);
433  }
Float_t x
Definition: compare.C:6
double fYInput
Max Y position.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:77
base_engine_t & getEngine() const
TString part[npart]
Definition: Style.C:32
double fSigmaT
Variation in t position (ns)
double fXHalfRange
Max X position.
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
std::vector< int > fPDG
PDG code of particles to generate.
int fTDist
How to distribute t (gaus, or uniform)
double fT0
Central t position (ns) in world coordinates.
std::pair< double, double > GetThetaAndEnergy(double rand1, double rand2)
double fZHalfRange
Max Z position.
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().

126 {
127  if (!moduleContext_)
128  return;
129 
130  // If none of the branches have missing consumes statements, exit early.
131  if (std::all_of(cbegin(missingConsumes_),
132  cend(missingConsumes_),
133  [](auto const& perBranch) { return perBranch.empty(); }))
134  return;
135 
136  constexpr cet::HorizontalRule rule{60};
137  mf::LogPrint log{"MTdiagnostics"};
138  log << '\n'
139  << rule('=') << '\n'
140  << "The following consumes (or mayConsume) statements are missing from\n"
141  << module_context(moduleDescription_) << '\n'
142  << rule('-') << '\n';
143 
144  cet::for_all_with_index(
145  missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
146  for (auto const& pi : perBranch) {
147  log << " "
148  << assemble_consumes_statement(static_cast<BranchType>(i), pi)
149  << '\n';
150  }
151  });
152  log << rule('=');
153 }
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
bool moduleContext_
Definition: Consumer.h:136
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

103 {
104  // Early exits if consumes tracking has been disabled or if the
105  // consumed product is an allowed consumable.
106  if (!moduleContext_)
107  return;
108 
109  if (cet::binary_search_all(consumables_[bt], pi))
110  return;
111 
112  if (requireConsumes_) {
114  "Consumer: an error occurred during validation of a "
115  "retrieved product\n\n")
116  << "The following consumes (or mayConsume) statement is missing from\n"
117  << module_context(moduleDescription_) << ":\n\n"
118  << " " << assemble_consumes_statement(bt, pi) << "\n\n";
119  }
120 
121  missingConsumes_[bt].insert(pi);
122 }
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
bool requireConsumes_
Definition: Consumer.h:137
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139

Member Data Documentation

double evgen::GaisserParam::DirCosineX
private

Definition at line 170 of file GaisserParam_module.cc.

double evgen::GaisserParam::DirCosineY
private

Definition at line 170 of file GaisserParam_module.cc.

double evgen::GaisserParam::DirCosineZ
private

Definition at line 170 of file GaisserParam_module.cc.

double evgen::GaisserParam::Energy
private

Definition at line 169 of file GaisserParam_module.cc.

double evgen::GaisserParam::fCenterX = 0
private

Definition at line 164 of file GaisserParam_module.cc.

double evgen::GaisserParam::fCenterZ = 0
private

Definition at line 165 of file GaisserParam_module.cc.

int evgen::GaisserParam::fCharge
private

Charge.

Definition at line 116 of file GaisserParam_module.cc.

double evgen::GaisserParam::fCryoBoundaries[6]
private

Definition at line 159 of file GaisserParam_module.cc.

int evgen::GaisserParam::fEBinsHigh
private

Number of high energy Bins.

Definition at line 123 of file GaisserParam_module.cc.

int evgen::GaisserParam::fEBinsLow
private

Number of low energy Bins.

Definition at line 122 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEmax
private

Maximum Kinetic Energy (GeV)

Definition at line 120 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEmid
private

Energy to go from low to high (GeV)

Definition at line 121 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEmin
private

Minimum Kinetic Energy (GeV)

Definition at line 119 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEpsilon
private

Minimum integration sum....

Definition at line 141 of file GaisserParam_module.cc.

std::string evgen::GaisserParam::fInputDir
private

Input Directory.

Definition at line 117 of file GaisserParam_module.cc.

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 108 of file GaisserParam_module.cc.

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 111 of file GaisserParam_module.cc.

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

PDG code of particles to generate.

Definition at line 115 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetParam
private

Which version of Gaissers Param.

Definition at line 137 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetRead
private

Whether to Read.

Definition at line 138 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetReWrite
private

Whether to ReWrite pdfs.

Definition at line 140 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetWrite
private

Whether to Write.

Definition at line 139 of file GaisserParam_module.cc.

double evgen::GaisserParam::fSigmaT
private

Variation in t position (ns)

Definition at line 134 of file GaisserParam_module.cc.

double evgen::GaisserParam::fT0
private

Central t position (ns) in world coordinates.

Definition at line 133 of file GaisserParam_module.cc.

int evgen::GaisserParam::fTDist
private

How to distribute t (gaus, or uniform)

Definition at line 135 of file GaisserParam_module.cc.

int evgen::GaisserParam::fThetaBins
private

Number of theta Bins.

Definition at line 127 of file GaisserParam_module.cc.

double evgen::GaisserParam::fThetamax
private

Maximum theta.

Definition at line 126 of file GaisserParam_module.cc.

double evgen::GaisserParam::fThetamin
private

Minimum theta.

Definition at line 125 of file GaisserParam_module.cc.

TTree* evgen::GaisserParam::fTree
private

Definition at line 168 of file GaisserParam_module.cc.

double evgen::GaisserParam::fXHalfRange
private

Max X position.

Definition at line 129 of file GaisserParam_module.cc.

double evgen::GaisserParam::fYInput
private

Max Y position.

Definition at line 130 of file GaisserParam_module.cc.

double evgen::GaisserParam::fZHalfRange
private

Max Z position.

Definition at line 131 of file GaisserParam_module.cc.

double evgen::GaisserParam::Gamma
private

Definition at line 169 of file GaisserParam_module.cc.

const int evgen::GaisserParam::kGAUS = 1
staticprivate

Definition at line 106 of file GaisserParam_module.cc.

double evgen::GaisserParam::KinEnergy
private

Definition at line 169 of file GaisserParam_module.cc.

TFile* evgen::GaisserParam::m_File

Definition at line 89 of file GaisserParam_module.cc.

dhist_Map* evgen::GaisserParam::m_PDFmap

Definition at line 90 of file GaisserParam_module.cc.

TH1* evgen::GaisserParam::m_thetaHist

Definition at line 91 of file GaisserParam_module.cc.

double evgen::GaisserParam::Momentum
private

Definition at line 169 of file GaisserParam_module.cc.

double evgen::GaisserParam::Phi
private

Definition at line 169 of file GaisserParam_module.cc.

double evgen::GaisserParam::pnorm
private

Definition at line 169 of file GaisserParam_module.cc.

double evgen::GaisserParam::Theta
private

Definition at line 169 of file GaisserParam_module.cc.

double evgen::GaisserParam::Time
private

Definition at line 169 of file GaisserParam_module.cc.

double evgen::GaisserParam::xNeg = 0
private

Definition at line 160 of file GaisserParam_module.cc.

double evgen::GaisserParam::xPos = 0
private

Definition at line 161 of file GaisserParam_module.cc.

double evgen::GaisserParam::XPosition
private

Definition at line 170 of file GaisserParam_module.cc.

double evgen::GaisserParam::YPosition
private

Definition at line 170 of file GaisserParam_module.cc.

double evgen::GaisserParam::zNeg = 0
private

Definition at line 162 of file GaisserParam_module.cc.

double evgen::GaisserParam::zPos = 0
private

Definition at line 163 of file GaisserParam_module.cc.

double evgen::GaisserParam::ZPosition
private

Definition at line 170 of file GaisserParam_module.cc.


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