gUpMuFluxGen.cxx File Reference

#include <cassert>
#include <cstdlib>
#include <cctype>
#include <string>
#include <vector>
#include <sstream>
#include <map>
#include <TRotation.h>
#include <TH1D.h>
#include <TH3D.h>
#include <TTree.h>
#include <TLorentzVector.h>
#include "Framework/Algorithm/Algorithm.h"
#include "Framework/Algorithm/AlgFactory.h"
#include "Framework/EventGen/XSecAlgorithmI.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/GFluxI.h"
#include "Framework/EventGen/GMCJDriver.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/UnitUtils.h"
#include "Framework/Utils/KineUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"

Go to the source code of this file.

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	PrintSyntax (void)
GFluxI *	GetFlux (void)
void	GenerateUpNu (GFluxI *flux_driver)
TH3D *	BuildEmuEnuCosThetaPdf (int nu_code)
TH1D *	GetEmuPdf (double Enu, double costheta, const TH3D *pdf3d)
TVector3	GetDetectorVertex (double CosTheta, double Enu)
double	GetCrossSection (int nu_code, double Enu, double Emu)
double	ProbabilityEmu (int nu_code, double Enu, double Emu)
int	main (int argc, char **argv)

Variables
Long_t	gOptRunNu
string	gOptFluxSim
map< int, string >	gOptFluxFiles
int	gOptNev = -1
double	gOptDetectorSide
long int	gOptRanSeed
string	gOptInpXSecFile
const double	a = 2e+6
const double	e = 500e+9
double	kDefOptDetectorSide = 1e+5

Function Documentation

TH3D * BuildEmuEnuCosThetaPdf

(

int

nu_code

)

Definition at line 444 of file gUpMuFluxGen.cxx.

{
// Set up a 3D histogram, with axes Emu, Enu, CosTheta.
// Bin convention is defined at the start.
// Content of each bin is given by the probability of getting a muon
// of energy Emu from a neutrino of energy Enu and zenith ange costheta

  // Bin convention for Enu, consistent with BGLRS
  const double Enumin            = 0.1;
  const int    nEnubinsPerDecade = 10;
  const int    nEnubins          = 31;

  // Bin convention for Emu, consistent with BGLRS
  const double Emumin            = 0.1;
  const int    nEmubinsPerDecade = 10;
  const int    nEmubins          = 31;

  // Bin convention for CosTheta, consistent with BGLRS
  const double costheta_min = -1;
  const double costheta_max = +1;
  const int    ncostheta    = 20;

  // Set up an array of CosTheta bins.
  double costhetabinwidth = ((costheta_max-costheta_min)/ncostheta);
  double CosThetaBinEdges[ncostheta+1];
  for (int i=0; i<=ncostheta; i++)
    {
      CosThetaBinEdges[i] = costheta_min + i*costhetabinwidth;
    }

  // Set up an array of Enu bins.
  Double_t MinLogEnu = log(Enumin);
  Double_t MaxLogEnu = log(10*Enumin);
  Double_t LogBinWidthEnu = ((MaxLogEnu-MinLogEnu)/nEnubinsPerDecade);
  Double_t EnuBinEdges[nEnubins+1];
  for (int i=0; i<=nEnubins; i++)
    {
      EnuBinEdges[i] = exp(MinLogEnu + i*LogBinWidthEnu);
    }

  // Set up an array of Emu bins.
  Double_t MinLogEmu = log(Emumin);
  Double_t MaxLogEmu = log(10*Emumin);
  Double_t LogBinWidthEmu = ((MaxLogEmu-MinLogEmu)/nEmubinsPerDecade);
  Double_t EmuBinEdges[nEmubins+2];
  for (int i=0; i<=nEmubins+1; i++)
    {
      EmuBinEdges[i] = exp(MinLogEmu + i*LogBinWidthEmu);
    }

  // Create 3D histogram. X-axis: Emu; Y-axis: Enu; Z-axis: CosTheta.
  TH3D *h3 = new TH3D("h3","",nEmubins+1,EmuBinEdges,nEnubins,EnuBinEdges,ncostheta,CosThetaBinEdges);

  // Draw histogram.
  //h3->Draw();

  // Calculate Probability for each triplet.
  for (int i=1; i< h3->GetNbinsX(); i++)
    {
      double Emu = h3->GetXaxis()->GetBinCenter(i);
      for (int j=1; j<= h3->GetNbinsY(); j++)
        {
          double Enu = h3->GetBinCenter(j);
          double Int = ProbabilityEmu(nu_code,Enu,Emu);
          for (int k=1; k<= h3->GetNbinsZ(); k++)
            {
              h3->SetBinContent(i,j,k,Int); // Bin Content will is Int.
              //h3->SetBinContent(i,j,k,1); // Bin content constant (to test)
            }
        }
    }
  return h3;
}

