Public Member Functions | Private Types | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
evgen::GaisserParam Class Reference

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

Inheritance diagram for evgen::GaisserParam:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 GaisserParam (fhicl::ParameterSet const &pset)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private 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_Mapm_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
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

module to produce single or multiple specified particles in the detector

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

Constructor & Destructor Documentation

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

Definition at line 162 of file GaisserParam_module.cc.

163  : art::EDProducer{pset}
164  // create a default random engine; obtain the random seed from NuRandomService,
165  // unless overridden in configuration with key "Seed"
166  , fEngine(art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*this, pset, "Seed"))
167  , fMode {pset.get< int >("ParticleSelectionMode")}
168  , fPadOutVectors{pset.get< bool >("PadOutVectors")}
169  , fPDG {pset.get< std::vector<int> >("PDG")}
170  , fCharge {pset.get< int >("Charge")}
171  , fInputDir {pset.get< std::string >("InputDir")}
172  , fEmin {pset.get< double >("Emin")}
173  , fEmax {pset.get< double >("Emax")}
174  , fEmid {pset.get< double >("Emid")}
175  , fEBinsLow {pset.get< int >("EBinsLow")}
176  , fEBinsHigh {pset.get< int >("EBinsHigh")}
177  , fThetamin {pset.get< double >("Thetamin")}
178  , fThetamax {pset.get< double >("Thetamax")}
179  , fThetaBins {pset.get< int >("ThetaBins")}
180  , fXHalfRange {pset.get< double >("XHalfRange")}
181  , fYInput {pset.get< double >("YInput")}
182  , fZHalfRange {pset.get< double >("ZHalfRange")}
183  , fT0 {pset.get< double >("T0")}
184  , fSigmaT {pset.get< double >("SigmaT")}
185  , fTDist {pset.get< int >("TDist")}
186  , fSetParam {pset.get< bool >("SetParam")}
187  , fSetRead {pset.get< bool >("SetRead")}
188  , fSetWrite {pset.get< bool >("SetWrite")}
189  , fSetReWrite {pset.get< bool >("SetReWrite")}
190  , fEpsilon {pset.get< double >("Epsilon")}
191  {
192  produces< std::vector<simb::MCTruth> >();
193  produces< sumdata::RunData, art::InRun >();
194  }
double fEpsilon
Minimum integration sum....
double fThetamax
Maximum theta.
double fYInput
Max Y position.
std::string string
Definition: nybbler.cc:12
int fThetaBins
Number of theta Bins.
double fEmid
Energy to go from low to high (GeV)
bool fSetReWrite
Whether to ReWrite pdfs.
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
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.

Member Function Documentation

void evgen::GaisserParam::beginJob ( )
overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 197 of file GaisserParam_module.cc.

198  {
199  //Work out center of cryostat(s)
201  for (unsigned int i=0; i < geom->Ncryostats() ; i++ ) {
203  if ( xNeg > fCryoBoundaries[0] ) xNeg = fCryoBoundaries[0];
204  if ( xPos < fCryoBoundaries[1] ) xPos = fCryoBoundaries[1];
205  if ( zNeg > fCryoBoundaries[4] ) zNeg = fCryoBoundaries[4];
206  if ( zPos < fCryoBoundaries[5] ) zPos = fCryoBoundaries[5];
207  }
208  fCenterX = xNeg + (xPos-xNeg)/2;
209  fCenterZ = zNeg + (zPos-zNeg)/2;
210 
211  // Make the Histograms....
213  /*
214  fPositionX = tfs->make<TH1D>("fPositionX" ,"Position (cm)" ,500,fCenterX-(fXHalfRange+10) ,fCenterX+(fXHalfRange+10));
215  fPositionY = tfs->make<TH1D>("fPositionY" ,"Position (cm)" ,500,-(fYInput+10),(fYInput+10));
216  fPositionZ = tfs->make<TH1D>("fPositionZ" ,"Position (cm)" ,500,fCenterZ-(fZHalfRange+10) ,fCenterZ+(fZHalfRange+10));
217  fTime = tfs->make<TH1D>("fTime" ,"Time (s)" ,500,0,1e6);
218  fMomentumHigh = tfs->make<TH1D>("fMomentumHigh","Momentum (GeV)",500,0,fEmax);
219  fMomentum = tfs->make<TH1D>("fMomentum" ,"Momentum (GeV)",500,0,100);
220  fEnergy = tfs->make<TH1D>("fEnergy" ,"Energy (GeV)" ,500,0,fEmax);
221 
222  fDirCosineX = tfs->make<TH1D>("fDirCosineX","Normalised Direction cosine",500,-1,1);
223  fDirCosineY = tfs->make<TH1D>("fDirCosineY","Normalised Direction cosine",500,-1,1);
224  fDirCosineZ = tfs->make<TH1D>("fDirCosineZ","Normalised Direction cosine",500,-1,1);
225 
226  fTheta = tfs->make<TH1D>("fTheta" ,"Angle (radians)",500,-365,365);
227  fPhi = tfs->make<TH1D>("fPhi" ,"Angle (radians)",500,-365,365);
228  */
229  fTree = tfs->make<TTree>("Generator","analysis tree");
230  fTree->Branch("XPosition" ,&XPosition ,"XPosition/D" );
231  fTree->Branch("YPosition" ,&YPosition ,"YPosition/D" );
232  fTree->Branch("ZPosition" ,&ZPosition ,"ZPosition/D" );
233  fTree->Branch("Time" ,&Time ,"Time/D" );
234  fTree->Branch("Momentum" ,&Momentum ,"Momentum/D" );
235  fTree->Branch("KinEnergy" ,&KinEnergy ,"KinEnergy/D" );
236  fTree->Branch("Energy" ,&Energy ,"Energy/D" );
237  fTree->Branch("DirCosineX",&DirCosineX,"DirCosineX/D");
238  fTree->Branch("DirCosineY",&DirCosineY,"DirCosineY/D");
239  fTree->Branch("DirCosineZ",&DirCosineZ,"DirCosineZ/D");
240  fTree->Branch("Theta" ,&Theta ,"Theta/D" );
241  fTree->Branch("Phi" ,&Phi ,"Phi/D" );
242  }
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)
overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 245 of file GaisserParam_module.cc.

