gSplineXml2Root.cxx File Reference

#include <cassert>
#include <string>
#include <sstream>
#include <vector>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TDirectory.h>
#include <TPostScript.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TGraph.h>
#include <TPaveText.h>
#include <TString.h>
#include <TH1F.h>
#include "Framework/ParticleData/BaryonResonance.h"
#include "Framework/ParticleData/BaryonResUtils.h"
#include "Framework/Conventions/XmlParserStatus.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/EventGen/InteractionList.h"
#include "Framework/EventGen/GEVGDriver.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/Spline.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Registry/Registry.h"
#include "Framework/Algorithm/AlgConfigPool.h"

Go to the source code of this file.

Functions
void	LoadSplines (void)
GEVGDriver	GetEventGenDriver (void)
void	SaveToPsFile (void)
void	SaveGraphsToRootFile (void)
void	SaveNtupleToRootFile (void)
void	GetCommandLineArgs (int argc, char **argv)
void	PrintSyntax (void)
PDGCodeList	GetPDGCodeListFromString (std::string s)
int	main (int argc, char **argv)
void	FormatXSecGraph (TGraph *g)

Variables
string	gOptXMLFilename
string	gOptROOTFilename
PDGCodeList	gOptProbePdgList
PDGCodeList	gOptTgtPdgList
int	gOptProbePdgCode
int	gOptTgtPdgCode
bool	gWriteOutPlots
double	gEmin
double	gEmax
bool	gInlogE
int	kNP = 300
int	kNSplineP = 1000
const int	kPsType = 111

Function Documentation

void FormatXSecGraph

(

TGraph *

g

)

Definition at line 516 of file gSplineXml2Root.cxx.

{
  g->SetTitle("GENIE cross section graph");
  g->GetXaxis()->SetTitle("Ev (GeV)");
  g->GetYaxis()->SetTitle("#sigma_{nuclear} (10^{-38} cm^{2})");
}

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 1356 of file gSplineXml2Root.cxx.

{
  LOG("gspl2root", 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);

  // input XML file name:
  if( parser.OptionExists('f') ) {
    LOG("gspl2root", pINFO) << "Reading input XML filename";
    gOptXMLFilename = parser.ArgAsString('f');
  } else {
    LOG("gspl2root", pFATAL) << "Unspecified input XML file!";
    PrintSyntax();
    exit(1);
  }

  // probe PDG code:
  if( parser.OptionExists('p') ) {
    LOG("gspl2root", pINFO) << "Reading probe PDG code";
    gOptProbePdgList = GetPDGCodeListFromString(parser.ArgAsString('p'));
  } else {
    LOG("gspl2root", pFATAL)
       << "Unspecified probe PDG code - Exiting";
    PrintSyntax();
    exit(1);
  }

  // target PDG code:
  if( parser.OptionExists('t') ) {
    LOG("gspl2root", pINFO) << "Reading target PDG code";
    gOptTgtPdgList = GetPDGCodeListFromString(parser.ArgAsString('t'));
  } else {
    LOG("gspl2root", pFATAL)
      << "Unspecified target PDG code - Exiting";
    PrintSyntax();
    exit(1);
  }

  // min,max neutrino energy
  if( parser.OptionExists('e') ) {
    LOG("gspl2root", pINFO) << "Reading neutrino energy";
    string nue = parser.ArgAsString('e');
    // is it just a value or a range (comma separated set of values)
    if(nue.find(",") != string::npos) {
       // split the comma separated list
       vector<string> nurange = utils::str::Split(nue, ",");
       assert(nurange.size() == 2);
       gEmin = atof(nurange[0].c_str());
       gEmax = atof(nurange[1].c_str());
    } else {
      const Registry * val_reg = AlgConfigPool::Instance() -> CommonList( "Param", "Validation" ) ;
      gEmin = val_reg -> GetDouble( "GVLD-Emin" ) ; 
      gEmax = atof(nue.c_str());
      LOG("gspl2root", pDEBUG)
        << "Unspecified Emin - Setting to " << gEmin << " GeV as per configuration";
    }
  } else {
    const Registry * val_reg = AlgConfigPool::Instance() -> CommonList("Param", "Validation" ) ;
    gEmin = val_reg -> GetDouble( "GVLD-Emin" ) ; 
    gEmax = 100;
    LOG("gspl2root", pDEBUG)
      << "Unspecified Emin,Emax - Setting to " << gEmin << ",100 GeV ";

  }
  assert(gEmin<gEmax);

  // output ROOT file name:
  if( parser.OptionExists('o') ) {
    LOG("gspl2root", pINFO) << "Reading output ROOT filename";
    gOptROOTFilename = parser.ArgAsString('o');
  } else {
    LOG("gspl2root", pDEBUG)
       << "Unspecified output ROOT file. Using default: gxsec.root";
    gOptROOTFilename = "gxsec.root";
  }

  // write-out a PS file with plots
  gWriteOutPlots = parser.OptionExists('w');

  // use same abscissa points as splines
  //not yet//gKeepSplineKnots = parser.OptionExists('k');

  gInlogE = parser.OptionExists('l');

  // print the options you got from command line arguments
  LOG("gspl2root", pINFO) << "Command line arguments:";
  LOG("gspl2root", pINFO) << "  Input XML file  = " << gOptXMLFilename;
  LOG("gspl2root", pINFO) << "  Probe PDG code  = " << gOptProbePdgCode;
  LOG("gspl2root", pINFO) << "  Target PDG code = " << gOptTgtPdgCode;
  LOG("gspl2root", pINFO) << "  Min neutrino E  = " << gEmin;
  LOG("gspl2root", pINFO) << "  Max neutrino E  = " << gEmax;
  LOG("gspl2root", pINFO) << "  In logE         = " << gInlogE;
  //not yet//LOG("gspl2root", pINFO) << "  Keep spline knots  = " << (gKeepSplineKnots?"true":"false");
}

