gSplineXml2Root.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*!

\program gspl2root

\brief   Utility to load a GENIE XML cross section spline file and convert
         splines into a ROOT format.

         Syntax :
           gspl2root -f xml_file -p nu -t tgt [-e emax]
                     [-o root_file] [-w] [-k] [-l]
                     [--message-thresholds xml_file]
                     [--event-generator-list list_name]

         Options :
           []  denotes an optional argument

           -f
              the input XML file containing the cross section spline data
           -p
              the neutrino pdg code
           -t
              the target pdg code (format: 10LZZZAAAI)
           -e
              the minimum and maximum energy (in generated plots -- use it to zoom at low E)
           -o
              output ROOT file name
           -w
              write out plots in a postscipt file
           -l
              energy bins in log10 scale
           -k
              keep spline knot points  (not yet implemented).
           --message-thresholds
              Allows users to customize the message stream thresholds.
           --event-generator-list
              List of event generators to load in event generation drivers.
              [default: "Default"].

         Example:

           Extract all numu+n, numu+p, numu+O16 cross section splines from the
           input XML file `mysplines.xml', convert splines into a ROOT format
           and save them into a single ROOT file.

           shell$ gspl2root -f mysplines.xml -p 14 -t 1000000010 -o xsec.root
           shell$ gspl2root -f mysplines.xml -p 14 -t 1000010010 -o xsec.root
           shell$ gspl2root -f mysplines.xml -p 14 -t 1000080160 -o xsec.root

           A large number of graphs (one per simulated process + appropriate
           totals) will be generated in each case. Each set of plots is saved
           into its own ROOT TDirectory named after the specified initial state.

         Important Note:

           The stored graphs can be used for cross section interpolation.
           Essentially, this _single_ ROOT file can contain _all_ the GENIE cross
           section `functions' you may need.
           For instance, the `xsec.root' file generated in the above example will
           contain a `nu_mu_O16' directory (generated by the last command)
           which will include cross section graphs for all numu + O16 processes.
           To extract the numu+O16 DIS CC cross section graph for hit u valence
           quarks in a bound proton and evaluate the cross section at energy E,
           then type:

           root[0] TDirectory * dir = (TDirectory*) file->Get("nu_mu_O16");
           root[1] TGraph * graph = (TGraph*) dir->Get("dis_cc_p_uval");
           root[2] cout << graph->Evaluate(E) << endl;

           See the User Manual for more details.

\author  Costas Andreopoulos <c.andreopoulos \at cern.ch>
 University of Liverpool

\created December 15, 2005

\cpright Copyright (c) 2003-2024, The GENIE Collaboration
         For the full text of the license visit http://copyright.genie-mc.org
*/
//____________________________________________________________________________

#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"


using std::string;
using std::vector;
using std::ostringstream;

using namespace genie;
using namespace genie::utils;

//Prototypes:
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);

//User-specified options:
string gOptXMLFilename;  // input XML filename
string gOptROOTFilename; // output ROOT filename
PDGCodeList gOptProbePdgList;  // list of probe PDG codes
PDGCodeList gOptTgtPdgList;    // list of target PDG codes
int    gOptProbePdgCode; // probe PDG code (currently being processed)
int    gOptTgtPdgCode;   // target PDG code
bool   gWriteOutPlots;   // write out a postscript file with plots
//bool   gKeepSplineKnots; // use spline abscissa points rather than equi-spaced

//Globals & constants
double gEmin;
double gEmax;
bool gInlogE;
int    kNP       = 300;
int    kNSplineP = 1000;
const int    kPsType   = 111;  // ps type: portrait

