GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions | Variables
gSplineXml2Root.cxx File Reference
#include <cassert>
#include <string>
#include <sstream>
#include <vector>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TDirectory.h>
#include <TPostScript.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TGraph.h>
#include <TPaveText.h>
#include <TString.h>
#include <TH1F.h>
#include "Framework/ParticleData/BaryonResonance.h"
#include "Framework/ParticleData/BaryonResUtils.h"
#include "Framework/Conventions/XmlParserStatus.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/EventGen/InteractionList.h"
#include "Framework/EventGen/GEVGDriver.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/Spline.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Registry/Registry.h"
#include "Framework/Algorithm/AlgConfigPool.h"
Include dependency graph for gSplineXml2Root.cxx:

Go to the source code of this file.

Functions

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

Variables

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

Function Documentation

void FormatXSecGraph ( TGraph *  g)

Definition at line 516 of file gSplineXml2Root.cxx.

Referenced by SaveGraphsToRootFile().

517 {
518  g->SetTitle("GENIE cross section graph");
519  g->GetXaxis()->SetTitle("Ev (GeV)");
520  g->GetYaxis()->SetTitle("#sigma_{nuclear} (10^{-38} cm^{2})");
521 }
genie::units::g
static constexpr double g
Definition: Units.h:144
void GetCommandLineArgs ( int  argc,
char **  argv 
)
GEVGDriver GetEventGenDriver ( void  )

Definition at line 200 of file gSplineXml2Root.cxx.

References genie::GEVGDriver::Configure(), genie::GEVGDriver::CreateSplines(), genie::GEVGDriver::CreateXSecSumSpline(), gEmax, gEmin, gOptProbePdgCode, gOptTgtPdgCode, genie::RunOpt::Instance(), and genie::GEVGDriver::SetEventGeneratorList().

Referenced by SaveGraphsToRootFile(), and SaveToPsFile().

201 {
202 // create an event genartion driver configured for the specified initial state
203 // (so that cross section splines will be accessed through that driver as in
204 // event generation mode)
205 
206  InitialState init_state(gOptTgtPdgCode, gOptProbePdgCode);
207 
208  GEVGDriver evg_driver;
209  evg_driver.SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
210  evg_driver.Configure(init_state);
211  evg_driver.CreateSplines();
212  evg_driver.CreateXSecSumSpline (100, gEmin, gEmax);
213 
214  return evg_driver;
215 }
genie::GEVGDriver::CreateXSecSumSpline
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
Definition: GEVGDriver.cxx:440
gEmin
double gEmin
Definition: gSplineXml2Root.cxx:150
genie::GEVGDriver::SetEventGeneratorList
void SetEventGeneratorList(string listname)
Definition: GEVGDriver.cxx:348
genie::GEVGDriver
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:54
gOptProbePdgCode
int gOptProbePdgCode
Definition: gEvGenHadronNucleus.cxx:146
gEmax
double gEmax
Definition: gSplineXml2Root.cxx:151
genie::GEVGDriver::Configure
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:137
genie::EventGeneratorList
A vector of EventGeneratorI objects.
Definition: EventGeneratorList.h:34
genie::GEVGDriver::CreateSplines
void CreateSplines(int nknots=-1, double emax=-1, bool inLogE=true)
Definition: GEVGDriver.cxx:577
genie::InitialState
Initial State information.
Definition: InitialState.h:48
gOptTgtPdgCode
int gOptTgtPdgCode
Definition: gEvGenHadronNucleus.cxx:147
PDGCodeList GetPDGCodeListFromString ( std::string  s)

Definition at line 1575 of file gSplineXml2Root.cxx.

References genie::PDGCodeList::push_back(), and genie::utils::str::Split().

1576 {
1577  vector<string> isvec = utils::str::Split(s,",");
1578 
1579  // fill in the PDG code list
1580  PDGCodeList list;
1581  vector<string>::const_iterator iter;
1582  for(iter = isvec.begin(); iter != isvec.end(); ++iter) {
1583  list.push_back( atoi(iter->c_str()) );
1584  }
1585  return list;
1586 
1587 }
genie::units::s
static constexpr double s
Definition: Units.h:95
genie::PDGCodeList
A list of PDG codes.
Definition: PDGCodeList.h:32
genie::utils::str::Split
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
void LoadSplines ( void  )

Definition at line 191 of file gSplineXml2Root.cxx.

References gOptXMLFilename, genie::XSecSplineList::Instance(), genie::kXmlOK, and genie::XSecSplineList::LoadFromXml().

Referenced by main().

192 {
193 // load the cross section splines specified at the cmd line
194 
195  XSecSplineList * splist = XSecSplineList::Instance();
196  XmlParserStatus_t ist = splist->LoadFromXml(gOptXMLFilename);
197  assert(ist == kXmlOK);
198 }
gOptXMLFilename
string gOptXMLFilename
Definition: gMaxPathLengths.cxx:103
genie::XmlParserStatus_t
enum genie::EXmlParseStatus XmlParserStatus_t
genie::XSecSplineList
List of cross section vs energy splines.
Definition: XSecSplineList.h:45
genie::XSecSplineList::LoadFromXml
XmlParserStatus_t LoadFromXml(const string &filename, bool keep=false)
Definition: XSecSplineList.cxx:402
int main ( int  argc,
char **  argv 
)

Definition at line 158 of file gSplineXml2Root.cxx.

References genie::RunOpt::BuildTune(), GetCommandLineArgs(), gOptProbePdgCode, gOptProbePdgList, gOptTgtPdgCode, gOptTgtPdgList, genie::RunOpt::Instance(), genie::XSecSplineList::Instance(), LoadSplines(), LOG, genie::utils::app_init::MesgThresholds(), pFATAL, pWARN, SaveGraphsToRootFile(), and SaveToPsFile().

159 {
160  GetCommandLineArgs(argc,argv);
161 
162  if ( ! RunOpt::Instance()->Tune() ) {
163  LOG("gslp2root", pFATAL) << " No TuneId in RunOption";
164  exit(-1);
165  }
166  RunOpt::Instance()->BuildTune();
167 
168  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
169 
170  // load the x-section splines xml file specified by the user
171  LoadSplines();
172 
173  if ( ! XSecSplineList::Instance() -> HasSplineFromTune(RunOpt::Instance() -> Tune() -> Name() ) ) {
174  LOG("gspl2root", pWARN) << "No splines loaded for tune " << RunOpt::Instance() -> Tune() -> Name() ;
175  }
176 
177  for (unsigned int indx_p = 0; indx_p < gOptProbePdgList.size(); ++indx_p ) {
178  for (unsigned int indx_t = 0; indx_t < gOptTgtPdgList.size(); ++indx_t ) {
179  gOptProbePdgCode = gOptProbePdgList[indx_p];
180  gOptTgtPdgCode = gOptTgtPdgList[indx_t];
181  // save the cross section plots in a postscript file
182  SaveToPsFile();
183  // save the cross section graphs at a root file
184  SaveGraphsToRootFile();
185  }
186  }
187 
188  return 0;
189 }
LoadSplines
void LoadSplines(void)
Definition: gSplineXml2Root.cxx:191
pFATAL
#define pFATAL
Definition: Messenger.h:56
SaveToPsFile
void SaveToPsFile(void)
Definition: gSplineXml2Root.cxx:217
gOptProbePdgList
PDGCodeList gOptProbePdgList
Definition: gSplineXml2Root.cxx:142
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
gOptProbePdgCode
int gOptProbePdgCode
Definition: gEvGenHadronNucleus.cxx:146
SaveGraphsToRootFile
void SaveGraphsToRootFile(void)
Definition: gSplineXml2Root.cxx:523
pWARN
#define pWARN
Definition: Messenger.h:60
genie::utils::app_init::MesgThresholds
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:99
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
gOptTgtPdgList
PDGCodeList gOptTgtPdgList
Definition: gSplineXml2Root.cxx:143
gOptTgtPdgCode
int gOptTgtPdgCode
Definition: gEvGenHadronNucleus.cxx:147
void PrintSyntax ( void  )
void SaveGraphsToRootFile ( void  )

Definition at line 523 of file gSplineXml2Root.cxx.