void GenerateUpNu

(

GFluxI *

flux_driver

)

Definition at line 416 of file gUpMuFluxGen.cxx.

{
// Generate a Neutrino. Keep generating neutrinos until an upgoing one is generated.

  while (1)
  {
    LOG("gevgen_upmu", pINFO) << "Pulling next neurtino from the flux driver...";
    flux_driver->GenerateNext();

    int nu_code = flux_driver->PdgCode();

    const TLorentzVector & p4 = flux_driver->Momentum();
    double costheta = -p4.Pz() / p4.Vect().Mag();

    bool keep = (costheta<0) && (nu_code==kPdgNuMu || nu_code==kPdgAntiNuMu);
    if(keep) return;
  }
}

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 643 of file gUpMuFluxGen.cxx.

{
// Get the command line arguments

  LOG("gevgen_upmu", pINFO) << "Parsing command line arguments";

  // Common run options. Set defaults and read.
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);

  // Parse run options for this app

  CmdLnArgParser parser(argc,argv);

  // help?
  bool help = parser.OptionExists('h');
  if(help) {
    PrintSyntax();
    exit(0);
  }

  //
  // *** run number
  //
  if( parser.OptionExists('r') ) {
    LOG("gevgen_upmu", pDEBUG) << "Reading MC run number";
    gOptRunNu = parser.ArgAsLong('r');
  } else {
    LOG("gevgen_upmu", pDEBUG) << "Unspecified run number - Using default";
    gOptRunNu = 1000;
  } //-r

  //
  // *** exposure
  //

  // in number of events
  bool have_required_statistics = false;
  if( parser.OptionExists('n') ) {
    LOG("gevgen_upmu", pDEBUG)
        << "Reading number of events to generate";
    gOptNev = parser.ArgAsInt('n');
    have_required_statistics = true;
  }//-n
  if(!have_required_statistics) {
    LOG("gevgen_upmu", pFATAL)
       << "You must request exposure either in terms of number of events"
       << "\nUse the -n option";
    PrintSyntax();
    gAbortingInErr = true;
    exit(1);
  }

  //
  // *** detector side length
  //

  if( parser.OptionExists('d') ) {
    LOG("gevgen_upmu", pDEBUG)
      << "Reading detector side length";
    gOptDetectorSide = parser.ArgAsDouble('d');
  } else {
    LOG("gevgen_upmu", pDEBUG)
      << "Unspecified detector side length - Using default";
    gOptDetectorSide = kDefOptDetectorSide;
  }//-d

  //
  // *** flux files
  //

  // syntax:
  // simulation:/path/file.data[neutrino_code],/path/file.data[neutrino_code],...
  //
  if( parser.OptionExists('f') ) {
    LOG("gevgen_upmu", pDEBUG) << "Getting input flux files";
    string flux = parser.ArgAsString('f');

    // get flux simulation info (FLUKA,BGLRS) so as to instantiate the
    // appropriate flux driver
    string::size_type jsimend = flux.find_first_of(":",0);
    if(jsimend==string::npos) {
       LOG("gevgen_upmu", pFATAL)
           << "You need to specify the flux file source";
       PrintSyntax();
       gAbortingInErr = true;
       exit(1);
    }
    gOptFluxSim = flux.substr(0,jsimend);
    for(string::size_type i=0; i<gOptFluxSim.size(); i++) {
       gOptFluxSim[i] = toupper(gOptFluxSim[i]);
    }
    if((gOptFluxSim != "FLUKA") && (gOptFluxSim != "BGLRS")) {
        LOG("gevgen_upmu", pFATAL)
             << "The flux file source needs to be one of <FLUKA,BGLRS>";
        PrintSyntax();
        gAbortingInErr = true;
        exit(1);
    }
    // now get the list of input files and the corresponding neutrino codes.
    flux.erase(0,jsimend+1);
    vector<string> fluxv = utils::str::Split(flux,",");
    vector<string>::const_iterator fluxiter = fluxv.begin();
    for( ; fluxiter != fluxv.end(); ++fluxiter) {
       string filename_and_pdg = *fluxiter;
       string::size_type open_bracket  = filename_and_pdg.find("[");
       string::size_type close_bracket = filename_and_pdg.find("]");
       if (open_bracket ==string::npos ||
           close_bracket==string::npos)
       {
           LOG("gevgen_upmu", pFATAL)
              << "You made an error in specifying the flux info";
           PrintSyntax();
           gAbortingInErr = true;
           exit(1);
       }
       string::size_type ibeg = 0;
       string::size_type iend = open_bracket;
       string::size_type jbeg = open_bracket+1;
       string::size_type jend = close_bracket;
       string flux_filename   = filename_and_pdg.substr(ibeg,iend-ibeg);
       string neutrino_pdg    = filename_and_pdg.substr(jbeg,jend-jbeg);
       gOptFluxFiles.insert(
          map<int,string>::value_type(atoi(neutrino_pdg.c_str()), flux_filename));
    }
    if(gOptFluxFiles.size() == 0) {
       LOG("gevgen_upmu", pFATAL)
          << "You must specify at least one flux file!";
       PrintSyntax();
       gAbortingInErr = true;
       exit(1);
    }

  } else {
    LOG("gevgen_upmu", pFATAL) << "No flux info was specified! Use the -f option.";
    PrintSyntax();
    gAbortingInErr = true;
    exit(1);
  }

  //
  // *** random number seed
  //
  if( parser.OptionExists("seed") ) {
    LOG("gevgen_upmu", pINFO) << "Reading random number seed";
    gOptRanSeed = parser.ArgAsLong("seed");
  } else {
    LOG("gevgen_upmu", pINFO) << "Unspecified random number seed - Using default";
    gOptRanSeed = -1;
  }

  //
  // *** input cross-section file
  //
  if( parser.OptionExists("cross-sections") ) {
    LOG("gevgen_upmu", pINFO) << "Reading cross-section file";
    gOptInpXSecFile = parser.ArgAsString("cross-sections");
  } else {
    LOG("gevgen_upmu", pINFO) << "Unspecified cross-section file";
    gOptInpXSecFile = "";
  }


  //
  // print-out summary
  //

  PDGLibrary * pdglib = PDGLibrary::Instance();

  ostringstream fluxinfo;
  fluxinfo << "Using " << gOptFluxSim << " flux files: ";
  map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
  for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
     int neutrino_code = file_iter->first;
     string filename   = file_iter->second;
     TParticlePDG * p = pdglib->Find(neutrino_code);
     if(p) {
        string name = p->GetName();
        fluxinfo << "(" << name << ") -> " << filename << " / ";
     }
  }

  ostringstream expinfo;
  if(gOptNev > 0) { expinfo << gOptNev << " events";   }

  LOG("gevgen_atmo", pNOTICE)
     << "\n\n"
     << utils::print::PrintFramedMesg("gevgen_upmu job configuration");

  LOG("gevgen_upmu", pNOTICE)
     << "\n"
     << "\n @@ Run number: " << gOptRunNu
     << "\n @@ Random number seed: " << gOptRanSeed
     << "\n @@ Using cross-section file: " << gOptInpXSecFile
     << "\n @@ Flux"
     << "\n\t" << fluxinfo.str()
     << "\n @@ Exposure"
     << "\n\t" << expinfo.str()
     << "\n\n";
}