246  {
247  // Check fcl parameters were set correctly
248  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";
249  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;
250  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;
251 
252  // grab the geometry object to see what geometry we are using
254  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
255  MakePDF ();
256  }
double fThetamax
Maximum theta.
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
#define M_PI
Definition: includeROOT.h:54
double fThetamin
Minimum theta.
LArSoft geometry interface.
Definition: ChannelGeo.h:16
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
double evgen::GaisserParam::GaisserFlux ( double  e,
double  theta 
)
private

Definition at line 661 of file GaisserParam_module.cc.

661  {
662 
663  double ct = cos(theta);
664  double di;
665  if(fSetParam){
666  // double gamma=2.77; // LVD spectrum: spectral index
667  // double A=1.84*0.14; // normalisation
668  double gamma = 2.7;
669  double A = 0.14;
670  double rc = 1.e-4; // fraction of prompt muons
671  double c1 = sqrt(1.-(1.-pow(ct,2.0))/pow( (1.+32./6370.) ,2.0)); // Earth curvature
672  double deltae = 2.06e-3 * (1030. / c1 - 120.); // muon energy loss in atmosphere
673  double em = e + deltae/2.;
674  double e1 = e + deltae;
675  double pdec = pow( (120. / (1030. / c1)), (1.04 / c1 / em)); // muon decay
676  di=A*pow(e1,-gamma)*(1./(1.+1.1*e1*c1/115.) + 0.054/(1.+1.1*e1*c1/850.) + rc)*pdec;
677  }
678  else{
679  double gamma=2.7; // spectral index
680  double A=0.14; // normalisation
681  double C = 3.64;
682  double gamma2 = 1.29;
683  double ct_star = sqrt( (pow(ct,2) + 0.102573*0.102573 - 0.068287*pow(ct,0.958633)
684  + 0.0407253*pow(ct,0.817285) )/(1.0 + 0.102573*0.102573 - 0.068287 + 0.0407253) );
685  double eMod = e*(1.0+C/(e*pow(ct_star,gamma2)));
686  di=A*pow(eMod,-gamma)*(1./(1.+1.1*e*ct_star/115.) + 0.054/(1.+1.1*e*ct_star/850.));
687  }
688 
689  return di;
690  } // GaisserFlux
constexpr T pow(T x)
Definition: pow.h:72
const double e
bool fSetParam
Which version of Gaissers Param.
double gamma(double KE, const simb::MCParticle *part)
double evgen::GaisserParam::GaisserMuonFlux_Integrand ( Double_t *  x,
Double_t *  par 
)
private

Definition at line 693 of file GaisserParam_module.cc.

693  {
694 
695  //---- calculate the flux
696  double flux = 2.0*M_PI*sin(x[1])*GaisserFlux(x[0],x[1]);
697 
698  return flux;
699  } // MuonFluxIntegrand
double GaisserFlux(double e, double theta)
#define M_PI
Definition: includeROOT.h:54
list x
Definition: train.py:276
std::vector< double > evgen::GaisserParam::GetBinning ( const TAxis *  axis,
bool  finalEdge = true 
)
private

Definition at line 702 of file GaisserParam_module.cc.

