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.

References h3, MECModelEnuComparisons::i, calib::j, test_ParserArtEvents::log, and ProbabilityEmu().

Referenced by main().

{
// 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.

References genie::GFluxI::GenerateNext(), sanity_check_grl::keep, genie::kPdgAntiNuMu, genie::kPdgNuMu, LOG, genie::GFluxI::Momentum(), make_associated_cosmic_defs::p4, genie::GFluxI::PdgCode(), and pINFO.

Referenced by main().

{
// 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.

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), exit(), shutoffs::filename, genie::PDGLibrary::Find(), flux, genie::gAbortingInErr, gOptDetectorSide, gOptFluxFiles, gOptFluxSim, gOptInpXSecFile, gOptNev, gOptRanSeed, gOptRunNu, plot_validation_datamc::help, MECModelEnuComparisons::i, genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), kDefOptDetectorSide, LOG, genie::CmdLnArgParser::OptionExists(), plot_validation_datamc::parser, pDEBUG, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), genie::RunOpt::ReadFromCommandLine(), and genie::utils::str::Split().

Referenced by main().

{
// 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.

References genie::AlgFactory::AdoptAlgorithm(), genie::units::cm2, genie::Interaction::DISCC(), dx, dy, MECModelEnuComparisons::i, genie::AlgFactory::Instance(), calib::j, genie::Interaction::KinePtr(), genie::constants::kNucleonMass, genie::kPdgNeutron, genie::kPdgProton, genie::kPdgTgtFreeN, genie::kPdgTgtFreeP, genie::kPSxyfE, genie::utils::kinematics::Q2(), genie::Kinematics::SetQ2(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), W, submit_syst::x, genie::XSecAlgorithmI::XSec(), genie::utils::kinematics::XYtoWQ2(), submit_syst::y, and ymax.

Referenced by main(), and ProbabilityEmu().

{
// 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.

References gOptDetectorSide, genie::RandomGen::Instance(), kPi, submit_hadd::l, make_associated_cosmic_defs::p3, Munits::rad, generate_hists::rnd, genie::RandomGen::RndFlux(), submit_syst::x, submit_syst::y, and test::z.

Referenced by main().

{
// 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.

Referenced by main().

{
  // 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.

References genie::flux::GAtmoFlux::AddFluxFile(), exit(), shutoffs::filename, genie::gAbortingInErr, genie::flux::GAtmoFlux::GenerateWeighted(), gOptFluxFiles, gOptFluxSim, genie::flux::GAtmoFlux::LoadFluxData(), LOG, pFATAL, and genie::flux::GAtmoFlux::SetRadii().

Referenced by main().

{
  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.

References BuildEmuEnuCosThetaPdf(), genie::RunOpt::BuildTune(), exit(), GenerateUpNu(), GetCommandLineArgs(), GetCrossSection(), GetDetectorVertex(), GetEmuPdf(), GetFlux(), gOptInpXSecFile, gOptNev, gOptRanSeed, gOptRunNu, genie::RunOpt::Instance(), genie::kPdgAntiNuMu, genie::kPdgNuMu, LOG, genie::utils::app_init::MesgThresholds(), genie::GFluxI::Momentum(), outf, genie::GFluxI::PdgCode(), pFATAL, pNOTICE, genie::utils::app_init::RandGen(), genie::GFluxI::Weight(), Write(), and genie::utils::app_init::XSecTable().

{
  // 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.

References LOG, and pFATAL.

Referenced by GetCommandLineArgs().

{
  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.

References a, e, GetCrossSection(), and genie::constants::kNA.

Referenced by BuildEmuEnuCosThetaPdf().

{
// 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.

Referenced by ana::_HistAxis< Var >::_HistAxis(), genie::Spline::Add(), evdb::View2D::AddArc(), ana::KrigeKernel::AddPoint(), ana::AlmostEqual(), slicer::S4DParamCalc::analyze(), bpfit::BPFdEdxHistoMaker::analyze(), slicer::SlicerAna::analyze(), bpfit::BPFTmvaTrainer::analyze(), air::AirKalmanAna::analyze(), novaddt::HotMapMaker::AppendCold(), novaddt::HotMapMaker::AppendHot(), novaddt::HotMapMaker::AppendRate(), novaddt::utils::DetectorUtils::AttenuationInFiber(), genie::KNOHadronization::AverageChMult(), gibuu::GiBUURegen::Basis(), calib::BetheBlochAnalytic::BetheBlochAnalytic(), genie::utils::nuclear::BindEnergy(), by_charge(), genie::QPMDISStrucFuncBase::Calculate(), beamlinereco::WCTrackAlg::calculateTheMomentum(), earms::GridSearch::ChooseBestProngDir(), genie::utils::kinematics::CohQ2Lim(), ana::CutBase::CombineExposures(), boost::python::map_indexing_suite< Container, NoProxy, final_map_derived_policies< Container, NoProxy > >::compare_index(), genie::RSHelicityAmplModelCC::Compute(), genie::RSHelicityAmplModelNCn::Compute(), genie::RSHelicityAmplModelNCp::Compute(), genie::alvarezruso::AREikonalSolution::Con(), CLHEP::MixMaxRng::convert1double(), genie::geometry::PlaneParam::ConvertMaster2Top(), genie::mueloss::BetheBlochModel::dE_dx(), fhicl::detail::decode(), fhicl::detail::decode_tuple(), calib::BetheBlochAnalytic::Delta(), genie::Spline::Divide(), genie::mueloss::gsl::BezrukovBugaevIntegrand::DoEval(), genie::mueloss::gsl::PetrukhinShestakovIntegrand::DoEval(), genie::utils::gsl::dXSec_Log_Wrapper::DoEval(), earms::ElasticArms::E(), earmsnt(), dt::Cand::EstimateStraightLine(), osc::_PMNS< T >::EvalEqn11(), osc::_PMNS< T >::EvalEqn2(), fiducial_mass(), art::View< T >::fill(), fhicl::detail::find_an_any(), util::find_median(), ana::GenieMultiverseSpectra::FindBandBoundaries(), lem::FindMatches(), cmf::SpectrumPredictionMaker::FindOscillationPoints(), fit_thresh_corrs(), fit_thresh_corrs_fb(), Flatplots(), g4lbneCrossChecks(), g4minervaCrossChecks(), GammaFit(), genie::BBA03ELFormFactorsModel::Gen(), ana::get3FlavorAna2020SmallXsecSysts(), fhicl::ParameterSet::get_one_(), ana::getAna2017SmallXsecSysts(), ana::getAna2018SmallXsecSysts(), genie::Spline::GetAsTSpline(), numusand::NumuSandFxs::getAveTrackdEdxLast4Cells(), numusand::NumuSandFxs::getAveTrackdEdxLast6Cells(), numusand::NumuSandFxs::getAveTrackdEdxLast8Cells(), trackinfo::TrackInfoFxs::getAveTrackdEdxLastCells(), HepGeom::BasicVector3D< float >::getEta(), genie::alvarezruso::AlvarezRusoCOHPiPDXSec::GetLeptonMass(), evg::MCTruthToDk2NuHackItr::GetMCTruthPtr(), ana::GenieMultiverseSpectra::GetSigmaHistograms(), LSTME::Png3dVars::getSortedIndices(), SliceLID::Png3dVars::getSortedIndices(), cheat::BackTracker::HitCollectionEfficiency(), genie::utils::kinematics::InelYLim_X(), genie::utils::kinematics::electromagnetic::InelYLim_X(), daqdataformats::VERSION_NAMESPACE::RawConfigurationBlock::init(), daqdataformats::VERSION_NAMESPACE::RawSummaryDroppedMicroblock::init(), daqdataformats::VERSION_NAMESPACE::RawTriggerTime::init(), rb::Track::InterpolateDir(), rb::Track::InterpolateXY(), fhicl::ParameterSet::is_key_to_atom(), isCloseEnough(), genie::geometry::PlaneParam::IsValid(), hough::IslandsTH2::LabelIslands(), ana::LabelsAndBins::LabelsAndBins(), earms::ElasticArms::LeadingArm(), novaddt::NNFastMMTrigger::LinFit(), zcl::FmmTrackerValidation::LinFit(), novaddt::UpMuTrigger::LinFit(), zcl::FastMMStudy::LinFit(), novaddt::UpMuAna::LinFit(), upmuana::UpMuAnalysis::LinFit(), upmuana::UpMuRecoAna::LinFit(), novaddt::UpMuTrigger::LinFitLLR(), novaddt::UpMuAna::LinFitLLR(), upmuana::LlrUpmu::LLR(), upmuana::UpMuAnalysis::LLR(), upmuana::UpMuRecoAna::LLR(), genie::NievesQELCCPXSec::LmunuAnumu(), nerd::NERDEval::load_anchors(), g4n::LoadNeutronFates(), genie::MECHadronTensor::LoadTensorTables(), make_cm_pullterm_pdf(), make_xsec_tuning_hists_mp5(), make_xsec_wgts_2018_hists(), nerd::NERDProng::MakeCluster(), novaddt::HotMapMaker::MakeMask(), MakeTestStandPlots_AllRuns(), MakeTestStandPlots_NewRuns(), MAX(), nl::TimeSeriesData::MeanValue(), dt::DiscreteTracker::MergeParasiteTracks(), ana::PredictionExtendToPeripheral::MergePeripheral(), min(), evd::RawDrawingOptions::ModuleList(), evd::RecoDrawingOptions::ModuleList(), hough::MultiHough2P::MultiMap(), genie::Spline::Multiply(), ana::TargetCount::NNuclei(), genie::geometry::PlaneParam::Normalize(), ana::Nue2018ThreeBinAxis(), Nue2018ThreeBinAxis_Reduced(), ana::Nue2018ThreeBinAxis_Reduced(), nus17_fiducial_mass(), OnMonPlotMaker(), ana::CutBase::operator!(), art::operator!(), ana::CutBase::operator&&(), art::operator&&(), order_cuts::operator()(), order_chns::operator()(), ana::CutBase::operator()(), std::less< art::InputTag >::operator()(), ana::VarBase::operator*(), HepGeom::BasicVector3D< T >::operator*(), HepGeom::BasicVector3D< float >::operator*=(), ana::VarBase::operator+(), ana::VarBase::operator-(), ana::VarBase::operator/(), HepGeom::BasicVector3D< T >::operator/(), HepGeom::BasicVector3D< float >::operator/=(), lem::LiteHit::operator<(), ana::operator<(), nova::database::Column::operator<(), nutools::dbi::Column::operator<(), nova::dbi::Column::operator<(), om::OStream::operator<<(), nova::dbi::operator<<(), nova::database::operator<<(), genie::flux::operator<<(), nova::database::Column::operator<=(), nutools::dbi::Column::operator<=(), nova::dbi::Column::operator<=(), art::operator==(), nova::database::Column::operator==(), ana::operator>(), nova::database::Column::operator>(), nutools::dbi::Column::operator>(), nova::dbi::Column::operator>(), nova::database::Column::operator>=(), nutools::dbi::Column::operator>=(), nova::dbi::Column::operator>=(), ana::CutBase::operator||(), art::operator||(), pie_arc(), genie::geometry::PlaneParam::PlaneParam(), plot_3NDvsFD_FHC(), plot_3NDvsFD_RHC(), PlotPionPlots(), sim::PhotonSignal::PoissonLambda(), g4n::NeutronSubstitutionProcess::PostStepDoIt(), pots(), genie::P33PaschosLalakulichPXSec::PPiStar(), evg::MCTruthToDk2NuHackItr::pretty_darn_close(), genie::geometry::PlaneParam::Print(), ProbabilityEmu(), genie::FGMBodekRitchie::ProbDistro(), genie::LocalFGM::ProbDistro(), cvntf::CVNCosmicTFFilt::produce(), cvn::CVNCosmicMapper::produce(), earms::ElasticArmsHS::produce(), bpfit::BPFEnergyEstimator::produce(), bpfit::BPFPIdMaker::produce(), genie::SmithMonizUtils::Q2lim1_SM(), genie::SmithMonizUtils::Q2lim2_SM(), ana::RefineSeeds(), HepGeom::Reflect3D::Reflect3D(), genie::NievesQELCCPXSec::relLindhardIm(), hough::MultiHough2P::ReweightHits(), genie::alvarezruso::integrationtools::RG202D(), genie::alvarezruso::integrationtools::RG482D(), art::right(), HepGeom::Rotate3D::Rotate3D(), genie::BYStrucFunc::ScalingVar(), setA(), earms::ElasticArms::SetArm(), ana::covmx::Sample::SetAxis(), SetGraphErrors(), ana::FitSinSq2Theta13::SetValue(), daqdataformats::VERSION_NAMESPACE::RawMicroBlockHeader::setVersion(), HepGeom::BasicVector3D< float >::setX(), HepGeom::BasicVector3D< float >::setY(), HepGeom::BasicVector3D< float >::setZ(), genie::alvarezruso::integrationtools::SG20R(), genie::alvarezruso::integrationtools::SG48R(), show_x_cut(), slid::ShowerLID::ShowerLID(), hough::Hough2P::SmoothMap(), hough::MultiHough2P::SmoothMap(), solve(), ana::MCMCSamples::SortedLLs(), evd::RawDataDrawer::StickFigure(), earms::ElasticArms::SumVi(), art::ValidHandle< T >::swap(), calib::ThresholdCorrMap::ThresholdCorrAt(), art::detail::throw_if_invalid(), lem::dec::Forest::TrainSingle(), rb::Track::Trajectory(), ncid::NCNNKeras::TransMomFraction(), earms::ElasticArms::UpdateDaa(), earms::ElasticArms::UpdateMia(), earms::ElasticArms::UpdateVia(), uptime(), osc::VacuumHamiltonian(), genie::geometry::PlaneParam::Vd(), genie::geometry::PlaneParam::Vn(), genie::SmithMonizUtils::vQES_SM_lim(), fhicl::detail::ParameterWalker< C >::walk_over(), genie::BSKLNBaseRESPXSec2014::XSec(), Y11(), Y12(), Y22(), rb::Track::ZipWith(), and calib::ThresholdCorrMap::~ThresholdCorrMap().

const double e = 500e+9

Definition at line 165 of file gUpMuFluxGen.cxx.

Referenced by ProbabilityEmu().

double gOptDetectorSide

Definition at line 159 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and GetDetectorVertex().

map<int,string> gOptFluxFiles

Definition at line 157 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

string gOptFluxSim

Definition at line 156 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and GetFlux().

string gOptInpXSecFile

Definition at line 161 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and main().

int gOptNev = -1

Definition at line 158 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and main().

long int gOptRanSeed

Definition at line 160 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and main().

Long_t gOptRunNu

Definition at line 155 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double kDefOptDetectorSide = 1e+5

Definition at line 169 of file gUpMuFluxGen.cxx.

Referenced by GetCommandLineArgs().