double GetCrossSection	(	int	nu_code,
		double	Enu,
		double	Emu
	)

Definition at line 538 of file gUpMuFluxGen.cxx.

{
// Get the interaction cross section for a neutrino of energy Enu
// to produce a muon of energy Emu.

  double dxsec_dxdy = 0;

  if ( Emu >= Enu ) {
    dxsec_dxdy = 0;
  }
  else {
    // get the algorithm factory & config pool
    AlgFactory * algf = AlgFactory::Instance();

    // instantiate some algorithms
    AlgId id("genie::QPMDISPXSec","Default");
    Algorithm * alg = algf->AdoptAlgorithm(id);
    XSecAlgorithmI * xsecalg = dynamic_cast<XSecAlgorithmI*>(alg);

    Interaction * vp = Interaction::DISCC(kPdgTgtFreeP, kPdgProton,  nu_code, Enu);
    Interaction * vn = Interaction::DISCC(kPdgTgtFreeN, kPdgNeutron, nu_code, Enu);

    // Integrate over x [0,1] and y [0,1-Emu/Enu], 100 steps for each.
    double ymax = 1 - Emu/Enu;

    double dy   = ymax/10;
    double dx   = 0.1;
    double x    = 0;
    double y    = 0;

    for (int i = 0; i<=10; i++){
      x = i*dx;
      for (int j = 0; j<=10; j++){
        y = j * dy;
        double W  = 0;
        double Q2 = 0;
        utils::kinematics::XYtoWQ2(Enu,kNucleonMass,W,Q2,x,y);
        vp->KinePtr()->Setx(x);
        vp->KinePtr()->Sety(y);
        vp->KinePtr()->SetQ2(Q2);
        vp->KinePtr()->SetW(W);
        vn->KinePtr()->Setx(x);
        vn->KinePtr()->Sety(y);
        vn->KinePtr()->SetQ2(Q2);
        vn->KinePtr()->SetW(W);

        dxsec_dxdy += 0.5*(xsecalg->XSec(vp,kPSxyfE) + xsecalg->XSec(vn,kPSxyfE)) / units::cm2;
        }
    }

    delete vp;
    delete vn;
  }

  return dxsec_dxdy;
}