702  {
703  std::vector<double> vbins;
704  for(int ibin=1; ibin<=axis->GetNbins()+1; ibin++){
705  if(ibin<=axis->GetNbins()) vbins.push_back(axis->GetBinLowEdge(ibin));
706  else if(finalEdge) vbins.push_back(axis->GetBinLowEdge(ibin) + axis->GetBinWidth(ibin));
707  }
708  return vbins;
709  } // Get Binning
std::pair< double, double > evgen::GaisserParam::GetThetaAndEnergy ( double  rand1,
double  rand2 
)
private

Definition at line 416 of file GaisserParam_module.cc.

417  {
418  if(rand1 < 0 || rand1 > 1) std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand1 << std::endl;
419  if(rand2 < 0 || rand2 > 1) std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand2 << std::endl;
420 
421  int thetaBin = 0;
422  m_thetaHist->GetBinWithContent(double(rand1),thetaBin,1,m_thetaHist->GetNbinsX(),1.0);
423  if(m_thetaHist->GetBinContent(thetaBin) < rand1) thetaBin += 1;
424  double drand1 = (m_thetaHist->GetBinContent(thetaBin) - rand1)/(m_thetaHist->GetBinContent(thetaBin) - m_thetaHist->GetBinContent(thetaBin-1));
425  double thetaLow = m_thetaHist->GetXaxis()->GetBinLowEdge(thetaBin);
426  double thetaUp = m_thetaHist->GetXaxis()->GetBinUpEdge(thetaBin);
427  double theta = drand1*(thetaUp-thetaLow) + thetaLow;
428 
429  int energyBin = 0;
430  TH1* energyHist = 0;
431  bool notFound = true;
432  for(dhist_Map_it mapit=m_PDFmap->begin(); mapit!=m_PDFmap->end(); mapit++){
433  if( fabs(mapit->first+thetaLow)<0.000001 ) {
434  energyHist = mapit->second;
435  notFound = false;
436  break;
437  }
438  }
439  if(notFound) std::cout << "GetThetaAndEnergy: ERROR:\tInvalid theta!" << std::endl;
440  // MSG("string = " << To_TString(thetaLow) << ", double = " << thetaLow );
441 
442  energyHist->GetBinWithContent(double(rand2),energyBin,1,energyHist->GetNbinsX(),1.0);
443  if(energyHist->GetBinContent(energyBin) < rand2) energyBin += 1;
444  double drand2 = (energyHist->GetBinContent(energyBin) - rand2)/(energyHist->GetBinContent(energyBin) - energyHist->GetBinContent(energyBin-1));
445  double energyLow = energyHist->GetXaxis()->GetBinLowEdge(energyBin);
446  double energyUp = energyHist->GetXaxis()->GetBinUpEdge(energyBin);
447  double energy = drand2*(energyUp-energyLow) + energyLow;
448  // MSG("MuFlux::GetThetaEnergy()\te = " << energy*1000. );
449 
450  return std::make_pair(theta,energy);
451  }
std::map< double, TH1 * >::iterator dhist_Map_it
QTextStream & endl(QTextStream &s)
void evgen::GaisserParam::MakePDF ( )
private

Definition at line 454 of file GaisserParam_module.cc.