References genie::utils::res::AsString(), genie::XclsTag::CharmHadronPdg(), genie::units::cm2, e, genie::Spline::Evaluate(), genie::Interaction::ExclTag(), genie::utils::str::FilterString(), genie::XclsTag::FinalLeptonPdg(), genie::XclsTag::FinalQuarkPdg(), genie::PDGLibrary::Find(), FormatXSecGraph(), gEmax, gEmin, GetEventGenDriver(), gInlogE, gOptProbePdgCode, gOptROOTFilename, gOptTgtPdgCode, genie::Target::HitNucIsSet(), genie::Target::HitNucPdg(), genie::Target::HitQrkIsSet(), genie::Target::HitQrkPdg(), genie::Target::HitSeaQrk(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::GEVGDriver::Interactions(), genie::pdg::Is2NucleonCluster(), genie::pdg::IsAntiBQuark(), genie::pdg::IsAntiCQuark(), genie::pdg::IsAntiDQuark(), genie::pdg::IsAntiSQuark(), genie::pdg::IsAntiTQuark(), genie::pdg::IsAntiUQuark(), genie::pdg::IsBQuark(), genie::pdg::IsChargedLepton(), genie::XclsTag::IsCharmEvent(), genie::ProcessInfo::IsCoherentElastic(), genie::ProcessInfo::IsCoherentProduction(), genie::pdg::IsCQuark(), genie::ProcessInfo::IsDarkNeutralCurrent(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::pdg::IsDQuark(), genie::pdg::IsElectron(), genie::ProcessInfo::IsEM(), genie::XclsTag::IsFinalLeptonEvent(), genie::XclsTag::IsFinalQuarkEvent(), genie::ProcessInfo::IsGlashowResonance(), genie::ProcessInfo::IsIMDAnnihilation(), genie::XclsTag::IsInclusiveCharm(), genie::XclsTag::IsInclusiveStrange(), genie::ProcessInfo::IsInverseMuDecay(), genie::ProcessInfo::IsMEC(), genie::pdg::IsMuon(), genie::pdg::IsNeutralLepton(), genie::pdg::IsNeutron(), genie::ProcessInfo::IsNuElectronElastic(), genie::ProcessInfo::IsPhotonCoherent(), genie::ProcessInfo::IsPhotonResonance(), genie::pdg::IsPion(), genie::pdg::IsProton(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::pdg::IsSQuark(), genie::XclsTag::IsStrangeEvent(), genie::pdg::IsTau(), genie::pdg::IsTQuark(), genie::pdg::IsUQuark(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakMix(), genie::ProcessInfo::IsWeakNC(), genie::XclsTag::KnownResonance(), kNSplineP, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, LOG, genie::XclsTag::NPi0(), genie::XclsTag::NPiMinus(), genie::XclsTag::NPions(), genie::XclsTag::NPiPlus(), genie::XclsTag::NRhos(), genie::XclsTag::NSingleGammas(), pINFO, genie::Interaction::ProcInfo(), pWARN, genie::XclsTag::Resonance(), genie::XclsTag::StrangeHadronPdg(), genie::InitialState::Tgt(), and genie::GEVGDriver::XSecSpline().

Referenced by main().

524 {
525  //-- get the event generation driver
526  GEVGDriver evg_driver = GetEventGenDriver();
527 
528  //-- get the list of interactions that can be simulated by the driver
529  const InteractionList * ilist = evg_driver.Interactions();
530 
531  //-- check whether the splines will be saved in a ROOT file - if not, exit now
532  bool save_in_root = gOptROOTFilename.size()>0;
533  if(!save_in_root) {
534 
535  LOG("gspl2root", pWARN) << "No Interaction List available" ;
536  return;
537  }
538 
539  //-- get pdglibrary for mapping pdg codes to names
540  PDGLibrary * pdglib = PDGLibrary::Instance();
541 
542  //-- check whether the requested filename exists
543  // if yes, then open the file in 'update' mode
544  bool exists = !(gSystem->AccessPathName(gOptROOTFilename.c_str()));
545 
546  TFile * froot = 0;
547  if(exists) froot = new TFile(gOptROOTFilename.c_str(), "UPDATE");
548  else froot = new TFile(gOptROOTFilename.c_str(), "RECREATE");
549  assert(froot);
550 
551  //-- create directory
552  ostringstream dptr;
553 
554  string probe_name = pdglib->Find(gOptProbePdgCode)->GetName();
555  string tgt_name = (gOptTgtPdgCode==1000000010) ?
556  "n" : pdglib->Find(gOptTgtPdgCode)->GetName();
557 
558  dptr << probe_name << "_" << tgt_name;
559  ostringstream dtitle;
560  dtitle << "Cross sections for: "
561  << pdglib->Find(gOptProbePdgCode)->GetName() << "+"
562  << pdglib->Find(gOptTgtPdgCode)->GetName();
563 
564  LOG("gspl2root", pINFO)
565  << "Will store graphs in root directory = " << dptr.str();
566  TDirectory * topdir =
567  dynamic_cast<TDirectory *> (froot->Get(dptr.str().c_str()));
568  if(topdir) {
569  LOG("gspl2root", pINFO)
570  << "Directory: " << dptr.str() << " already exists!! Exiting";
571  froot->Close();
572  delete froot;
573  return;
574  }
575 
576  topdir = froot->mkdir(dptr.str().c_str(),dtitle.str().c_str());
577  topdir->cd();
578 
579  double de = (gInlogE) ? (TMath::Log(gEmax)-TMath::Log(gEmin))/(kNSplineP-1) : (gEmax-gEmin)/(kNSplineP-1);
580  double * e = new double[kNSplineP];
581  for(int i=0; i<kNSplineP; i++) { e[i] = (gInlogE) ? TMath::Exp(TMath::Log(gEmin) + i*de) : gEmin + i*de; }
582 
583  double * xs = new double[kNSplineP];
584 
585  InteractionList::const_iterator ilistiter = ilist->begin();
586 
587  for(; ilistiter != ilist->end(); ++ilistiter) {
588 
589  const Interaction * interaction = *ilistiter;
590 
591  const ProcessInfo & proc = interaction->ProcInfo();
592  const XclsTag & xcls = interaction->ExclTag();
593  const InitialState & init = interaction->InitState();
594  const Target & tgt = init.Tgt();
595 
596  // graph title
597  ostringstream title;
598 
599  if (proc.IsQuasiElastic() ) { title << "qel"; }
600  else if (proc.IsMEC() ) { title << "mec"; }
601  else if (proc.IsResonant() ) {
602  title << "res";
603  if ( ! xcls.KnownResonance() )
604  title << "_single_pion" ;
605  }
606  else if (proc.IsDeepInelastic() ) { title << "dis"; }
607  else if (proc.IsDiffractive() ) { title << "dfr"; }
608  else if (proc.IsCoherentProduction() ) {
609  title << "coh";
610  if ( xcls.NSingleGammas() > 0 ) title << "_gamma" ;
611  else if ( xcls.NPions() > 0 ) title << "_pion" ;
612  else if ( xcls.NRhos() > 0 ) title << "_rho" ;
613  else title << "_other" ;
614  }
615  else if (proc.IsCoherentElastic() ) { title << "cevns"; }
616  else if (proc.IsInverseMuDecay() ) { title << "imd"; }
617  else if (proc.IsIMDAnnihilation() ) { title << "imdanh";}
618  else if (proc.IsNuElectronElastic()) { title << "ve"; }
619  else if (proc.IsGlashowResonance() ) { title << "glres"; }
620  else if (proc.IsPhotonResonance() ) { title << "phres"; }
621  else if (proc.IsPhotonCoherent() ) { title << "phcoh"; }
622  else {
623  LOG("gspl2root", pWARN) << "Process " << proc
624  << " scattering type not recognised: spline not added" ;
625  continue; }
626 
627  if (proc.IsWeakCC()) { title << "_cc"; }
628  else if (proc.IsWeakNC()) { title << "_nc"; }
629  else if (proc.IsWeakMix()) { title << "_ccncmix"; }
630  else if (proc.IsEM() ) { title << "_em"; }
631  else if (proc.IsDarkNeutralCurrent() ) { title << "_dark"; }
632  else {
633  LOG("gspl2root", pWARN) << "Process " << proc
634  << " interaction type has not recongnised: spline not added " ;
635  continue; }
636 
637  if(tgt.HitNucIsSet()) {
638  int hitnuc = tgt.HitNucPdg();
639  if ( pdg::IsProton (hitnuc) ) { title << "_p"; }
640  else if ( pdg::IsNeutron(hitnuc) ) { title << "_n"; }
641  else if ( pdg::Is2NucleonCluster(hitnuc) )
642  {
643  if (hitnuc == kPdgClusterNN) { title << "_nn"; }
644  else if (hitnuc == kPdgClusterNP) { title << "_np"; }
645  else if (hitnuc == kPdgClusterPP) { title << "_pp"; }
646  else {
647  LOG("gspl2root", pWARN) << "Can't handle hit 2-nucleon cluster PDG = " << hitnuc;
648  }
649  }
650  else {
651  LOG("gspl2root", pWARN) << "Can't handle hit nucleon PDG = " << hitnuc;
652  }
653 
654  if(tgt.HitQrkIsSet()) {
655  int qrkpdg = tgt.HitQrkPdg();
656  bool insea = tgt.HitSeaQrk();
657 
658  if ( pdg::IsUQuark(qrkpdg) ) { title << "_u"; }
659  else if ( pdg::IsDQuark(qrkpdg) ) { title << "_d"; }
660  else if ( pdg::IsSQuark(qrkpdg) ) { title << "_s"; }
661  else if ( pdg::IsCQuark(qrkpdg) ) { title << "_c"; }
662  else if ( pdg::IsBQuark(qrkpdg) ) { title << "_b"; }
663  else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
664  else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
665  else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
666  else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
667  else if ( pdg::IsAntiBQuark(qrkpdg) ) { title << "_bbar"; }
668 
669  if(insea) { title << "sea"; }
670  else { title << "val"; }
671  }
672  }
673  if(proc.IsResonant()) {
674  if ( xcls.KnownResonance() ) {
675  Resonance_t res = xcls.Resonance();
676  string resname = res::AsString(res);
677  resname = str::FilterString(")", resname);
678  resname = str::FilterString("(", resname);
679  title << "_" << resname.substr(3,4) << resname.substr(0,3);
680  }
681  else if ( xcls.NPions() == 1 ) {
682  // we are in the case of the single pion production
683  // since the hiht nucleon is known and we also know if
684  // it is CC or NC, the only missing information to identify the channel
685  // is only the flavour of the pion
686  string channel ;
687  if ( xcls.NPiPlus() == 1 ) {
688  channel = "pip" ;
689  }
690  else if ( xcls.NPiMinus() == 1 ) {
691  channel = "pim" ;
692  }
693  else if ( xcls.NPi0() == 1 ) {
694  channel = "pi0" ;
695  }
696 
697  title << '_' << channel ;
698  } // single resonsance or single pion
699  } // resonance case
700 
701  if(xcls.IsStrangeEvent()) {
702  title << "_strange";
703  if(!xcls.IsInclusiveStrange()) { title << xcls.StrangeHadronPdg(); }
704  }
705 
706  if(xcls.IsCharmEvent()) {
707  title << "_charm";
708  if(!xcls.IsInclusiveCharm()) { title << xcls.CharmHadronPdg(); }
709  }
710 
711  if(xcls.IsFinalQuarkEvent()) {
712  int qrkpdg = xcls.FinalQuarkPdg();
713  if ( pdg::IsUQuark(qrkpdg) ) { title << "_u"; }
714  else if ( pdg::IsDQuark(qrkpdg) ) { title << "_d"; }
715  else if ( pdg::IsSQuark(qrkpdg) ) { title << "_s"; }
716  else if ( pdg::IsCQuark(qrkpdg) ) { title << "_c"; }
717  else if ( pdg::IsBQuark(qrkpdg) ) { title << "_b"; }
718  else if ( pdg::IsTQuark(qrkpdg) ) { title << "_t"; }
719  else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
720  else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
721  else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
722  else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
723  else if ( pdg::IsAntiBQuark(qrkpdg) ) { title << "_bbar"; }
724  else if ( pdg::IsAntiTQuark(qrkpdg) ) { title << "_tbar"; }
725  }
726  if(xcls.IsFinalLeptonEvent()) {
727  int leppdg = TMath::Abs(xcls.FinalLeptonPdg());
728  if ( pdg::IsMuon(leppdg) ) { title << "_mu"; }
729  else if ( pdg::IsElectron(leppdg) ) { title << "_e"; }
730  else if ( pdg::IsTau(leppdg) ) { title << "_tau"; }
731  else if ( pdg::IsPion(leppdg) ) { title << "_had"; }
732  }
733 
734  const Spline * spl = evg_driver.XSecSpline(interaction);
735  for(int i=0; i<kNSplineP; i++) {
736  xs[i] = spl->Evaluate(e[i]) * (1E+38/units::cm2);
737  }
738 
739  TGraph * gr = new TGraph(kNSplineP, e, xs);
740  gr->SetName(title.str().c_str());
741  FormatXSecGraph(gr);
742  gr->SetTitle(spl->GetName());
743 
744  topdir->Add(gr);
745  }
746 
747 
748  //
749  // totals for (anti-)neutrino scattering
750  //
751 
752  bool is_neutrino = pdg::IsNeutralLepton(gOptProbePdgCode);
753 
754  if(is_neutrino) {
755 
756  //
757  // add-up all res channels
758  //
759 
760  double * xsresccp = new double[kNSplineP];
761  double * xsresccn = new double[kNSplineP];
762  double * xsresncp = new double[kNSplineP];
763  double * xsresncn = new double[kNSplineP];
764  for(int i=0; i<kNSplineP; i++) {
765  xsresccp[i] = 0;
766  xsresccn[i] = 0;
767  xsresncp[i] = 0;
768  xsresncn[i] = 0;
769  }
770 
771  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
772  const Interaction * interaction = *ilistiter;
773  const ProcessInfo & proc = interaction->ProcInfo();
774  const InitialState & init = interaction->InitState();
775  const Target & tgt = init.Tgt();
776 
777  const Spline * spl = evg_driver.XSecSpline(interaction);
778 
779  if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
780  for(int i=0; i<kNSplineP; i++) {
781  xsresccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
782  }
783  }
784  if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
785  for(int i=0; i<kNSplineP; i++) {
786  xsresccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
787  }
788  }
789  if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
790  for(int i=0; i<kNSplineP; i++) {
791  xsresncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
792  }
793  }
794  if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
795  for(int i=0; i<kNSplineP; i++) {
796  xsresncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
797  }
798  }
799  }
800 
801  TGraph * gr_resccp = new TGraph(kNSplineP, e, xsresccp);
802  gr_resccp->SetName("res_cc_p");
803  FormatXSecGraph(gr_resccp);
804  topdir->Add(gr_resccp);
805  TGraph * gr_resccn = new TGraph(kNSplineP, e, xsresccn);
806  gr_resccn->SetName("res_cc_n");
807  FormatXSecGraph(gr_resccn);
808  topdir->Add(gr_resccn);
809  TGraph * gr_resncp = new TGraph(kNSplineP, e, xsresncp);
810  gr_resncp->SetName("res_nc_p");
811  FormatXSecGraph(gr_resncp);
812  topdir->Add(gr_resncp);
813  TGraph * gr_resncn = new TGraph(kNSplineP, e, xsresncn);
814  gr_resncn->SetName("res_nc_n");
815  FormatXSecGraph(gr_resncn);
816  topdir->Add(gr_resncn);
817 
818  //
819  // add-up all dis channels
820  //
821 
822  double * xsdiscc = new double[kNSplineP];
823  double * xsdisnc = new double[kNSplineP];
824  double * xsdisccp = new double[kNSplineP];
825  double * xsdisccn = new double[kNSplineP];
826  double * xsdisncp = new double[kNSplineP];
827  double * xsdisncn = new double[kNSplineP];
828  for(int i=0; i<kNSplineP; i++) {
829  xsdiscc[i] = 0;
830  xsdisnc[i] = 0;
831  xsdisccp[i] = 0;
832  xsdisccn[i] = 0;
833  xsdisncp[i] = 0;
834  xsdisncn[i] = 0;
835  }
836  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
837  const Interaction * interaction = *ilistiter;
838  const ProcessInfo & proc = interaction->ProcInfo();
839  const XclsTag & xcls = interaction->ExclTag();
840  const InitialState & init = interaction->InitState();
841  const Target & tgt = init.Tgt();
842 
843  const Spline * spl = evg_driver.XSecSpline(interaction);
844 
845  if(xcls.IsCharmEvent()) continue;
846 
847  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
848  for(int i=0; i<kNSplineP; i++) {
849  xsdiscc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
850  xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
851  }
852  }
853  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
854  for(int i=0; i<kNSplineP; i++) {
855  xsdiscc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
856  xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
857  }
858  }
859  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
860  for(int i=0; i<kNSplineP; i++) {
861  xsdisnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
862  xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
863  }
864  }
865  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
866  for(int i=0; i<kNSplineP; i++) {
867  xsdisnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
868  xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
869  }
870  }
871  }
872  TGraph * gr_discc = new TGraph(kNSplineP, e, xsdiscc);
873  gr_discc->SetName("dis_cc");
874  FormatXSecGraph(gr_discc);
875  topdir->Add(gr_discc);
876  TGraph * gr_disnc = new TGraph(kNSplineP, e, xsdisnc);
877  gr_disnc->SetName("dis_nc");
878  FormatXSecGraph(gr_disnc);
879  topdir->Add(gr_disnc);
880  TGraph * gr_disccp = new TGraph(kNSplineP, e, xsdisccp);
881  gr_disccp->SetName("dis_cc_p");
882  FormatXSecGraph(gr_disccp);
883  topdir->Add(gr_disccp);
884  TGraph * gr_disccn = new TGraph(kNSplineP, e, xsdisccn);
885  gr_disccn->SetName("dis_cc_n");
886  FormatXSecGraph(gr_disccn);
887  topdir->Add(gr_disccn);
888  TGraph * gr_disncp = new TGraph(kNSplineP, e, xsdisncp);
889  gr_disncp->SetName("dis_nc_p");
890  FormatXSecGraph(gr_disncp);
891  topdir->Add(gr_disncp);
892  TGraph * gr_disncn = new TGraph(kNSplineP, e, xsdisncn);
893  gr_disncn->SetName("dis_nc_n");
894  FormatXSecGraph(gr_disncn);
895  topdir->Add(gr_disncn);
896 
897  //
898  // add-up all charm dis channels
899  //
900 
901  for(int i=0; i<kNSplineP; i++) {
902  xsdiscc[i] = 0;
903  xsdisnc[i] = 0;
904  xsdisccp[i] = 0;
905  xsdisccn[i] = 0;
906  xsdisncp[i] = 0;
907  xsdisncn[i] = 0;
908  }
909  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
910  const Interaction * interaction = *ilistiter;
911  const ProcessInfo & proc = interaction->ProcInfo();
912  const XclsTag & xcls = interaction->ExclTag();
913  const InitialState & init = interaction->InitState();
914  const Target & tgt = init.Tgt();
915 
916  const Spline * spl = evg_driver.XSecSpline(interaction);
917 
918  if(!xcls.IsCharmEvent()) continue;
919 
920  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
921  for(int i=0; i<kNSplineP; i++) {
922  xsdiscc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
923  xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
924  }
925  }
926  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
927  for(int i=0; i<kNSplineP; i++) {
928  xsdiscc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
929  xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
930  }
931  }
932  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
933  for(int i=0; i<kNSplineP; i++) {
934  xsdisnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
935  xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
936  }
937  }
938  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
939  for(int i=0; i<kNSplineP; i++) {
940  xsdisnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
941  xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
942  }
943  }
944  }
945  TGraph * gr_discc_charm = new TGraph(kNSplineP, e, xsdiscc);
946  gr_discc_charm->SetName("dis_cc_charm");
947  FormatXSecGraph(gr_discc_charm);
948  topdir->Add(gr_discc_charm);
949  TGraph * gr_disnc_charm = new TGraph(kNSplineP, e, xsdisnc);
950  gr_disnc_charm->SetName("dis_nc_charm");
951  FormatXSecGraph(gr_disnc_charm);
952  topdir->Add(gr_disnc_charm);
953  TGraph * gr_disccp_charm = new TGraph(kNSplineP, e, xsdisccp);
954  gr_disccp_charm->SetName("dis_cc_p_charm");
955  FormatXSecGraph(gr_disccp_charm);
956  topdir->Add(gr_disccp_charm);
957  TGraph * gr_disccn_charm = new TGraph(kNSplineP, e, xsdisccn);
958  gr_disccn_charm->SetName("dis_cc_n_charm");
959  FormatXSecGraph(gr_disccn_charm);
960  topdir->Add(gr_disccn_charm);
961  TGraph * gr_disncp_charm = new TGraph(kNSplineP, e, xsdisncp);
962  gr_disncp_charm->SetName("dis_nc_p_charm");
963  FormatXSecGraph(gr_disncp_charm);
964  topdir->Add(gr_disncp_charm);
965  TGraph * gr_disncn_charm = new TGraph(kNSplineP, e, xsdisncn);
966  gr_disncn_charm->SetName("dis_nc_n_charm");
967  FormatXSecGraph(gr_disncn_charm);
968  topdir->Add(gr_disncn_charm);
969 
970  //
971  // add-up all mec channels
972  //
973 
974  double * xsmeccc = new double[kNSplineP];
975  double * xsmecnc = new double[kNSplineP];
976  for(int i=0; i<kNSplineP; i++) {
977  xsmeccc[i] = 0;
978  xsmecnc[i] = 0;
979  }
980 
981  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
982  const Interaction * interaction = *ilistiter;
983  const ProcessInfo & proc = interaction->ProcInfo();
984 
985  const Spline * spl = evg_driver.XSecSpline(interaction);
986 
987  if (proc.IsMEC() && proc.IsWeakCC()) {
988  for(int i=0; i<kNSplineP; i++) {
989  xsmeccc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
990  }
991  }
992  if (proc.IsMEC() && proc.IsWeakNC()) {
993  for(int i=0; i<kNSplineP; i++) {
994  xsmecnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
995  }
996  }
997  }
998 
999  TGraph * gr_meccc = new TGraph(kNSplineP, e, xsmeccc);
1000  gr_meccc->SetName("mec_cc");
1001  FormatXSecGraph(gr_meccc);
1002  topdir->Add(gr_meccc);
1003  TGraph * gr_mecnc = new TGraph(kNSplineP, e, xsmecnc);
1004  gr_mecnc->SetName("mec_nc");
1005  FormatXSecGraph(gr_mecnc);
1006  topdir->Add(gr_mecnc);
1007 
1008  //
1009  // add-up all COH channels
1010  //
1011 
1012  double * xscohcc = new double[kNSplineP];
1013  double * xscohnc = new double[kNSplineP];
1014  double * xscohtot = new double[kNSplineP];
1015  for(int i=0; i<kNSplineP; i++) {
1016  xscohcc[i] = 0;
1017  xscohnc[i] = 0;
1018  xscohtot[i] = 0;
1019  }
1020 
1021  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1022 
1023  const Interaction * interaction = *ilistiter;
1024  const ProcessInfo & proc = interaction->ProcInfo();
1025 
1026  const Spline * spl = evg_driver.XSecSpline(interaction);
1027 
1028  if (proc.IsCoherentProduction() && proc.IsWeakCC()) {
1029  for(int i=0; i<kNSplineP; i++) {
1030  xscohcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1031  }
1032  }
1033  if (proc.IsCoherentProduction() && proc.IsWeakNC()) {
1034  for(int i=0; i<kNSplineP; i++) {
1035  xscohnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1036  }
1037  }
1038  if ( proc.IsCoherentProduction() ) {
1039  for(int i=0; i<kNSplineP; i++) {
1040  xscohtot[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1041  }
1042  }
1043 
1044  }
1045 
1046  TGraph * gr_cohcc = new TGraph(kNSplineP, e, xscohcc);
1047  gr_cohcc->SetName("coh_cc");
1048  FormatXSecGraph(gr_cohcc);
1049  topdir->Add(gr_cohcc);
1050 
1051  TGraph * gr_cohnc = new TGraph(kNSplineP, e, xscohnc);
1052  gr_cohnc->SetName("coh_nc");
1053  FormatXSecGraph(gr_cohnc);
1054  topdir->Add(gr_cohnc);
1055 
1056  TGraph * gr_cohtot = new TGraph(kNSplineP, e, xscohtot);
1057  gr_cohtot->SetName("coh");
1058  FormatXSecGraph(gr_cohtot);
1059  topdir->Add(gr_cohtot);
1060 
1061  //
1062  // add-up all glres and photon-res channels
1063  //
1064 
1065  double * xsglrescc = new double[kNSplineP];
1066  double * xsglresnc = new double[kNSplineP];
1067  double * xsphrescc = new double[kNSplineP];
1068  for(int i=0; i<kNSplineP; i++) {
1069  xsglrescc[i] = 0;
1070  xsglresnc[i] = 0;
1071  xsphrescc[i] = 0;
1072  }
1073  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1074  const Interaction * interaction = *ilistiter;
1075  const ProcessInfo & proc = interaction->ProcInfo();
1076 
1077  const Spline * spl = evg_driver.XSecSpline(interaction);
1078 
1079  if (proc.IsGlashowResonance() && proc.IsWeakCC()) {
1080  for(int i=0; i<kNSplineP; i++) {
1081  xsglrescc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1082  }
1083  }
1084  if (proc.IsGlashowResonance() && proc.IsWeakNC()) {
1085  for(int i=0; i<kNSplineP; i++) {
1086  xsglresnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1087  }
1088  }
1089  if (proc.IsPhotonResonance()) {
1090  for(int i=0; i<kNSplineP; i++) {
1091  xsphrescc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1092  }
1093 
1094  }
1095  }
1096  TGraph * gr_glrescc = new TGraph(kNSplineP, e, xsglrescc);
1097  gr_glrescc->SetName("glres_cc");
1098  FormatXSecGraph(gr_glrescc);
1099  topdir->Add(gr_glrescc);
1100  TGraph * gr_glresnc = new TGraph(kNSplineP, e, xsglresnc);
1101  gr_glresnc->SetName("glres_nc");
1102  FormatXSecGraph(gr_glresnc);
1103  topdir->Add(gr_glresnc);
1104  TGraph * gr_phrescc = new TGraph(kNSplineP, e, xsphrescc);
1105  gr_phrescc->SetName("phres_cc");
1106  FormatXSecGraph(gr_phrescc);
1107  topdir->Add(gr_phrescc);
1108 
1109 
1110  //
1111  // total cross sections
1112  //
1113  double * xstotcc = new double[kNSplineP];
1114  double * xstotccp = new double[kNSplineP];
1115  double * xstotccn = new double[kNSplineP];
1116  double * xstotnc = new double[kNSplineP];
1117  double * xstotncp = new double[kNSplineP];
1118  double * xstotncn = new double[kNSplineP];
1119  for(int i=0; i<kNSplineP; i++) {
1120  xstotcc [i] = 0;
1121  xstotccp[i] = 0;
1122  xstotccn[i] = 0;
1123  xstotnc [i] = 0;
1124  xstotncp[i] = 0;
1125  xstotncn[i] = 0;
1126  }
1127  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1128  const Interaction * interaction = *ilistiter;
1129  const ProcessInfo & proc = interaction->ProcInfo();
1130  const InitialState & init = interaction->InitState();
1131  const Target & tgt = init.Tgt();
1132 
1133  const Spline * spl = evg_driver.XSecSpline(interaction);
1134 
1135  bool iscc = proc.IsWeakCC();
1136  bool isnc = proc.IsWeakNC();
1137  bool offp = pdg::IsProton (tgt.HitNucPdg());
1138  bool offn = pdg::IsNeutron(tgt.HitNucPdg());
1139 
1140  if (iscc && offp) {
1141  for(int i=0; i<kNSplineP; i++) {
1142  xstotccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1143  }
1144  }
1145  if (iscc && offn) {
1146  for(int i=0; i<kNSplineP; i++) {
1147  xstotccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1148  }
1149  }
1150  if (isnc && offp) {
1151  for(int i=0; i<kNSplineP; i++) {
1152  xstotncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1153  }
1154  }
1155  if (isnc && offn) {
1156  for(int i=0; i<kNSplineP; i++) {
1157  xstotncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1158  }
1159  }
1160 
1161  if (iscc) {
1162  for(int i=0; i<kNSplineP; i++) {
1163  xstotcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1164  }
1165  }
1166  if (isnc) {
1167  for(int i=0; i<kNSplineP; i++) {
1168  xstotnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1169  }
1170  }
1171  }
1172 
1173  TGraph * gr_totcc = new TGraph(kNSplineP, e, xstotcc);
1174  gr_totcc->SetName("tot_cc");
1175  FormatXSecGraph(gr_totcc);
1176  topdir->Add(gr_totcc);
1177  TGraph * gr_totccp = new TGraph(kNSplineP, e, xstotccp);
1178  gr_totccp->SetName("tot_cc_p");
1179  FormatXSecGraph(gr_totccp);
1180  topdir->Add(gr_totccp);
1181  TGraph * gr_totccn = new TGraph(kNSplineP, e, xstotccn);
1182  gr_totccn->SetName("tot_cc_n");
1183  FormatXSecGraph(gr_totccn);
1184  topdir->Add(gr_totccn);
1185  TGraph * gr_totnc = new TGraph(kNSplineP, e, xstotnc);
1186  gr_totnc->SetName("tot_nc");
1187  FormatXSecGraph(gr_totnc);
1188  topdir->Add(gr_totnc);
1189  TGraph * gr_totncp = new TGraph(kNSplineP, e, xstotncp);
1190  gr_totncp->SetName("tot_nc_p");
1191  FormatXSecGraph(gr_totncp);
1192  topdir->Add(gr_totncp);
1193  TGraph * gr_totncn = new TGraph(kNSplineP, e, xstotncn);
1194  gr_totncn->SetName("tot_nc_n");
1195  FormatXSecGraph(gr_totncn);
1196  topdir->Add(gr_totncn);
1197 
1198  delete [] e;
1199  delete [] xs;
1200  delete [] xsresccp;
1201  delete [] xsresccn;
1202  delete [] xsresncp;
1203  delete [] xsresncn;
1204  delete [] xsdisccp;
1205  delete [] xsdisccn;
1206  delete [] xsdisncp;
1207  delete [] xsdisncn;
1208  delete [] xsdiscc;
1209  delete [] xsdisnc;
1210  delete [] xscohcc;
1211  delete [] xscohnc;
1212  delete [] xscohtot;
1213  delete [] xsglrescc;
1214  delete [] xsglresnc;
1215  delete [] xsphrescc;
1216  delete [] xstotcc;
1217  delete [] xstotccp;
1218  delete [] xstotccn;
1219  delete [] xstotnc;
1220  delete [] xstotncp;
1221  delete [] xstotncn;
1222 
1223  }// neutrinos
1224 
1225 
1226  //
1227  // totals for charged lepton scattering
1228  //
1229 
1230  bool is_charged_lepton = pdg::IsChargedLepton(gOptProbePdgCode);
1231 
1232  if(is_charged_lepton) {
1233 
1234  //
1235  // add-up all res channels
1236  //
1237 
1238  double * xsresemp = new double[kNSplineP];
1239  double * xsresemn = new double[kNSplineP];
1240  for(int i=0; i<kNSplineP; i++) {
1241  xsresemp[i] = 0;
1242  xsresemn[i] = 0;
1243  }
1244 
1245  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1246  const Interaction * interaction = *ilistiter;
1247  const ProcessInfo & proc = interaction->ProcInfo();
1248  const InitialState & init = interaction->InitState();
1249  const Target & tgt = init.Tgt();
1250 
1251  const Spline * spl = evg_driver.XSecSpline(interaction);
1252 
1253  if (proc.IsResonant() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1254  for(int i=0; i<kNSplineP; i++) {
1255  xsresemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1256  }
1257  }
1258  if (proc.IsResonant() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1259  for(int i=0; i<kNSplineP; i++) {
1260  xsresemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1261  }
1262  }
1263  }
1264 
1265  TGraph * gr_resemp = new TGraph(kNSplineP, e, xsresemp);
1266  gr_resemp->SetName("res_em_p");
1267  FormatXSecGraph(gr_resemp);
1268  topdir->Add(gr_resemp);
1269  TGraph * gr_resemn = new TGraph(kNSplineP, e, xsresemn);
1270  gr_resemn->SetName("res_em_n");
1271  FormatXSecGraph(gr_resemn);
1272  topdir->Add(gr_resemn);
1273 
1274  //
1275  // add-up all dis channels
1276  //
1277 
1278  double * xsdisemp = new double[kNSplineP];
1279  double * xsdisemn = new double[kNSplineP];
1280  for(int i=0; i<kNSplineP; i++) {
1281  xsdisemp[i] = 0;
1282  xsdisemn[i] = 0;
1283  }
1284  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1285  const Interaction * interaction = *ilistiter;
1286  const ProcessInfo & proc = interaction->ProcInfo();
1287  const XclsTag & xcls = interaction->ExclTag();
1288  const InitialState & init = interaction->InitState();
1289  const Target & tgt = init.Tgt();
1290 
1291  const Spline * spl = evg_driver.XSecSpline(interaction);
1292 
1293  if(xcls.IsCharmEvent()) continue;
1294 
1295  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1296  for(int i=0; i<kNSplineP; i++) {
1297  xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1298  }
1299  }
1300  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1301  for(int i=0; i<kNSplineP; i++) {
1302  xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1303  }
1304  }
1305  }
1306  TGraph * gr_disemp = new TGraph(kNSplineP, e, xsdisemp);
1307  gr_disemp->SetName("dis_em_p");
1308  FormatXSecGraph(gr_disemp);
1309  topdir->Add(gr_disemp);
1310  TGraph * gr_disemn = new TGraph(kNSplineP, e, xsdisemn);
1311  gr_disemn->SetName("dis_em_n");
1312  FormatXSecGraph(gr_disemn);
1313  topdir->Add(gr_disemn);
1314 
1315  //
1316  // add-up all charm dis channels
1317  //
1318 
1319  for(int i=0; i<kNSplineP; i++) {
1320  xsdisemp[i] = 0;
1321  xsdisemn[i] = 0;
1322  }
1323  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1324  const Interaction * interaction = *ilistiter;
1325  const ProcessInfo & proc = interaction->ProcInfo();
1326  const XclsTag & xcls = interaction->ExclTag();
1327  const InitialState & init = interaction->InitState();
1328  const Target & tgt = init.Tgt();
1329 
1330  const Spline * spl = evg_driver.XSecSpline(interaction);
1331 
1332  if(!xcls.IsCharmEvent()) continue;
1333 
1334  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1335  for(int i=0; i<kNSplineP; i++) {
1336  xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1337  }
1338  }
1339  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1340  for(int i=0; i<kNSplineP; i++) {
1341  xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1342  }
1343  }
1344  }
1345  TGraph * gr_disemp_charm = new TGraph(kNSplineP, e, xsdisemp);
1346  gr_disemp_charm->SetName("dis_em_p_charm");
1347  FormatXSecGraph(gr_disemp_charm);
1348  topdir->Add(gr_disemp_charm);
1349  TGraph * gr_disemn_charm = new TGraph(kNSplineP, e, xsdisemn);
1350  gr_disemn_charm->SetName("dis_em_n_charm");
1351  FormatXSecGraph(gr_disemn_charm);
1352  topdir->Add(gr_disemn_charm);
1353 
1354  //
1355  // total cross sections
1356  //
1357  double * xstotem = new double[kNSplineP];
1358  double * xstotemp = new double[kNSplineP];
1359  double * xstotemn = new double[kNSplineP];
1360  for(int i=0; i<kNSplineP; i++) {
1361  xstotem [i] = 0;
1362  xstotemp[i] = 0;
1363  xstotemn[i] = 0;
1364  }
1365  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1366  const Interaction * interaction = *ilistiter;
1367  const ProcessInfo & proc = interaction->ProcInfo();
1368  const InitialState & init = interaction->InitState();
1369  const Target & tgt = init.Tgt();
1370 
1371  const Spline * spl = evg_driver.XSecSpline(interaction);
1372 
1373  bool isem = proc.IsEM();
1374  bool offp = pdg::IsProton (tgt.HitNucPdg());
1375  bool offn = pdg::IsNeutron(tgt.HitNucPdg());
1376 
1377  if (isem && offp) {
1378  for(int i=0; i<kNSplineP; i++) {
1379  xstotemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1380  }
1381  }
1382  if (isem && offn) {
1383  for(int i=0; i<kNSplineP; i++) {
1384  xstotemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1385  }
1386  }
1387  if (isem) {
1388  for(int i=0; i<kNSplineP; i++) {
1389  xstotem[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1390  }
1391  }
1392  }
1393 
1394  TGraph * gr_totem = new TGraph(kNSplineP, e, xstotem);
1395  gr_totem->SetName("tot_em");
1396  FormatXSecGraph(gr_totem);
1397  topdir->Add(gr_totem);
1398  TGraph * gr_totemp = new TGraph(kNSplineP, e, xstotemp);
1399  gr_totemp->SetName("tot_em_p");
1400  FormatXSecGraph(gr_totemp);
1401  topdir->Add(gr_totemp);
1402  TGraph * gr_totemn = new TGraph(kNSplineP, e, xstotemn);
1403  gr_totemn->SetName("tot_em_n");
1404  FormatXSecGraph(gr_totemn);
1405  topdir->Add(gr_totemn);
1406 
1407  delete [] e;
1408  delete [] xs;
1409  delete [] xsresemp;
1410  delete [] xsresemn;
1411  delete [] xsdisemp;
1412  delete [] xsdisemn;
1413  delete [] xstotem;
1414  delete [] xstotemp;
1415  delete [] xstotemn;
1416 
1417  }// charged leptons
1418 
1419  topdir->Write();
1420 
1421  if(froot) {
1422  froot->Close();
1423  delete froot;
1424  }
1425 }
genie::ProcessInfo::IsResonant
bool IsResonant(void) const
Definition: ProcessInfo.cxx:99
genie::pdg::IsBQuark
bool IsBQuark(int pdgc)
Definition: PDGUtils.cxx:286
genie::ProcessInfo::IsPhotonResonance
bool IsPhotonResonance(void) const
Definition: ProcessInfo.cxx:166
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
genie::XclsTag::IsFinalLeptonEvent
bool IsFinalLeptonEvent(void) const
Definition: XclsTag.h:73
genie::ProcessInfo::IsWeakMix
bool IsWeakMix(void) const
Definition: ProcessInfo.cxx:218
genie::Target::HitSeaQrk
bool HitSeaQrk(void) const
Definition: Target.cxx:299
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::XclsTag::IsFinalQuarkEvent
bool IsFinalQuarkEvent(void) const
Definition: XclsTag.h:71
genie::XclsTag::FinalLeptonPdg
int FinalLeptonPdg(void) const
Definition: XclsTag.h:74
genie::pdg::IsUQuark
bool IsUQuark(int pdgc)
Definition: PDGUtils.cxx:266
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::ProcessInfo::IsDarkNeutralCurrent
bool IsDarkNeutralCurrent(void) const
Definition: ProcessInfo.cxx:146
genie::Target::HitQrkPdg
int HitQrkPdg(void) const
Definition: Target.cxx:242
genie::XclsTag::NPions
int NPions(void) const
Definition: XclsTag.h:62
genie::ProcessInfo::IsInverseMuDecay
bool IsInverseMuDecay(void) const
Definition: ProcessInfo.cxx:126
genie::kPdgClusterNP
const int kPdgClusterNP
Definition: PDGCodes.h:215
genie::ProcessInfo::IsQuasiElastic
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
FormatXSecGraph
void FormatXSecGraph(TGraph *g)
Definition: gSplineXml2Root.cxx:516
genie::XclsTag::KnownResonance
bool KnownResonance(void) const
Definition: XclsTag.h:68
genie::Spline
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
genie::XclsTag::CharmHadronPdg
int CharmHadronPdg(void) const
Definition: XclsTag.h:52
genie::XclsTag::IsStrangeEvent
bool IsStrangeEvent(void) const
Definition: XclsTag.h:53
genie::XclsTag::NPiMinus
int NPiMinus(void) const
Definition: XclsTag.h:60
genie::pdg::IsChargedLepton
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
gInlogE
bool gInlogE
Definition: gSplineXml2Root.cxx:152
genie::pdg::IsTQuark
bool IsTQuark(int pdgc)
Definition: PDGUtils.cxx:291
genie::pdg::IsSQuark
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:276
genie::ProcessInfo::IsDiffractive
bool IsDiffractive(void) const
Definition: ProcessInfo.cxx:181
genie::ProcessInfo::IsCoherentProduction
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:104
genie::ProcessInfo::IsIMDAnnihilation
bool IsIMDAnnihilation(void) const
Definition: ProcessInfo.cxx:131
genie::pdg::IsAntiSQuark
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:306
genie::Spline::Evaluate
double Evaluate(double x) const
Definition: Spline.cxx:363
genie::Resonance_t
enum genie::EResonance Resonance_t
genie::pdg::IsAntiDQuark
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:301
gEmin
double gEmin
Definition: gSplineXml2Root.cxx:150
genie::XclsTag
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
genie::kPdgClusterNN
const int kPdgClusterNN
Definition: PDGCodes.h:214
genie::XclsTag::IsCharmEvent
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
gOptROOTFilename
string gOptROOTFilename
Definition: gSplineXml2Root.cxx:141
genie::XclsTag::NSingleGammas
int NSingleGammas(void) const
Definition: XclsTag.h:64
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition: ProcessInfo.cxx:208
e
const double e
Definition: gUpMuFluxGen.cxx:165
genie::XclsTag::StrangeHadronPdg
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::ProcessInfo::IsCoherentElastic
bool IsCoherentElastic(void) const
Definition: ProcessInfo.cxx:109
genie::ProcessInfo::IsNuElectronElastic
bool IsNuElectronElastic(void) const
Definition: ProcessInfo.cxx:121
genie::GEVGDriver::Interactions
const InteractionList * Interactions(void) const
Definition: GEVGDriver.cxx:334
genie::units::cm2
static constexpr double cm2
Definition: Units.h:69
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::pdg::IsAntiTQuark
bool IsAntiTQuark(int pdgc)
Definition: PDGUtils.cxx:321
genie::XclsTag::NPiPlus
int NPiPlus(void) const
Definition: XclsTag.h:59
genie::pdg::IsTau
bool IsTau(int pdgc)
Definition: PDGUtils.cxx:208
genie::pdg::IsAntiBQuark
bool IsAntiBQuark(int pdgc)
Definition: PDGUtils.cxx:316
genie::GEVGDriver
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:54
pINFO
#define pINFO
Definition: Messenger.h:62
gOptProbePdgCode
int gOptProbePdgCode
Definition: gEvGenHadronNucleus.cxx:146
genie::XclsTag::NPi0
int NPi0(void) const
Definition: XclsTag.h:58
genie::pdg::IsMuon
bool IsMuon(int pdgc)
Definition: PDGUtils.cxx:198
GetEventGenDriver
GEVGDriver GetEventGenDriver(void)
Definition: gSplineXml2Root.cxx:200
genie::XclsTag::Resonance
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
genie::ProcessInfo::IsPhotonCoherent
bool IsPhotonCoherent(void) const
Definition: ProcessInfo.cxx:161
pWARN
#define pWARN
Definition: Messenger.h:60
gEmax
double gEmax
Definition: gSplineXml2Root.cxx:151
genie::ProcessInfo::IsMEC
bool IsMEC(void) const
Definition: ProcessInfo.cxx:176
genie::ProcessInfo::IsEM
bool IsEM(void) const
Definition: ProcessInfo.cxx:193
genie::XclsTag::FinalQuarkPdg
int FinalQuarkPdg(void) const
Definition: XclsTag.h:72
genie::pdg::IsNeutralLepton
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:95
genie::ProcessInfo::IsDeepInelastic
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
genie::pdg::Is2NucleonCluster
bool Is2NucleonCluster(int pdgc)
Definition: PDGUtils.cxx:402
genie::utils::str::FilterString
string FilterString(string filt, string input)
Definition: StringUtils.cxx:79
genie::XclsTag::IsInclusiveStrange
bool IsInclusiveStrange(void) const
Definition: XclsTag.cxx:71
genie::pdg::IsCQuark
bool IsCQuark(int pdgc)
Definition: PDGUtils.cxx:281
genie::GEVGDriver::XSecSpline
const Spline * XSecSpline(const Interaction *interaction) const
Definition: GEVGDriver.cxx:488
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::Target::HitNucIsSet
bool HitNucIsSet(void) const
Definition: Target.cxx:283
genie::Target::HitQrkIsSet
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
genie::pdg::IsAntiCQuark
bool IsAntiCQuark(int pdgc)
Definition: PDGUtils.cxx:311
genie::Interaction::ExclTag
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
genie::XclsTag::IsInclusiveCharm
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:54
genie::pdg::IsDQuark
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:271
genie::InteractionList
A vector of Interaction objects.
Definition: InteractionList.h:34
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition: Interaction.h:69
genie::utils::res::AsString
const char * AsString(Resonance_t res)
resonance id -&gt; string
Definition: BaryonResUtils.cxx:35
genie::Interaction::ProcInfo
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::ProcessInfo::IsGlashowResonance
bool IsGlashowResonance(void) const
Definition: ProcessInfo.cxx:156
genie::XclsTag::NRhos
int NRhos(void) const
Definition: XclsTag.h:63
genie::pdg::IsElectron
bool IsElectron(int pdgc)
Definition: PDGUtils.cxx:188
kNSplineP
int kNSplineP
Definition: gSplineXml2Root.cxx:154
genie::pdg::IsAntiUQuark
bool IsAntiUQuark(int pdgc)
Definition: PDGUtils.cxx:296
genie::InitialState
Initial State information.
Definition: InitialState.h:48
gOptTgtPdgCode
int gOptTgtPdgCode
Definition: gEvGenHadronNucleus.cxx:147
genie::kPdgClusterPP
const int kPdgClusterPP
Definition: PDGCodes.h:216
void SaveNtupleToRootFile ( void  )