TVector3 GetDetectorVertex	(	double	CosTheta,
		double	Enu
	)

Definition at line 288 of file gUpMuFluxGen.cxx.

{
// Find the point at which the muon crosses the detector. Returns 0,0,0 if the muon misses.
// Detector is a cube of side length l (=gOptDetectorSide).

  // Get a RandomGen instance
  RandomGen * rnd = RandomGen::Instance();

  // Generate random phi.
  double phi = 2.*kPi* rnd->RndFlux().Rndm();

  // Set distance at which incoming muon is considered
  double rad = 0.87*gOptDetectorSide;

  // Set transverse radius of a circle
  double rad_trans = 0.87*gOptDetectorSide;

  // Get necessary trig
  double sintheta = TMath::Sqrt(1-costheta*costheta);
  double cosphi   = TMath::Cos(phi);
  double sinphi   = TMath::Sin(phi);

  // Compute the muon position at distance Rad.
  double z = rad * costheta;
  double y = rad * sintheta * cosphi;
  double x = rad * sintheta * sinphi;

  // Displace muon randomly on a circular surface of radius rad_trans,
  // perpendicular to a sphere radius rad at that position [x,y,z].
  TVector3 vec(x,y,z);                // vector towards selected point
  TVector3 dvec1 = vec.Orthogonal();  // orthogonal vector
  TVector3 dvec2 = dvec1;             // second orthogonal vector
  dvec2.Rotate(-kPi/2.0,vec);         // rotate second vector by 90deg -> Cartesian coords

  // Select a random point on the transverse surface, within radius rad_trans
  double psi = 2 * kPi * rnd->RndFlux().Rndm();      // rndm angle [0,2pi]
  double random = rnd->RndFlux().Rndm();             // rndm number [0,1]
  dvec1.SetMag(TMath::Sqrt(random)*rad_trans*TMath::Cos(psi));
  dvec2.SetMag(TMath::Sqrt(random)*rad_trans*TMath::Sin(psi));
  x += dvec1.X() + dvec2.X();    // x-coord of a point the muon passes through
  y += dvec1.Y() + dvec2.Y();    // y-coord of a point the muon passes through
  z += dvec1.Z() + dvec2.Z();    // z-coord of a point the muon passes through

  // Find out if the muon passes through any side of the detector.

  // Get momentum vector
  double pz = Enu * costheta;
  double py = Enu * sintheta * sinphi;
  double px = Enu * sintheta * cosphi;
  TVector3 p3(-px,-py,-pz);

  // Set up other vectors needed for line-box intersection.
  TVector3 x3(x,y,z);
  TVector3 temp3(x3);
  TVector3 Hit3(0,0,0);

  // Find out if the line of the muon intersects the detector.
  bool HitFound = false;
  double l = gOptDetectorSide;
  TVector3 unit(0,0,0);
  unit = p3.Unit();
  double unitx = unit.X();
  double unity = unit.Y();
  double unitz = unit.Z();

  // Check to see if muon intersects with z=-l/2 surface.
  if (x3.Z() < -l/2 && unitz > 0 && !HitFound){
    while (temp3.Z() < -l/2){
      temp3 += unit;
    }
    if ( -l/2<temp3.X() && l/2>temp3.X() && -l/2<temp3.Y() && l/2>temp3.Y() ){
      Hit3.SetXYZ(temp3.X(),temp3.Y(),-l/2);
      HitFound = true;
    }
    else temp3 = x3;
  }

  // Check to see if muon intersects with x=l/2 surface.
  if (x3.X() > l/2 && unitx < 0 && !HitFound){
    while (temp3.X() > l/2){
      temp3 += p3.Unit();
    }
    if ( -l/2<temp3.Y() && l/2>temp3.Y() && -l/2<temp3.Z() && l/2>temp3.Z() ){
      Hit3.SetXYZ(l/2,temp3.Y(),temp3.Z());
      HitFound = true;
    }
    else temp3 = x3;
  }

  // Check to see if muon intersects with y=l/2 surface.
  if (x3.Y() > l/2 && unity < 0 && !HitFound){
    while (temp3.Y() > l/2){
      temp3 += p3.Unit();
    }
    if ( -l/2<temp3.X() && l/2>temp3.X() && -l/2<temp3.Z() && l/2>temp3.Z() ){
      Hit3.SetXYZ(temp3.X(),l/2,temp3.Z());
      HitFound = true;
    }
    else temp3 = x3;
  }

  // Check to see if muon intersects with x=-l/2 surface.
  if (x3.X() < -l/2 && unitx > 0 && !HitFound){
    while (temp3.X() < -l/2){
      temp3 += p3.Unit();
    }
    if ( -l/2<temp3.Y() && l/2>temp3.Y() && -l/2<temp3.Z() && l/2>temp3.Z() ){
      Hit3.SetXYZ(-l/2,temp3.Y(),temp3.Z());
      HitFound = true;
    }
    else temp3 = x3;
  }

  // Check to see if muon intersects with y=-l/2 surface.
  if (x3.Y() < -l/2 && unity > 0 && !HitFound){
    while (temp3.Y() < -l/2){
      temp3 += p3.Unit();
    }
    if ( -l/2<temp3.X() && l/2>temp3.X() && -l/2<temp3.Z() && l/2>temp3.Z() ){
      Hit3.SetXYZ(temp3.X(),-l/2,temp3.Z());
      HitFound = true;
    }
    else temp3 = x3;
  }

  return Hit3;
}