455  {
456  std::cout << "In my function MakePDF" << std::endl;
457  m_File=0;
458  m_PDFmap=0;
459  m_thetaHist=0;
460  double TotalMuonFlux=0;
461 
462  if(m_PDFmap){
463  std::cout << "PDFMAP" << std::endl;
464  if(m_PDFmap->size()>0 && !fSetReWrite){
465  std::cout << "MakePDF: Map has already been initialised. " << std::endl;
466  std::cout << "Do fSetReWrite - true if you really want to override this map." << std::endl;
467  return;
468  }
469  std::cout << fSetReWrite << std::endl;
470  if(fSetReWrite) ResetMap();
471  }
472  else{
473  m_PDFmap = new dhist_Map;
474  std::cout << "Making new dhist_Map called m_PDFmap....." << std::endl;
475  }
476 
477  TF2* muonSpec = new TF2("muonSpec", this,
480  "GaisserParam", "GaisserMuonFlux_Integrand"
481  );
482  //--------------------------------------------
483  //------------ Compute the pdfs
484 
485  //---- compute pdf for the theta
486  TotalMuonFlux = muonSpec->Integral(fEmin, fEmax, fThetamin, fThetamax, fEpsilon ); // Work out the muon flux at the surface
487  std::cout << "Surface flux of muons = " << TotalMuonFlux << " cm-2 s-1" << std::endl;
488 
489  //---- work out if we're reading a file, writing to file, or neither
490  std::ostringstream pdfFile;
491  pdfFile << "GaisserPDF_"<<fEmin<<"-"<<fEmid<<"-"<<fEmax<<"-"<< fEBinsLow<<"-"<<fEBinsHigh<<"-"<<fThetamin<<"-"<<fThetamax<<"-"<<fThetaBins<<".root";
492  std::string tmpfileName = pdfFile.str();
493  std::replace(tmpfileName.begin(),tmpfileName.end(),'+','0');
494  if (tmpfileName == "GaisserPDF_0-100-100100-1000-10000-0-1.5708-100.root") tmpfileName = "GaisserPDF_DefaultBins.root";
495  else if (tmpfileName == "GaisserPDF_0-100-4000-1000-1000-0-1.5708-100.root") tmpfileName = "GaisserPDF_LowEnergy.root";
496  else if (tmpfileName == "GaisserPDF_4000-10000-100000-1000-10000-0-1.5708-100.root") tmpfileName = "GaisserPDF_MidEnergy.root";
497  else if (tmpfileName == "GaisserPDF_100000-500000-1e007-10000-100000-0-1.5708-100.root") tmpfileName = "GaisserPDF_HighEnergy.root";
498 
499  std::ostringstream pdfFilePath;
500  pdfFilePath << fInputDir << tmpfileName;
501  std::string fileName = pdfFilePath.str();
502 
503  // fTemplateFile = pset.get< std::string >("TemplateFile");
504  // //fCalorimetryModuleLabel = pset.get< std::string >("CalorimetryModuleLabel");
505 
506  cet::search_path sp("FW_SEARCH_PATH");
507  std::string fROOTfile; //return /lbne/data/0-100-1.57.root
508  if( sp.find_file(tmpfileName, fROOTfile) ) fileName = fROOTfile;
509 
510  std::cout << "File path; " << fileName << std::endl;
511 
512  if(fSetRead){
513  struct stat buffer;
514  fSetRead = stat(fileName.c_str(), &buffer) == 0; // Check if file exists already
515  if(!fSetRead) std::cout << "WARNING- "+fileName+" does not exist." << std::endl;
516  else{
517  std::cout << "Reading PDF from file "+fileName << std::endl;
518  m_File = new TFile(fileName.c_str()); // Open file
519  if(m_File->IsZombie() || m_File->TestBit(TFile::kRecovered)){ // Check that file is not corrupted
520  std::cout << "WARNING- "+fileName+" is corrupted or cannot be read." << std::endl;
521  fSetRead = false;
522  }
523  }
524  }
525 
526  if(fSetRead){ // If the file exists then read it....
527  std::cout << "Now going to read file...." << std::endl;
528  fSetWrite = false; // Don't want to write as already exists
529  double thetalow = fThetamin;
530  m_thetaHist = (TH1D*) m_File->Get("pdf_theta");
531  for(int i=1; i<=fThetaBins; i++){
532  std::ostringstream pdfEnergyHist;
533  pdfEnergyHist << "pdf_energy_"<<i;
534  std::string pdfEnergyHiststr = pdfEnergyHist.str();
535 
536  TH1* pdf_hist = (TH1D*) m_File->Get( pdfEnergyHiststr.c_str() ); // Get the bin
537  m_PDFmap->insert(std::make_pair(-thetalow,pdf_hist)); //---- -ve theta for quicker sorting
538  thetalow = double(i)*(fThetamax)/double(fThetaBins); //---- increment the value of theta
539  }
540  } // ------------------------------------------
541  else { // File doesn't exist so want to make it.....
542  // ------------------------------------------
543  std::cout << "Generating a new muon flux PDF" << std::endl;
544  if(fSetWrite){
545  std::cout << "Writing to PDF to file "+fileName << std::endl;
546  m_File = new TFile(fileName.c_str(),"RECREATE");
547  }
548 
549  double dnbins_theta = double(fThetaBins);
550  m_thetaHist = new TH1D("pdf_theta", "pdf_theta", fThetaBins, fThetamin, fThetamax);
551  for(int i=1; i<=fThetaBins; i++){
552  double di = i==0 ? 0.1 : double(i);
553  double theta = di*(fThetamax)/dnbins_theta;
554  double int_i = muonSpec->Integral(fEmin, fEmax, fThetamin, theta, fEpsilon);
555  m_thetaHist->SetBinContent(i, int_i/TotalMuonFlux);
556  }
557  if(fSetWrite) m_thetaHist->Write();
558  std::cout << "theta PDF complete... now making the energy PDFs (this will take a little longer)... " << std::endl;
559 
560  //---- now compute the energy pdf
561  double thetalow = fThetamin;
562  for(int i=1; i<=fThetaBins; i++){
563 
564  double di = double(i);
565  double theta = di*(fThetamax)/fThetaBins;
566 
567  //---- compute the total integral
568  double int_tot = muonSpec->Integral(fEmin, fEmax, thetalow, theta, fEpsilon);
569 
570  //---- compute pdf for the low energy
571  int nbins = fEBinsLow;
572  TH1* pdf_lowenergy = new TH1D("pdf_lowenergy", "pdf_lowenergy", nbins, fEmin, fEmid);
573  double dnbins = double(nbins);
574  for(int j=1; j<=nbins; j++){
575  double dj = double(j);
576  double int_j = muonSpec->Integral(fEmin, fEmin + dj*(fEmid-fEmin)/dnbins, thetalow, theta, fEpsilon);
577  // std::std::cout << j << "(" << m_emin << " --> " << m_emin + dj*m_emid/dnbins << ") = " << int_j/int_tot << std::std::endl;
578  // std::std::cout << j << "(" << m_emin << " --> " << m_emin + dj*(m_emid-m_emin)/dnbins << ") = " << int_j/int_tot << std::std::endl;
579  pdf_lowenergy->SetBinContent(j, int_j/int_tot);
580  }
581 
582  //---- compute pdf for the high energy
583  nbins = fEBinsHigh;
584  dnbins=double(nbins);
585  TH1* pdf_highenergy = new TH1D("pdf_highenergy", "pdf_highenergy", nbins, fEmid, fEmax);
586  for(int j=1; j<=nbins; j++){
587  double dj = double(j);
588  double int_j = muonSpec->Integral(fEmin, fEmid + dj*(fEmax-fEmid)/dnbins, thetalow, theta, fEpsilon);
589  // std::cout << j << "(" << m_emin << " --> " << m_emid + dj*(m_emax-m_emid)/dnbins << ") = " << int_j/int_tot << std::endl;
590  pdf_highenergy->SetBinContent(j, int_j/int_tot);
591  }
592 
593  //---- now combine the two energy hists
594  std::vector<double> vxbins = GetBinning(pdf_lowenergy->GetXaxis(),false);
595  std::vector<double> vxbins2 = GetBinning(pdf_highenergy->GetXaxis());
596  vxbins.insert(vxbins.end(), vxbins2.begin(), vxbins2.end());
597 
598  int ibin = 0;
599  double* xbins = new double[vxbins.size()];
600  for(std::vector<double>::const_iterator binit=vxbins.begin(); binit!=vxbins.end(); binit++, ibin++) xbins[ibin]=(*binit);
601  TH1* pdf_energy = new TH1D("pdf_energy", "pdf_energy", vxbins.size()-1, xbins);
602  int ibin2 = 1;
603  for(ibin = 1; ibin<=pdf_lowenergy->GetNbinsX(); ibin++, ibin2++){
604  double content = pdf_lowenergy->GetBinContent(ibin);
605  if(ibin == 1) content = content - 0.00001;
606  pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), content);
607  }
608  for(ibin = 1; ibin<=pdf_highenergy->GetNbinsX(); ibin++, ibin2++){
609  pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), pdf_highenergy->GetBinContent(ibin));
610  }
611 
612  //---- and remove any negative bins
613  std::ostringstream Clonestr;
614  Clonestr << "pdf_energy_"<<i;
615  TH1* pdf_energy_noneg = (TH1D*) pdf_energy->Clone( Clonestr.str().c_str() );
616  pdf_energy_noneg->Reset();
617 
618  double PDF = 0.0;
619  double lastPD = 0.0;
620  int nSkip=0;
621  for(ibin = 1; ibin<=pdf_energy->GetNbinsX(); ibin++){
622  double newPD = pdf_energy->GetBinContent(ibin);
623  double probDiff = newPD - lastPD;
624  if(probDiff<0){
625  if(ibin!=pdf_energy->GetNbinsX()){
626  nSkip++;
627  continue;
628  }
629  else probDiff = 0;
630  }
631  else PDF += probDiff;
632  if(PDF > 1) PDF = 1;
633  for(int iskip=0; iskip <= nSkip; iskip++) pdf_energy_noneg->Fill(pdf_energy_noneg->GetBinCenter(ibin-iskip), PDF);
634  nSkip=0;
635  lastPD = newPD;
636  }
637 
638 
639  //---- add this hist, increment thetalow and delete unwanted TH1s
640  if(fSetWrite) pdf_energy_noneg->Write();
641  m_PDFmap->insert(std::make_pair(-thetalow,pdf_energy_noneg)); //---- -ve theta for quicker sorting
642 
643  //---- increment the value of theta
644  thetalow = theta;
645 
646  //---- free up memory from unwanted hists
647  delete pdf_lowenergy;
648  delete pdf_highenergy;
649  delete pdf_energy;
650 
651  std::cout << "\r===> " << 100.0*double(i)/double(fThetaBins) << "% complete... " << std::endl;
652  } // ThetaBins
653  std::cout << "finished the energy pdfs." << std::endl;
654  }//---- if(!m_doRead)
655 
656  delete muonSpec;
657  return;
658  } // Make PDF
double fEpsilon
Minimum integration sum....
double fThetamax
Maximum theta.
std::string string
Definition: nybbler.cc:12
int fThetaBins
Number of theta Bins.
double fEmid
Energy to go from low to high (GeV)
bool fSetReWrite
Whether to ReWrite pdfs.
intermediate_table::const_iterator const_iterator
std::string fInputDir
Input Directory.
int fEBinsLow
Number of low energy Bins.
static QChar PDF((ushort) 0x202c)
fileName
Definition: dumpTree.py:9
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.
QTextStream & endl(QTextStream &s)
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 259 of file GaisserParam_module.cc.