Definition at line 1427 of file gSplineXml2Root.cxx.

1428 {
1429 /*
1430  //-- create cross section ntuple
1431  //
1432  double brXSec;
1433  double brEv;
1434  bool brIsQel;
1435  bool brIsRes;
1436  bool brIsDis;
1437  bool brIsCoh;
1438  bool brIsImd;
1439  bool brIsNuEEl;
1440  bool brIsCC;
1441  bool brIsNC;
1442  int brNucleon;
1443  int brQrk;
1444  bool brIsSeaQrk;
1445  int brRes;
1446  bool brCharm;
1447  TTree * xsecnt = new TTree("xsecnt",dtitle.str().c_str());
1448  xsecnt->Branch("xsec", &brXSec, "xsec/D");
1449  xsecnt->Branch("e", &brEv, "e/D");
1450  xsecnt->Branch("qel", &brIsQel, "qel/O");
1451  xsecnt->Branch("res", &brIsRes, "res/O");
1452  xsecnt->Branch("dis", &brIsDis, "dis/O");
1453  xsecnt->Branch("coh", &brIsCoh, "coh/O");
1454  xsecnt->Branch("imd", &brIsImd, "imd/O");
1455  xsecnt->Branch("veel", &brIsNuEEl, "veel/O");
1456  xsecnt->Branch("cc", &brIsCC, "cc/O");
1457  xsecnt->Branch("nc", &brIsNC, "nc/O");
1458  xsecnt->Branch("nuc", &brNucleon, "nuc/I");
1459  xsecnt->Branch("qrk", &brQrk, "qrk/I");
1460  xsecnt->Branch("sea", &brIsSeaQrk, "sea/O");
1461  xsecnt->Branch("res", &brRes, "res/I");
1462  xsecnt->Branch("charm", &brCharm, "charm/O");
1463 */
1464 }
void SaveToPsFile ( void  )