TH1D * GetEmuPdf	(	double	Enu,
		double	costheta,
		const TH3D *	pdf3d
	)

Definition at line 518 of file gUpMuFluxGen.cxx.

{
  // Build 1-D Emu pdf
  int Emu_nbins = pdf3d->GetXaxis()->GetNbins();
  const double * Emu_bins = pdf3d->GetXaxis()->GetXbins()->GetArray();
  TH1D * pdf_Emu = new TH1D("pdf_Emu","",Emu_nbins,Emu_bins);
  pdf_Emu->SetDirectory(0);

  // Find appropriate bins for Enu and CosTheta.
  int Enu_bin      = pdf3d->GetYaxis()->FindBin(Enu);
  int costheta_bin = pdf3d->GetZaxis()->FindBin(costheta);

  // For those Enu and costheta bins, build Emu pdf
  for (int Emu_bin = 1; Emu_bin < pdf_Emu->GetNbinsX(); Emu_bin++) {
    pdf_Emu->SetBinContent(Emu_bin, pdf3d->GetBinContent(Emu_bin,Enu_bin,costheta_bin));
  }

  return pdf_Emu;
}

GFluxI* GetFlux

(

void

)

Definition at line 595 of file gUpMuFluxGen.cxx.

{
  GFluxI * flux_driver = 0;

#ifdef __GENIE_FLUX_DRIVERS_ENABLED__

  // Instantiate appropriate concrete flux driver
  GAtmoFlux * atmo_flux_driver = 0;
  if(gOptFluxSim == "FLUKA") {
     GFLUKAAtmoFlux * fluka_flux = new GFLUKAAtmoFlux;
     atmo_flux_driver = dynamic_cast<GAtmoFlux *>(fluka_flux);
  } else
  if(gOptFluxSim == "BGLRS") {
     GBGLRSAtmoFlux * bartol_flux = new GBGLRSAtmoFlux;
     atmo_flux_driver = dynamic_cast<GAtmoFlux *>(bartol_flux);
  } else {
     LOG("gevgen_upmu", pFATAL) << "Uknonwn flux simulation: " << gOptFluxSim;
     gAbortingInErr = true;
     exit(1);
  }


  atmo_flux_driver->GenerateWeighted(true);

  // Configure GAtmoFlux options (common to all concrete atmospheric flux drivers)
  // set flux files:
  map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
  for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
    int neutrino_code = file_iter->first;
    string filename   = file_iter->second;
    atmo_flux_driver->AddFluxFile(neutrino_code, filename);
  }
  atmo_flux_driver->LoadFluxData();
  // configure flux generation surface:
  atmo_flux_driver->SetRadii(1, 1);
  // Cast to GFluxI, the generic flux driver interface
  flux_driver = dynamic_cast<GFluxI *>(atmo_flux_driver);