GEVGDriver GetEventGenDriver

(

void

)

Definition at line 200 of file gSplineXml2Root.cxx.

{
// create an event genartion driver configured for the specified initial state
// (so that cross section splines will be accessed through that driver as in
// event generation mode)

  InitialState init_state(gOptTgtPdgCode, gOptProbePdgCode);

  GEVGDriver evg_driver;
  evg_driver.SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
  evg_driver.Configure(init_state);
  evg_driver.CreateSplines();
  evg_driver.CreateXSecSumSpline (100, gEmin, gEmax);

  return evg_driver;
}

PDGCodeList GetPDGCodeListFromString

(

std::string

s

)

Definition at line 1465 of file gSplineXml2Root.cxx.

{
  vector<string> isvec = utils::str::Split(s,",");

  // fill in the PDG code list
  PDGCodeList list;
  vector<string>::const_iterator iter;
  for(iter = isvec.begin(); iter != isvec.end(); ++iter) {
    list.push_back( atoi(iter->c_str()) );
  }
  return list;

}

void LoadSplines

(

void

)

Definition at line 191 of file gSplineXml2Root.cxx.

{
// load the cross section splines specified at the cmd line

  XSecSplineList * splist = XSecSplineList::Instance();
  XmlParserStatus_t ist = splist->LoadFromXml(gOptXMLFilename);
  assert(ist == kXmlOK);
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 158 of file gSplineXml2Root.cxx.

{
  GetCommandLineArgs(argc,argv);

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

  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());

  // load the x-section splines xml file specified by the user
  LoadSplines();

  if ( ! XSecSplineList::Instance() -> HasSplineFromTune(RunOpt::Instance() -> Tune() -> Name() ) ) {
    LOG("gspl2root", pWARN) << "No splines loaded for tune " << RunOpt::Instance() -> Tune() -> Name() ;
  }

  for (unsigned int indx_p = 0; indx_p < gOptProbePdgList.size(); ++indx_p ) {
    for (unsigned int indx_t = 0; indx_t < gOptTgtPdgList.size(); ++indx_t ) {
      gOptProbePdgCode = gOptProbePdgList[indx_p];
      gOptTgtPdgCode   = gOptTgtPdgList[indx_t];
      // save the cross section plots in a postscript file
      SaveToPsFile();
      // save the cross section graphs at a root file
      SaveGraphsToRootFile();
    }
  }

  return 0;
}

void PrintSyntax

(

void

)

Definition at line 1456 of file gSplineXml2Root.cxx.

{
  LOG("gspl2root", pNOTICE)
      << "\n\n" << "Syntax:" << "\n"
      << "   gspl2root -f xml_file -p probe_pdg -t target_pdg"
      << "            [-e emin,emax] [-o output_root_file] [-w] [-l]\n"
      << "            [--message-thresholds xml_file]\n";
}

void SaveGraphsToRootFile

(

void

)

Definition at line 523 of file gSplineXml2Root.cxx.

