gEvComp.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*!

\program gevcomp

\brief   A GENIE neutrino event sample comparison utility.
         Reads-in a GENIE GHEP event tree and generates a postscript file
         with MC truth quantity plots (& comparisons with a reference event 
         sample is specified)

         Syntax :
           gevcomp -f sample [-r reference_sample]

         Options:
           [] Denotes an optional argument
           -f Specifies the GENIE/ROOT file with the generated event sample
           -r Specifies another GENIE/ROOT event sample file for comparison 
           -n Specifies how many events to analyze [default: all]

         Notes:
           The input ROOT files are the gst summary ntuples generated by 
           GENIE's gntpc utility
                      
         Example:
           gevcomp -f /path/gntp.1.gst.root -r /path/gntp.2.gst.root
                      
\author  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
 University of Liverpool & STFC Rutherford Appleton Laboratory

\created June 06, 2008 

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

#include <cassert>
#include <sstream>
#include <string>

#include <TSystem.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TTree.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <TPostScript.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TH3D.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TPavesText.h>
#include <TText.h>
#include <TStyle.h>
#include <TLegend.h>

#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/Style.h"

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

using namespace genie;

// function prototypes
void   GetCommandLineArgs   (int argc, char ** argv);
void   PrintSyntax          (void);
bool   CheckRootFilename    (string filename);
string OutputFileName       (string input_file_name);
void   CreatePlots          (string filename, string filename_ref);

// command-line arguments
string   gOptInpFile     = ""; // (-f) input GENIE event sample file
string   gOptInpFileRef  = ""; // (-r) input GENIE event sample file (reference)

