evgen::GaisserParam Class Reference

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

Inheritance diagram for evgen::GaisserParam:

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

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

Private Member Functions
void	produce (art::Event &evt) override
void	beginJob () override
void	beginRun (art::Run &run) override
void	SampleOne (unsigned int i, simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
void	Sample (simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
std::pair< double, double >	GetThetaAndEnergy (double rand1, double rand2)
void	MakePDF ()
void	ResetMap ()
double	GaisserMuonFlux_Integrand (Double_t *x, Double_t *par)
double	GaisserFlux (double e, double theta)
std::vector< double >	GetBinning (const TAxis *axis, bool finalEdge=true)

Private Attributes
TFile *	m_File
dhist_Map *	m_PDFmap
TH1 *	m_thetaHist
CLHEP::HepRandomEngine &	fEngine
	art-managed random-number engine More...
int	fMode
bool	fPadOutVectors
std::vector< int >	fPDG
	PDG code of particles to generate. More...
int	fCharge
	Charge. More...
std::string	fInputDir
	Input Directory. More...
double	fEmin
	Minimum Kinetic Energy (GeV) More...
double	fEmax
	Maximum Kinetic Energy (GeV) More...
double	fEmid
	Energy to go from low to high (GeV) More...
int	fEBinsLow
	Number of low energy Bins. More...
int	fEBinsHigh
	Number of high energy Bins. More...
double	fThetamin
	Minimum theta. More...
double	fThetamax
	Maximum theta. More...
int	fThetaBins
	Number of theta Bins. More...
double	fXHalfRange
	Max X position. More...
double	fYInput
	Max Y position. More...
double	fZHalfRange
	Max Z position. More...
double	fT0
	Central t position (ns) in world coordinates. More...
double	fSigmaT
	Variation in t position (ns) More...
int	fTDist
	How to distribute t (gaus, or uniform) More...
bool	fSetParam
	Which version of Gaissers Param. More...
bool	fSetRead
	Whether to Read. More...
bool	fSetWrite
	Whether to Write. More...
bool	fSetReWrite
	Whether to ReWrite pdfs. More...
double	fEpsilon
	Minimum integration sum.... More...
double	fCryoBoundaries [6]
double	xNeg = 0
double	xPos = 0
double	zNeg = 0
double	zPos = 0
double	fCenterX = 0
double	fCenterZ = 0
TTree *	fTree
double	Time
double	Momentum
double	KinEnergy
double	Gamma
double	Energy
double	Theta
double	Phi
double	pnorm
double	XPosition
double	YPosition
double	ZPosition
double	DirCosineX
double	DirCosineY
double	DirCosineZ

Static Private Attributes
static constexpr int	kGAUS = 1

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

Detailed Description

module to 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.

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

Member Function Documentation

void evgen::GaisserParam::beginJob ( )

overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 197 of file GaisserParam_module.cc.

  {
    //Work out center of cryostat(s)
    art::ServiceHandle<geo::Geometry const> geom;
    for (unsigned int i=0; i < geom->Ncryostats() ; i++ ) {
      geom->CryostatBoundaries(fCryoBoundaries, i);
      if ( xNeg > fCryoBoundaries[0] ) xNeg = fCryoBoundaries[0];
      if ( xPos < fCryoBoundaries[1] ) xPos = fCryoBoundaries[1];
      if ( zNeg > fCryoBoundaries[4] ) zNeg = fCryoBoundaries[4];
      if ( zPos < fCryoBoundaries[5] ) zPos = fCryoBoundaries[5];
      }
    fCenterX = xNeg + (xPos-xNeg)/2;
    fCenterZ = zNeg + (zPos-zNeg)/2;

    // Make the Histograms....
    art::ServiceHandle<art::TFileService const> tfs;
    /*
    fPositionX    = tfs->make<TH1D>("fPositionX"   ,"Position (cm)" ,500,fCenterX-(fXHalfRange+10) ,fCenterX+(fXHalfRange+10));
    fPositionY    = tfs->make<TH1D>("fPositionY"   ,"Position (cm)" ,500,-(fYInput+10),(fYInput+10));
    fPositionZ    = tfs->make<TH1D>("fPositionZ"   ,"Position (cm)" ,500,fCenterZ-(fZHalfRange+10) ,fCenterZ+(fZHalfRange+10));
    fTime         = tfs->make<TH1D>("fTime"        ,"Time (s)"      ,500,0,1e6);
    fMomentumHigh = tfs->make<TH1D>("fMomentumHigh","Momentum (GeV)",500,0,fEmax);
    fMomentum     = tfs->make<TH1D>("fMomentum"    ,"Momentum (GeV)",500,0,100);
    fEnergy       = tfs->make<TH1D>("fEnergy"      ,"Energy (GeV)"  ,500,0,fEmax);

    fDirCosineX = tfs->make<TH1D>("fDirCosineX","Normalised Direction cosine",500,-1,1);
    fDirCosineY = tfs->make<TH1D>("fDirCosineY","Normalised Direction cosine",500,-1,1);
    fDirCosineZ = tfs->make<TH1D>("fDirCosineZ","Normalised Direction cosine",500,-1,1);

    fTheta      = tfs->make<TH1D>("fTheta"     ,"Angle (radians)",500,-365,365);
    fPhi        = tfs->make<TH1D>("fPhi"       ,"Angle (radians)",500,-365,365);
    */
    fTree = tfs->make<TTree>("Generator","analysis tree");
    fTree->Branch("XPosition" ,&XPosition ,"XPosition/D" );
    fTree->Branch("YPosition" ,&YPosition ,"YPosition/D" );
    fTree->Branch("ZPosition" ,&ZPosition ,"ZPosition/D" );
    fTree->Branch("Time"      ,&Time      ,"Time/D"      );
    fTree->Branch("Momentum"  ,&Momentum  ,"Momentum/D"  );
    fTree->Branch("KinEnergy" ,&KinEnergy ,"KinEnergy/D" );
    fTree->Branch("Energy"    ,&Energy    ,"Energy/D"    );
    fTree->Branch("DirCosineX",&DirCosineX,"DirCosineX/D");
    fTree->Branch("DirCosineY",&DirCosineY,"DirCosineY/D");
    fTree->Branch("DirCosineZ",&DirCosineZ,"DirCosineZ/D");
    fTree->Branch("Theta"     ,&Theta     ,"Theta/D"     );
    fTree->Branch("Phi"       ,&Phi       ,"Phi/D"       );
  }

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

overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 245 of file GaisserParam_module.cc.

  {
    // Check fcl parameters were set correctly
    if ( fThetamax > M_PI/2 + 0.01 ) throw cet::exception("GaisserParam")<< "\nThetamax has to be less than " << M_PI/2 << ", but was entered as " << fThetamax << ", this cause an error so leaving program now...\n\n";
    if ( fThetamin < 0         )     throw cet::exception("GaisserParam")<< "\nThetamin has to be more than 0, but was entered as " << fThetamin << ", this cause an error so leaving program now...\n\n" << std::endl;
    if ( fThetamax < fThetamin )     throw cet::exception("GaisserParam")<< "\nMinimum angle is bigger than maximum angle....causes an error so leaving program now....\n\n" << std::endl;

    // grab the geometry object to see what geometry we are using
    art::ServiceHandle<geo::Geometry const> geo;
    run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
    MakePDF ();
  }

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

private

Definition at line 661 of file GaisserParam_module.cc.

                                                        {

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

    return di;
  } // GaisserFlux

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

private

Definition at line 693 of file GaisserParam_module.cc.

                                                                      {

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

    return flux;
  } // MuonFluxIntegrand

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

private

Definition at line 702 of file GaisserParam_module.cc.

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

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

private

Definition at line 416 of file GaisserParam_module.cc.

  {
    if(rand1 < 0 || rand1 > 1) std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand1 << std::endl;
    if(rand2 < 0 || rand2 > 1) std::cerr << "GetThetaAndEnergy:\tInvalid random number " << rand2 << std::endl;

    int thetaBin = 0;
    m_thetaHist->GetBinWithContent(double(rand1),thetaBin,1,m_thetaHist->GetNbinsX(),1.0);
    if(m_thetaHist->GetBinContent(thetaBin) < rand1) thetaBin += 1;
    double drand1 = (m_thetaHist->GetBinContent(thetaBin) - rand1)/(m_thetaHist->GetBinContent(thetaBin) - m_thetaHist->GetBinContent(thetaBin-1));
    double thetaLow = m_thetaHist->GetXaxis()->GetBinLowEdge(thetaBin);
    double thetaUp  = m_thetaHist->GetXaxis()->GetBinUpEdge(thetaBin);
    double theta = drand1*(thetaUp-thetaLow) + thetaLow;

    int energyBin = 0;
    TH1* energyHist = 0;
    bool notFound = true;
    for(dhist_Map_it mapit=m_PDFmap->begin(); mapit!=m_PDFmap->end(); mapit++){
      if( fabs(mapit->first+thetaLow)<0.000001 ) {
        energyHist = mapit->second;
        notFound = false;
        break;
      }
    }
    if(notFound) std::cout << "GetThetaAndEnergy: ERROR:\tInvalid theta!" << std::endl;
    //  MSG("string = " << To_TString(thetaLow) << ", double = " << thetaLow );

    energyHist->GetBinWithContent(double(rand2),energyBin,1,energyHist->GetNbinsX(),1.0);
    if(energyHist->GetBinContent(energyBin) < rand2) energyBin += 1;
    double drand2 = (energyHist->GetBinContent(energyBin) - rand2)/(energyHist->GetBinContent(energyBin) - energyHist->GetBinContent(energyBin-1));
    double energyLow = energyHist->GetXaxis()->GetBinLowEdge(energyBin);
    double energyUp  = energyHist->GetXaxis()->GetBinUpEdge(energyBin);
    double energy = drand2*(energyUp-energyLow) + energyLow;
    //  MSG("MuFlux::GetThetaEnergy()\te = " << energy*1000. );

    return std::make_pair(theta,energy);
  }

void evgen::GaisserParam::MakePDF ( )

private

Definition at line 454 of file GaisserParam_module.cc.

  {
    std::cout << "In my function MakePDF" << std::endl;
    m_File=0;
    m_PDFmap=0;
    m_thetaHist=0;
    double TotalMuonFlux=0;

    if(m_PDFmap){
      std::cout << "PDFMAP" << std::endl;
      if(m_PDFmap->size()>0 && !fSetReWrite){
        std::cout << "MakePDF: Map has already been initialised. " << std::endl;
        std::cout << "Do fSetReWrite - true if you really want to override this map." << std::endl;
        return;
        }
      std::cout << fSetReWrite << std::endl;
      if(fSetReWrite) ResetMap();
    }
    else{
      m_PDFmap = new dhist_Map;
      std::cout << "Making new dhist_Map called m_PDFmap....." << std::endl;
    }

    TF2* muonSpec = new TF2("muonSpec", this,
                            &evgen::GaisserParam::GaisserMuonFlux_Integrand,
                            fEmin, fEmax, fThetamin, fThetamax, 0,
                            "GaisserParam", "GaisserMuonFlux_Integrand"
                            );
    //--------------------------------------------
    //------------ Compute the pdfs

    //---- compute pdf for the theta
    TotalMuonFlux = muonSpec->Integral(fEmin, fEmax, fThetamin, fThetamax, fEpsilon ); // Work out the muon flux at the surface
    std::cout << "Surface flux of muons = " << TotalMuonFlux << " cm-2 s-1" << std::endl;

    //---- work out if we're reading a file, writing to file, or neither
    std::ostringstream pdfFile;
    pdfFile << "GaisserPDF_"<<fEmin<<"-"<<fEmid<<"-"<<fEmax<<"-"<< fEBinsLow<<"-"<<fEBinsHigh<<"-"<<fThetamin<<"-"<<fThetamax<<"-"<<fThetaBins<<".root";
    std::string tmpfileName = pdfFile.str();
    std::replace(tmpfileName.begin(),tmpfileName.end(),'+','0');
    if      (tmpfileName == "GaisserPDF_0-100-100100-1000-10000-0-1.5708-100.root")          tmpfileName = "GaisserPDF_DefaultBins.root";
    else if (tmpfileName == "GaisserPDF_0-100-4000-1000-1000-0-1.5708-100.root")             tmpfileName = "GaisserPDF_LowEnergy.root";
    else if (tmpfileName == "GaisserPDF_4000-10000-100000-1000-10000-0-1.5708-100.root")     tmpfileName = "GaisserPDF_MidEnergy.root";
    else if (tmpfileName == "GaisserPDF_100000-500000-1e007-10000-100000-0-1.5708-100.root") tmpfileName = "GaisserPDF_HighEnergy.root";

    std::ostringstream pdfFilePath;
    pdfFilePath << fInputDir << tmpfileName;
    std::string fileName = pdfFilePath.str();

  // fTemplateFile           = pset.get< std::string >("TemplateFile");
  // //fCalorimetryModuleLabel = pset.get< std::string >("CalorimetryModuleLabel");

   cet::search_path sp("FW_SEARCH_PATH");
   std::string fROOTfile; //return /lbne/data/0-100-1.57.root
   if( sp.find_file(tmpfileName, fROOTfile) ) fileName = fROOTfile;

   std::cout << "File path; " << fileName << std::endl;

    if(fSetRead){
      struct stat buffer;
      fSetRead = stat(fileName.c_str(), &buffer) == 0; // Check if file exists already
      if(!fSetRead) std::cout << "WARNING- "+fileName+" does not exist." << std::endl;
      else{
        std::cout << "Reading PDF from file "+fileName << std::endl;
        m_File = new TFile(fileName.c_str()); // Open file
        if(m_File->IsZombie() || m_File->TestBit(TFile::kRecovered)){ // Check that file is not corrupted
          std::cout << "WARNING- "+fileName+" is corrupted or cannot be read." << std::endl;
          fSetRead = false;
        }
      }
    }

    if(fSetRead){ // If the file exists then read it....
      std::cout << "Now going to read file...." << std::endl;
      fSetWrite = false; // Don't want to write as already exists
      double thetalow  = fThetamin;
      m_thetaHist = (TH1D*) m_File->Get("pdf_theta");
      for(int i=1; i<=fThetaBins; i++){
        std::ostringstream pdfEnergyHist;
        pdfEnergyHist << "pdf_energy_"<<i;
        std::string pdfEnergyHiststr = pdfEnergyHist.str();

        TH1* pdf_hist = (TH1D*) m_File->Get( pdfEnergyHiststr.c_str() ); // Get the bin
        m_PDFmap->insert(std::make_pair(-thetalow,pdf_hist));         //---- -ve theta for quicker sorting
        thetalow = double(i)*(fThetamax)/double(fThetaBins); //---- increment the value of theta
      }
    }      // ------------------------------------------
    else { // File doesn't exist so want to make it.....
           // ------------------------------------------
      std::cout << "Generating a new muon flux PDF" << std::endl;
      if(fSetWrite){
        std::cout << "Writing to PDF to file "+fileName << std::endl;
        m_File = new TFile(fileName.c_str(),"RECREATE");
      }

      double dnbins_theta = double(fThetaBins);
      m_thetaHist = new TH1D("pdf_theta", "pdf_theta", fThetaBins, fThetamin, fThetamax);
      for(int i=1; i<=fThetaBins; i++){
        double di = i==0 ? 0.1 : double(i);
        double theta = di*(fThetamax)/dnbins_theta;
        double int_i = muonSpec->Integral(fEmin, fEmax, fThetamin, theta, fEpsilon);
        m_thetaHist->SetBinContent(i, int_i/TotalMuonFlux);
      }
      if(fSetWrite) m_thetaHist->Write();
      std::cout << "theta PDF complete... now making the energy PDFs (this will take a little longer)... " << std::endl;

      //---- now compute the energy pdf
      double thetalow  = fThetamin;
      for(int i=1; i<=fThetaBins; i++){

        double di = double(i);
        double theta = di*(fThetamax)/fThetaBins;

        //---- compute the total integral
        double int_tot = muonSpec->Integral(fEmin, fEmax, thetalow, theta, fEpsilon);

        //---- compute pdf for the low energy
        int nbins = fEBinsLow;
        TH1* pdf_lowenergy = new TH1D("pdf_lowenergy", "pdf_lowenergy", nbins, fEmin, fEmid);
        double dnbins = double(nbins);
        for(int j=1; j<=nbins; j++){
          double dj = double(j);
          double int_j = muonSpec->Integral(fEmin, fEmin + dj*(fEmid-fEmin)/dnbins, thetalow, theta, fEpsilon);
          //    std::std::cout << j << "(" << m_emin << " --> " << m_emin + dj*m_emid/dnbins << ") = " << int_j/int_tot << std::std::endl;
          //    std::std::cout << j << "(" << m_emin << " --> " << m_emin + dj*(m_emid-m_emin)/dnbins << ") = " << int_j/int_tot << std::std::endl;
          pdf_lowenergy->SetBinContent(j, int_j/int_tot);
        }

        //---- compute pdf for the high energy
        nbins = fEBinsHigh;
        dnbins=double(nbins);
        TH1* pdf_highenergy = new TH1D("pdf_highenergy", "pdf_highenergy", nbins, fEmid, fEmax);
        for(int j=1; j<=nbins; j++){
          double dj = double(j);
          double int_j = muonSpec->Integral(fEmin, fEmid + dj*(fEmax-fEmid)/dnbins, thetalow, theta, fEpsilon);
          //      std::cout << j << "(" << m_emin << " --> " << m_emid + dj*(m_emax-m_emid)/dnbins << ") = " << int_j/int_tot << std::endl;
          pdf_highenergy->SetBinContent(j, int_j/int_tot);
        }

        //---- now combine the two energy hists
        std::vector<double> vxbins  = GetBinning(pdf_lowenergy->GetXaxis(),false);
        std::vector<double> vxbins2 = GetBinning(pdf_highenergy->GetXaxis());
        vxbins.insert(vxbins.end(), vxbins2.begin(), vxbins2.end());

        int ibin = 0;
        double* xbins = new double[vxbins.size()];
        for(std::vector<double>::const_iterator binit=vxbins.begin(); binit!=vxbins.end(); binit++, ibin++) xbins[ibin]=(*binit);
        TH1* pdf_energy = new TH1D("pdf_energy", "pdf_energy", vxbins.size()-1, xbins);
        int ibin2 = 1;
        for(ibin = 1; ibin<=pdf_lowenergy->GetNbinsX();   ibin++, ibin2++){
          double content = pdf_lowenergy->GetBinContent(ibin);
          if(ibin == 1) content = content - 0.00001;
          pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), content);
        }
        for(ibin = 1; ibin<=pdf_highenergy->GetNbinsX();  ibin++, ibin2++){
          pdf_energy->Fill(pdf_energy->GetBinCenter(ibin2), pdf_highenergy->GetBinContent(ibin));
        }

        //---- and remove any negative bins
        std::ostringstream Clonestr;
        Clonestr << "pdf_energy_"<<i;
        TH1* pdf_energy_noneg = (TH1D*) pdf_energy->Clone( Clonestr.str().c_str() );
        pdf_energy_noneg->Reset();

        double PDF = 0.0;
        double lastPD = 0.0;
        int nSkip=0;
        for(ibin = 1; ibin<=pdf_energy->GetNbinsX(); ibin++){
          double newPD    = pdf_energy->GetBinContent(ibin);
          double probDiff = newPD - lastPD;
          if(probDiff<0){
            if(ibin!=pdf_energy->GetNbinsX()){
              nSkip++;
              continue;
            }
            else probDiff = 0;
          }
          else PDF += probDiff;
          if(PDF > 1) PDF = 1;
          for(int iskip=0; iskip <= nSkip; iskip++) pdf_energy_noneg->Fill(pdf_energy_noneg->GetBinCenter(ibin-iskip), PDF);
          nSkip=0;
          lastPD = newPD;
        }


        //---- add this hist, increment thetalow and delete unwanted TH1s
        if(fSetWrite) pdf_energy_noneg->Write();
        m_PDFmap->insert(std::make_pair(-thetalow,pdf_energy_noneg)); //---- -ve theta for quicker sorting

        //---- increment the value of theta
        thetalow = theta;

        //---- free up memory from unwanted hists
        delete pdf_lowenergy;
        delete pdf_highenergy;
        delete pdf_energy;

        std::cout << "\r===> " << 100.0*double(i)/double(fThetaBins) << "% complete... " << std::endl;
      } // ThetaBins
      std::cout << "finished the energy pdfs." << std::endl;
    }//---- if(!m_doRead)

    delete muonSpec;
    return;
  } // Make PDF

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.

  {
    ///unique_ptr allows ownership to be transferred to the art::Event after the put statement
    std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);

    simb::MCTruth truth;
    truth.SetOrigin(simb::kSingleParticle);

    Sample(truth, fEngine);

    MF_LOG_DEBUG("GaisserParam") << truth;

    truthcol->push_back(truth);

    evt.put(std::move(truthcol));
  }