Definition at line 217 of file gSplineXml2Root.cxx.

References genie::Interaction::AsString(), genie::units::cm2, genie::Interaction::ExclTag(), genie::PDGLibrary::Find(), gEmax, gEmin, genie::Spline::GetAsTGraph(), GetEventGenDriver(), gOptProbePdgCode, gOptTgtPdgCode, gWriteOutPlots, genie::PDGLibrary::Instance(), genie::GEVGDriver::Interactions(), genie::XclsTag::IsCharmEvent(), genie::ProcessInfo::IsCoherentProduction(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsIMDAnnihilation(), genie::ProcessInfo::IsInverseMuDecay(), genie::ProcessInfo::IsNuElectronElastic(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), kNP, kPsType, LOG, pINFO, genie::Interaction::ProcInfo(), genie::units::ps, pWARN, genie::GEVGDriver::XSecSpline(), and genie::GEVGDriver::XSecSumSpline().

Referenced by main().

218 {
219  if(!gWriteOutPlots) return;
220 
221  //-- get the event generation driver
222  GEVGDriver evg_driver = GetEventGenDriver();
223 
224  //-- define some marker styles / colors
225  const unsigned int kNMarkers = 5;
226  const unsigned int kNColors = 6;
227  unsigned int markers[kNMarkers] = {20, 28, 29, 27, 3};
228  unsigned int colors [kNColors] = {1, 2, 4, 6, 8, 28};
229 
230  //-- create a postscript document for saving cross section plpots
231 
232  TCanvas * c = new TCanvas("c","",20,20,500,850);
233  c->SetBorderMode(0);
234  c->SetFillColor(0);
235  TLegend * legend = new TLegend(0.01,0.01,0.99,0.99);
236  legend->SetFillColor(0);
237  legend->SetBorderSize(0);
238 
239  //-- get pdglibrary for mapping pdg codes to names
240  PDGLibrary * pdglib = PDGLibrary::Instance();
241  ostringstream filename;
242  filename << "xsec-splines-"
243  << pdglib->Find(gOptProbePdgCode)->GetName() << "-"
244  << pdglib->Find(gOptTgtPdgCode)->GetName() << ".ps";
245  TPostScript * ps = new TPostScript(filename.str().c_str(), kPsType);
246 
247  //-- get the list of interactions that can be simulated by the driver
248  const InteractionList * ilist = evg_driver.Interactions();
249  unsigned int nspl = ilist->size();
250 
251  //-- book enough space for xsec plots (last one is the sum)
252  TGraph * gr[nspl+1];
253 
254  //-- loop over all the simulated interactions & create the cross section graphs
255  InteractionList::const_iterator ilistiter = ilist->begin();
256  unsigned int i=0;
257  for(; ilistiter != ilist->end(); ++ilistiter) {
258 
259  const Interaction * interaction = *ilistiter;
260  LOG("gspl2root", pINFO)
261  << "Current interaction: " << interaction->AsString();
262 
263 
264  //-- access the cross section spline
265  const Spline * spl = evg_driver.XSecSpline(interaction);
266  if(!spl) {
267  LOG("gspl2root", pWARN)
268  << "Can't get spline for: " << interaction->AsString();
269  exit(2);
270  }
271 
272  //-- set graph color/style
273  int icol = TMath::Min( i % kNColors, kNColors-1 );
274  int isty = TMath::Min( i / kNMarkers, kNMarkers-1 );
275  int col = colors[icol];
276  int sty = markers[isty];
277 
278  LOG("gspl2root", pINFO)
279  << "color = " << col << ", marker = " << sty;
280 
281  //-- export Spline as TGraph / set color & stype
282  gr[i] = spl->GetAsTGraph(kNP,true,true,1.,1./units::cm2);
283  gr[i]->SetLineColor(col);
284  gr[i]->SetMarkerColor(col);
285  gr[i]->SetMarkerStyle(sty);
286  gr[i]->SetMarkerSize(0.5);
287 
288  i++;
289  }
290 
291  //-- now, get the sum...
292  const Spline * splsum = evg_driver.XSecSumSpline();
293  if(!splsum) {
294  LOG("gspl2root", pWARN)
295  << "Can't get the cross section sum spline";
296  exit(2);
297  }
298  gr[nspl] = splsum->GetAsTGraph(kNP,true,true,1.,1./units::cm2);
299 
300  //-- figure out the minimum / maximum xsec in plotted range
301  double XSmax = -9999;
302  double XSmin = 9999;
303  double x=0, y=0;
304  for(int j=0; j<kNP; j++) {
305  gr[nspl]->GetPoint(j,x,y);
306  XSmax = TMath::Max(XSmax,y);
307  }
308  XSmin = XSmax/100000.;
309  XSmax = XSmax*1.2;
310 
311  LOG("gspl2root", pINFO) << "Drawing frame: E = (" << gEmin << ", " << gEmax << ")";
312  LOG("gspl2root", pINFO) << "Drawing frame: XSec = (" << XSmin << ", " << XSmax << ")";
313 
314  //-- ps output: add the 1st page with _all_ xsec spline plots
315  //
316  //c->Draw();
317  TH1F * h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
318  for(unsigned int ispl = 0; ispl <= nspl; ispl++) if(gr[ispl]) { gr[ispl]->Draw("LP"); }
319  h->GetXaxis()->SetTitle("Ev (GeV)");
320  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
321  c->SetLogx();
322  c->SetLogy();
323  c->SetGridx();
324  c->SetGridy();
325  c->Update();
326 
327  //-- plot QEL xsecs only
328  //
329  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
330  i=0;
331  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
332  const Interaction * interaction = *ilistiter;
333  if(interaction->ProcInfo().IsQuasiElastic() && ! interaction->ExclTag().IsCharmEvent()) {
334  gr[i]->Draw("LP");
335  TString spltitle(interaction->AsString());
336  spltitle = spltitle.ReplaceAll(";",1," ",1);
337  legend->AddEntry(gr[i], spltitle.Data(),"LP");
338  }
339  i++;
340  }
341  legend->SetHeader("QEL Cross Sections");
342  gr[nspl]->Draw("LP");
343  legend->AddEntry(gr[nspl], "sum","LP");
344  h->GetXaxis()->SetTitle("Ev (GeV)");
345  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
346  c->SetLogx();
347  c->SetLogy();
348  c->SetGridx();
349  c->SetGridy();
350  c->Update();
351  c->Clear();
352  c->Range(0,0,1,1);
353  legend->Draw();
354  c->Update();
355 
356  //-- plot RES xsecs only
357  //
358  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
359  legend->Clear();
360  i=0;
361  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
362  const Interaction * interaction = *ilistiter;
363  if(interaction->ProcInfo().IsResonant()) {
364  gr[i]->Draw("LP");
365  TString spltitle(interaction->AsString());
366  spltitle = spltitle.ReplaceAll(";",1," ",1);
367  legend->AddEntry(gr[i], spltitle.Data(),"LP");
368  }
369  i++;
370  }
371  legend->SetHeader("RES Cross Sections");
372  gr[nspl]->Draw("LP");
373  legend->AddEntry(gr[nspl], "sum","LP");
374  h->GetXaxis()->SetTitle("Ev (GeV)");
375  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
376  c->SetLogx();
377  c->SetLogy();
378  c->SetGridx();
379  c->SetGridy();
380  c->Update();
381  c->Clear();
382  c->Range(0,0,1,1);
383  legend->Draw();
384  c->Update();
385 
386  //-- plot DIS xsecs only
387  //
388  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
389  legend->Clear();
390  i=0;
391  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
392  const Interaction * interaction = *ilistiter;
393  if(interaction->ProcInfo().IsDeepInelastic()) {
394  gr[i]->Draw("LP");
395  TString spltitle(interaction->AsString());
396  spltitle = spltitle.ReplaceAll(";",1," ",1);
397  legend->AddEntry(gr[i], spltitle.Data(),"LP");
398  }
399  i++;
400  }
401  legend->SetHeader("DIS Cross Sections");
402  gr[nspl]->Draw("LP");
403  legend->AddEntry(gr[nspl], "sum","LP");
404  h->GetXaxis()->SetTitle("Ev (GeV)");
405  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
406  c->SetLogx();
407  c->SetLogy();
408  c->SetGridx();
409  c->SetGridy();
410  c->Update();
411  c->Clear();
412  c->Range(0,0,1,1);
413  legend->Draw();
414  c->Update();
415 
416  //-- plot COH xsecs only
417  //
418  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
419  legend->Clear();
420  i=0;
421  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
422  const Interaction * interaction = *ilistiter;
423  if(interaction->ProcInfo().IsCoherentProduction()) {
424  gr[i]->Draw("LP");
425  TString spltitle(interaction->AsString());
426  spltitle = spltitle.ReplaceAll(";",1," ",1);
427  legend->AddEntry(gr[i], spltitle.Data(),"LP");
428  }
429  i++;
430  }
431  legend->SetHeader("COH Cross Sections");
432  gr[nspl]->Draw("LP");
433  legend->AddEntry(gr[nspl], "sum","LP");
434  h->GetXaxis()->SetTitle("Ev (GeV)");
435  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
436  c->SetLogx();
437  c->SetLogy();
438  c->SetGridx();
439  c->SetGridy();
440  c->Update();
441  c->Clear();
442  c->Range(0,0,1,1);
443  legend->Draw();
444  c->Update();
445 
446  //-- plot charm xsecs only
447  //
448  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
449  i=0;
450  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
451  const Interaction * interaction = *ilistiter;
452  if(interaction->ExclTag().IsCharmEvent()) {
453  gr[i]->Draw("LP");
454  TString spltitle(interaction->AsString());
455  spltitle = spltitle.ReplaceAll(";",1," ",1);
456  legend->AddEntry(gr[i], spltitle.Data(),"LP");
457  }
458  i++;
459  }
460  legend->SetHeader("Charm Prod. Cross Sections");
461  //gr[nspl]->Draw("LP");
462  legend->AddEntry(gr[nspl], "sum","LP");
463  h->GetXaxis()->SetTitle("Ev (GeV)");
464  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
465  c->SetLogx();
466  c->SetLogy();
467  c->SetGridx();
468  c->SetGridy();
469  c->Update();
470  c->Clear();
471  c->Range(0,0,1,1);
472  legend->Draw();
473  c->Update();
474 
475  //-- plot ve xsecs only
476  //
477  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
478  legend->Clear();
479  i=0;
480  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
481  const Interaction * interaction = *ilistiter;
482  if(interaction->ProcInfo().IsInverseMuDecay() ||
483  interaction->ProcInfo().IsIMDAnnihilation() ||
484  interaction->ProcInfo().IsNuElectronElastic()) {
485  gr[i]->Draw("LP");
486  TString spltitle(interaction->AsString());
487  spltitle = spltitle.ReplaceAll(";",1," ",1);
488  legend->AddEntry(gr[i], spltitle.Data(),"LP");
489  }
490  i++;
491  }
492  legend->SetHeader("IMD and ve Elastic Cross Sections");
493  gr[nspl]->Draw("LP");
494  legend->AddEntry(gr[nspl], "sum","LP");
495  h->GetXaxis()->SetTitle("Ev (GeV)");
496  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
497  c->SetLogx();
498  c->SetLogy();
499  c->SetGridx();
500  c->SetGridy();
501  c->Update();
502  c->Clear();
503  c->Range(0,0,1,1);
504  legend->Draw();
505  c->Update();
506 
507  //-- close the postscript document
508  ps->Close();
509 
510  //-- clean-up
511  for(unsigned int j=0; j<=nspl; j++) { if(gr[j]) delete gr[j]; }
512  delete c;
513  delete ps;
514 }
genie::ProcessInfo::IsResonant
bool IsResonant(void) const
Definition: ProcessInfo.cxx:99
kNP
int kNP
Definition: gSplineXml2Root.cxx:153
genie::Spline::GetAsTGraph
TGraph * GetAsTGraph(int np=500, bool xscaling=false, bool inlog=false, double fx=1., double fy=1.) const
Definition: Spline.cxx:496
genie::ProcessInfo::IsInverseMuDecay
bool IsInverseMuDecay(void) const
Definition: ProcessInfo.cxx:126
genie::ProcessInfo::IsQuasiElastic
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
kPsType
const int kPsType
Definition: gSplineXml2Root.cxx:155
genie::Spline
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
genie::ProcessInfo::IsCoherentProduction
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:104
genie::ProcessInfo::IsIMDAnnihilation
bool IsIMDAnnihilation(void) const
Definition: ProcessInfo.cxx:131
genie::GEVGDriver::XSecSumSpline
const Spline * XSecSumSpline(void) const
Definition: GEVGDriver.h:87
genie::Interaction::AsString
string AsString(void) const
Definition: Interaction.cxx:246
gEmin
double gEmin
Definition: gSplineXml2Root.cxx:150
genie::XclsTag::IsCharmEvent
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::ProcessInfo::IsNuElectronElastic
bool IsNuElectronElastic(void) const
Definition: ProcessInfo.cxx:121
genie::GEVGDriver::Interactions
const InteractionList * Interactions(void) const
Definition: GEVGDriver.cxx:334
genie::units::cm2
static constexpr double cm2
Definition: Units.h:69
genie::GEVGDriver
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:54
pINFO
#define pINFO
Definition: Messenger.h:62
gOptProbePdgCode
int gOptProbePdgCode
Definition: gEvGenHadronNucleus.cxx:146
genie::units::ps
static constexpr double ps
Definition: Units.h:99
GetEventGenDriver
GEVGDriver GetEventGenDriver(void)
Definition: gSplineXml2Root.cxx:200
pWARN
#define pWARN
Definition: Messenger.h:60
gEmax
double gEmax
Definition: gSplineXml2Root.cxx:151
genie::ProcessInfo::IsDeepInelastic
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
genie::GEVGDriver::XSecSpline
const Spline * XSecSpline(const Interaction *interaction) const
Definition: GEVGDriver.cxx:488
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::Interaction::ExclTag
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
genie::InteractionList
A vector of Interaction objects.
Definition: InteractionList.h:34
genie::Interaction::ProcInfo
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
gWriteOutPlots
bool gWriteOutPlots
Definition: gSplineXml2Root.cxx:146
gOptTgtPdgCode
int gOptTgtPdgCode
Definition: gEvGenHadronNucleus.cxx:147