#else
  LOG("gevgen_upmu", pFATAL) << "You need to enable the GENIE flux drivers first!";
  LOG("gevgen_upmu", pFATAL) << "Use --enable-flux-drivers at the configuration step.";
  gAbortingInErr = true;
  exit(1);
#endif

  return flux_driver;
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 172 of file gUpMuFluxGen.cxx.

{
  // Parse command line arguments
  GetCommandLineArgs(argc,argv);

  if ( ! RunOpt::Instance()->Tune() ) {
    LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
    exit(-1);
  }
  RunOpt::Instance()->BuildTune();

  // Init
  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
  utils::app_init::RandGen(gOptRanSeed);
  utils::app_init::XSecTable(gOptInpXSecFile, true);


  // Get requested flux driver
  GFluxI * flux_driver = GetFlux();

  // Create output tree to store generated up-going muons
  TTree * ntupmuflux = new TTree("ntupmuflux","GENIE Upgoing Muon Event Tree");
  // Tree branches
  int    brIev        = 0;      // Event number
  int    brNuCode     = 0;      // Neutrino PDG code
  double brEmu        = 0;      // Muon energy (GeV)
  double brEnu        = 0;      // Neutrino energy (GeV)
  double brCosTheta   = 0;      // Neutrino cos(zenith angle), muon assumed to be collinear
  double brWghtFlxNu  = 0;      // Weight associated with the current flux neutrino
  double brWghtEmuPdf = 0;      // Weight associated with the Emu pdf for current Enu, costheta bins
  double brVx         = 0;      // Muon x (mm) - intersection with box surrounding the detector volume
  double brVy         = 0;      // Muon y (mm) - intersection with box surrounding the detector volume
  double brVz         = 0;      // Muon z (mm) - intersection with box surrounding the detector volume
  double brXSec       = 0;      //
  ntupmuflux->Branch("iev",          &brIev,         "iev/I"         );
  ntupmuflux->Branch("nu_code",      &brNuCode,      "nu_code/I"     );
  ntupmuflux->Branch("Emu",          &brEmu,         "Emu/D"         );
  ntupmuflux->Branch("Enu",          &brEnu,         "Enu/D"         );
  ntupmuflux->Branch("costheta",     &brCosTheta,    "costheta/D"    );
  ntupmuflux->Branch("wght_fluxnu",  &brWghtFlxNu,   "wght_fluxnu/D" );
  ntupmuflux->Branch("wght_emupdf",  &brWghtEmuPdf,  "wght_emupdf/D" );
  ntupmuflux->Branch("vx",           &brVx,          "vx/D"          );
  ntupmuflux->Branch("vy",           &brVy,          "vy/D"          );
  ntupmuflux->Branch("vz",           &brVz,          "vz/D"          );
  ntupmuflux->Branch("xsec",         &brXSec,        "xsec/D"        );

  // Build 3-D pdfs describing the the probability of a muon neutrino (or anti-neutrino)
  // of energy Enu and zenith angle costheta producing a mu- (or mu+) of energy E_mu
  TH3D * pdf3d_numu    = BuildEmuEnuCosThetaPdf (kPdgNuMu    );
  TH3D * pdf3d_numubar = BuildEmuEnuCosThetaPdf (kPdgAntiNuMu);

  // Up-going muon event loop
  for(brIev = 0; brIev < gOptNev; brIev++) {

    // Loop until upgoing nu is generated.
    GenerateUpNu(flux_driver);

    // Get neutrino code, Enu, costheta and weight for the generated neutrino
    brNuCode    = flux_driver->PdgCode();
    brEnu       = flux_driver->Momentum().E();
    brCosTheta  = -1. * flux_driver->Momentum().Pz() / flux_driver->Momentum().Vect().Mag();
    brWghtFlxNu = flux_driver->Weight();

    LOG("gevgen_upmu", pNOTICE)
        << "Generated flux neutrino: code = " << brNuCode
        << ", Ev = " << brEnu << " GeV"
        << ", cos(theta) = " << brCosTheta
        << ", weight = " << brWghtFlxNu;

    // Get Emu pdf my slicing the 3-D Enu,Emu,costheta pdf.
    TH3D * pdf3d  = (brNuCode == kPdgNuMu) ? pdf3d_numu : pdf3d_numubar;
    TH1D * pdfEmu = GetEmuPdf(brEnu,brCosTheta,pdf3d);

    // Get a random Emu from the pdf, and get the weight for that Emu.
    brEmu        = pdfEmu->GetRandom();
    brWghtEmuPdf = pdfEmu->Integral("width");
    LOG("gevgen_upmu", pNOTICE)
        << "Selected muon has energy Emu = " << brEmu
        << " and Emu pdg weight = " << brWghtEmuPdf;

    // Randomly select the point at which the neutrino crosses the detector.
    // assumes a simple cube detector of side length gOptDetectorSide.
    TVector3 Vertex = GetDetectorVertex(brCosTheta,brEnu);
    brVx = Vertex.X();
    brVy = Vertex.Y();
    brVz = Vertex.Z();

    LOG("gevgen_upmu", pNOTICE)
      << "Generated muon position: (" << brVx << ", " << brVy << ", " << brVz << ") m";

    // Save the cross section for the interaction generated.
    brXSec = GetCrossSection(brNuCode, brEnu, brEmu);

    // save all relevant values to the ntuple
    ntupmuflux->Fill();

    // Clean-up
    delete pdfEmu;
  }

  // Save the muon ntuple and calculate 3-D pdfs

  ostringstream outfname;
  outfname << "./genie." << gOptRunNu << ".upmu.root";
  TFile outf(outfname.str().c_str(), "recreate");
  ntupmuflux    -> Write("ntupmu");
  pdf3d_numu    -> Write("pdf3d_numu");
  pdf3d_numubar -> Write("pdf3d_numubar");
  outf.Close();

  // Clean-up
  delete flux_driver;

  return 0;
}