void evgen::GaisserParam::ResetMap ( )

private

Definition at line 712 of file GaisserParam_module.cc.

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

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

private

Definition at line 388 of file GaisserParam_module.cc.

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

        unsigned int i=flat.fireInt(fPDG.size());
        SampleOne(i, mct, engine);
      }
      break;
    default:
      mf::LogWarning("UnrecognizeOption") << "GaisserParam does not recognize ParticleSelectionMode "
                                          << fMode;
      break;
    } // switch on fMode

    return;
  }

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

private

Definition at line 279 of file GaisserParam_module.cc.

                                                                                            {

    CLHEP::RandFlat   flat(engine);
    CLHEP::RandGaussQ gauss(engine);

    Time = Momentum = KinEnergy = Gamma = Energy = Theta = Phi = pnorm = 0.0;
    XPosition = YPosition = ZPosition = DirCosineX = DirCosineY = DirCosineZ = 0.0;

    TVector3 x;
    TVector3 pmom;

    // set track id to -i as these are all primary particles and have id <= 0
    int trackid = -1*(i+1);
    std::string primary("primary");

    // Work out whether particle/antiparticle, and mass...
    double m =0.0;
    fPDG[i] = 13;
    if (fCharge == 0 ) {
      if(1.0-2.0*flat.fire() > 0) fPDG[i]=-fPDG[i];
    }
    else if ( fCharge == 1 ) fPDG[i] = 13;
    else if ( fCharge == 2 ) fPDG[i] = -13;
    static TDatabasePDG  pdgt;
    TParticlePDG* pdgp = pdgt.GetParticle(fPDG[i]);
    if (pdgp) m = pdgp->Mass();

    // Work out T0...
    if(fTDist==kGAUS){
      Time = gauss.fire(fT0, fSigmaT);
    }
    else {
      Time = fT0 + fSigmaT*(2.0*flat.fire()-1.0);
    }

    // Work out Positioning.....
    x[0] = (1.0 - 2.0*flat.fire())*fXHalfRange + fCenterX;
    x[1] = fYInput;
    x[2] = (1.0 - 2.0*flat.fire())*fZHalfRange + fCenterZ;

    // Make Lorentz vector for x and t....
    TLorentzVector pos(x[0], x[1], x[2], Time);

    // Access the pdf map which has been loaded.....
    if(m_PDFmap) {

      //---- get the muon theta and energy from histograms using 2 random numbers
      std::pair<double,double> theta_energy; //---- muon theta and energy
      theta_energy = GetThetaAndEnergy(flat.fire(),flat.fire());

      //---- Set theta, phi
      Theta = theta_energy.first;  // Angle returned by GetThetaAndEnergy between 0 and pi/2
      Phi   = M_PI*( 1.0-2.0*flat.fire() ); // Randomly generated angle between -pi and pi

      //---- Set KE, E, p
      KinEnergy = theta_energy.second; // Energy returned by GetThetaAndEnergy
      Gamma         = 1 + (KinEnergy/m);
      Energy        = Gamma * m;
      Momentum      = std::sqrt(Energy*Energy-m*m); // Get momentum

      pmom[0] = Momentum*std::sin(Theta)*std::cos(Phi);
      pmom[1] = -Momentum*std::cos(Theta);
      pmom[2] = Momentum*std::sin(Theta)*std::sin(Phi);

      pnorm   = std::sqrt( pmom[0]*pmom[0] + pmom[1]*pmom[1] + pmom[2]*pmom[2] );
      DirCosineX = pmom[0] / pnorm;
      DirCosineY = pmom[1] / pnorm;
      DirCosineZ = pmom[2] / pnorm;
    }
    else {
      std::cout << "MuFlux map hasn't been initialised, aborting...." << std::endl;
      return;
    }

    TLorentzVector pvec(pmom[0], pmom[1], pmom[2], Energy );

    simb::MCParticle part(trackid, fPDG[i], primary);
    part.AddTrajectoryPoint(pos, pvec);
    /*
    fTree->Branch();
    fTree->Branch();
    fTree->Branch();
    fPositionX    ->Fill (x[0]);
    fPositionY    ->Fill (x[1]);
    fPositionZ    ->Fill (x[2]);
    fTime         ->Fill (Time);
    fMomentumHigh ->Fill (Momentum);
    fMomentum     ->Fill (Momentum);
    fEnergy       ->Fill (Energy);
    fDirCosineX   ->Fill (DirCosineX);
    fDirCosineY   ->Fill (DirCosineY);
    fDirCosineZ   ->Fill (DirCosineZ);
    fTheta        ->Fill (Theta*180/M_PI);
    fPhi          ->Fill (Phi  *180/M_PI);
    */
    XPosition = x[0];
    YPosition = x[1];
    ZPosition = x[2];

    std::cout << "Theta: " << Theta << " Phi: " << Phi << " KineticEnergy: " << KinEnergy << " Energy: " << Energy << " Momentum: " << Momentum << std::endl;
    std::cout << "x: " <<  pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
    std::cout << "Px: " <<  pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
    std::cout << "Normalised..." << DirCosineX << " " << DirCosineY << " " << DirCosineZ << std::endl;

    fTree->Fill();
    mct.Add(part);
  }

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:

• larsim/larsim/EventGenerator/MuonPropagation/GaisserParam_module.cc