//____________________________________________________________________________
int main(int argc, char ** argv)
{
  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 LoadSplines(void)
{
// load the cross section splines specified at the cmd line

  XSecSplineList * splist = XSecSplineList::Instance();
  XmlParserStatus_t ist = splist->LoadFromXml(gOptXMLFilename);
  assert(ist == kXmlOK);
}
//____________________________________________________________________________
GEVGDriver GetEventGenDriver(void)
{
// 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;
}
//____________________________________________________________________________
void SaveToPsFile(void)
{
  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;
}
//____________________________________________________________________________
void FormatXSecGraph(TGraph * g)
{
  g->SetTitle("GENIE cross section graph");
  g->GetXaxis()->SetTitle("Ev (GeV)");
  g->GetYaxis()->SetTitle("#sigma_{nuclear} (10^{-38} cm^{2})");
}
//____________________________________________________________________________
void SaveGraphsToRootFile(void)
{
  //-- 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";   
      if ( ! xcls.KnownResonance() ) 
        title << "_single_pion" ;
    }
    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 if (proc.IsPhotonResonance() )  { title << "phres"; }
    else if (proc.IsPhotonCoherent() )   { title << "phcoh"; }
    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()) {
      if ( xcls.KnownResonance() ) {
        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);
      }
      else if ( xcls.NPions() == 1 ) {
        // we are in the case of the single pion production
        // since the hiht nucleon is known and we also know if 
        // it is CC or NC, the only missing information to identify the channel
        // is only the flavour of the pion
        string channel ;
        if ( xcls.NPiPlus() == 1 ) {
          channel = "pip" ;
        }
        else if ( xcls.NPiMinus() == 1 ) {
          channel = "pim" ;
        }
        else if ( xcls.NPi0() == 1 ) {
          channel = "pi0" ;
        }
        
        title << '_' << channel ;
      }  // single resonsance or single pion
    }  // resonance case

    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 = TMath::Abs(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 * xsdiscc  = new double[kNSplineP];
    double * xsdisnc  = new double[kNSplineP];
    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++) {
       xsdiscc[i]  = 0;
       xsdisnc[i]  = 0;
       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++) {
             xsdiscc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             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++) {
             xsdiscc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             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++) {
             xsdisnc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             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++) {
             xsdisnc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }
    TGraph * gr_discc = new TGraph(kNSplineP, e, xsdiscc);
    gr_discc->SetName("dis_cc");
    FormatXSecGraph(gr_discc);
    topdir->Add(gr_discc);
    TGraph * gr_disnc = new TGraph(kNSplineP, e, xsdisnc);
    gr_disnc->SetName("dis_nc");
    FormatXSecGraph(gr_disnc);
    topdir->Add(gr_disnc);
    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++) {
       xsdiscc[i]  = 0;
       xsdisnc[i]  = 0;
       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++) {
             xsdiscc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             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++) {
             xsdiscc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             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++) {
             xsdisnc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             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++) {
             xsdisnc[i]  += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
             xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
         }
       }
    }
    TGraph * gr_discc_charm = new TGraph(kNSplineP, e, xsdiscc);
    gr_discc_charm->SetName("dis_cc_charm");
    FormatXSecGraph(gr_discc_charm);
    topdir->Add(gr_discc_charm);
    TGraph * gr_disnc_charm = new TGraph(kNSplineP, e, xsdisnc);
    gr_disnc_charm->SetName("dis_nc_charm");
    FormatXSecGraph(gr_disnc_charm);
    topdir->Add(gr_disnc_charm);
    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);

    //
    // add-up all glres and photon-res channels
    //

    double * xsglrescc = new double[kNSplineP];
    double * xsglresnc = new double[kNSplineP];
    double * xsphrescc = new double[kNSplineP];
    for(int i=0; i<kNSplineP; i++) {
       xsglrescc[i] = 0;
       xsglresnc[i] = 0;
       xsphrescc[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.IsGlashowResonance() && proc.IsWeakCC()) {
         for(int i=0; i<kNSplineP; i++) { 
             xsglrescc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); 
         }        
       } 
       if (proc.IsGlashowResonance() && proc.IsWeakNC()) {
         for(int i=0; i<kNSplineP; i++) { 
             xsglresnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); 
         }        
       } 
       if (proc.IsPhotonResonance()) {
         for(int i=0; i<kNSplineP; i++) { 
             xsphrescc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); 
         }    
             
       } 
    }
    TGraph * gr_glrescc = new TGraph(kNSplineP, e, xsglrescc);
    gr_glrescc->SetName("glres_cc");
    FormatXSecGraph(gr_glrescc);
    topdir->Add(gr_glrescc);
    TGraph * gr_glresnc = new TGraph(kNSplineP, e, xsglresnc);
    gr_glresnc->SetName("glres_nc");
    FormatXSecGraph(gr_glresnc);
    topdir->Add(gr_glresnc);
    TGraph * gr_phrescc = new TGraph(kNSplineP, e, xsphrescc);
    gr_phrescc->SetName("phres_cc");
    FormatXSecGraph(gr_phrescc);
    topdir->Add(gr_phrescc);


    //
    // 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 [] xsdiscc;
    delete [] xsdisnc;
    delete [] xscohcc;
    delete [] xscohnc;
    delete [] xscohtot;
    delete [] xsglrescc;
    delete [] xsglresnc;
    delete [] xsphrescc;
    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)
{
/*
  //-- 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 GetCommandLineArgs(int argc, char ** argv)
{
  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");
}
//____________________________________________________________________________
void PrintSyntax(void)
{
  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";
}
//____________________________________________________________________________
PDGCodeList GetPDGCodeListFromString(std::string s)
{
  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;

}
//____________________________________________________________________________