Variable Documentation

double gEmax

Definition at line 151 of file gSplineXml2Root.cxx.

Referenced by GetEventGenDriver(), SaveGraphsToRootFile(), and SaveToPsFile().

double gEmin

Definition at line 150 of file gSplineXml2Root.cxx.

Referenced by GetEventGenDriver(), SaveGraphsToRootFile(), and SaveToPsFile().

bool gInlogE

Definition at line 152 of file gSplineXml2Root.cxx.

Referenced by SaveGraphsToRootFile().

int gOptProbePdgCode

Definition at line 144 of file gSplineXml2Root.cxx.

PDGCodeList gOptProbePdgList

Definition at line 142 of file gSplineXml2Root.cxx.

Referenced by main().

string gOptROOTFilename

Definition at line 141 of file gSplineXml2Root.cxx.

Referenced by SaveGraphsToRootFile().

int gOptTgtPdgCode

Definition at line 145 of file gSplineXml2Root.cxx.

PDGCodeList gOptTgtPdgList

Definition at line 143 of file gSplineXml2Root.cxx.

Referenced by main().

string gOptXMLFilename

Definition at line 140 of file gSplineXml2Root.cxx.

bool gWriteOutPlots

Definition at line 146 of file gSplineXml2Root.cxx.

Referenced by SaveToPsFile().

int kNP = 300

Definition at line 153 of file gSplineXml2Root.cxx.

Referenced by SaveToPsFile().

int kNSplineP = 1000

Definition at line 154 of file gSplineXml2Root.cxx.

Referenced by SaveGraphsToRootFile().

const int kPsType = 111

Definition at line 155 of file gSplineXml2Root.cxx.

Referenced by SaveToPsFile().

Generated on Wed Jul 17 2024 14:46:48 for GENIEGenerator by   doxygen 1.8.5