{
  //-- get the event generation driver
  GEVGDriver evg_driver = GetEventGenDriver();

  //-- get the list of interactions that can be simulated by the driver
  const InteractionList * ilist = evg_driver.Interactions();

  //-- check whether the splines will be saved in a ROOT file - if not, exit now
  bool save_in_root = gOptROOTFilename.size()>0;
  if(!save_in_root) {

    LOG("gspl2root", pWARN) << "No Interaction  List available" ;
    return;
  }

  //-- get pdglibrary for mapping pdg codes to names
  PDGLibrary * pdglib = PDGLibrary::Instance();

  //-- check whether the requested filename exists
  //   if yes, then open the file in 'update' mode
  bool exists = !(gSystem->AccessPathName(gOptROOTFilename.c_str()));

  TFile * froot = 0;
  if(exists) froot = new TFile(gOptROOTFilename.c_str(), "UPDATE");
  else       froot = new TFile(gOptROOTFilename.c_str(), "RECREATE");
  assert(froot);

  //-- create directory
  ostringstream dptr;

  string probe_name = pdglib->Find(gOptProbePdgCode)->GetName();
  string tgt_name   = (gOptTgtPdgCode==1000000010) ?
                      "n" : pdglib->Find(gOptTgtPdgCode)->GetName();

  dptr << probe_name << "_" << tgt_name;
  ostringstream dtitle;
  dtitle << "Cross sections for: "
         << pdglib->Find(gOptProbePdgCode)->GetName() << "+"
         << pdglib->Find(gOptTgtPdgCode)->GetName();

  LOG("gspl2root", pINFO)
           << "Will store graphs in root directory = " << dptr.str();
  TDirectory * topdir =
         dynamic_cast<TDirectory *> (froot->Get(dptr.str().c_str()));
  if(topdir) {
     LOG("gspl2root", pINFO)
       << "Directory: " << dptr.str() << " already exists!! Exiting";
     froot->Close();
     delete froot;
     return;
  }

  topdir = froot->mkdir(dptr.str().c_str(),dtitle.str().c_str());
  topdir->cd();

  double   de = (gInlogE) ? (TMath::Log(gEmax)-TMath::Log(gEmin))/(kNSplineP-1) : (gEmax-gEmin)/(kNSplineP-1);
  double * e  = new double[kNSplineP];
  for(int i=0; i<kNSplineP; i++) {  e[i]  = (gInlogE) ? TMath::Exp(TMath::Log(gEmin) + i*de) : gEmin + i*de; }

  double * xs = new double[kNSplineP];

  InteractionList::const_iterator ilistiter = ilist->begin();

  for(; ilistiter != ilist->end(); ++ilistiter) {

    const Interaction * interaction = *ilistiter;

    const ProcessInfo &  proc = interaction->ProcInfo();
    const XclsTag &      xcls = interaction->ExclTag();
    const InitialState & init = interaction->InitState();
    const Target &       tgt  = init.Tgt();

    // graph title
    ostringstream title;

    if      (proc.IsQuasiElastic()     ) { title << "qel";   }
    else if (proc.IsMEC()              ) { title << "mec";   }
    else if (proc.IsResonant()         ) { title << "res";   }
    else if (proc.IsDeepInelastic()    ) { title << "dis";   }
    else if (proc.IsDiffractive()      ) { title << "dfr";   }
    else if (proc.IsCoherentProduction() ) {
      title << "coh";
      if      ( xcls.NSingleGammas() > 0 ) title << "_gamma" ;
      else if ( xcls.NPions() > 0 )        title << "_pion"  ;
      else if ( xcls.NRhos() > 0 )         title << "_rho"   ;
      else                                 title << "_other" ;
    }
    else if (proc.IsCoherentElastic()  ) { title << "cevns"; }
    else if (proc.IsInverseMuDecay()   ) { title << "imd";   }
    else if (proc.IsIMDAnnihilation()  ) { title << "imdanh";}
    else if (proc.IsNuElectronElastic()) { title << "ve";    }
    else if (proc.IsGlashowResonance() ) { title << "glres"; }
    else                                 { 
      LOG("gspl2root", pWARN) << "Process " << proc
                              << " scattering type not recognised: spline not added" ;
      continue;         }
    
    if      (proc.IsWeakCC())  { title << "_cc";      }
    else if (proc.IsWeakNC())  { title << "_nc";      }
    else if (proc.IsWeakMix()) { title << "_ccncmix"; }
    else if (proc.IsEM()    )  { title << "_em";      }
    else if (proc.IsDarkNeutralCurrent() )  { title << "_dark";  }
    else                       {
      LOG("gspl2root", pWARN) << "Process " << proc
                              << " interaction type has not recongnised: spline not added " ;
      continue;            }

    if(tgt.HitNucIsSet()) {
      int hitnuc = tgt.HitNucPdg();
      if      ( pdg::IsProton (hitnuc) ) { title << "_p"; }
      else if ( pdg::IsNeutron(hitnuc) ) { title << "_n"; }
      else if ( pdg::Is2NucleonCluster(hitnuc) )
      {
        if      (hitnuc == kPdgClusterNN) { title << "_nn"; }
        else if (hitnuc == kPdgClusterNP) { title << "_np"; }
        else if (hitnuc == kPdgClusterPP) { title << "_pp"; }
        else {
           LOG("gspl2root", pWARN) << "Can't handle hit 2-nucleon cluster PDG = " << hitnuc;
        }
      }
      else {
        LOG("gspl2root", pWARN) << "Can't handle hit nucleon PDG = " << hitnuc;
      }

      if(tgt.HitQrkIsSet()) {
        int  qrkpdg = tgt.HitQrkPdg();
        bool insea  = tgt.HitSeaQrk();

        if      ( pdg::IsUQuark(qrkpdg)     ) { title << "_u";    }
        else if ( pdg::IsDQuark(qrkpdg)     ) { title << "_d";    }
        else if ( pdg::IsSQuark(qrkpdg)     ) { title << "_s";    }
        else if ( pdg::IsCQuark(qrkpdg)     ) { title << "_c";    }
        else if ( pdg::IsBQuark(qrkpdg)     ) { title << "_b";    }
        else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
        else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
        else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
        else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
        else if ( pdg::IsAntiBQuark(qrkpdg) ) { title << "_bbar"; }

        if(insea) { title << "sea"; }
        else      { title << "val"; }
      }
    }
    if(proc.IsResonant()) {
       Resonance_t res = xcls.Resonance();
       string resname = res::AsString(res);
       resname = str::FilterString(")", resname);
       resname = str::FilterString("(", resname);
       title << "_" << resname.substr(3,4) << resname.substr(0,3);
    }

    if(xcls.IsStrangeEvent()) {
            title << "_strange";
            if(!xcls.IsInclusiveStrange()) { title << xcls.StrangeHadronPdg(); }
    }

    if(xcls.IsCharmEvent()) {
        title << "_charm";
        if(!xcls.IsInclusiveCharm()) { title << xcls.CharmHadronPdg(); }
    }

    if(xcls.IsFinalQuarkEvent()) {
        int  qrkpdg = xcls.FinalQuarkPdg();
        if      ( pdg::IsUQuark(qrkpdg)     ) { title << "_u";    }
        else if ( pdg::IsDQuark(qrkpdg)     ) { title << "_d";    }
        else if ( pdg::IsSQuark(qrkpdg)     ) { title << "_s";    }
        else if ( pdg::IsCQuark(qrkpdg)     ) { title << "_c";    }
        else if ( pdg::IsBQuark(qrkpdg)     ) { title << "_b";    }
        else if ( pdg::IsTQuark(qrkpdg)     ) { title << "_t";    }
        else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
        else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
        else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
        else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
        else if ( pdg::IsAntiBQuark(qrkpdg) ) { title << "_bbar"; }
        else if ( pdg::IsAntiTQuark(qrkpdg) ) { title << "_tbar"; }
    }
    if(xcls.IsFinalLeptonEvent()) {
        int  leppdg = xcls.FinalLeptonPdg();
        if      ( pdg::IsMuon(leppdg)     ) { title << "_mu";     }
        else if ( pdg::IsElectron(leppdg) ) { title << "_e";      }
        else if ( pdg::IsTau(leppdg)      ) { title << "_tau";    }
        else if ( pdg::IsPion(leppdg)     ) { title << "_had";    }
    }

    const Spline * spl = evg_driver.XSecSpline(interaction);
    for(int i=0; i<kNSplineP; i++) {
      xs[i] = spl->Evaluate(e[i]) * (1E+38/units::cm2);
    }

    TGraph * gr = new TGraph(kNSplineP, e, xs);
    gr->SetName(title.str().c_str());
    FormatXSecGraph(gr);
    gr->SetTitle(spl->GetName());

    topdir->Add(gr);
  }


  //
  // totals for (anti-)neutrino scattering
  //

  bool is_neutrino = pdg::IsNeutralLepton(gOptProbePdgCode);

  if(is_neutrino) {

    //
    // add-up all res channels
    //

    double * xsresccp = new double[kNSplineP];
    double * xsresccn = new double[kNSplineP];
    double * xsresncp = new double[kNSplineP];
    double * xsresncn = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
       xsresccp[i] = 0;
       xsresccn[i] = 0;
       xsresncp[i] = 0;
       xsresncn[i] = 0;
    }

    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
       const Interaction * interaction = *ilistiter;
       const ProcessInfo &  proc = interaction->ProcInfo();
       const InitialState & init = interaction->InitState();
       const Target &       tgt  = init.Tgt();

       const Spline * spl = evg_driver.XSecSpline(interaction);

       if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsresccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsresccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsresncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsresncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }

    TGraph * gr_resccp = new TGraph(kNSplineP, e, xsresccp);
    gr_resccp->SetName("res_cc_p");
    FormatXSecGraph(gr_resccp);
    topdir->Add(gr_resccp);
    TGraph * gr_resccn = new TGraph(kNSplineP, e, xsresccn);
    gr_resccn->SetName("res_cc_n");
    FormatXSecGraph(gr_resccn);
    topdir->Add(gr_resccn);
    TGraph * gr_resncp = new TGraph(kNSplineP, e, xsresncp);
    gr_resncp->SetName("res_nc_p");
    FormatXSecGraph(gr_resncp);
    topdir->Add(gr_resncp);
    TGraph * gr_resncn = new TGraph(kNSplineP, e, xsresncn);
    gr_resncn->SetName("res_nc_n");
    FormatXSecGraph(gr_resncn);
    topdir->Add(gr_resncn);

    //
    // add-up all dis channels
    //

    double * xsdisccp = new double[kNSplineP];
    double * xsdisccn = new double[kNSplineP];
    double * xsdisncp = new double[kNSplineP];
    double * xsdisncn = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
       xsdisccp[i] = 0;
       xsdisccn[i] = 0;
       xsdisncp[i] = 0;
       xsdisncn[i] = 0;
    }
    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
       const Interaction * interaction = *ilistiter;
       const ProcessInfo &  proc = interaction->ProcInfo();
       const XclsTag &      xcls = interaction->ExclTag();
       const InitialState & init = interaction->InitState();
       const Target &       tgt  = init.Tgt();

       const Spline * spl = evg_driver.XSecSpline(interaction);

       if(xcls.IsCharmEvent()) continue;

       if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }
    TGraph * gr_disccp = new TGraph(kNSplineP, e, xsdisccp);
    gr_disccp->SetName("dis_cc_p");
    FormatXSecGraph(gr_disccp);
    topdir->Add(gr_disccp);
    TGraph * gr_disccn = new TGraph(kNSplineP, e, xsdisccn);
    gr_disccn->SetName("dis_cc_n");
    FormatXSecGraph(gr_disccn);
    topdir->Add(gr_disccn);
    TGraph * gr_disncp = new TGraph(kNSplineP, e, xsdisncp);
    gr_disncp->SetName("dis_nc_p");
    FormatXSecGraph(gr_disncp);
    topdir->Add(gr_disncp);
    TGraph * gr_disncn = new TGraph(kNSplineP, e, xsdisncn);
    gr_disncn->SetName("dis_nc_n");
    FormatXSecGraph(gr_disncn);
    topdir->Add(gr_disncn);

    //
    // add-up all charm dis channels
    //

    for(int i=0; i<kNSplineP; i++) {
      xsdisccp[i] = 0;
      xsdisccn[i] = 0;
      xsdisncp[i] = 0;
      xsdisncn[i] = 0;
    }
    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
      const Interaction * interaction = *ilistiter;
      const ProcessInfo &  proc = interaction->ProcInfo();
      const XclsTag &      xcls = interaction->ExclTag();
      const InitialState & init = interaction->InitState();
      const Target &       tgt  = init.Tgt();

      const Spline * spl = evg_driver.XSecSpline(interaction);

      if(!xcls.IsCharmEvent()) continue;

      if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
        for(int i=0; i<kNSplineP; i++) {
            xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
        for(int i=0; i<kNSplineP; i++) {
            xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
        for(int i=0; i<kNSplineP; i++) {
            xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
        for(int i=0; i<kNSplineP; i++) {
            xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
    }
    TGraph * gr_disccp_charm = new TGraph(kNSplineP, e, xsdisccp);
    gr_disccp_charm->SetName("dis_cc_p_charm");
    FormatXSecGraph(gr_disccp_charm);
    topdir->Add(gr_disccp_charm);
    TGraph * gr_disccn_charm = new TGraph(kNSplineP, e, xsdisccn);
    gr_disccn_charm->SetName("dis_cc_n_charm");
    FormatXSecGraph(gr_disccn_charm);
    topdir->Add(gr_disccn_charm);
    TGraph * gr_disncp_charm = new TGraph(kNSplineP, e, xsdisncp);
    gr_disncp_charm->SetName("dis_nc_p_charm");
    FormatXSecGraph(gr_disncp_charm);
    topdir->Add(gr_disncp_charm);
    TGraph * gr_disncn_charm = new TGraph(kNSplineP, e, xsdisncn);
    gr_disncn_charm->SetName("dis_nc_n_charm");
    FormatXSecGraph(gr_disncn_charm);
    topdir->Add(gr_disncn_charm);

    //
    // add-up all mec channels
    //

    double * xsmeccc = new double[kNSplineP];
    double * xsmecnc = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
       xsmeccc[i] = 0;
       xsmecnc[i] = 0;
    }

    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
       const Interaction * interaction = *ilistiter;
       const ProcessInfo &  proc = interaction->ProcInfo();

       const Spline * spl = evg_driver.XSecSpline(interaction);

       if (proc.IsMEC() && proc.IsWeakCC()) {
         for(int i=0; i<kNSplineP; i++) {
             xsmeccc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsMEC() && proc.IsWeakNC()) {
         for(int i=0; i<kNSplineP; i++) {
             xsmecnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }

    TGraph * gr_meccc = new TGraph(kNSplineP, e, xsmeccc);
    gr_meccc->SetName("mec_cc");
    FormatXSecGraph(gr_meccc);
    topdir->Add(gr_meccc);
    TGraph * gr_mecnc = new TGraph(kNSplineP, e, xsmecnc);
    gr_mecnc->SetName("mec_nc");
    FormatXSecGraph(gr_mecnc);
    topdir->Add(gr_mecnc);

    //
    // add-up all COH channels
    //

    double * xscohcc = new double[kNSplineP];
    double * xscohnc = new double[kNSplineP];
    double * xscohtot = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
      xscohcc[i] = 0;
      xscohnc[i] = 0;
      xscohtot[i] = 0;
    }

    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {

      const Interaction * interaction = *ilistiter;
      const ProcessInfo &  proc = interaction->ProcInfo();
      
      const Spline * spl = evg_driver.XSecSpline(interaction);
      
      if (proc.IsCoherentProduction() && proc.IsWeakCC()) {
        for(int i=0; i<kNSplineP; i++) {
          xscohcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (proc.IsCoherentProduction() && proc.IsWeakNC()) {
        for(int i=0; i<kNSplineP; i++) {
          xscohnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if ( proc.IsCoherentProduction() ) {
        for(int i=0; i<kNSplineP; i++) {
          xscohtot[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }

    }

    TGraph * gr_cohcc = new TGraph(kNSplineP, e, xscohcc);
    gr_cohcc->SetName("coh_cc");
    FormatXSecGraph(gr_cohcc);
    topdir->Add(gr_cohcc);

    TGraph * gr_cohnc = new TGraph(kNSplineP, e, xscohnc);
    gr_cohnc->SetName("coh_nc");
    FormatXSecGraph(gr_cohnc);
    topdir->Add(gr_cohnc);

    TGraph * gr_cohtot = new TGraph(kNSplineP, e, xscohtot);
    gr_cohtot->SetName("coh");
    FormatXSecGraph(gr_cohtot);
    topdir->Add(gr_cohtot);

    //
    // total cross sections
    //
    double * xstotcc  = new double[kNSplineP];
    double * xstotccp = new double[kNSplineP];
    double * xstotccn = new double[kNSplineP];
    double * xstotnc  = new double[kNSplineP];
    double * xstotncp = new double[kNSplineP];
    double * xstotncn = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
      xstotcc [i] = 0;
      xstotccp[i] = 0;
      xstotccn[i] = 0;
      xstotnc [i] = 0;
      xstotncp[i] = 0;
      xstotncn[i] = 0;
    }
    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
      const Interaction * interaction = *ilistiter;
      const ProcessInfo &  proc = interaction->ProcInfo();
      const InitialState & init = interaction->InitState();
      const Target &       tgt  = init.Tgt();

      const Spline * spl = evg_driver.XSecSpline(interaction);

      bool iscc = proc.IsWeakCC();
      bool isnc = proc.IsWeakNC();
      bool offp = pdg::IsProton (tgt.HitNucPdg());
      bool offn = pdg::IsNeutron(tgt.HitNucPdg());

      if (iscc && offp) {
        for(int i=0; i<kNSplineP; i++) {
            xstotccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (iscc && offn) {
        for(int i=0; i<kNSplineP; i++) {
            xstotccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (isnc && offp) {
        for(int i=0; i<kNSplineP; i++) {
            xstotncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (isnc && offn) {
        for(int i=0; i<kNSplineP; i++) {
            xstotncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }

      if (iscc) {
        for(int i=0; i<kNSplineP; i++) {
            xstotcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (isnc) {
        for(int i=0; i<kNSplineP; i++) {
            xstotnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
    }

    TGraph * gr_totcc = new TGraph(kNSplineP, e, xstotcc);
    gr_totcc->SetName("tot_cc");
    FormatXSecGraph(gr_totcc);
    topdir->Add(gr_totcc);
    TGraph * gr_totccp = new TGraph(kNSplineP, e, xstotccp);
    gr_totccp->SetName("tot_cc_p");
    FormatXSecGraph(gr_totccp);
    topdir->Add(gr_totccp);
    TGraph * gr_totccn = new TGraph(kNSplineP, e, xstotccn);
    gr_totccn->SetName("tot_cc_n");
    FormatXSecGraph(gr_totccn);
    topdir->Add(gr_totccn);
    TGraph * gr_totnc = new TGraph(kNSplineP, e, xstotnc);
    gr_totnc->SetName("tot_nc");
    FormatXSecGraph(gr_totnc);
    topdir->Add(gr_totnc);
    TGraph * gr_totncp = new TGraph(kNSplineP, e, xstotncp);
    gr_totncp->SetName("tot_nc_p");
    FormatXSecGraph(gr_totncp);
    topdir->Add(gr_totncp);
    TGraph * gr_totncn = new TGraph(kNSplineP, e, xstotncn);
    gr_totncn->SetName("tot_nc_n");
    FormatXSecGraph(gr_totncn);
    topdir->Add(gr_totncn);

    delete [] e;
    delete [] xs;
    delete [] xsresccp;
    delete [] xsresccn;
    delete [] xsresncp;
    delete [] xsresncn;
    delete [] xsdisccp;
    delete [] xsdisccn;
    delete [] xsdisncp;
    delete [] xsdisncn;
    delete [] xscohcc;
    delete [] xscohnc;
    delete [] xscohtot;
    delete [] xstotcc;
    delete [] xstotccp;
    delete [] xstotccn;
    delete [] xstotnc;
    delete [] xstotncp;
    delete [] xstotncn;

  }// neutrinos


  //
  // totals for charged lepton scattering
  //

  bool is_charged_lepton = pdg::IsChargedLepton(gOptProbePdgCode);

  if(is_charged_lepton) {

    //
    // add-up all res channels
    //

    double * xsresemp = new double[kNSplineP];
    double * xsresemn = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
       xsresemp[i] = 0;
       xsresemn[i] = 0;
    }

    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
       const Interaction * interaction = *ilistiter;
       const ProcessInfo &  proc = interaction->ProcInfo();
       const InitialState & init = interaction->InitState();
       const Target &       tgt  = init.Tgt();

       const Spline * spl = evg_driver.XSecSpline(interaction);

       if (proc.IsResonant() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsresemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsResonant() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsresemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }

    TGraph * gr_resemp = new TGraph(kNSplineP, e, xsresemp);
    gr_resemp->SetName("res_em_p");
    FormatXSecGraph(gr_resemp);
    topdir->Add(gr_resemp);
    TGraph * gr_resemn = new TGraph(kNSplineP, e, xsresemn);
    gr_resemn->SetName("res_em_n");
    FormatXSecGraph(gr_resemn);
    topdir->Add(gr_resemn);

    //
    // add-up all dis channels
    //

    double * xsdisemp = new double[kNSplineP];
    double * xsdisemn = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
       xsdisemp[i] = 0;
       xsdisemn[i] = 0;
    }
    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
       const Interaction * interaction = *ilistiter;
       const ProcessInfo &  proc = interaction->ProcInfo();
       const XclsTag &      xcls = interaction->ExclTag();
       const InitialState & init = interaction->InitState();
       const Target &       tgt  = init.Tgt();

       const Spline * spl = evg_driver.XSecSpline(interaction);

       if(xcls.IsCharmEvent()) continue;

       if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
       if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
         for(int i=0; i<kNSplineP; i++) {
             xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }
    TGraph * gr_disemp = new TGraph(kNSplineP, e, xsdisemp);
    gr_disemp->SetName("dis_em_p");
    FormatXSecGraph(gr_disemp);
    topdir->Add(gr_disemp);
    TGraph * gr_disemn = new TGraph(kNSplineP, e, xsdisemn);
    gr_disemn->SetName("dis_em_n");
    FormatXSecGraph(gr_disemn);
    topdir->Add(gr_disemn);

    //
    // add-up all charm dis channels
    //

    for(int i=0; i<kNSplineP; i++) {
      xsdisemp[i] = 0;
      xsdisemn[i] = 0;
    }
    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
      const Interaction * interaction = *ilistiter;
      const ProcessInfo &  proc = interaction->ProcInfo();
      const XclsTag &      xcls = interaction->ExclTag();
      const InitialState & init = interaction->InitState();
      const Target &       tgt  = init.Tgt();

      const Spline * spl = evg_driver.XSecSpline(interaction);

      if(!xcls.IsCharmEvent()) continue;

      if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
        for(int i=0; i<kNSplineP; i++) {
            xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
        for(int i=0; i<kNSplineP; i++) {
            xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
    }
    TGraph * gr_disemp_charm = new TGraph(kNSplineP, e, xsdisemp);
    gr_disemp_charm->SetName("dis_em_p_charm");
    FormatXSecGraph(gr_disemp_charm);
    topdir->Add(gr_disemp_charm);
    TGraph * gr_disemn_charm = new TGraph(kNSplineP, e, xsdisemn);
    gr_disemn_charm->SetName("dis_em_n_charm");
    FormatXSecGraph(gr_disemn_charm);
    topdir->Add(gr_disemn_charm);

    //
    // total cross sections
    //
    double * xstotem  = new double[kNSplineP];
    double * xstotemp = new double[kNSplineP];
    double * xstotemn = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
      xstotem [i] = 0;
      xstotemp[i] = 0;
      xstotemn[i] = 0;
    }
    for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
      const Interaction * interaction = *ilistiter;
      const ProcessInfo &  proc = interaction->ProcInfo();
      const InitialState & init = interaction->InitState();
      const Target &       tgt  = init.Tgt();

      const Spline * spl = evg_driver.XSecSpline(interaction);

      bool isem = proc.IsEM();
      bool offp = pdg::IsProton (tgt.HitNucPdg());
      bool offn = pdg::IsNeutron(tgt.HitNucPdg());

      if (isem && offp) {
        for(int i=0; i<kNSplineP; i++) {
            xstotemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (isem && offn) {
        for(int i=0; i<kNSplineP; i++) {
            xstotemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
      if (isem) {
        for(int i=0; i<kNSplineP; i++) {
            xstotem[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
        }
      }
    }

    TGraph * gr_totem = new TGraph(kNSplineP, e, xstotem);
    gr_totem->SetName("tot_em");
    FormatXSecGraph(gr_totem);
    topdir->Add(gr_totem);
    TGraph * gr_totemp = new TGraph(kNSplineP, e, xstotemp);
    gr_totemp->SetName("tot_em_p");
    FormatXSecGraph(gr_totemp);
    topdir->Add(gr_totemp);
    TGraph * gr_totemn = new TGraph(kNSplineP, e, xstotemn);
    gr_totemn->SetName("tot_em_n");
    FormatXSecGraph(gr_totemn);
    topdir->Add(gr_totemn);

    delete [] e;
    delete [] xs;
    delete [] xsresemp;
    delete [] xsresemn;
    delete [] xsdisemp;
    delete [] xsdisemn;
    delete [] xstotem;
    delete [] xstotemp;
    delete [] xstotemn;

  }// charged leptons

  topdir->Write();

  if(froot) {
    froot->Close();
    delete froot;
  }
}

void SaveNtupleToRootFile

(

void

)

Definition at line 1317 of file gSplineXml2Root.cxx.

{
/*
  //-- create cross section ntuple
  //
  double brXSec;
  double brEv;
  bool   brIsQel;
  bool   brIsRes;
  bool   brIsDis;
  bool   brIsCoh;
  bool   brIsImd;
  bool   brIsNuEEl;
  bool   brIsCC;
  bool   brIsNC;
  int    brNucleon;
  int    brQrk;
  bool   brIsSeaQrk;
  int    brRes;
  bool   brCharm;
  TTree * xsecnt = new TTree("xsecnt",dtitle.str().c_str());
  xsecnt->Branch("xsec",  &brXSec,     "xsec/D");
  xsecnt->Branch("e",     &brEv,       "e/D");
  xsecnt->Branch("qel",   &brIsQel,    "qel/O");
  xsecnt->Branch("res",   &brIsRes,    "res/O");
  xsecnt->Branch("dis",   &brIsDis,    "dis/O");
  xsecnt->Branch("coh",   &brIsCoh,    "coh/O");
  xsecnt->Branch("imd",   &brIsImd,    "imd/O");
  xsecnt->Branch("veel",  &brIsNuEEl,  "veel/O");
  xsecnt->Branch("cc",    &brIsCC,     "cc/O");
  xsecnt->Branch("nc",    &brIsNC,     "nc/O");
  xsecnt->Branch("nuc",   &brNucleon,  "nuc/I");
  xsecnt->Branch("qrk",   &brQrk,      "qrk/I");
  xsecnt->Branch("sea",   &brIsSeaQrk, "sea/O");
  xsecnt->Branch("res",   &brRes,      "res/I");
  xsecnt->Branch("charm", &brCharm,    "charm/O");
*/
}

void SaveToPsFile

(

void

)

Definition at line 217 of file gSplineXml2Root.cxx.

{
  if(!gWriteOutPlots) return;

  //-- get the event generation driver
  GEVGDriver evg_driver = GetEventGenDriver();

  //-- define some marker styles / colors
  const unsigned int kNMarkers = 5;
  const unsigned int kNColors  = 6;
  unsigned int markers[kNMarkers] = {20, 28, 29, 27, 3};
  unsigned int colors [kNColors]  = {1, 2, 4, 6, 8, 28};

  //-- create a postscript document for saving cross section plpots

  TCanvas * c = new TCanvas("c","",20,20,500,850);
  c->SetBorderMode(0);
  c->SetFillColor(0);
  TLegend * legend = new TLegend(0.01,0.01,0.99,0.99);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);

  //-- get pdglibrary for mapping pdg codes to names
  PDGLibrary * pdglib = PDGLibrary::Instance();
  ostringstream filename;
  filename << "xsec-splines-"
          <<  pdglib->Find(gOptProbePdgCode)->GetName()  << "-"
          <<  pdglib->Find(gOptTgtPdgCode)->GetName() << ".ps";
  TPostScript * ps = new TPostScript(filename.str().c_str(), kPsType);

  //-- get the list of interactions that can be simulated by the driver
  const InteractionList * ilist = evg_driver.Interactions();
  unsigned int nspl = ilist->size();

  //-- book enough space for xsec plots (last one is the sum)
  TGraph * gr[nspl+1];

  //-- loop over all the simulated interactions & create the cross section graphs
  InteractionList::const_iterator ilistiter = ilist->begin();
  unsigned int i=0;
  for(; ilistiter != ilist->end(); ++ilistiter) {

    const Interaction * interaction = *ilistiter;
    LOG("gspl2root", pINFO)
       << "Current interaction: " << interaction->AsString();


    //-- access the cross section spline
    const Spline * spl = evg_driver.XSecSpline(interaction);
    if(!spl) {
      LOG("gspl2root", pWARN)
         << "Can't get spline for: " << interaction->AsString();
      exit(2);
    }

    //-- set graph color/style
    int icol = TMath::Min( i % kNColors,  kNColors-1  );
    int isty = TMath::Min( i / kNMarkers, kNMarkers-1 );
    int col = colors[icol];
    int sty = markers[isty];

    LOG("gspl2root", pINFO)
      << "color = " << col << ", marker = " << sty;

    //-- export Spline as TGraph / set color & stype
    gr[i] = spl->GetAsTGraph(kNP,true,true,1.,1./units::cm2);
    gr[i]->SetLineColor(col);
    gr[i]->SetMarkerColor(col);
    gr[i]->SetMarkerStyle(sty);
    gr[i]->SetMarkerSize(0.5);

    i++;
  }

  //-- now, get the sum...
  const Spline * splsum = evg_driver.XSecSumSpline();
  if(!splsum) {
      LOG("gspl2root", pWARN)
          << "Can't get the cross section sum spline";
      exit(2);
  }
  gr[nspl] = splsum->GetAsTGraph(kNP,true,true,1.,1./units::cm2);

  //-- figure out the minimum / maximum xsec in plotted range
  double XSmax = -9999;
  double XSmin =  9999;
  double x=0, y=0;
  for(int j=0; j<kNP; j++) {
    gr[nspl]->GetPoint(j,x,y);
    XSmax = TMath::Max(XSmax,y);
  }
  XSmin = XSmax/100000.;
  XSmax = XSmax*1.2;

  LOG("gspl2root", pINFO) << "Drawing frame: E    = (" << gEmin  << ", " << gEmax << ")";
  LOG("gspl2root", pINFO) << "Drawing frame: XSec = (" << XSmin  << ", " << XSmax << ")";

  //-- ps output: add the 1st page with _all_ xsec spline plots
  //
  //c->Draw();
  TH1F * h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  for(unsigned int ispl = 0; ispl <= nspl; ispl++) if(gr[ispl]) { gr[ispl]->Draw("LP"); }
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();

  //-- plot QEL xsecs only
  //
  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  i=0;
  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
    const Interaction * interaction = *ilistiter;
    if(interaction->ProcInfo().IsQuasiElastic() && ! interaction->ExclTag().IsCharmEvent()) {
        gr[i]->Draw("LP");
        TString spltitle(interaction->AsString());
        spltitle = spltitle.ReplaceAll(";",1," ",1);
        legend->AddEntry(gr[i], spltitle.Data(),"LP");
    }
    i++;
  }
  legend->SetHeader("QEL Cross Sections");
  gr[nspl]->Draw("LP");
  legend->AddEntry(gr[nspl], "sum","LP");
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();
  c->Clear();
  c->Range(0,0,1,1);
  legend->Draw();
  c->Update();

  //-- plot RES xsecs only
  //
  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  legend->Clear();
  i=0;
  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
    const Interaction * interaction = *ilistiter;
    if(interaction->ProcInfo().IsResonant()) {
        gr[i]->Draw("LP");
        TString spltitle(interaction->AsString());
        spltitle = spltitle.ReplaceAll(";",1," ",1);
        legend->AddEntry(gr[i], spltitle.Data(),"LP");
    }
    i++;
  }
  legend->SetHeader("RES Cross Sections");
  gr[nspl]->Draw("LP");
  legend->AddEntry(gr[nspl], "sum","LP");
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();
  c->Clear();
  c->Range(0,0,1,1);
  legend->Draw();
  c->Update();

  //-- plot DIS xsecs only
  //
  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  legend->Clear();
  i=0;
  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
    const Interaction * interaction = *ilistiter;
    if(interaction->ProcInfo().IsDeepInelastic()) {
        gr[i]->Draw("LP");
        TString spltitle(interaction->AsString());
        spltitle = spltitle.ReplaceAll(";",1," ",1);
        legend->AddEntry(gr[i], spltitle.Data(),"LP");
    }
    i++;
  }
  legend->SetHeader("DIS Cross Sections");
  gr[nspl]->Draw("LP");
  legend->AddEntry(gr[nspl], "sum","LP");
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();
  c->Clear();
  c->Range(0,0,1,1);
  legend->Draw();
  c->Update();

  //-- plot COH xsecs only
  //
  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  legend->Clear();
  i=0;
  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
    const Interaction * interaction = *ilistiter;
    if(interaction->ProcInfo().IsCoherentProduction()) {
        gr[i]->Draw("LP");
        TString spltitle(interaction->AsString());
        spltitle = spltitle.ReplaceAll(";",1," ",1);
        legend->AddEntry(gr[i], spltitle.Data(),"LP");
    }
    i++;
  }
  legend->SetHeader("COH Cross Sections");
  gr[nspl]->Draw("LP");
  legend->AddEntry(gr[nspl], "sum","LP");
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();
  c->Clear();
  c->Range(0,0,1,1);
  legend->Draw();
  c->Update();

  //-- plot charm xsecs only
  //
  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  i=0;
  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
    const Interaction * interaction = *ilistiter;
    if(interaction->ExclTag().IsCharmEvent()) {
        gr[i]->Draw("LP");
        TString spltitle(interaction->AsString());
        spltitle = spltitle.ReplaceAll(";",1," ",1);
        legend->AddEntry(gr[i], spltitle.Data(),"LP");
    }
    i++;
  }
  legend->SetHeader("Charm Prod. Cross Sections");
  //gr[nspl]->Draw("LP");
  legend->AddEntry(gr[nspl], "sum","LP");
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();
  c->Clear();
  c->Range(0,0,1,1);
  legend->Draw();
  c->Update();

  //-- plot ve xsecs only
  //
  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
  legend->Clear();
  i=0;
  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
    const Interaction * interaction = *ilistiter;
    if(interaction->ProcInfo().IsInverseMuDecay() ||
       interaction->ProcInfo().IsIMDAnnihilation() ||
       interaction->ProcInfo().IsNuElectronElastic()) {
        gr[i]->Draw("LP");
        TString spltitle(interaction->AsString());
        spltitle = spltitle.ReplaceAll(";",1," ",1);
        legend->AddEntry(gr[i], spltitle.Data(),"LP");
    }
    i++;
  }
  legend->SetHeader("IMD and ve Elastic Cross Sections");
  gr[nspl]->Draw("LP");
  legend->AddEntry(gr[nspl], "sum","LP");
  h->GetXaxis()->SetTitle("Ev (GeV)");
  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
  c->SetLogx();
  c->SetLogy();
  c->SetGridx();
  c->SetGridy();
  c->Update();
  c->Clear();
  c->Range(0,0,1,1);
  legend->Draw();
  c->Update();

  //-- close the postscript document
  ps->Close();

  //-- clean-up
  for(unsigned int j=0; j<=nspl; j++) { if(gr[j]) delete gr[j]; }
  delete c;
  delete ps;
}

Variable Documentation

double gEmax

Definition at line 151 of file gSplineXml2Root.cxx.

double gEmin

Definition at line 150 of file gSplineXml2Root.cxx.

bool gInlogE

Definition at line 152 of file gSplineXml2Root.cxx.

int gOptProbePdgCode

Definition at line 144 of file gSplineXml2Root.cxx.

PDGCodeList gOptProbePdgList

Definition at line 142 of file gSplineXml2Root.cxx.

string gOptROOTFilename

Definition at line 141 of file gSplineXml2Root.cxx.

int gOptTgtPdgCode

Definition at line 145 of file gSplineXml2Root.cxx.

PDGCodeList gOptTgtPdgList

Definition at line 143 of file gSplineXml2Root.cxx.

string gOptXMLFilename

Definition at line 140 of file gSplineXml2Root.cxx.

bool gWriteOutPlots

Definition at line 146 of file gSplineXml2Root.cxx.

int kNP = 300

Definition at line 153 of file gSplineXml2Root.cxx.

int kNSplineP = 1000

Definition at line 154 of file gSplineXml2Root.cxx.

const int kPsType = 111

Definition at line 155 of file gSplineXml2Root.cxx.