gSplineXml2Root.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gspl2root
5 
6 \brief Utility to load a GENIE XML cross section spline file and convert
7  splines into a ROOT format.
8 
9  Syntax :
10  gspl2root -f xml_file -p nu -t tgt [-e emax]
11  [-o root_file] [-w] [-k]
12  [--message-thresholds xml_file]
13  [--event-generator-list list_name]
14 
15  Options :
16  [] denotes an optional argument
17 
18  -f
19  the input XML file containing the cross section spline data
20  -p
21  the neutrino pdg code
22  -t
23  the target pdg code (format: 10LZZZAAAI)
24  -e
25  the maximum energy (in generated plots -- use it to zoom at low E)
26  -o
27  output ROOT file name
28  -w
29  write out plots in a postscipt file
30  -k
31  keep spline knot points (not yet implemented).
32  --message-thresholds
33  Allows users to customize the message stream thresholds.
34  --event-generator-list
35  List of event generators to load in event generation drivers.
36  [default: "Default"].
37 
38  Example:
39 
40  Extract all numu+n, numu+p, numu+O16 cross section splines from the
41  input XML file `mysplines.xml', convert splines into a ROOT format
42  and save them into a single ROOT file.
43 
44  shell$ gspl2root -f mysplines.xml -p 14 -t 1000000010 -o xsec.root
45  shell$ gspl2root -f mysplines.xml -p 14 -t 1000010010 -o xsec.root
46  shell$ gspl2root -f mysplines.xml -p 14 -t 1000080160 -o xsec.root
47 
48  A large number of graphs (one per simulated process + appropriate
49  totals) will be generated in each case. Each set of plots is saved
50  into its own ROOT TDirectory named after the specified initial state.
51 
52  Important Note:
53 
54  The stored graphs can be used for cross section interpolation.
55  Essentially, this _single_ ROOT file can contain _all_ the GENIE cross
56  section `functions' you may need.
57  For instance, the `xsec.root' file generated in the above example will
58  contain a `nu_mu_O16' directory (generated by the last command)
59  which will include cross section graphs for all numu + O16 processes.
60  To extract the numu+O16 DIS CC cross section graph for hit u valence
61  quarks in a bound proton and evaluate the cross section at energy E,
62  then type:
63 
64  root[0] TDirectory * dir = (TDirectory*) file->Get("nu_mu_O16");
65  root[1] TGraph * graph = (TGraph*) dir->Get("dis_cc_p_uval");
66  root[2] cout << graph->Evaluate(E) << endl;
67 
68  See the User Manual for more details.
69 
70 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
71  University of Liverpool & STFC Rutherford Appleton Lab
72 
73 \created December 15, 2005
74 
75 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
76  For the full text of the license visit http://copyright.genie-mc.org
77  or see $GENIE/LICENSE
78 */
79 //____________________________________________________________________________
80 
81 #include <cassert>
82 #include <string>
83 #include <sstream>
84 #include <vector>
85 
86 #include <TSystem.h>
87 #include <TFile.h>
88 #include <TTree.h>
89 #include <TDirectory.h>
90 #include <TPostScript.h>
91 #include <TCanvas.h>
92 #include <TLegend.h>
93 #include <TGraph.h>
94 #include <TPaveText.h>
95 #include <TString.h>
96 #include <TH1F.h>
97 
111 #include "Framework/Utils/AppInit.h"
112 #include "Framework/Utils/RunOpt.h"
116 
117 using std::string;
118 using std::vector;
119 using std::ostringstream;
120 
121 using namespace genie;
122 using namespace genie::utils;
123 
124 //Prototypes:
125 void LoadSplines (void);
127 void SaveToPsFile (void);
128 void SaveGraphsToRootFile (void);
129 void SaveNtupleToRootFile (void);
130 void GetCommandLineArgs (int argc, char ** argv);
131 void PrintSyntax (void);
133 
134 //User-specified options:
135 string gOptXMLFilename; // input XML filename
136 string gOptROOTFilename; // output ROOT filename
137 double gOptNuEnergy; // Ev(max)
138 PDGCodeList gOptProbePdgList; // list of probe PDG codes
139 PDGCodeList gOptTgtPdgList; // list of target PDG codes
140 int gOptProbePdgCode; // probe PDG code (currently being processed)
141 int gOptTgtPdgCode; // target PDG code
142 bool gWriteOutPlots; // write out a postscript file with plots
143 //bool gKeepSplineKnots; // use spline abscissa points rather than equi-spaced
144 
145 //Globals & constants
146 double gEmin;
147 double gEmax;
148 int kNP = 300;
149 int kNSplineP = 1000;
150 const int kPsType = 111; // ps type: portrait
151 const double kEmin = 0.01; // minimum energy in plots (GeV)
152 
153 //____________________________________________________________________________
154 int main(int argc, char ** argv)
155 {
156  GetCommandLineArgs(argc,argv);
157 
158  if ( ! RunOpt::Instance()->Tune() ) {
159  LOG("gslp2root", pFATAL) << " No TuneId in RunOption";
160  exit(-1);
161  }
163 
164  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
165 
166  // load the x-section splines xml file specified by the user
167  LoadSplines();
168 
169  if ( ! XSecSplineList::Instance() -> HasSplineFromTune(RunOpt::Instance() -> Tune() -> Name() ) ) {
170  LOG("gspl2root", pWARN) << "No splines loaded for tune " << RunOpt::Instance() -> Tune() -> Name() ;
171  }
172 
173  for (unsigned int indx_p = 0; indx_p < gOptProbePdgList.size(); ++indx_p ) {
174  for (unsigned int indx_t = 0; indx_t < gOptTgtPdgList.size(); ++indx_t ) {
175  gOptProbePdgCode = gOptProbePdgList[indx_p];
176  gOptTgtPdgCode = gOptTgtPdgList[indx_t];
177  // save the cross section plots in a postscript file
178  SaveToPsFile();
179  // save the cross section graphs at a root file
181  }
182  }
183 
184  return 0;
185 }
186 //____________________________________________________________________________
187 void LoadSplines(void)
188 {
189 // load the cross section splines specified at the cmd line
190 
193  assert(ist == kXmlOK);
194 }
195 //____________________________________________________________________________
197 {
198 // create an event genartion driver configured for the specified initial state
199 // (so that cross section splines will be accessed through that driver as in
200 // event generation mode)
201 
203 
204  GEVGDriver evg_driver;
206  evg_driver.Configure(init_state);
207  evg_driver.CreateSplines();
208  evg_driver.CreateXSecSumSpline (100, gEmin, gEmax);
209 
210  return evg_driver;
211 }
212 //____________________________________________________________________________
213 void SaveToPsFile(void)
214 {
215  if(!gWriteOutPlots) return;
216 
217  //-- get the event generation driver
218  GEVGDriver evg_driver = GetEventGenDriver();
219 
220  //-- define some marker styles / colors
221  const unsigned int kNMarkers = 5;
222  const unsigned int kNColors = 6;
223  unsigned int markers[kNMarkers] = {20, 28, 29, 27, 3};
224  unsigned int colors [kNColors] = {1, 2, 4, 6, 8, 28};
225 
226  //-- create a postscript document for saving cross section plpots
227 
228  TCanvas * c = new TCanvas("c","",20,20,500,850);
229  c->SetBorderMode(0);
230  c->SetFillColor(0);
231  TLegend * legend = new TLegend(0.01,0.01,0.99,0.99);
232  legend->SetFillColor(0);
233  legend->SetBorderSize(0);
234 
235  //-- get pdglibrary for mapping pdg codes to names
236  PDGLibrary * pdglib = PDGLibrary::Instance();
237  ostringstream filename;
238  filename << "xsec-splines-"
239  << pdglib->Find(gOptProbePdgCode)->GetName() << "-"
240  << pdglib->Find(gOptTgtPdgCode)->GetName() << ".ps";
241  TPostScript * ps = new TPostScript(filename.str().c_str(), kPsType);
242 
243  //-- get the list of interactions that can be simulated by the driver
244  const InteractionList * ilist = evg_driver.Interactions();
245  unsigned int nspl = ilist->size();
246 
247  //-- book enough space for xsec plots (last one is the sum)
248  TGraph * gr[nspl+1];
249 
250  //-- loop over all the simulated interactions & create the cross section graphs
251  InteractionList::const_iterator ilistiter = ilist->begin();
252  unsigned int i=0;
253  for(; ilistiter != ilist->end(); ++ilistiter) {
254 
255  const Interaction * interaction = *ilistiter;
256  LOG("gspl2root", pINFO)
257  << "Current interaction: " << interaction->AsString();
258 
259 
260  //-- access the cross section spline
261  const Spline * spl = evg_driver.XSecSpline(interaction);
262  if(!spl) {
263  LOG("gspl2root", pWARN)
264  << "Can't get spline for: " << interaction->AsString();
265  exit(2);
266  }
267 
268  //-- set graph color/style
269  int icol = TMath::Min( i % kNColors, kNColors-1 );
270  int isty = TMath::Min( i / kNMarkers, kNMarkers-1 );
271  int col = colors[icol];
272  int sty = markers[isty];
273 
274  LOG("gspl2root", pINFO)
275  << "color = " << col << ", marker = " << sty;
276 
277  //-- export Spline as TGraph / set color & stype
278  gr[i] = spl->GetAsTGraph(kNP,true,true,1.,1./units::cm2);
279  gr[i]->SetLineColor(col);
280  gr[i]->SetMarkerColor(col);
281  gr[i]->SetMarkerStyle(sty);
282  gr[i]->SetMarkerSize(0.5);
283 
284  i++;
285  }
286 
287  //-- now, get the sum...
288  const Spline * splsum = evg_driver.XSecSumSpline();
289  if(!splsum) {
290  LOG("gspl2root", pWARN)
291  << "Can't get the cross section sum spline";
292  exit(2);
293  }
294  gr[nspl] = splsum->GetAsTGraph(kNP,true,true,1.,1./units::cm2);
295 
296  //-- figure out the minimum / maximum xsec in plotted range
297  double XSmax = -9999;
298  double XSmin = 9999;
299  double x=0, y=0;
300  for(int j=0; j<kNP; j++) {
301  gr[nspl]->GetPoint(j,x,y);
302  XSmax = TMath::Max(XSmax,y);
303  }
304  XSmin = XSmax/100000.;
305  XSmax = XSmax*1.2;
306 
307  LOG("gspl2root", pINFO) << "Drawing frame: E = (" << gEmin << ", " << gEmax << ")";
308  LOG("gspl2root", pINFO) << "Drawing frame: XSec = (" << XSmin << ", " << XSmax << ")";
309 
310  //-- ps output: add the 1st page with _all_ xsec spline plots
311  //
312  //c->Draw();
313  TH1F * h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
314  for(unsigned int ispl = 0; ispl <= nspl; ispl++) if(gr[ispl]) { gr[ispl]->Draw("LP"); }
315  h->GetXaxis()->SetTitle("Ev (GeV)");
316  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
317  c->SetLogx();
318  c->SetLogy();
319  c->SetGridx();
320  c->SetGridy();
321  c->Update();
322 
323  //-- plot QEL xsecs only
324  //
325  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
326  i=0;
327  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
328  const Interaction * interaction = *ilistiter;
329  if(interaction->ProcInfo().IsQuasiElastic() && ! interaction->ExclTag().IsCharmEvent()) {
330  gr[i]->Draw("LP");
331  TString spltitle(interaction->AsString());
332  spltitle = spltitle.ReplaceAll(";",1," ",1);
333  legend->AddEntry(gr[i], spltitle.Data(),"LP");
334  }
335  i++;
336  }
337  legend->SetHeader("QEL Cross Sections");
338  gr[nspl]->Draw("LP");
339  legend->AddEntry(gr[nspl], "sum","LP");
340  h->GetXaxis()->SetTitle("Ev (GeV)");
341  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
342  c->SetLogx();
343  c->SetLogy();
344  c->SetGridx();
345  c->SetGridy();
346  c->Update();
347  c->Clear();
348  c->Range(0,0,1,1);
349  legend->Draw();
350  c->Update();
351 
352  //-- plot RES xsecs only
353  //
354  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
355  legend->Clear();
356  i=0;
357  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
358  const Interaction * interaction = *ilistiter;
359  if(interaction->ProcInfo().IsResonant()) {
360  gr[i]->Draw("LP");
361  TString spltitle(interaction->AsString());
362  spltitle = spltitle.ReplaceAll(";",1," ",1);
363  legend->AddEntry(gr[i], spltitle.Data(),"LP");
364  }
365  i++;
366  }
367  legend->SetHeader("RES Cross Sections");
368  gr[nspl]->Draw("LP");
369  legend->AddEntry(gr[nspl], "sum","LP");
370  h->GetXaxis()->SetTitle("Ev (GeV)");
371  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
372  c->SetLogx();
373  c->SetLogy();
374  c->SetGridx();
375  c->SetGridy();
376  c->Update();
377  c->Clear();
378  c->Range(0,0,1,1);
379  legend->Draw();
380  c->Update();
381 
382  //-- plot DIS xsecs only
383  //
384  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
385  legend->Clear();
386  i=0;
387  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
388  const Interaction * interaction = *ilistiter;
389  if(interaction->ProcInfo().IsDeepInelastic()) {
390  gr[i]->Draw("LP");
391  TString spltitle(interaction->AsString());
392  spltitle = spltitle.ReplaceAll(";",1," ",1);
393  legend->AddEntry(gr[i], spltitle.Data(),"LP");
394  }
395  i++;
396  }
397  legend->SetHeader("DIS Cross Sections");
398  gr[nspl]->Draw("LP");
399  legend->AddEntry(gr[nspl], "sum","LP");
400  h->GetXaxis()->SetTitle("Ev (GeV)");
401  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
402  c->SetLogx();
403  c->SetLogy();
404  c->SetGridx();
405  c->SetGridy();
406  c->Update();
407  c->Clear();
408  c->Range(0,0,1,1);
409  legend->Draw();
410  c->Update();
411 
412  //-- plot COH xsecs only
413  //
414  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
415  legend->Clear();
416  i=0;
417  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
418  const Interaction * interaction = *ilistiter;
419  if(interaction->ProcInfo().IsCoherent()) {
420  gr[i]->Draw("LP");
421  TString spltitle(interaction->AsString());
422  spltitle = spltitle.ReplaceAll(";",1," ",1);
423  legend->AddEntry(gr[i], spltitle.Data(),"LP");
424  }
425  i++;
426  }
427  legend->SetHeader("COH Cross Sections");
428  gr[nspl]->Draw("LP");
429  legend->AddEntry(gr[nspl], "sum","LP");
430  h->GetXaxis()->SetTitle("Ev (GeV)");
431  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
432  c->SetLogx();
433  c->SetLogy();
434  c->SetGridx();
435  c->SetGridy();
436  c->Update();
437  c->Clear();
438  c->Range(0,0,1,1);
439  legend->Draw();
440  c->Update();
441 
442  //-- plot charm xsecs only
443  //
444  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
445  i=0;
446  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
447  const Interaction * interaction = *ilistiter;
448  if(interaction->ExclTag().IsCharmEvent()) {
449  gr[i]->Draw("LP");
450  TString spltitle(interaction->AsString());
451  spltitle = spltitle.ReplaceAll(";",1," ",1);
452  legend->AddEntry(gr[i], spltitle.Data(),"LP");
453  }
454  i++;
455  }
456  legend->SetHeader("Charm Prod. Cross Sections");
457  //gr[nspl]->Draw("LP");
458  legend->AddEntry(gr[nspl], "sum","LP");
459  h->GetXaxis()->SetTitle("Ev (GeV)");
460  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
461  c->SetLogx();
462  c->SetLogy();
463  c->SetGridx();
464  c->SetGridy();
465  c->Update();
466  c->Clear();
467  c->Range(0,0,1,1);
468  legend->Draw();
469  c->Update();
470 
471  //-- plot ve xsecs only
472  //
473  h = (TH1F*) c->DrawFrame(gEmin, XSmin, gEmax, XSmax);
474  legend->Clear();
475  i=0;
476  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
477  const Interaction * interaction = *ilistiter;
478  if(interaction->ProcInfo().IsInverseMuDecay() ||
479  interaction->ProcInfo().IsIMDAnnihilation() ||
480  interaction->ProcInfo().IsNuElectronElastic()) {
481  gr[i]->Draw("LP");
482  TString spltitle(interaction->AsString());
483  spltitle = spltitle.ReplaceAll(";",1," ",1);
484  legend->AddEntry(gr[i], spltitle.Data(),"LP");
485  }
486  i++;
487  }
488  legend->SetHeader("IMD and ve Elastic Cross Sections");
489  gr[nspl]->Draw("LP");
490  legend->AddEntry(gr[nspl], "sum","LP");
491  h->GetXaxis()->SetTitle("Ev (GeV)");
492  h->GetYaxis()->SetTitle("#sigma_{nuclear}/Ev (cm^{2}/GeV)");
493  c->SetLogx();
494  c->SetLogy();
495  c->SetGridx();
496  c->SetGridy();
497  c->Update();
498  c->Clear();
499  c->Range(0,0,1,1);
500  legend->Draw();
501  c->Update();
502 
503  //-- close the postscript document
504  ps->Close();
505 
506  //-- clean-up
507  for(unsigned int j=0; j<=nspl; j++) { if(gr[j]) delete gr[j]; }
508  delete c;
509  delete ps;
510 }
511 //____________________________________________________________________________
512 void FormatXSecGraph(TGraph * g)
513 {
514  g->SetTitle("GENIE cross section graph");
515  g->GetXaxis()->SetTitle("Ev (GeV)");
516  g->GetYaxis()->SetTitle("#sigma_{nuclear} (10^{-38} cm^{2})");
517 }
518 //____________________________________________________________________________
520 {
521  //-- get the event generation driver
522  GEVGDriver evg_driver = GetEventGenDriver();
523 
524  //-- get the list of interactions that can be simulated by the driver
525  const InteractionList * ilist = evg_driver.Interactions();
526 
527  //-- check whether the splines will be saved in a ROOT file - if not, exit now
528  bool save_in_root = gOptROOTFilename.size()>0;
529  if(!save_in_root) return;
530 
531  //-- get pdglibrary for mapping pdg codes to names
532  PDGLibrary * pdglib = PDGLibrary::Instance();
533 
534  //-- check whether the requested filename exists
535  // if yes, then open the file in 'update' mode
536  bool exists = !(gSystem->AccessPathName(gOptROOTFilename.c_str()));
537 
538  TFile * froot = 0;
539  if(exists) froot = new TFile(gOptROOTFilename.c_str(), "UPDATE");
540  else froot = new TFile(gOptROOTFilename.c_str(), "RECREATE");
541  assert(froot);
542 
543  //-- create directory
544  ostringstream dptr;
545 
546  string probe_name = pdglib->Find(gOptProbePdgCode)->GetName();
547  string tgt_name = (gOptTgtPdgCode==1000000010) ?
548  "n" : pdglib->Find(gOptTgtPdgCode)->GetName();
549 
550  dptr << probe_name << "_" << tgt_name;
551  ostringstream dtitle;
552  dtitle << "Cross sections for: "
553  << pdglib->Find(gOptProbePdgCode)->GetName() << "+"
554  << pdglib->Find(gOptTgtPdgCode)->GetName();
555 
556  LOG("gspl2root", pINFO)
557  << "Will store graphs in root directory = " << dptr.str();
558  TDirectory * topdir =
559  dynamic_cast<TDirectory *> (froot->Get(dptr.str().c_str()));
560  if(topdir) {
561  LOG("gspl2root", pINFO)
562  << "Directory: " << dptr.str() << " already exists!! Exiting";
563  froot->Close();
564  delete froot;
565  return;
566  }
567 
568  topdir = froot->mkdir(dptr.str().c_str(),dtitle.str().c_str());
569  topdir->cd();
570 
571  double de = (gEmax-gEmin)/(kNSplineP-1);
572  double * e = new double[kNSplineP];
573  for(int i=0; i<kNSplineP; i++) { e[i] = gEmin + i*de; }
574 
575  double * xs = new double[kNSplineP];
576 
577  InteractionList::const_iterator ilistiter = ilist->begin();
578 
579  for(; ilistiter != ilist->end(); ++ilistiter) {
580 
581  const Interaction * interaction = *ilistiter;
582 
583  const ProcessInfo & proc = interaction->ProcInfo();
584  const XclsTag & xcls = interaction->ExclTag();
585  const InitialState & init = interaction->InitState();
586  const Target & tgt = init.Tgt();
587 
588  // graph title
589  ostringstream title;
590 
591  if (proc.IsQuasiElastic() ) { title << "qel"; }
592  else if (proc.IsMEC() ) { title << "mec"; }
593  else if (proc.IsResonant() ) { title << "res"; }
594  else if (proc.IsDeepInelastic() ) { title << "dis"; }
595  else if (proc.IsDiffractive() ) { title << "dfr"; }
596  else if (proc.IsCoherent() ) { title << "coh"; }
597  else if (proc.IsCoherentElas() ) { title << "cohel"; }
598  else if (proc.IsInverseMuDecay() ) { title << "imd"; }
599  else if (proc.IsIMDAnnihilation() ) { title << "imdanh";}
600  else if (proc.IsNuElectronElastic()) { title << "ve"; }
601  else { continue; }
602 
603  if (proc.IsWeakCC()) { title << "_cc"; }
604  else if (proc.IsWeakNC()) { title << "_nc"; }
605  else if (proc.IsWeakMix()) { title << "_ccncmix"; }
606  else if (proc.IsEM() ) { title << "_em"; }
607  else { continue; }
608 
609  if(tgt.HitNucIsSet()) {
610  int hitnuc = tgt.HitNucPdg();
611  if ( pdg::IsProton (hitnuc) ) { title << "_p"; }
612  else if ( pdg::IsNeutron(hitnuc) ) { title << "_n"; }
613  else if ( pdg::Is2NucleonCluster(hitnuc) )
614  {
615  if (hitnuc == kPdgClusterNN) { title << "_nn"; }
616  else if (hitnuc == kPdgClusterNP) { title << "_np"; }
617  else if (hitnuc == kPdgClusterPP) { title << "_pp"; }
618  else {
619  LOG("gspl2root", pWARN) << "Can't handle hit 2-nucleon cluster PDG = " << hitnuc;
620  }
621  }
622  else {
623  LOG("gspl2root", pWARN) << "Can't handle hit nucleon PDG = " << hitnuc;
624  }
625 
626  if(tgt.HitQrkIsSet()) {
627  int qrkpdg = tgt.HitQrkPdg();
628  bool insea = tgt.HitSeaQrk();
629 
630  if ( pdg::IsUQuark(qrkpdg) ) { title << "_u"; }
631  else if ( pdg::IsDQuark(qrkpdg) ) { title << "_d"; }
632  else if ( pdg::IsSQuark(qrkpdg) ) { title << "_s"; }
633  else if ( pdg::IsCQuark(qrkpdg) ) { title << "_c"; }
634  else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
635  else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
636  else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
637  else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
638 
639  if(insea) { title << "sea"; }
640  else { title << "val"; }
641  }
642  }
643  if(proc.IsResonant()) {
644  Resonance_t res = xcls.Resonance();
645  string resname = res::AsString(res);
646  resname = str::FilterString(")", resname);
647  resname = str::FilterString("(", resname);
648  title << "_" << resname.substr(3,4) << resname.substr(0,3);
649  }
650 
651  if(xcls.IsStrangeEvent()) {
652  title << "_strange";
653  if(!xcls.IsInclusiveStrange()) { title << xcls.StrangeHadronPdg(); }
654  }
655 
656  if(xcls.IsCharmEvent()) {
657  title << "_charm";
658  if(!xcls.IsInclusiveCharm()) { title << xcls.CharmHadronPdg(); }
659  }
660 
661  const Spline * spl = evg_driver.XSecSpline(interaction);
662  for(int i=0; i<kNSplineP; i++) {
663  xs[i] = spl->Evaluate(e[i]) * (1E+38/units::cm2);
664  }
665 
666  TGraph * gr = new TGraph(kNSplineP, e, xs);
667  gr->SetName(title.str().c_str());
668  FormatXSecGraph(gr);
669  gr->SetTitle(spl->GetName());
670 
671  topdir->Add(gr);
672  }
673 
674 
675  //
676  // totals for (anti-)neutrino scattering
677  //
678 
679  bool is_neutrino = pdg::IsNeutralLepton(gOptProbePdgCode);
680 
681  if(is_neutrino) {
682 
683  //
684  // add-up all res channels
685  //
686 
687  double * xsresccp = new double[kNSplineP];
688  double * xsresccn = new double[kNSplineP];
689  double * xsresncp = new double[kNSplineP];
690  double * xsresncn = new double[kNSplineP];
691  for(int i=0; i<kNSplineP; i++) {
692  xsresccp[i] = 0;
693  xsresccn[i] = 0;
694  xsresncp[i] = 0;
695  xsresncn[i] = 0;
696  }
697 
698  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
699  const Interaction * interaction = *ilistiter;
700  const ProcessInfo & proc = interaction->ProcInfo();
701  const InitialState & init = interaction->InitState();
702  const Target & tgt = init.Tgt();
703 
704  const Spline * spl = evg_driver.XSecSpline(interaction);
705 
706  if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
707  for(int i=0; i<kNSplineP; i++) {
708  xsresccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
709  }
710  }
711  if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
712  for(int i=0; i<kNSplineP; i++) {
713  xsresccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
714  }
715  }
716  if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
717  for(int i=0; i<kNSplineP; i++) {
718  xsresncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
719  }
720  }
721  if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
722  for(int i=0; i<kNSplineP; i++) {
723  xsresncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
724  }
725  }
726  }
727 
728  TGraph * gr_resccp = new TGraph(kNSplineP, e, xsresccp);
729  gr_resccp->SetName("res_cc_p");
730  FormatXSecGraph(gr_resccp);
731  topdir->Add(gr_resccp);
732  TGraph * gr_resccn = new TGraph(kNSplineP, e, xsresccn);
733  gr_resccn->SetName("res_cc_n");
734  FormatXSecGraph(gr_resccn);
735  topdir->Add(gr_resccn);
736  TGraph * gr_resncp = new TGraph(kNSplineP, e, xsresncp);
737  gr_resncp->SetName("res_nc_p");
738  FormatXSecGraph(gr_resncp);
739  topdir->Add(gr_resncp);
740  TGraph * gr_resncn = new TGraph(kNSplineP, e, xsresncn);
741  gr_resncn->SetName("res_nc_n");
742  FormatXSecGraph(gr_resncn);
743  topdir->Add(gr_resncn);
744 
745  //
746  // add-up all dis channels
747  //
748 
749  double * xsdisccp = new double[kNSplineP];
750  double * xsdisccn = new double[kNSplineP];
751  double * xsdisncp = new double[kNSplineP];
752  double * xsdisncn = new double[kNSplineP];
753  for(int i=0; i<kNSplineP; i++) {
754  xsdisccp[i] = 0;
755  xsdisccn[i] = 0;
756  xsdisncp[i] = 0;
757  xsdisncn[i] = 0;
758  }
759  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
760  const Interaction * interaction = *ilistiter;
761  const ProcessInfo & proc = interaction->ProcInfo();
762  const XclsTag & xcls = interaction->ExclTag();
763  const InitialState & init = interaction->InitState();
764  const Target & tgt = init.Tgt();
765 
766  const Spline * spl = evg_driver.XSecSpline(interaction);
767 
768  if(xcls.IsCharmEvent()) continue;
769 
770  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
771  for(int i=0; i<kNSplineP; i++) {
772  xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
773  }
774  }
775  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
776  for(int i=0; i<kNSplineP; i++) {
777  xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
778  }
779  }
780  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
781  for(int i=0; i<kNSplineP; i++) {
782  xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
783  }
784  }
785  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
786  for(int i=0; i<kNSplineP; i++) {
787  xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
788  }
789  }
790  }
791  TGraph * gr_disccp = new TGraph(kNSplineP, e, xsdisccp);
792  gr_disccp->SetName("dis_cc_p");
793  FormatXSecGraph(gr_disccp);
794  topdir->Add(gr_disccp);
795  TGraph * gr_disccn = new TGraph(kNSplineP, e, xsdisccn);
796  gr_disccn->SetName("dis_cc_n");
797  FormatXSecGraph(gr_disccn);
798  topdir->Add(gr_disccn);
799  TGraph * gr_disncp = new TGraph(kNSplineP, e, xsdisncp);
800  gr_disncp->SetName("dis_nc_p");
801  FormatXSecGraph(gr_disncp);
802  topdir->Add(gr_disncp);
803  TGraph * gr_disncn = new TGraph(kNSplineP, e, xsdisncn);
804  gr_disncn->SetName("dis_nc_n");
805  FormatXSecGraph(gr_disncn);
806  topdir->Add(gr_disncn);
807 
808  //
809  // add-up all charm dis channels
810  //
811 
812  for(int i=0; i<kNSplineP; i++) {
813  xsdisccp[i] = 0;
814  xsdisccn[i] = 0;
815  xsdisncp[i] = 0;
816  xsdisncn[i] = 0;
817  }
818  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
819  const Interaction * interaction = *ilistiter;
820  const ProcessInfo & proc = interaction->ProcInfo();
821  const XclsTag & xcls = interaction->ExclTag();
822  const InitialState & init = interaction->InitState();
823  const Target & tgt = init.Tgt();
824 
825  const Spline * spl = evg_driver.XSecSpline(interaction);
826 
827  if(!xcls.IsCharmEvent()) continue;
828 
829  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
830  for(int i=0; i<kNSplineP; i++) {
831  xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
832  }
833  }
834  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
835  for(int i=0; i<kNSplineP; i++) {
836  xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
837  }
838  }
839  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
840  for(int i=0; i<kNSplineP; i++) {
841  xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
842  }
843  }
844  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
845  for(int i=0; i<kNSplineP; i++) {
846  xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
847  }
848  }
849  }
850  TGraph * gr_disccp_charm = new TGraph(kNSplineP, e, xsdisccp);
851  gr_disccp_charm->SetName("dis_cc_p_charm");
852  FormatXSecGraph(gr_disccp_charm);
853  topdir->Add(gr_disccp_charm);
854  TGraph * gr_disccn_charm = new TGraph(kNSplineP, e, xsdisccn);
855  gr_disccn_charm->SetName("dis_cc_n_charm");
856  FormatXSecGraph(gr_disccn_charm);
857  topdir->Add(gr_disccn_charm);
858  TGraph * gr_disncp_charm = new TGraph(kNSplineP, e, xsdisncp);
859  gr_disncp_charm->SetName("dis_nc_p_charm");
860  FormatXSecGraph(gr_disncp_charm);
861  topdir->Add(gr_disncp_charm);
862  TGraph * gr_disncn_charm = new TGraph(kNSplineP, e, xsdisncn);
863  gr_disncn_charm->SetName("dis_nc_n_charm");
864  FormatXSecGraph(gr_disncn_charm);
865  topdir->Add(gr_disncn_charm);
866 
867  //
868  // add-up all mec channels
869  //
870 
871  double * xsmeccc = new double[kNSplineP];
872  double * xsmecnc = new double[kNSplineP];
873  for(int i=0; i<kNSplineP; i++) {
874  xsmeccc[i] = 0;
875  xsmecnc[i] = 0;
876  }
877 
878  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
879  const Interaction * interaction = *ilistiter;
880  const ProcessInfo & proc = interaction->ProcInfo();
881 
882  const Spline * spl = evg_driver.XSecSpline(interaction);
883 
884  if (proc.IsMEC() && proc.IsWeakCC()) {
885  for(int i=0; i<kNSplineP; i++) {
886  xsmeccc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
887  }
888  }
889  if (proc.IsMEC() && proc.IsWeakNC()) {
890  for(int i=0; i<kNSplineP; i++) {
891  xsmecnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
892  }
893  }
894  }
895 
896  TGraph * gr_meccc = new TGraph(kNSplineP, e, xsmeccc);
897  gr_meccc->SetName("mec_cc");
898  FormatXSecGraph(gr_meccc);
899  topdir->Add(gr_meccc);
900  TGraph * gr_mecnc = new TGraph(kNSplineP, e, xsmecnc);
901  gr_mecnc->SetName("mec_nc");
902  FormatXSecGraph(gr_mecnc);
903  topdir->Add(gr_mecnc);
904 
905  //
906  // total cross sections
907  //
908  double * xstotcc = new double[kNSplineP];
909  double * xstotccp = new double[kNSplineP];
910  double * xstotccn = new double[kNSplineP];
911  double * xstotnc = new double[kNSplineP];
912  double * xstotncp = new double[kNSplineP];
913  double * xstotncn = new double[kNSplineP];
914  for(int i=0; i<kNSplineP; i++) {
915  xstotcc [i] = 0;
916  xstotccp[i] = 0;
917  xstotccn[i] = 0;
918  xstotnc [i] = 0;
919  xstotncp[i] = 0;
920  xstotncn[i] = 0;
921  }
922  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
923  const Interaction * interaction = *ilistiter;
924  const ProcessInfo & proc = interaction->ProcInfo();
925  const InitialState & init = interaction->InitState();
926  const Target & tgt = init.Tgt();
927 
928  const Spline * spl = evg_driver.XSecSpline(interaction);
929 
930  bool iscc = proc.IsWeakCC();
931  bool isnc = proc.IsWeakNC();
932  bool offp = pdg::IsProton (tgt.HitNucPdg());
933  bool offn = pdg::IsNeutron(tgt.HitNucPdg());
934 
935  if (iscc && offp) {
936  for(int i=0; i<kNSplineP; i++) {
937  xstotccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
938  }
939  }
940  if (iscc && offn) {
941  for(int i=0; i<kNSplineP; i++) {
942  xstotccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
943  }
944  }
945  if (isnc && offp) {
946  for(int i=0; i<kNSplineP; i++) {
947  xstotncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
948  }
949  }
950  if (isnc && offn) {
951  for(int i=0; i<kNSplineP; i++) {
952  xstotncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
953  }
954  }
955 
956  if (iscc) {
957  for(int i=0; i<kNSplineP; i++) {
958  xstotcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
959  }
960  }
961  if (isnc) {
962  for(int i=0; i<kNSplineP; i++) {
963  xstotnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
964  }
965  }
966  }
967 
968  TGraph * gr_totcc = new TGraph(kNSplineP, e, xstotcc);
969  gr_totcc->SetName("tot_cc");
970  FormatXSecGraph(gr_totcc);
971  topdir->Add(gr_totcc);
972  TGraph * gr_totccp = new TGraph(kNSplineP, e, xstotccp);
973  gr_totccp->SetName("tot_cc_p");
974  FormatXSecGraph(gr_totccp);
975  topdir->Add(gr_totccp);
976  TGraph * gr_totccn = new TGraph(kNSplineP, e, xstotccn);
977  gr_totccn->SetName("tot_cc_n");
978  FormatXSecGraph(gr_totccn);
979  topdir->Add(gr_totccn);
980  TGraph * gr_totnc = new TGraph(kNSplineP, e, xstotnc);
981  gr_totnc->SetName("tot_nc");
982  FormatXSecGraph(gr_totnc);
983  topdir->Add(gr_totnc);
984  TGraph * gr_totncp = new TGraph(kNSplineP, e, xstotncp);
985  gr_totncp->SetName("tot_nc_p");
986  FormatXSecGraph(gr_totncp);
987  topdir->Add(gr_totncp);
988  TGraph * gr_totncn = new TGraph(kNSplineP, e, xstotncn);
989  gr_totncn->SetName("tot_nc_n");
990  FormatXSecGraph(gr_totncn);
991  topdir->Add(gr_totncn);
992 
993  delete [] e;
994  delete [] xs;
995  delete [] xsresccp;
996  delete [] xsresccn;
997  delete [] xsresncp;
998  delete [] xsresncn;
999  delete [] xsdisccp;
1000  delete [] xsdisccn;
1001  delete [] xsdisncp;
1002  delete [] xsdisncn;
1003  delete [] xstotcc;
1004  delete [] xstotccp;
1005  delete [] xstotccn;
1006  delete [] xstotnc;
1007  delete [] xstotncp;
1008  delete [] xstotncn;
1009 
1010  }// neutrinos
1011 
1012 
1013  //
1014  // totals for charged lepton scattering
1015  //
1016 
1017  bool is_charged_lepton = pdg::IsChargedLepton(gOptProbePdgCode);
1018 
1019  if(is_charged_lepton) {
1020 
1021  //
1022  // add-up all res channels
1023  //
1024 
1025  double * xsresemp = new double[kNSplineP];
1026  double * xsresemn = new double[kNSplineP];
1027  for(int i=0; i<kNSplineP; i++) {
1028  xsresemp[i] = 0;
1029  xsresemn[i] = 0;
1030  }
1031 
1032  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1033  const Interaction * interaction = *ilistiter;
1034  const ProcessInfo & proc = interaction->ProcInfo();
1035  const InitialState & init = interaction->InitState();
1036  const Target & tgt = init.Tgt();
1037 
1038  const Spline * spl = evg_driver.XSecSpline(interaction);
1039 
1040  if (proc.IsResonant() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1041  for(int i=0; i<kNSplineP; i++) {
1042  xsresemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1043  }
1044  }
1045  if (proc.IsResonant() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1046  for(int i=0; i<kNSplineP; i++) {
1047  xsresemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1048  }
1049  }
1050  }
1051 
1052  TGraph * gr_resemp = new TGraph(kNSplineP, e, xsresemp);
1053  gr_resemp->SetName("res_em_p");
1054  FormatXSecGraph(gr_resemp);
1055  topdir->Add(gr_resemp);
1056  TGraph * gr_resemn = new TGraph(kNSplineP, e, xsresemn);
1057  gr_resemn->SetName("res_em_n");
1058  FormatXSecGraph(gr_resemn);
1059  topdir->Add(gr_resemn);
1060 
1061  //
1062  // add-up all dis channels
1063  //
1064 
1065  double * xsdisemp = new double[kNSplineP];
1066  double * xsdisemn = new double[kNSplineP];
1067  for(int i=0; i<kNSplineP; i++) {
1068  xsdisemp[i] = 0;
1069  xsdisemn[i] = 0;
1070  }
1071  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1072  const Interaction * interaction = *ilistiter;
1073  const ProcessInfo & proc = interaction->ProcInfo();
1074  const XclsTag & xcls = interaction->ExclTag();
1075  const InitialState & init = interaction->InitState();
1076  const Target & tgt = init.Tgt();
1077 
1078  const Spline * spl = evg_driver.XSecSpline(interaction);
1079 
1080  if(xcls.IsCharmEvent()) continue;
1081 
1082  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1083  for(int i=0; i<kNSplineP; i++) {
1084  xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1085  }
1086  }
1087  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1088  for(int i=0; i<kNSplineP; i++) {
1089  xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1090  }
1091  }
1092  }
1093  TGraph * gr_disemp = new TGraph(kNSplineP, e, xsdisemp);
1094  gr_disemp->SetName("dis_em_p");
1095  FormatXSecGraph(gr_disemp);
1096  topdir->Add(gr_disemp);
1097  TGraph * gr_disemn = new TGraph(kNSplineP, e, xsdisemn);
1098  gr_disemn->SetName("dis_em_n");
1099  FormatXSecGraph(gr_disemn);
1100  topdir->Add(gr_disemn);
1101 
1102  //
1103  // add-up all charm dis channels
1104  //
1105 
1106  for(int i=0; i<kNSplineP; i++) {
1107  xsdisemp[i] = 0;
1108  xsdisemn[i] = 0;
1109  }
1110  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1111  const Interaction * interaction = *ilistiter;
1112  const ProcessInfo & proc = interaction->ProcInfo();
1113  const XclsTag & xcls = interaction->ExclTag();
1114  const InitialState & init = interaction->InitState();
1115  const Target & tgt = init.Tgt();
1116 
1117  const Spline * spl = evg_driver.XSecSpline(interaction);
1118 
1119  if(!xcls.IsCharmEvent()) continue;
1120 
1121  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1122  for(int i=0; i<kNSplineP; i++) {
1123  xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1124  }
1125  }
1126  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1127  for(int i=0; i<kNSplineP; i++) {
1128  xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1129  }
1130  }
1131  }
1132  TGraph * gr_disemp_charm = new TGraph(kNSplineP, e, xsdisemp);
1133  gr_disemp_charm->SetName("dis_em_p_charm");
1134  FormatXSecGraph(gr_disemp_charm);
1135  topdir->Add(gr_disemp_charm);
1136  TGraph * gr_disemn_charm = new TGraph(kNSplineP, e, xsdisemn);
1137  gr_disemn_charm->SetName("dis_em_n_charm");
1138  FormatXSecGraph(gr_disemn_charm);
1139  topdir->Add(gr_disemn_charm);
1140 
1141  //
1142  // total cross sections
1143  //
1144  double * xstotem = new double[kNSplineP];
1145  double * xstotemp = new double[kNSplineP];
1146  double * xstotemn = new double[kNSplineP];
1147  for(int i=0; i<kNSplineP; i++) {
1148  xstotem [i] = 0;
1149  xstotemp[i] = 0;
1150  xstotemn[i] = 0;
1151  }
1152  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1153  const Interaction * interaction = *ilistiter;
1154  const ProcessInfo & proc = interaction->ProcInfo();
1155  const InitialState & init = interaction->InitState();
1156  const Target & tgt = init.Tgt();
1157 
1158  const Spline * spl = evg_driver.XSecSpline(interaction);
1159 
1160  bool isem = proc.IsEM();
1161  bool offp = pdg::IsProton (tgt.HitNucPdg());
1162  bool offn = pdg::IsNeutron(tgt.HitNucPdg());
1163 
1164  if (isem && offp) {
1165  for(int i=0; i<kNSplineP; i++) {
1166  xstotemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1167  }
1168  }
1169  if (isem && offn) {
1170  for(int i=0; i<kNSplineP; i++) {
1171  xstotemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1172  }
1173  }
1174  if (isem) {
1175  for(int i=0; i<kNSplineP; i++) {
1176  xstotem[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1177  }
1178  }
1179  }
1180 
1181  TGraph * gr_totem = new TGraph(kNSplineP, e, xstotem);
1182  gr_totem->SetName("tot_em");
1183  FormatXSecGraph(gr_totem);
1184  topdir->Add(gr_totem);
1185  TGraph * gr_totemp = new TGraph(kNSplineP, e, xstotemp);
1186  gr_totemp->SetName("tot_em_p");
1187  FormatXSecGraph(gr_totemp);
1188  topdir->Add(gr_totemp);
1189  TGraph * gr_totemn = new TGraph(kNSplineP, e, xstotemn);
1190  gr_totemn->SetName("tot_em_n");
1191  FormatXSecGraph(gr_totemn);
1192  topdir->Add(gr_totemn);
1193 
1194  delete [] e;
1195  delete [] xs;
1196  delete [] xsresemp;
1197  delete [] xsresemn;
1198  delete [] xsdisemp;
1199  delete [] xsdisemn;
1200  delete [] xstotem;
1201  delete [] xstotemp;
1202  delete [] xstotemn;
1203 
1204  }// charged leptons
1205 
1206  topdir->Write();
1207 
1208  if(froot) {
1209  froot->Close();
1210  delete froot;
1211  }
1212 }
1213 //____________________________________________________________________________
1215 {
1216 /*
1217  //-- create cross section ntuple
1218  //
1219  double brXSec;
1220  double brEv;
1221  bool brIsQel;
1222  bool brIsRes;
1223  bool brIsDis;
1224  bool brIsCoh;
1225  bool brIsImd;
1226  bool brIsNuEEl;
1227  bool brIsCC;
1228  bool brIsNC;
1229  int brNucleon;
1230  int brQrk;
1231  bool brIsSeaQrk;
1232  int brRes;
1233  bool brCharm;
1234  TTree * xsecnt = new TTree("xsecnt",dtitle.str().c_str());
1235  xsecnt->Branch("xsec", &brXSec, "xsec/D");
1236  xsecnt->Branch("e", &brEv, "e/D");
1237  xsecnt->Branch("qel", &brIsQel, "qel/O");
1238  xsecnt->Branch("res", &brIsRes, "res/O");
1239  xsecnt->Branch("dis", &brIsDis, "dis/O");
1240  xsecnt->Branch("coh", &brIsCoh, "coh/O");
1241  xsecnt->Branch("imd", &brIsImd, "imd/O");
1242  xsecnt->Branch("veel", &brIsNuEEl, "veel/O");
1243  xsecnt->Branch("cc", &brIsCC, "cc/O");
1244  xsecnt->Branch("nc", &brIsNC, "nc/O");
1245  xsecnt->Branch("nuc", &brNucleon, "nuc/I");
1246  xsecnt->Branch("qrk", &brQrk, "qrk/I");
1247  xsecnt->Branch("sea", &brIsSeaQrk, "sea/O");
1248  xsecnt->Branch("res", &brRes, "res/I");
1249  xsecnt->Branch("charm", &brCharm, "charm/O");
1250 */
1251 }
1252 //____________________________________________________________________________
1253 void GetCommandLineArgs(int argc, char ** argv)
1254 {
1255  LOG("gspl2root", pINFO) << "Parsing command line arguments";
1256 
1257  // Common run options. Set defaults and read.
1258  RunOpt::Instance()->ReadFromCommandLine(argc,argv);
1259 
1260  // Parse run options for this app
1261 
1262  CmdLnArgParser parser(argc,argv);
1263 
1264  // input XML file name:
1265  if( parser.OptionExists('f') ) {
1266  LOG("gspl2root", pINFO) << "Reading input XML filename";
1267  gOptXMLFilename = parser.ArgAsString('f');
1268  } else {
1269  LOG("gspl2root", pFATAL) << "Unspecified input XML file!";
1270  PrintSyntax();
1271  exit(1);
1272  }
1273 
1274  // probe PDG code:
1275  if( parser.OptionExists('p') ) {
1276  LOG("gspl2root", pINFO) << "Reading probe PDG code";
1277  gOptProbePdgList = GetPDGCodeListFromString(parser.ArgAsString('p'));
1278  } else {
1279  LOG("gspl2root", pFATAL)
1280  << "Unspecified probe PDG code - Exiting";
1281  PrintSyntax();
1282  exit(1);
1283  }
1284 
1285  // target PDG code:
1286  if( parser.OptionExists('t') ) {
1287  LOG("gspl2root", pINFO) << "Reading target PDG code";
1288  gOptTgtPdgList = GetPDGCodeListFromString(parser.ArgAsString('t'));
1289  } else {
1290  LOG("gspl2root", pFATAL)
1291  << "Unspecified target PDG code - Exiting";
1292  PrintSyntax();
1293  exit(1);
1294  }
1295 
1296  // max neutrino energy
1297  if( parser.OptionExists('e') ) {
1298  LOG("gspl2root", pINFO) << "Reading neutrino energy";
1299  gOptNuEnergy = parser.ArgAsDouble('e');
1300  } else {
1301  LOG("gspl2root", pDEBUG)
1302  << "Unspecified Emax - Setting to 100 GeV";
1303  gOptNuEnergy = 100;
1304  }
1305 
1306  // output ROOT file name:
1307  if( parser.OptionExists('o') ) {
1308  LOG("gspl2root", pINFO) << "Reading output ROOT filename";
1309  gOptROOTFilename = parser.ArgAsString('o');
1310  } else {
1311  LOG("gspl2root", pDEBUG)
1312  << "Unspecified output ROOT file. Using default: gxsec.root";
1313  gOptROOTFilename = "gxsec.root";
1314  }
1315 
1316  // write-out a PS file with plots
1317  gWriteOutPlots = parser.OptionExists('w');
1318 
1319  // use same abscissa points as splines
1320  //not yet//gKeepSplineKnots = parser.OptionExists('k');
1321 
1322 
1323  gEmin = kEmin;
1324  gEmax = gOptNuEnergy;
1325  assert(gEmin<gEmax);
1326 
1327  // print the options you got from command line arguments
1328  LOG("gspl2root", pINFO) << "Command line arguments:";
1329  LOG("gspl2root", pINFO) << " Input XML file = " << gOptXMLFilename;
1330  LOG("gspl2root", pINFO) << " Probe PDG code = " << gOptProbePdgCode;
1331  LOG("gspl2root", pINFO) << " Target PDG code = " << gOptTgtPdgCode;
1332  LOG("gspl2root", pINFO) << " Max neutrino E = " << gOptNuEnergy;
1333  //not yet//LOG("gspl2root", pINFO) << " Keep spline knots = " << (gKeepSplineKnots?"true":"false");
1334 }
1335 //____________________________________________________________________________
1336 void PrintSyntax(void)
1337 {
1338  LOG("gspl2root", pNOTICE)
1339  << "\n\n" << "Syntax:" << "\n"
1340  << " gspl2root -f xml_file -p probe_pdg -t target_pdg"
1341  << " [-e emax] [-o output_root_file] [-w]\n"
1342  << " [--message-thresholds xml_file]\n";
1343 }
1344 //____________________________________________________________________________
1346 {
1347  vector<string> isvec = utils::str::Split(s,",");
1348 
1349  // fill in the PDG code list
1350  PDGCodeList list;
1351  vector<string>::const_iterator iter;
1352  for(iter = isvec.begin(); iter != isvec.end(); ++iter) {
1353  list.push_back( atoi(iter->c_str()) );
1354  }
1355  return list;
1356 
1357 }
1358 //____________________________________________________________________________
bool IsResonant(void) const
Definition: ProcessInfo.cxx:92
int kNP
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
Definition: GEVGDriver.cxx:441
bool IsWeakMix(void) const
double ArgAsDouble(char opt)
bool HitSeaQrk(void) const
Definition: Target.cxx:316
bool IsWeakCC(void) const
TGraph * GetAsTGraph(int np=500, bool xscaling=false, bool inlog=false, double fx=1., double fy=1.) const
Definition: Spline.cxx:490
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
bool IsUQuark(int pdgc)
Definition: PDGUtils.cxx:249
int HitNucPdg(void) const
Definition: Target.cxx:321
int HitQrkPdg(void) const
Definition: Target.cxx:259
bool IsInverseMuDecay(void) const
const int kPdgClusterNP
Definition: PDGCodes.h:192
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
void LoadSplines(void)
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
void FormatXSecGraph(TGraph *g)
const int kPsType
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
int CharmHadronPdg(void) const
Definition: XclsTag.h:50
#define pFATAL
Definition: Messenger.h:57
bool IsStrangeEvent(void) const
Definition: XclsTag.h:51
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:99
string gOptXMLFilename
void SaveToPsFile(void)
string filename
Definition: shutoffs.py:106
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:259
bool IsDiffractive(void) const
bool IsIMDAnnihilation(void) const
const Spline * XSecSumSpline(void) const
Definition: GEVGDriver.h:88
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:279
static XSecSplineList * Instance()
PDGCodeList gOptProbePdgList
double Evaluate(double x) const
Definition: Spline.cxx:362
static const double cm2
Definition: Units.h:77
enum genie::EResonance Resonance_t
void PrintSyntax(void)
string AsString(void) const
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:274
list markers
Definition: styles.py:7
double gEmin
const double kEmin
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgClusterNN
Definition: PDGCodes.h:191
bool IsCharmEvent(void) const
Definition: XclsTag.h:48
string gOptROOTFilename
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
const XML_Char * s
Definition: expat.h:262
Summary information for an interaction.
Definition: Interaction.h:56
int gOptProbePdgCode
int colors[6]
Definition: tools.h:1
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
bool IsWeakNC(void) const
int StrangeHadronPdg(void) const
Definition: XclsTag.h:53
int gOptTgtPdgCode
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void SetEventGeneratorList(string listname)
Definition: GEVGDriver.cxx:349
bool IsNuElectronElastic(void) const
const InteractionList * Interactions(void) const
Definition: GEVGDriver.cxx:335
Float_t E
Definition: plot.C:20
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
Int_t col[ntarg]
Definition: Style.C:29
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
const double j
Definition: BetheBloch.cxx:29
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:55
#define pINFO
Definition: Messenger.h:63
GEVGDriver GetEventGenDriver(void)
Resonance_t Resonance(void) const
Definition: XclsTag.h:62
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:100
void SaveGraphsToRootFile(void)
#define pWARN
Definition: Messenger.h:61
double gEmax
bool IsMEC(void) const
bool IsEM(void) const
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:93
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:82
bool Is2NucleonCluster(int pdgc)
Definition: PDGUtils.cxx:365
PDGCodeList GetPDGCodeListFromString(std::string s)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
string FilterString(string filt, string input)
Definition: StringUtils.cxx:85
bool IsInclusiveStrange(void) const
Definition: XclsTag.cxx:80
int main(int argc, char **argv)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static RunOpt * Instance(void)
Definition: RunOpt.cxx:62
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
double gOptNuEnergy
bool IsCQuark(int pdgc)
Definition: PDGUtils.cxx:264
const Spline * XSecSpline(const Interaction *interaction) const
Definition: GEVGDriver.cxx:489
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
bool HitNucIsSet(void) const
Definition: Target.cxx:300
bool HitQrkIsSet(void) const
Definition: Target.cxx:309
bool IsAntiCQuark(int pdgc)
Definition: PDGUtils.cxx:284
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:63
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:254
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:138
assert(nhit_max >=nhit_nbins)
A vector of Interaction objects.
const InitialState & InitState(void) const
Definition: Interaction.h:69
const char * AsString(Resonance_t res)
resonance id -> string
A vector of EventGeneratorI objects.
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void CreateSplines(int nknots=-1, double emax=-1, bool inLogE=true)
Definition: GEVGDriver.cxx:578
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
void GetCommandLineArgs(int argc, char **argv)
void SaveNtupleToRootFile(void)
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
const Target & Tgt(void) const
Definition: InitialState.h:67
bool IsCoherentElas(void) const
Float_t e
Definition: plot.C:35
bool gWriteOutPlots
enum genie::EXmlParseStatus XmlParserStatus_t
PDGCodeList gOptTgtPdgList
List of cross section vs energy splines.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
XmlParserStatus_t LoadFromXml(const string &filename, bool keep=false)
bool OptionExists(char opt)
was option set?
Root of GENIE utility namespaces.
int kNSplineP
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
bool IsAntiUQuark(int pdgc)
Definition: PDGUtils.cxx:269
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
const int kPdgClusterPP
Definition: PDGCodes.h:193