260  {
261  ///unique_ptr allows ownership to be transferred to the art::Event after the put statement
262  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
263 
264  simb::MCTruth truth;
266 
267  Sample(truth, fEngine);
268 
269  MF_LOG_DEBUG("GaisserParam") << truth;
270 
271  truthcol->push_back(truth);
272 
273  evt.put(std::move(truthcol));
274  }
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
single particles thrown at the detector
Definition: MCTruth.h:26
def move(depos, offset)
Definition: depos.py:107
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
void Sample(simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
#define MF_LOG_DEBUG(id)
Event generator information.
Definition: MCTruth.h:32
void evgen::GaisserParam::ResetMap ( )
private

Definition at line 712 of file GaisserParam_module.cc.

712  {
713  if(m_thetaHist) delete m_thetaHist;
714  if(m_PDFmap){
715  for(dhist_Map_it mapit=m_PDFmap->begin(); mapit!=m_PDFmap->end(); mapit++){
716  if(mapit->second) delete mapit->second;
717  }
718  m_PDFmap->clear();
719  delete m_PDFmap;
720  std::cout << "Reset PDFmap and thetaHist..." << std::endl;
721  }
722  } // ResetMap
std::map< double, TH1 * >::iterator dhist_Map_it
QTextStream & endl(QTextStream &s)
void evgen::GaisserParam::Sample ( simb::MCTruth mct,
CLHEP::HepRandomEngine &  engine 
)
private

Definition at line 388 of file GaisserParam_module.cc.

389  {
390  switch (fMode) {
391  case 0: // List generation mode: every event will have one of each
392  // particle species in the fPDG array
393  for (unsigned int i=0; i<fPDG.size(); ++i) {
394  SampleOne(i, mct, engine);
395  }//end loop over particles
396  break;
397  case 1: // Random selection mode: every event will exactly one particle
398  // selected randomly from the fPDG array
399  {
400  CLHEP::RandFlat flat(engine);
401 
402  unsigned int i=flat.fireInt(fPDG.size());
403  SampleOne(i, mct, engine);
404  }
405  break;
406  default:
407  mf::LogWarning("UnrecognizeOption") << "GaisserParam does not recognize ParticleSelectionMode "
408  << fMode;
409  break;
410  } // switch on fMode
411 
412  return;
413  }
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< int > fPDG
PDG code of particles to generate.
void SampleOne(unsigned int i, simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
void evgen::GaisserParam::SampleOne ( unsigned int  i,
simb::MCTruth mct,
CLHEP::HepRandomEngine &  engine 
)
private

Definition at line 279 of file GaisserParam_module.cc.

279  {
280 
281  CLHEP::RandFlat flat(engine);
282  CLHEP::RandGaussQ gauss(engine);
283 
284  Time = Momentum = KinEnergy = Gamma = Energy = Theta = Phi = pnorm = 0.0;
286 
287  TVector3 x;
288  TVector3 pmom;
289 
290  // set track id to -i as these are all primary particles and have id <= 0
291  int trackid = -1*(i+1);
292  std::string primary("primary");
293 
294  // Work out whether particle/antiparticle, and mass...
295  double m =0.0;
296  fPDG[i] = 13;
297  if (fCharge == 0 ) {
298  if(1.0-2.0*flat.fire() > 0) fPDG[i]=-fPDG[i];
299  }
300  else if ( fCharge == 1 ) fPDG[i] = 13;
301  else if ( fCharge == 2 ) fPDG[i] = -13;
302  static TDatabasePDG pdgt;
303  TParticlePDG* pdgp = pdgt.GetParticle(fPDG[i]);
304  if (pdgp) m = pdgp->Mass();
305 
306  // Work out T0...
307  if(fTDist==kGAUS){
308  Time = gauss.fire(fT0, fSigmaT);
309  }
310  else {
311  Time = fT0 + fSigmaT*(2.0*flat.fire()-1.0);
312  }
313 
314  // Work out Positioning.....
315  x[0] = (1.0 - 2.0*flat.fire())*fXHalfRange + fCenterX;
316  x[1] = fYInput;
317  x[2] = (1.0 - 2.0*flat.fire())*fZHalfRange + fCenterZ;
318 
319  // Make Lorentz vector for x and t....
320  TLorentzVector pos(x[0], x[1], x[2], Time);
321 
322  // Access the pdf map which has been loaded.....
323  if(m_PDFmap) {
324 
325  //---- get the muon theta and energy from histograms using 2 random numbers
326  std::pair<double,double> theta_energy; //---- muon theta and energy
327  theta_energy = GetThetaAndEnergy(flat.fire(),flat.fire());
328 
329  //---- Set theta, phi
330  Theta = theta_energy.first; // Angle returned by GetThetaAndEnergy between 0 and pi/2
331  Phi = M_PI*( 1.0-2.0*flat.fire() ); // Randomly generated angle between -pi and pi
332 
333  //---- Set KE, E, p
334  KinEnergy = theta_energy.second; // Energy returned by GetThetaAndEnergy
335  Gamma = 1 + (KinEnergy/m);
336  Energy = Gamma * m;
337  Momentum = std::sqrt(Energy*Energy-m*m); // Get momentum
338 
339  pmom[0] = Momentum*std::sin(Theta)*std::cos(Phi);
340  pmom[1] = -Momentum*std::cos(Theta);
341  pmom[2] = Momentum*std::sin(Theta)*std::sin(Phi);
342 
343  pnorm = std::sqrt( pmom[0]*pmom[0] + pmom[1]*pmom[1] + pmom[2]*pmom[2] );
344  DirCosineX = pmom[0] / pnorm;
345  DirCosineY = pmom[1] / pnorm;
346  DirCosineZ = pmom[2] / pnorm;
347  }
348  else {
349  std::cout << "MuFlux map hasn't been initialised, aborting...." << std::endl;
350  return;
351  }
352 
353  TLorentzVector pvec(pmom[0], pmom[1], pmom[2], Energy );
354 
355  simb::MCParticle part(trackid, fPDG[i], primary);
356  part.AddTrajectoryPoint(pos, pvec);
357  /*
358  fTree->Branch();
359  fTree->Branch();
360  fTree->Branch();
361  fPositionX ->Fill (x[0]);
362  fPositionY ->Fill (x[1]);
363  fPositionZ ->Fill (x[2]);
364  fTime ->Fill (Time);
365  fMomentumHigh ->Fill (Momentum);
366  fMomentum ->Fill (Momentum);
367  fEnergy ->Fill (Energy);
368  fDirCosineX ->Fill (DirCosineX);
369  fDirCosineY ->Fill (DirCosineY);
370  fDirCosineZ ->Fill (DirCosineZ);
371  fTheta ->Fill (Theta*180/M_PI);
372  fPhi ->Fill (Phi *180/M_PI);
373  */
374  XPosition = x[0];
375  YPosition = x[1];
376  ZPosition = x[2];
377 
378  std::cout << "Theta: " << Theta << " Phi: " << Phi << " KineticEnergy: " << KinEnergy << " Energy: " << Energy << " Momentum: " << Momentum << std::endl;
379  std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
380  std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
381  std::cout << "Normalised..." << DirCosineX << " " << DirCosineY << " " << DirCosineZ << std::endl;
382 
383  fTree->Fill();
384  mct.Add(part);
385  }
static constexpr int kGAUS
double fYInput
Max Y position.
std::string string
Definition: nybbler.cc:12
double fSigmaT
Variation in t position (ns)
double fXHalfRange
Max X position.
#define M_PI
Definition: includeROOT.h:54
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
std::vector< int > fPDG
PDG code of particles to generate.
list x
Definition: train.py:276
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)
QTextStream & endl(QTextStream &s)
double fZHalfRange
Max Z position.

Member Data Documentation

double evgen::GaisserParam::DirCosineX
private

Definition at line 155 of file GaisserParam_module.cc.

double evgen::GaisserParam::DirCosineY
private

Definition at line 155 of file GaisserParam_module.cc.

double evgen::GaisserParam::DirCosineZ
private

Definition at line 155 of file GaisserParam_module.cc.

double evgen::GaisserParam::Energy
private

Definition at line 154 of file GaisserParam_module.cc.

double evgen::GaisserParam::fCenterX = 0
private

Definition at line 149 of file GaisserParam_module.cc.

double evgen::GaisserParam::fCenterZ = 0
private

Definition at line 150 of file GaisserParam_module.cc.

int evgen::GaisserParam::fCharge
private

Charge.

Definition at line 101 of file GaisserParam_module.cc.

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

Definition at line 144 of file GaisserParam_module.cc.

int evgen::GaisserParam::fEBinsHigh
private

Number of high energy Bins.

Definition at line 108 of file GaisserParam_module.cc.

int evgen::GaisserParam::fEBinsLow
private

Number of low energy Bins.

Definition at line 107 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEmax
private

Maximum Kinetic Energy (GeV)

Definition at line 105 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEmid
private

Energy to go from low to high (GeV)

Definition at line 106 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEmin
private

Minimum Kinetic Energy (GeV)

Definition at line 104 of file GaisserParam_module.cc.

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

art-managed random-number engine

Definition at line 89 of file GaisserParam_module.cc.

double evgen::GaisserParam::fEpsilon
private

Minimum integration sum....

Definition at line 126 of file GaisserParam_module.cc.

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

Input Directory.

Definition at line 102 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 93 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 96 of file GaisserParam_module.cc.

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

PDG code of particles to generate.

Definition at line 100 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetParam
private

Which version of Gaissers Param.

Definition at line 122 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetRead
private

Whether to Read.

Definition at line 123 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetReWrite
private

Whether to ReWrite pdfs.

Definition at line 125 of file GaisserParam_module.cc.

bool evgen::GaisserParam::fSetWrite
private

Whether to Write.

Definition at line 124 of file GaisserParam_module.cc.

double evgen::GaisserParam::fSigmaT
private

Variation in t position (ns)

Definition at line 119 of file GaisserParam_module.cc.

double evgen::GaisserParam::fT0
private

Central t position (ns) in world coordinates.

Definition at line 118 of file GaisserParam_module.cc.

int evgen::GaisserParam::fTDist
private

How to distribute t (gaus, or uniform)

Definition at line 120 of file GaisserParam_module.cc.

int evgen::GaisserParam::fThetaBins
private

Number of theta Bins.

Definition at line 112 of file GaisserParam_module.cc.

double evgen::GaisserParam::fThetamax
private

Maximum theta.

Definition at line 111 of file GaisserParam_module.cc.

double evgen::GaisserParam::fThetamin
private

Minimum theta.

Definition at line 110 of file GaisserParam_module.cc.

TTree* evgen::GaisserParam::fTree
private

Definition at line 153 of file GaisserParam_module.cc.

double evgen::GaisserParam::fXHalfRange
private

Max X position.

Definition at line 114 of file GaisserParam_module.cc.

double evgen::GaisserParam::fYInput
private

Max Y position.

Definition at line 115 of file GaisserParam_module.cc.

double evgen::GaisserParam::fZHalfRange
private

Max Z position.

Definition at line 116 of file GaisserParam_module.cc.

double evgen::GaisserParam::Gamma
private

Definition at line 154 of file GaisserParam_module.cc.

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

Definition at line 91 of file GaisserParam_module.cc.

double evgen::GaisserParam::KinEnergy
private

Definition at line 154 of file GaisserParam_module.cc.

TFile* evgen::GaisserParam::m_File
private

Definition at line 75 of file GaisserParam_module.cc.

dhist_Map* evgen::GaisserParam::m_PDFmap
private

Definition at line 76 of file GaisserParam_module.cc.

TH1* evgen::GaisserParam::m_thetaHist
private

Definition at line 77 of file GaisserParam_module.cc.

double evgen::GaisserParam::Momentum
private

Definition at line 154 of file GaisserParam_module.cc.

double evgen::GaisserParam::Phi
private

Definition at line 154 of file GaisserParam_module.cc.

double evgen::GaisserParam::pnorm
private

Definition at line 154 of file GaisserParam_module.cc.

double evgen::GaisserParam::Theta
private

Definition at line 154 of file GaisserParam_module.cc.

double evgen::GaisserParam::Time
private

Definition at line 154 of file GaisserParam_module.cc.

double evgen::GaisserParam::xNeg = 0
private

Definition at line 145 of file GaisserParam_module.cc.

double evgen::GaisserParam::xPos = 0
private

Definition at line 146 of file GaisserParam_module.cc.

double evgen::GaisserParam::XPosition
private

Definition at line 155 of file GaisserParam_module.cc.

double evgen::GaisserParam::YPosition
private

Definition at line 155 of file GaisserParam_module.cc.

double evgen::GaisserParam::zNeg = 0
private

Definition at line 147 of file GaisserParam_module.cc.

double evgen::GaisserParam::zPos = 0
private

Definition at line 148 of file GaisserParam_module.cc.

double evgen::GaisserParam::ZPosition
private

Definition at line 155 of file GaisserParam_module.cc.


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