void PrintSyntax

(

void

)

Definition at line 843 of file gUpMuFluxGen.cxx.

{
  LOG("gevgen_upmu", pFATAL)
   << "\n **Syntax**"
   << "\n gevgen_upmu [-h]"
   << "\n             [-r run#]"
   << "\n              -f simulation:flux_file[neutrino_code],..."
   << "\n              -n n_of_events,"
   << "\n             [-d detector side length (mm)]"
   << "\n             [--seed random_number_seed]"
   << "\n              --cross-sections xml_file"
   << "\n            [--message-thresholds xml_file]"
   << "\n"
   << " Please also read the detailed documentation at http://www.genie-mc.org"
   << "\n";
}

double ProbabilityEmu	(	int	nu_code,
		double	Enu,
		double	Emu
	)

Definition at line 435 of file gUpMuFluxGen.cxx.

{
// Calculate the probability of an incoming neutrino of energy Enu
// generating a muon of energy Emu.
  double dxsec_dxdy = GetCrossSection(nu_code,Enu,Emu);
  double Int = e * constants::kNA * dxsec_dxdy / (a * (1 + Emu/e));
  return Int;
}

Variable Documentation

const double a = 2e+6

Definition at line 164 of file gUpMuFluxGen.cxx.

const double e = 500e+9

Definition at line 165 of file gUpMuFluxGen.cxx.

double gOptDetectorSide

Definition at line 159 of file gUpMuFluxGen.cxx.

map<int,string> gOptFluxFiles

Definition at line 157 of file gUpMuFluxGen.cxx.

string gOptFluxSim

Definition at line 156 of file gUpMuFluxGen.cxx.

string gOptInpXSecFile

Definition at line 161 of file gUpMuFluxGen.cxx.

int gOptNev = -1

Definition at line 158 of file gUpMuFluxGen.cxx.

long int gOptRanSeed

Definition at line 160 of file gUpMuFluxGen.cxx.

Long_t gOptRunNu

Definition at line 155 of file gUpMuFluxGen.cxx.

double kDefOptDetectorSide = 1e+5

Definition at line 169 of file gUpMuFluxGen.cxx.