//_________________________________________________________________________________
int main(int argc, char ** argv)
{
  GetCommandLineArgs(argc,argv);
  utils::style::SetDefaultStyle();
  CreatePlots(gOptInpFile,gOptInpFileRef);
  
  LOG("gevcomp", pINFO)  << "Done!";
  return 0;
}
//_________________________________________________________________________________
void CreatePlots(string inp_filename, string inp_filename_ref)
{
  TFile * fin_0 = 0;
  TFile * fin_1 = 0;
  TTree * gst_0 = 0;
  TTree * gst_1 = 0;

  if(!CheckRootFilename(inp_filename)) {
    LOG("gevcomp", pERROR) << "Input file: " << inp_filename << " doesn't exist";
    return;
  }
  fin_0 = new TFile(inp_filename.c_str(),"READ");
  gst_0 = (TTree *) fin_0->Get("gst");
  assert(gst_0);
  
  if(CheckRootFilename(inp_filename_ref)) {
     fin_1 = new TFile(inp_filename_ref.c_str(),"READ");
     gst_1 = (TTree *) fin_1->Get("gst");
     assert(gst_1);
  }
  
  gst_0->SetLineColor(kBlack);
  gst_0->SetLineWidth(3);
  if(gst_1) {
    gst_1->SetLineColor(kRed);
    gst_1->SetMarkerColor(kRed);
    gst_1->SetLineWidth(2);
    gst_1->SetMarkerStyle(20);
    gst_1->SetMarkerSize(1);
  }
  
  // Set global plot style
  //
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetHistTopMargin(0.33);
  gStyle->SetHistMinimumZero(true);
  
  // Plotting options
  //
  bool monoenergetic_sample = true;
  bool show_coh_plots       = true;
  bool show_calc_kinematics = true;
  bool show_mult_per_proc   = true;
  bool show_primary_hadsyst = true;
  
  gst_0->Draw("1","tgt>1000010010","GOFF");
  show_coh_plots = (gst_0->GetSelectedRows() > 0); 

  TCanvas * c = new TCanvas("c","",20,20,500,650);
  c->SetBorderMode(0);
  c->SetFillColor(0);
  c->SetGridx();
  c->SetGridy();

  TLegend * ls = new TLegend(0.20,0.94,0.99,0.99);
  ls->SetFillColor(0);
  ls->SetBorderSize(0);

  string ps_filename = OutputFileName(inp_filename);
  TPostScript * ps = new TPostScript(ps_filename.c_str(), 111);
  
  //
  // SECTION: PS File Header
  //

  ps->NewPage();
  c->Range(0,0,100,100);
  TPavesText hdr(10,40,90,70,3,"tr");
  hdr.AddText("GENIE Event Sample Comparisons");
  hdr.AddText(" ");
  hdr.AddText(" ");
  hdr.AddText("Event Sample:");
  //hdr.AddText(sample);
  hdr.AddText(" ");
  hdr.AddText("Notes:");
  hdr.AddText(" ");
  hdr.AddText(" ");
  hdr.Draw();
  c->Update();

  //
  // SECTION: Event Numbers
  //

  TH1F * h0num = new TH1F("h0num","", 2, 0., 2.);
  TH1F * h1num = new TH1F("h1num","", 2, 0., 2.);
  float n0=0, n1=0;
  ps->NewPage();
  c->Range(0,0,100,100);
  TPavesText evn(10,10,90,90,3,"tr");
  evn.AddText("Event Numbers:");
  evn.AddText("  ");
 
            gst_0->Draw("1>>h0num","1","goff");
  if(gst_1) gst_1->Draw("1>>h1num","1","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("ALL    : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","qel","goff");
  if(gst_1) gst_1->Draw("1>>h1num","qel","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("QEL    : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","qel&&cc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","qel&&cc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("QEL-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","qel&&nc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","qel&&nc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("QEL-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","res","goff");
  if(gst_1) gst_1->Draw("1>>h1num","res","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("RES    : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","res&&cc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","res&&cc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("RES-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","res&&nc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","res&&nc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("RES-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","dis","goff");
  if(gst_1) gst_1->Draw("1>>h1num","dis","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("DIS    : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","dis&&cc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","dis&&cc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("DIS-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","dis&&nc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","dis&&nc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("DIS-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","coh","goff");
  if(gst_1) gst_1->Draw("1>>h1num","coh","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("COH      : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","coh&&cc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","coh&&cc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("COH-CC   : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","coh&&nc","goff");
  if(gst_1) gst_1->Draw("1>>h1num","coh&&nc","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("COH-NC   : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","imd","goff");
  if(gst_1) gst_1->Draw("1>>h1num","imd","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("IMD    : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","nuel","goff");
  if(gst_1) gst_1->Draw("1>>h1num","nuel","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("NuE-EL : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","dis&&cc&&charm","goff");
  if(gst_1) gst_1->Draw("1>>h1num","dis&&cc&&charm","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("DIS-CHARM: %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

            gst_0->Draw("1>>h0num","qel&&cc&&charm","goff");
  if(gst_1) gst_1->Draw("1>>h1num","qel&&cc&&charm","goff");
  n0 =           h0num->GetEntries();
  n1 = (gst_1) ? h1num->GetEntries() : 0;
  evn.AddText( Form("QEL-CHARM: %7.0f [test sample], %7.0f [ref sample]", n0, n1) );

  evn.Draw();
  c->Update();

  if(!monoenergetic_sample) {
              gst_0->Draw("Ev","","");
    if(gst_1) gst_1->Draw("Ev","","perrsame");
    ls->Clear();
    ls->SetHeader("Neutrino Energy Spectrum");
    ls->Draw();
    c->Update();  
  }

  //
  // SECTION: Kinematics
  //
  ps->NewPage();
  c->Clear();
  c->Range(0,0,100,100);
  TPavesText hdrk(10,40,90,70,3,"tr");
  hdrk.AddText("Selected Kinematical Quantities");
  hdrk.AddText(" ");
  hdrk.Draw();
  c->Update();

  //------ selected Q2 for all events
  ps->NewPage();
  gst_0->Draw("Q2s","Q2s>0","");
  if(gst_1) gst_1->Draw("Q2s","Q2s>0","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for all events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for QEL
  ps->NewPage();
  gst_0->Draw("Q2s","qel&&!charm","");
  if(gst_1) gst_1->Draw("Q2s","qel&&!charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for QEL events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for QEL CC
  ps->NewPage();
  gst_0->Draw("Q2s","qel&&cc&&!charm","");
  if(gst_1) gst_1->Draw("Q2s","qel&&cc&&!charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for QEL CC events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for QEL NC
  ps->NewPage();
  gst_0->Draw("Q2s","qel&&nc&&!charm","");
  if(gst_1) gst_1->Draw("Q2s","qel&&nc&&!charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for QEL NC events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for RES
  ps->NewPage();
  gst_0->Draw("Q2s","res","");
  if(gst_1) gst_1->Draw("Q2s","res","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for RES events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for RES CC
  ps->NewPage();
  gst_0->Draw("Q2s","res&&cc","");
  if(gst_1) gst_1->Draw("Q2s","res&&cc","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for RES CC events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for RES NC
  ps->NewPage();
  gst_0->Draw("Q2s","res&&nc","");
  if(gst_1) gst_1->Draw("Q2s","res&&nc","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for RES NC events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for DIS
  ps->NewPage();
  gst_0->Draw("Q2s","dis","");
  if(gst_1) gst_1->Draw("Q2s","dis","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for DIS events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for DIS CC
  ps->NewPage();
  gst_0->Draw("Q2s","dis&&cc","");
  if(gst_1) gst_1->Draw("Q2s","dis&&cc","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for DIS CC events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for DIS NC
  ps->NewPage();
  gst_0->Draw("Q2s","dis&&nc","");
  if(gst_1) gst_1->Draw("Q2s","dis&&nc","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for DIS NC events");
  ls->Draw();
  c->Update();

  //------ selected Q2 for Charm/DIS
  ps->NewPage();
  gst_0->Draw("Q2s","dis&&charm","");
  if(gst_1) gst_1->Draw("Q2s","dis&&charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected Q2 for Charm/DIS events");
  ls->Draw();
  c->Update();

  if(show_coh_plots) {
     //------ selected Q2 for COH
     ps->NewPage();
     gst_0->Draw("Q2s","coh","");
     if(gst_1) gst_1->Draw("Q2s","coh","perrsame");
     ls->Clear();
     ls->SetHeader("selected Q2 for COH events");
     ls->Draw();
     c->Update();

     //------ selected Q2 for COH CC
     ps->NewPage();
     gst_0->Draw("Q2s","coh&&cc","");
     if(gst_1) gst_1->Draw("Q2s","coh&&cc","perrsame");
     ls->Clear();
     ls->SetHeader("selected Q2 for COH CC events");
     ls->Draw();
     c->Update();

     //------ selected Q2 for COH NC
     ps->NewPage();
     gst_0->Draw("Q2s","coh&&nc","");
     if(gst_1) gst_1->Draw("Q2s","coh&&nc","perrsame");
     ls->Clear();
     ls->SetHeader("selected Q2 for COH NC events");
     ls->Draw();
     c->Update();
  }
  
  //------ selected W for all events
  ps->NewPage();
  gst_0->Draw("Ws","Ws>0","");
  if(gst_1) gst_1->Draw("Ws","Ws>0","perrsame");
  ls->Clear();
  ls->SetHeader("selected W for all events");
  ls->Draw();
  c->Update();

  //------ selected W for QEL
  ps->NewPage();
  gst_0->Draw("Ws","qel&&!charm","");
  if(gst_1) gst_1->Draw("Ws","qel&&!charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected W for QEL events");
  ls->Draw();
  c->Update();

  //------ selected W for RES
  ps->NewPage();
  gst_0->Draw("Ws","res","");
  if(gst_1) gst_1->Draw("Ws","res","perrsame");
  ls->Clear();
  ls->SetHeader("selected W for RES events");
  ls->Draw();
  c->Update();

  //------ selected W for DIS
  ps->NewPage();
  gst_0->Draw("Ws","dis","");
  if(gst_1) gst_1->Draw("Ws","dis","perrsame");
  ls->Clear();
  ls->SetHeader("selected W for DIS events");
  ls->Draw();
  c->Update();

  //------ selected W for DIS CC
  ps->NewPage();
  gst_0->Draw("Ws","dis&&cc","");
  if(gst_1) gst_1->Draw("Ws","dis&&cc","perrsame");
  ls->Clear();
  ls->SetHeader("selected W for DIS CC events");
  ls->Draw();
  c->Update();

  //------ selected W for DIS NC
  ps->NewPage();
  gst_0->Draw("Ws","dis&&nc","");
  if(gst_1) gst_1->Draw("Ws","dis&&nc","perrsame");
  ls->Clear();
  ls->SetHeader("selected W for DIS NC events");
  ls->Draw();
  c->Update();

  //------ selected x for all events
  ps->NewPage();
  gst_0->Draw("xs","","");
  if(gst_1) gst_1->Draw("xs","","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for all events");
  ls->Draw();
  c->Update();

  //------ selected x for QEL
  ps->NewPage();
  gst_0->Draw("xs","qel&&!charm","");
  if(gst_1) gst_1->Draw("xs","qel&&!charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for QEL events");
  ls->Draw();
  c->Update();

  //------ selected x for RES
  ps->NewPage();
  gst_0->Draw("xs","res","");
  if(gst_1) gst_1->Draw("xs","res","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for RES events");
  ls->Draw();
  c->Update();

  //------ selected x for DIS
  ps->NewPage();
  gst_0->Draw("xs","dis","");
  if(gst_1) gst_1->Draw("xs","dis","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for DIS events");
  ls->Draw();
  c->Update();

  //------ selected x for DIS CC
  ps->NewPage();
  gst_0->Draw("xs","dis&&cc","");
  if(gst_1) gst_1->Draw("xs","dis&&cc","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for DIS CC events");
  ls->Draw();
  c->Update();

  //------ selected x for DIS NC
  ps->NewPage();
  gst_0->Draw("xs","dis&&nc","");
  if(gst_1) gst_1->Draw("xs","dis&&nc","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for DIS NC events");
  ls->Draw();
  c->Update();

  //------ selected x for Charm/DIS 
  ps->NewPage();
  gst_0->Draw("xs","dis&&charm","");
  if(gst_1) gst_1->Draw("xs","dis&&charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected x for Charm/DIS events");
  ls->Draw();
  c->Update();

  if(show_coh_plots) {
     //------ selected x for COH
     ps->NewPage();
     gst_0->Draw("xs","coh","");
     if(gst_1) gst_1->Draw("xs","coh","perrsame");
     ls->Clear();
     ls->SetHeader("selected x for COH events");
     ls->Draw();
     c->Update();

     //------ selected x for COH CC
     ps->NewPage();
     gst_0->Draw("xs","coh&&cc","");
     if(gst_1) gst_1->Draw("xs","coh&&cc","perrsame");
     ls->Clear();
     ls->SetHeader("selected x for COH CC events");
     ls->Draw();
     c->Update();

     //------ selected x for COH NC
     ps->NewPage();
     gst_0->Draw("xs","coh&&nc","");
     if(gst_1) gst_1->Draw("xs","coh&&nc","perrsame");
     ls->Clear();
     ls->SetHeader("selected x for COH NC events");
     ls->Draw();
     c->Update();
  }

  //------ selected y for all events
  ps->NewPage();
  gst_0->Draw("ys","","");
  if(gst_1) gst_1->Draw("ys","","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for all events");
  ls->Draw();
  c->Update();

  //------ selected y for QEL
  ps->NewPage();
  gst_0->Draw("ys","qel&&!charm","");
  if(gst_1) gst_1->Draw("ys","qel&&!charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for QEL events");
  ls->Draw();
  c->Update();

  //------ selected y for RES
  ps->NewPage();
  gst_0->Draw("ys","res","");
  if(gst_1) gst_1->Draw("ys","res","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for RES events");
  ls->Draw();
  c->Update();

  //------ selected y for DIS
  ps->NewPage();
  gst_0->Draw("ys","dis","");
  if(gst_1) gst_1->Draw("ys","dis","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for DIS events");
  ls->Draw();
  c->Update();

  //------ selected y for DIS CC
  ps->NewPage();
  gst_0->Draw("ys","dis&&cc","");
  if(gst_1) gst_1->Draw("ys","dis&&cc","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for DIS CC events");
  ls->Draw();
  c->Update();

  //------ selected y for DIS NC
  ps->NewPage();
  gst_0->Draw("ys","dis&&nc","");
  if(gst_1) gst_1->Draw("ys","dis&&nc","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for DIS NC events");
  ls->Draw();
  c->Update();

  //------ selected y for Charm/DIS 
  ps->NewPage();
  gst_0->Draw("ys","dis&&charm","");
  if(gst_1) gst_1->Draw("ys","dis&&charm","perrsame");
  ls->Clear();
  ls->SetHeader("selected y for Charm/DIS events");
  ls->Draw();
  c->Update();

  if(show_coh_plots) {
     //------ selected y for COH
     ps->NewPage();
     gst_0->Draw("ys","coh","");
     if(gst_1) gst_1->Draw("ys","coh","perrsame");
     ls->Clear();
     ls->SetHeader("selected y for COH events");
     ls->Draw();
     c->Update();

     //------ selected y for COH CC
     ps->NewPage();
     gst_0->Draw("ys","coh&&cc","");
     if(gst_1) gst_1->Draw("ys","coh&&cc","perrsame");
     ls->Clear();
     ls->SetHeader("selected y for COH CC events");
     ls->Draw();
     c->Update();

     //------ selected y for COH NC
     ps->NewPage();
     gst_0->Draw("ys","coh&&nc","");
     if(gst_1) gst_1->Draw("ys","coh&&nc","perrsame");
     ls->Clear();
     ls->SetHeader("selected y for COH NC events");
     ls->Draw();
     c->Update();

     //------ selected t for COH
     ps->NewPage();
     gst_0->Draw("ts","coh","");
     if(gst_1) gst_1->Draw("ts","coh","perrsame");
     ls->Clear();
     ls->SetHeader("selected t for COH events");
     ls->Draw();
     c->Update();
  }

  if(show_calc_kinematics) {

     //
     // SECTION: Computed Kinematics 
     //
     ps->NewPage();
     c->Clear();
     c->Range(0,0,100,100);
     TPavesText hdrck(10,40,90,70,3,"tr");
     hdrck.AddText("Kinematical Quantities");
     hdrck.AddText(" ");
     hdrck.AddText(" ");
     hdrck.AddText("Similar to the previous set of plots but");
     hdrck.AddText("showing 'computed' rather than 'selected' variables");
     hdrck.Draw();
     c->Update();

     //------ Q2 for all events
     ps->NewPage();
               gst_0->Draw("Q2","","");
     if(gst_1) gst_1->Draw("Q2","","perrsame");
     ls->Clear();
     ls->SetHeader("computed Q2 for all events");
     ls->Draw();
     c->Update();

     //------ Q2 for QEL
     ps->NewPage();
               gst_0->Draw("Q2","qel&&!charm","");
     if(gst_1) gst_1->Draw("Q2","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("computed Q2 for QEL events");
     ls->Draw();
     c->Update();

     //------ Q2 for RES
     ps->NewPage();
               gst_0->Draw("Q2","res","");
     if(gst_1) gst_1->Draw("Q2","res","perrsame");
     ls->Clear();
     ls->SetHeader("computed Q2 for RES events");
     ls->Draw();
     c->Update();

     //------ Q2 for DIS
     ps->NewPage();
               gst_0->Draw("Q2","dis","");
     if(gst_1) gst_1->Draw("Q2","dis","perrsame");
     ls->Clear();
     ls->SetHeader("computed Q2 for DIS events");
     ls->Draw();
     c->Update();

     //------ x for all events
     ps->NewPage();
     gst_0->Draw("x","","");
     if(gst_1) gst_1->Draw("x","","perrsame");
     ls->Clear();
     ls->SetHeader("computed x for all events");
     ls->Draw();
     c->Update();

     //------ x for QEL
     ps->NewPage();
     gst_0->Draw("x","qel&&!charm","");
     if(gst_1) gst_1->Draw("x","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("computed x for QEL events");
     ls->Draw();
     c->Update();

     //------ x for RES
     ps->NewPage();
     gst_0->Draw("x","res","");
     if(gst_1) gst_1->Draw("x","res","perrsame");
     ls->Clear();
     ls->SetHeader("computed x for RES events");
     ls->Draw();
     c->Update();

     //------ x for DIS
     ps->NewPage();
     gst_0->Draw("x","dis","");
     if(gst_1) gst_1->Draw("x","dis","perrsame");
     ls->Clear();
     ls->SetHeader("computed x for DIS events");
     ls->Draw();
     c->Update();

     //------ y for all events
     ps->NewPage();
     gst_0->Draw("y","","");
     if(gst_1) gst_1->Draw("y","","perrsame");
     ls->Clear();
     ls->SetHeader("computed y for all events");
     ls->Draw();
     c->Update();

     //------ y for QEL
     ps->NewPage();
     gst_0->Draw("y","qel&&!charm","");
     if(gst_1) gst_1->Draw("y","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("computed y for QEL events");
     ls->Draw();
     c->Update();

     //------ y for RES
     ps->NewPage();
     gst_0->Draw("y","res","");
     if(gst_1) gst_1->Draw("y","res","perrsame");
     ls->Clear();
     ls->SetHeader("computed y for RES events");
     ls->Draw();
     c->Update();

     //------ y for DIS
     ps->NewPage();
     gst_0->Draw("y","dis","");
     if(gst_1) gst_1->Draw("y","dis","perrsame");
     ls->Clear();
     ls->SetHeader("computed y for DIS events");
     ls->Draw();
     c->Update();

  }//show?


  //
  // SECTION: Initial State nucleon
  //
  ps->NewPage();
  c->Clear();
  c->Range(0,0,100,100);
  TPavesText hdrinuc(10,40,90,70,3,"tr");
  hdrinuc.AddText("Initial state nucleon 4-Momentum");
  hdrinuc.Draw();
  c->Update();

  //------ selected hit nucleon px
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxn","","");
  if(gst_1) gst_1->Draw("pxn","","perrsame");
  c->cd(2);
  gst_0->Draw("pyn","","");
  if(gst_1) gst_1->Draw("pyn","","perrsame");
  c->cd(3);
  gst_0->Draw("pzn","","");
  if(gst_1) gst_1->Draw("pzn","","perrsame");
  c->cd(4);
  gst_0->Draw("En","En>.2","");
  if(gst_1) gst_1->Draw("En","En>.2","perrsame");
  c->Update();

  //
  // SECTION: Final State Primary Lepton
  //
  ps->NewPage();
  c->Clear();
  c->Range(0,0,100,100);
  TPavesText hdrfsl(10,40,90,70,3,"tr");
  hdrfsl.AddText("Final State Primary Lepton 4-Momentum");
  hdrfsl.Draw();
  c->Update();

  //------ f/s primary lepton : all events
  ps->NewPage();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxl","","");
  if(gst_1) gst_1->Draw("pxl","","perrsame");
  c->cd(2);
  gst_0->Draw("pyl","","");
  if(gst_1) gst_1->Draw("pyl","","perrsame");
  c->cd(3);
  gst_0->Draw("pzl","","");
  if(gst_1) gst_1->Draw("pzl","","perrsame");
  c->cd(4);
  gst_0->Draw("El","","");
  if(gst_1) gst_1->Draw("El","","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state primary lepton 4-p: All events");
  ls->Draw();
  c->Update();

  //------ f/s primary lepton : all CC events
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxl","cc","");
  if(gst_1) gst_1->Draw("pxl","cc","perrsame");
  c->cd(2);
  gst_0->Draw("pyl","cc","");
  if(gst_1) gst_1->Draw("pyl","cc","perrsame");
  c->cd(3);
  gst_0->Draw("pzl","cc","");
  if(gst_1) gst_1->Draw("pzl","cc","perrsame");
  c->cd(4);
  gst_0->Draw("El","cc","");
  if(gst_1) gst_1->Draw("El","cc","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state primary lepton 4-p: All CC events");
  ls->Draw();
  c->Update();

  //------ f/s primary lepton : all NC events
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxl","nc","");
  if(gst_1) gst_1->Draw("pxl","nc","perrsame");
  c->cd(2);
  gst_0->Draw("pyl","nc","");
  if(gst_1) gst_1->Draw("pyl","nc","perrsame");
  c->cd(3);
  gst_0->Draw("pzl","nc","");
  if(gst_1) gst_1->Draw("pzl","nc","perrsame");
  c->cd(4);
  gst_0->Draw("El","nc","");
  if(gst_1) gst_1->Draw("El","nc","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state primary lepton 4-p: All NC events");
  ls->Draw();
  c->Update();

  //------ f/s primary lepton : QEL events
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxl","qel&&!charm","");
  if(gst_1) gst_1->Draw("pxl","qel&&!charm","perrsame");
  c->cd(2);
  gst_0->Draw("pyl","qel&&!charm","");
  if(gst_1) gst_1->Draw("pyl","qel&&!charm","perrsame");
  c->cd(3);
  gst_0->Draw("pzl","qel&&!charm","");
  if(gst_1) gst_1->Draw("pzl","qel&&!charm","perrsame");
  c->cd(4);
  gst_0->Draw("El","qel&&!charm","");
  if(gst_1) gst_1->Draw("El","qel&&!charm","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state primary lepton 4-p: QEL events");
  ls->Draw();
  c->Update();

  //------ f/s primary lepton : RES events
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxl","res","");
  if(gst_1) gst_1->Draw("pxl","res","perrsame");
  c->cd(2);
  gst_0->Draw("pyl","res","");
  if(gst_1) gst_1->Draw("pyl","res","perrsame");
  c->cd(3);
  gst_0->Draw("pzl","res","");
  if(gst_1) gst_1->Draw("pzl","res","perrsame");
  c->cd(4);
  gst_0->Draw("El","res","");
  if(gst_1) gst_1->Draw("El","res","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state primary lepton 4-p: RES events");
  ls->Draw();
  c->Update();

  //------ f/s primary lepton : DIS events
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxl","dis","");
  if(gst_1) gst_1->Draw("pxl","dis","perrsame");
  c->cd(2);
  gst_0->Draw("pyl","dis","");
  if(gst_1) gst_1->Draw("pyl","dis","perrsame");
  c->cd(3);
  gst_0->Draw("pzl","dis","");
  if(gst_1) gst_1->Draw("pzl","dis","perrsame");
  c->cd(4);
  gst_0->Draw("El","dis","");
  if(gst_1) gst_1->Draw("El","dis","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state primary lepton 4-p: All DIS events");
  ls->Draw();
  c->Update();

  if(show_coh_plots) {
     //------ f/s primary lepton : COH events
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxl","coh","");
     if(gst_1) gst_1->Draw("pxl","coh","perrsame");
     c->cd(2);
     gst_0->Draw("pyl","coh","");
     if(gst_1) gst_1->Draw("pyl","coh","perrsame");
     c->cd(3);
     gst_0->Draw("pzl","coh","");
     if(gst_1) gst_1->Draw("pzl","coh","perrsame");
     c->cd(4);
     gst_0->Draw("El","coh","");
     if(gst_1) gst_1->Draw("El","coh","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state primary lepton 4-p: COH events");
     ls->Draw();
     c->Update();
  }

  //
  // SECTION: Final State Hadronic System Multiplicities & 4P
  //
  ps->NewPage();
  c->Clear();
  c->Range(0,0,100,100);
  TPavesText hdrfhad(10,40,90,70,3,"tr");
  hdrfhad.AddText("Final State Hadronic System");
  hdrfhad.AddText("Multiplicities and 4-Momenta");
  hdrfhad.AddText(" ");
  hdrfhad.AddText(" ");
  hdrfhad.AddText(" ");
  hdrfhad.AddText(" ");
  hdrfhad.AddText("Note:");
  hdrfhad.AddText("For nuclear targets these plots include the effect");
  hdrfhad.AddText("of intranuclear hadron transport / rescattering");
  hdrfhad.Draw();
  c->Update();

  //------ number of final state p
  ps->NewPage();
  gst_0->Draw("nfp","","");
  if(gst_1) gst_1->Draw("nfp","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state protons");
  ls->Draw();
  c->Update();

  //------ number of final state n
  ps->NewPage();
  gst_0->Draw("nfn","","");
  if(gst_1) gst_1->Draw("nfn","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state neutrons");
  ls->Draw();
  c->Update();

  //------ number of final state pi+
  ps->NewPage();
  gst_0->Draw("nfpip","","");
  if(gst_1) gst_1->Draw("nfpip","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state pi+");
  ls->Draw();
  c->Update();

  //------ number of final state pi-
  ps->NewPage();
  gst_0->Draw("nfpim","","");
  if(gst_1) gst_1->Draw("nfpim","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state pi-");
  ls->Draw();
  c->Update();

  //------ number of final state pi0
  ps->NewPage();
  gst_0->Draw("nfpi0","","");
  if(gst_1) gst_1->Draw("nfpi0","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state pi0");
  ls->Draw();
  c->Update();

  //------ number of final state K+
  ps->NewPage();
  gst_0->Draw("nfkp","","");
  if(gst_1) gst_1->Draw("nfkp","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state K+");
  ls->Draw();
  c->Update();

  //------ number of final state K-
  ps->NewPage();
  gst_0->Draw("nfkm","","");
  if(gst_1) gst_1->Draw("nfkm","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state K-");
  ls->Draw();
  c->Update();

  //------ number of final state K0
  ps->NewPage();
  gst_0->Draw("nfk0","","");
  if(gst_1) gst_1->Draw("nfk0","","perrsame");
  ls->Clear();
  ls->SetHeader("Number of final state K0");
  ls->Draw();
  c->Update();

  //------ momentum of final state p
  ps->NewPage();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxf","pdgf==2212","");
  if(gst_1) gst_1->Draw("pxf","pdgf==2212","perrsame");
  c->cd(2);
  gst_0->Draw("pyf","pdgf==2212","");
  if(gst_1) gst_1->Draw("pyf","pdgf==2212","perrsame");
  c->cd(3);
  gst_0->Draw("pzf","pdgf==2212","");
  if(gst_1) gst_1->Draw("pzf","pdgf==2212","perrsame"); 
  c->cd(4);
  gst_0->Draw("Ef","pdgf==2212","");
  if(gst_1) gst_1->Draw("Ef","pdgf==2212","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state protons 4-momentum");
  ls->Draw();
  c->Update();

  //------ momentum of final state n
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxf","pdgf==2112","");
  if(gst_1) gst_1->Draw("pxf","pdgf==2112","perrsame");
  c->cd(2);
  gst_0->Draw("pyf","pdgf==2112","");
  if(gst_1) gst_1->Draw("pyf","pdgf==2112","perrsame");
  c->cd(3);
  gst_0->Draw("pzf","pdgf==2112","");
  if(gst_1) gst_1->Draw("pzf","pdgf==2112","perrsame");
  c->cd(4);
  gst_0->Draw("Ef","pdgf==2112","");
  if(gst_1) gst_1->Draw("Ef","pdgf==2112","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state neutrons 4-momentum");
  ls->Draw();
  c->Update();

  //------ momentum of final state pi0
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxf","pdgf==111","");
  if(gst_1) gst_1->Draw("pxf","pdgf==111","perrsame");
  c->cd(2);
  gst_0->Draw("pyf","pdgf==111","");
  if(gst_1) gst_1->Draw("pyf","pdgf==111","perrsame");
  c->cd(3);
  gst_0->Draw("pzf","pdgf==111","");
  if(gst_1) gst_1->Draw("pzf","pdgf==111","perrsame");
  c->cd(4);
  gst_0->Draw("Ef","pdgf==111","");
  if(gst_1) gst_1->Draw("Ef","pdgf==111","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state pi0's 4-momentum");
  ls->Draw();
  c->Update();

  //------ momentum of final state pi+
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxf","pdgf==211","");
  if(gst_1) gst_1->Draw("pxf","pdgf==211","perrsame");
  c->cd(2);
  gst_0->Draw("pyf","pdgf==211","");
  if(gst_1) gst_1->Draw("pyf","pdgf==211","perrsame");
  c->cd(3);
  gst_0->Draw("pzf","pdgf==211","");
  if(gst_1) gst_1->Draw("pzf","pdgf==211","perrsame");
  c->cd(4);
  gst_0->Draw("Ef","pdgf==211","");
  if(gst_1) gst_1->Draw("Ef","pdgf==211","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state pi+'s 4-momentum");
  ls->Draw();
  c->Update();

  //------ momentum of final state pi+
  ps->NewPage();
  c->Clear();
  c->Divide(2,2);
  c->cd(1);
  gst_0->Draw("pxf","pdgf==-211","");
  if(gst_1) gst_1->Draw("pxf","pdgf==-211","perrsame");
  c->cd(2);
  gst_0->Draw("pyf","pdgf==-211","");
  if(gst_1) gst_1->Draw("pyf","pdgf==-211","perrsame");
  c->cd(3);
  gst_0->Draw("pzf","pdgf==-211","");
  if(gst_1) gst_1->Draw("pzf","pdgf==-211","perrsame");
  c->cd(4);
  gst_0->Draw("Ef","pdgf==-211","");
  if(gst_1) gst_1->Draw("Ef","pdgf==-211","perrsame");
  c->cd();
  ls->Clear();
  ls->SetHeader("Final state pi-'s 4-momentum");
  ls->Draw();
  c->Update();

  if(show_mult_per_proc) {

     //
     // similarly but for QEL events only
     //

     //------ number of final state p /QEL
     ps->NewPage();
     if(gst_1) gst_1->Draw("nfp","qel&&!charm","");
     gst_0->Draw("nfp","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state protons / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state n /QEL
     ps->NewPage();
     if(gst_1) gst_1->Draw("nfn","qel&&!charm","");
     gst_0->Draw("nfn","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state neutrons / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state pi+ /QEL
     ps->NewPage();
     gst_0->Draw("nfpip","qel&&!charm","");
     if(gst_1) gst_1->Draw("nfpip","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi+ / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state pi- /QEL
     ps->NewPage();
     gst_0->Draw("nfpim","qel&&!charm","");
     if(gst_1) gst_1->Draw("nfpim","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi- / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state pi0 /QEL
     ps->NewPage();
     gst_0->Draw("nfpi0","qel&&!charm","");
     if(gst_1) gst_1->Draw("nfpi0","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi0 / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state K+ /QEL
     ps->NewPage();
     gst_0->Draw("nfkp","qel&&!charm","");
     if(gst_1) gst_1->Draw("nfkp","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K+ / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state K- /QEL
     ps->NewPage();
     gst_0->Draw("nfkm","qel&&!charm","");
     if(gst_1) gst_1->Draw("nfkm","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K- / QEL only");
     ls->Draw();
     c->Update();

     //------ number of final state K0 /QEL
     ps->NewPage();
     gst_0->Draw("nfk0","qel&&!charm","");
     if(gst_1) gst_1->Draw("nfk0","qel&&!charm","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K0 / QEL only");
     ls->Draw();
     c->Update();

     //------ momentum of final state p /QEL
     ps->NewPage();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","qel&&!charm&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==2212","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","qel&&!charm&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==2212","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","qel&&!charm&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==2212","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","qel&&!charm&&pdgf==2212","");
     if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==2212","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state protons 4-momentum / QEL only");
     ls->Draw();
     c->Update();

     //------ momentum of final state n /QEL
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","qel&&!charm&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==2112","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","qel&&!charm&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==2112","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","qel&&!charm&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==2112","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","qel&&!charm&&pdgf==2112","");
     if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==2112","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state neutrons 4-momentum / QEL only");
     ls->Draw();
     c->Update();

     //------ momentum of final state pi0 /QEL
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","qel&&!charm&&pdgf==111","");
     if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==111","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","qel&&!charm&&pdgf==111","");
     if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==111","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","qel&&!charm&&pdgf==111","");
     if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==111","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","qel&&!charm&&pdgf==111","");
     if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==111","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi0's 4-momentum / QEL only");
     ls->Draw();
     c->Update();

     //------ momentum of final state pi+ /QEL
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","qel&&!charm&&pdgf==211","");
     if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==211","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","qel&&!charm&&pdgf==211","");
     if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==211","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","qel&&!charm&&pdgf==211","");
     if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==211","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","qel&&!charm&&pdgf==211","");
     if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi+'s 4-momentum / QEL only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state pi+ /QEL
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","qel&&!charm&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==-211","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","qel&&!charm&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==-211","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","qel&&!charm&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==-211","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","qel&&!charm&&pdgf==-211","");
     if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==-211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi-'s 4-momentum/ QEL only");
     ls->Draw();
     c->Update();
          
     //
     // similarly but for RES events only
     //
     
     //------ number of final state p /RES
     ps->NewPage();
               gst_0->Draw("nfp","res","");
     if(gst_1) gst_1->Draw("nfp","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state protons / RES only");
     ls->Draw();
     c->Update();
     
     //------ number of final state n /RES
     ps->NewPage();
               gst_0->Draw("nfn","res","");
     if(gst_1) gst_1->Draw("nfn","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state neutrons / RES only");
     ls->Draw();
     c->Update();
     
     //------ number of final state pi+ /RES
     ps->NewPage();
               gst_0->Draw("nfpip","res","");
     if(gst_1) gst_1->Draw("nfpip","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi+ / RES only");
     ls->Draw();
     c->Update();
     
     //------ number of final state pi- /RES
     ps->NewPage();
               gst_0->Draw("nfpim","res","");
     if(gst_1) gst_1->Draw("nfpim","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi- / RES only");
     ls->Draw();
     c->Update();

     //------ number of final state pi0 /RES
     ps->NewPage();
               gst_0->Draw("nfpi0","res","");
     if(gst_1) gst_1->Draw("nfpi0","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi0 / RES only");
     ls->Draw();
     c->Update();
     
     //------ number of final state K+ /RES
     ps->NewPage();
     gst_0->Draw("nfkp","res","");
     if(gst_1) gst_1->Draw("nfkp","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K+ / RES only");
     ls->Draw();
     c->Update();
     
     //------ number of final state K- /RES
     ps->NewPage();
     gst_0->Draw("nfkm","res","");
     if(gst_1) gst_1->Draw("nfkm","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K- / RES only");
     ls->Draw();
     c->Update();
     
     //------ number of final state K0 /RES
     ps->NewPage();
     gst_0->Draw("nfk0","res","");
     if(gst_1) gst_1->Draw("nfk0","res","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K0 / RES only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state p /RES
     ps->NewPage();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","res&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pxf","res&&pdgf==2212","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","res&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pyf","res&&pdgf==2212","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","res&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pzf","res&&pdgf==2212","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","res&&pdgf==2212","");
     if(gst_1) gst_1->Draw("Ef","res&&pdgf==2212","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state protons 4-momentum / RES only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state n /RES
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","res&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pxf","res&&pdgf==2112","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","res&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pyf","res&&pdgf==2112","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","res&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pzf","res&&pdgf==2112","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","res&&pdgf==2112","");
     if(gst_1) gst_1->Draw("Ef","res&&pdgf==2112","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state neutrons 4-momentum / RES only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state pi0 /RES
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","res&&pdgf==111","");
     if(gst_1) gst_1->Draw("pxf","res&&pdgf==111","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","res&&pdgf==111","");
     if(gst_1) gst_1->Draw("pyf","res&&pdgf==111","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","res&&pdgf==111","");
     if(gst_1) gst_1->Draw("pzf","res&&pdgf==111","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","res&&pdgf==111","");
     if(gst_1) gst_1->Draw("Ef","res&&pdgf==111","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi0's 4-momentum / RES only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state pi+ /RES
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","res&&pdgf==211","");
     if(gst_1) gst_1->Draw("pxf","res&&pdgf==211","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","res&&pdgf==211","");
     if(gst_1) gst_1->Draw("pyf","res&&pdgf==211","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","res&&pdgf==211","");
     if(gst_1) gst_1->Draw("pzf","res&&pdgf==211","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","res&&pdgf==211","");
     if(gst_1) gst_1->Draw("Ef","res&&pdgf==211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi+'s 4-momentum / RES only");
     ls->Draw();
     c->Update();

     //------ momentum of final state pi+ /RES
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","res&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pxf","res&&pdgf==-211","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","res&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pyf","res&&pdgf==-211","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","res&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pzf","res&&pdgf==-211","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","res&&pdgf==-211","");
     if(gst_1) gst_1->Draw("Ef","res&&pdgf==-211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi-'s 4-momentum/ RES only");
     ls->Draw();
     c->Update();
     
     //
     // similarly but for DIS events only
     //
     
     //------ number of final state p /DIS
     ps->NewPage();
     gst_0->Draw("nfp","dis","");
     if(gst_1) gst_1->Draw("nfp","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state protons / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state n /DIS
     ps->NewPage();
     gst_0->Draw("nfn","dis","");
     if(gst_1) gst_1->Draw("nfn","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state neutrons / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state pi+ /DIS
     ps->NewPage();
     gst_0->Draw("nfpip","dis","");
     if(gst_1) gst_1->Draw("nfpip","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi+ / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state pi- /DIS
     ps->NewPage();
     gst_0->Draw("nfpim","dis","");
     if(gst_1) gst_1->Draw("nfpim","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi- / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state pi0 /DIS
     ps->NewPage();
     gst_0->Draw("nfpi0","dis","");
     if(gst_1) gst_1->Draw("nfpi0","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state pi0 / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state K+ /DIS
     ps->NewPage();
     gst_0->Draw("nfkp","dis","");
     if(gst_1) gst_1->Draw("nfkp","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K+ / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state K- /DIS
     ps->NewPage();
     gst_0->Draw("nfkm","dis","");
     if(gst_1) gst_1->Draw("nfkm","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K- / DIS only");
     ls->Draw();
     c->Update();
     
     //------ number of final state K0 /DIS
     ps->NewPage();
     gst_0->Draw("nfk0","dis","");
     if(gst_1) gst_1->Draw("nfk0","dis","perrsame");
     ls->Clear();
     ls->SetHeader("Number of final state K0 / DIS only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state p /DIS
     ps->NewPage();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","dis&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pxf","dis&&pdgf==2212","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","dis&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pyf","dis&&pdgf==2212","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","dis&&pdgf==2212","");
     if(gst_1) gst_1->Draw("pzf","dis&&pdgf==2212","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","dis&&pdgf==2212","");
     if(gst_1) gst_1->Draw("Ef","dis&&pdgf==2212","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state protons 4-momentum / DIS only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state n /DIS
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","dis&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pxf","dis&&pdgf==2112","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","dis&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pyf","dis&&pdgf==2112","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","dis&&pdgf==2112","");
     if(gst_1) gst_1->Draw("pzf","dis&&pdgf==2112","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","dis&&pdgf==2112","");
     if(gst_1) gst_1->Draw("Ef","dis&&pdgf==2112","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state neutrons 4-momentum / DIS only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state pi0 /DIS
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","dis&&pdgf==111","");
     if(gst_1) gst_1->Draw("pxf","dis&&pdgf==111","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","dis&&pdgf==111","");
     if(gst_1) gst_1->Draw("pyf","dis&&pdgf==111","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","dis&&pdgf==111","");
     if(gst_1) gst_1->Draw("pzf","dis&&pdgf==111","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","dis&&pdgf==111","");
     if(gst_1) gst_1->Draw("Ef","dis&&pdgf==111","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi0's 4-momentum / DIS only");
     ls->Draw();
     c->Update();

     //------ momentum of final state pi+ /DIS
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","dis&&pdgf==211","");
     if(gst_1) gst_1->Draw("pxf","dis&&pdgf==211","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","dis&&pdgf==211","");
     if(gst_1) gst_1->Draw("pyf","dis&&pdgf==211","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","dis&&pdgf==211","");
     if(gst_1) gst_1->Draw("pzf","dis&&pdgf==211","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","dis&&pdgf==211","");
     if(gst_1) gst_1->Draw("Ef","dis&&pdgf==211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi+'s 4-momentum / DIS only");
     ls->Draw();
     c->Update();
     
     //------ momentum of final state pi+ /DIS
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxf","dis&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pxf","dis&&pdgf==-211","perrsame");
     c->cd(2);
     gst_0->Draw("pyf","dis&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pyf","dis&&pdgf==-211","perrsame");
     c->cd(3);
     gst_0->Draw("pzf","dis&&pdgf==-211","");
     if(gst_1) gst_1->Draw("pzf","dis&&pdgf==-211","perrsame");
     c->cd(4);
     gst_0->Draw("Ef","dis&&pdgf==-211","");
     if(gst_1) gst_1->Draw("Ef","dis&&pdgf==-211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Final state pi-'s 4-momentum/ DIS only");
     ls->Draw();
     c->Update();
     
  } // per-proc

  //
  // SECTION: Primary Hadronic System Multiplicities & 4P
  //
  if(show_primary_hadsyst) {
     
     ps->NewPage();
     c->Clear();
     c->Range(0,0,100,100);
     TPavesText hdrihad(10,40,90,70,3,"tr");
     hdrihad.AddText("Parimary Hadronic System");
     hdrihad.AddText("Multiplicities and 4-Momenta");
     hdrihad.AddText(" ");
     hdrihad.AddText(" ");
     hdrihad.AddText(" ");
     hdrihad.AddText(" ");
     hdrihad.AddText("Note:");
     hdrihad.AddText("For nuclear targets these plots show the hadronic system");
     hdrihad.AddText("BEFORE any intranuclear hadron transport / rescattering");
     hdrihad.Draw();
     c->Update();
     
     //------ number of prim p
     ps->NewPage();
     gst_0->Draw("nip","","");
     if(gst_1) gst_1->Draw("nip","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of protons");
     ls->Draw();
     c->Update();
     
     //------ number of prim n
     ps->NewPage();
     gst_0->Draw("nin","","");
     if(gst_1) gst_1->Draw("nin","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of neutrons");
     ls->Draw();
     c->Update();
     
     //------ number of prim pi+
     ps->NewPage();
     gst_0->Draw("nipip","","");
     if(gst_1) gst_1->Draw("nipip","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of pi+");
     ls->Draw();
     c->Update();
     
     //------ number of prim pi-
     ps->NewPage();
     gst_0->Draw("nipim","","");
     if(gst_1) gst_1->Draw("nipim","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of pi-");
     ls->Draw();
     c->Update();
     
     //------ number of prim pi0
     ps->NewPage();
     gst_0->Draw("nipi0","","");
     if(gst_1) gst_1->Draw("nipi0","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of pi0");
     ls->Draw();
     c->Update();
     
     //------ number of prim K+
     ps->NewPage();
     gst_0->Draw("nikp","","");
     if(gst_1) gst_1->Draw("nikp","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of K+");
     ls->Draw();
     c->Update();
     
     //------ number of prim K-
     ps->NewPage();
     gst_0->Draw("nikm","","");
     if(gst_1) gst_1->Draw("nikm","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of K-");
     ls->Draw();
     c->Update();
     
     //------ number of prim K0
     ps->NewPage();
     gst_0->Draw("nik0","","");
     if(gst_1) gst_1->Draw("nik0","","perrsame");
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: Number of K0");
     ls->Draw();
     c->Update();
     
     //------ momentum of prim, p
     ps->NewPage();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxi","pdgi==2212","");
     if(gst_1) gst_1->Draw("pxi","pdgi==2212","perrsame");
     c->cd(2);
     gst_0->Draw("pyi","pdgi==2212","");
     if(gst_1) gst_1->Draw("pyi","pdgi==2212","perrsame");
     c->cd(3);
     gst_0->Draw("pzi","pdgi==2212","");
     if(gst_1) gst_1->Draw("pzi","pdgi==2212","perrsame");
     c->cd(4);
     gst_0->Draw("Ei","pdgi==2212","");
     if(gst_1) gst_1->Draw("Ei","pdgi==2212","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: proton 4-momentum");
     ls->Draw();
     c->Update();
     
     //------ momentum of prim. n
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxi","pdgi==2112","");
     if(gst_1) gst_1->Draw("pxi","pdgi==2112","perrsame");
     c->cd(2);
     gst_0->Draw("pyi","pdgi==2112","");
     if(gst_1) gst_1->Draw("pyi","pdgi==2112","perrsame");
     c->cd(3);
     gst_0->Draw("pzi","pdgi==2112","");
     if(gst_1) gst_1->Draw("pzi","pdgi==2112","perrsame");
     c->cd(4);
     gst_0->Draw("Ei","pdgi==2112","");
     if(gst_1) gst_1->Draw("Ei","pdgi==2112","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: neutron 4-momentum");
     ls->Draw();
     c->Update();
     
     //------ momentum of prim. pi0
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxi","pdgi==111","");
     if(gst_1) gst_1->Draw("pxi","pdgi==111","perrsame");
     c->cd(2);
     gst_0->Draw("pyi","pdgi==111","");
     if(gst_1) gst_1->Draw("pyi","pdgi==111","perrsame");
     c->cd(3);
     gst_0->Draw("pzi","pdgi==111","");
     if(gst_1) gst_1->Draw("pzi","pdgi==111","perrsame");
     c->cd(4);
     gst_0->Draw("Ei","pdgi==111","");
     if(gst_1) gst_1->Draw("Ei","pdgi==111","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Primary Hadronic System: pi0's 4-momentum");
     ls->Draw();
     c->Update();

     //------ momentum of prim pi+
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxi","pdgi==211","");
     if(gst_1) gst_1->Draw("pxi","pdgi==211","perrsame");
     c->cd(2);
     gst_0->Draw("pyi","pdgi==211","");
     if(gst_1) gst_1->Draw("pyi","pdgi==211","perrsame");
     c->cd(3);
     gst_0->Draw("pzi","pdgi==211","");
     if(gst_1) gst_1->Draw("pzi","pdgi==211","perrsame");
     c->cd(4);
     gst_0->Draw("Ei","pdgi==211","");
     if(gst_1) gst_1->Draw("Ei","pdgi==211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Primary Hadronic System:pi+'s 4-momentum");
     ls->Draw();
     c->Update();
     
     //------ momentum of prim. pi+
     ps->NewPage();
     c->Clear();
     c->Divide(2,2);
     c->cd(1);
     gst_0->Draw("pxi","pdgi==-211","");
     if(gst_1) gst_1->Draw("pxi","pdgi==-211","perrsame");
     c->cd(2);
     gst_0->Draw("pyi","pdgi==-211","");
     if(gst_1) gst_1->Draw("pyi","pdgi==-211","perrsame");
     c->cd(3);
     gst_0->Draw("pzi","pdgi==-211","");
     if(gst_1) gst_1->Draw("pzi","pdgi==-211","perrsame");
     c->cd(4);
     gst_0->Draw("Ei","pdgi==-211","");
     if(gst_1) gst_1->Draw("Ei","pdgi==-211","perrsame");
     c->cd();
     ls->Clear();
     ls->SetHeader("Primary Hadronic System:  pi-'s 4-momentum");
     ls->Draw();
     c->Update();
  }//show?
       
  ps->Close();

  fin_0->Close();
  delete fin_0;
  if(fin_1) {
     fin_1->Close();
     delete fin_1;
  }
}
//_________________________________________________________________________________
string OutputFileName(string inpname)
{
// Builds the output filename based on the name of the input filename
// Perfors the following conversion: name.root -> name.nusample_test.ps

  unsigned int L = inpname.length();

  // if the last 4 characters are "root" (ROOT file extension) then
  // remove them
  if(inpname.substr(L-4, L).find("root") != string::npos) {
    inpname.erase(L-4, L);
  }

  ostringstream name;
  name << inpname << "sample_test.ps";

  return gSystem->BaseName(name.str().c_str());
}
//_________________________________________________________________________________
void GetCommandLineArgs(int argc, char ** argv)
{
  LOG("gevcomp", pNOTICE) << "*** Parsing command line arguments";

  CmdLnArgParser parser(argc,argv);

  // get GENIE summary ntuple
  if( parser.OptionExists('f') ) {
    LOG("gevcomp", pINFO) << "Reading filename for tested event sample";
    gOptInpFile = parser.ArgAsString('f');
  } else {
    LOG("gevcomp", pFATAL) << "Unspecified input filename - Exiting";
    PrintSyntax();
    exit(1);
  }

  // get another (reference) GENIE summary ntuple
  if( parser.OptionExists('r') ) {
    LOG("gevcomp", pINFO) << "Reading filename for reference event sample";
    gOptInpFileRef = parser.ArgAsString('r');
  } else {
    LOG("gevcomp", pNOTICE) << "Unspecified 'reference' event sample";
  }
}
//_________________________________________________________________________________
void PrintSyntax(void)
{
  LOG("gevcomp", pNOTICE)
    << "\n\n" << "Syntax:" << "\n"
    << " gevcomp -f sample.root [-n nev] [-r reference_sample.root]\n";
}
//_________________________________________________________________________________
bool CheckRootFilename(string filename)
{
  if(filename.size() == 0) return false;

  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));
  if (!is_accessible) {
   LOG("gevcomp", pERROR)
       << "The input ROOT file [" << filename << "] is not accessible";
   return false;
  }
  return true;
}
//_________________________________________________________________________________