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] [-l]
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 minimum and 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  -l
31  energy bins in log10 scale
32  -k
33  keep spline knot points (not yet implemented).
34  --message-thresholds
35  Allows users to customize the message stream thresholds.
36  --event-generator-list
37  List of event generators to load in event generation drivers.
38  [default: "Default"].
39 
40  Example:
41 
42  Extract all numu+n, numu+p, numu+O16 cross section splines from the
43  input XML file `mysplines.xml', convert splines into a ROOT format
44  and save them into a single ROOT file.
45 
46  shell$ gspl2root -f mysplines.xml -p 14 -t 1000000010 -o xsec.root
47  shell$ gspl2root -f mysplines.xml -p 14 -t 1000010010 -o xsec.root
48  shell$ gspl2root -f mysplines.xml -p 14 -t 1000080160 -o xsec.root
49 
50  A large number of graphs (one per simulated process + appropriate
51  totals) will be generated in each case. Each set of plots is saved
52  into its own ROOT TDirectory named after the specified initial state.
53 
54  Important Note:
55 
56  The stored graphs can be used for cross section interpolation.
57  Essentially, this _single_ ROOT file can contain _all_ the GENIE cross
58  section `functions' you may need.
59  For instance, the `xsec.root' file generated in the above example will
60  contain a `nu_mu_O16' directory (generated by the last command)
61  which will include cross section graphs for all numu + O16 processes.
62  To extract the numu+O16 DIS CC cross section graph for hit u valence
63  quarks in a bound proton and evaluate the cross section at energy E,
64  then type:
65 
66  root[0] TDirectory * dir = (TDirectory*) file->Get("nu_mu_O16");
67  root[1] TGraph * graph = (TGraph*) dir->Get("dis_cc_p_uval");
68  root[2] cout << graph->Evaluate(E) << endl;
69 
70  See the User Manual for more details.
71 
72 \author Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
73  University of Liverpool & STFC Rutherford Appleton Laboratory
74 
75 \created December 15, 2005
76 
77 \cpright Copyright (c) 2003-2020, The GENIE Collaboration
78  For the full text of the license visit http://copyright.genie-mc.org
79 */
80 //____________________________________________________________________________
81 
82 #include <cassert>
83 #include <string>
84 #include <sstream>
85 #include <vector>
86 
87 #include <TSystem.h>
88 #include <TFile.h>
89 #include <TTree.h>
90 #include <TDirectory.h>
91 #include <TPostScript.h>
92 #include <TCanvas.h>
93 #include <TLegend.h>
94 #include <TGraph.h>
95 #include <TPaveText.h>
96 #include <TString.h>
97 #include <TH1F.h>
98 
113 #include "Framework/Utils/AppInit.h"
114 #include "Framework/Utils/RunOpt.h"
120 
121 
122 using std::string;
123 using std::vector;
124 using std::ostringstream;
125 
126 using namespace genie;
127 using namespace genie::utils;
128 
129 //Prototypes:
130 void LoadSplines (void);
132 void SaveToPsFile (void);
133 void SaveGraphsToRootFile (void);
134 void SaveNtupleToRootFile (void);
135 void GetCommandLineArgs (int argc, char ** argv);
136 void PrintSyntax (void);
138 
139 //User-specified options:
140 string gOptXMLFilename; // input XML filename
141 string gOptROOTFilename; // output ROOT filename
142 PDGCodeList gOptProbePdgList; // list of probe PDG codes
143 PDGCodeList gOptTgtPdgList; // list of target PDG codes
144 int gOptProbePdgCode; // probe PDG code (currently being processed)
145 int gOptTgtPdgCode; // target PDG code
146 bool gWriteOutPlots; // write out a postscript file with plots
147 //bool gKeepSplineKnots; // use spline abscissa points rather than equi-spaced
148 
149 //Globals & constants
150 double gEmin;
151 double gEmax;
152 bool gInlogE;
153 int kNP = 300;
154 int kNSplineP = 1000;
155 const int kPsType = 111; // ps type: portrait
156 
157 //____________________________________________________________________________
158 int main(int argc, char ** argv)
159 {
160  GetCommandLineArgs(argc,argv);
161 
162  if ( ! RunOpt::Instance()->Tune() ) {
163  LOG("gslp2root", pFATAL) << " No TuneId in RunOption";
164  exit(-1);
165  }
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
185  }
186  }
187 
188  return 0;
189 }
190 //____________________________________________________________________________
191 void LoadSplines(void)
192 {
193 // load the cross section splines specified at the cmd line
194 
197  assert(ist == kXmlOK);
198 }
199 //____________________________________________________________________________
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 
207 
208  GEVGDriver evg_driver;
210  evg_driver.Configure(init_state);
211  evg_driver.CreateSplines();
212  evg_driver.CreateXSecSumSpline (100, gEmin, gEmax);
213 
214  return evg_driver;
215 }
216 //____________________________________________________________________________
217 void SaveToPsFile(void)
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 }
515 //____________________________________________________________________________
516 void FormatXSecGraph(TGraph * g)
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 }
522 //____________________________________________________________________________
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() ) { title << "res"; }
602  else if (proc.IsDeepInelastic() ) { title << "dis"; }
603  else if (proc.IsDiffractive() ) { title << "dfr"; }
604  else if (proc.IsCoherentProduction() ) {
605  title << "coh";
606  if ( xcls.NSingleGammas() > 0 ) title << "_gamma" ;
607  else if ( xcls.NPions() > 0 ) title << "_pion" ;
608  else if ( xcls.NRhos() > 0 ) title << "_rho" ;
609  else title << "_other" ;
610  }
611  else if (proc.IsCoherentElastic() ) { title << "cevns"; }
612  else if (proc.IsInverseMuDecay() ) { title << "imd"; }
613  else if (proc.IsIMDAnnihilation() ) { title << "imdanh";}
614  else if (proc.IsNuElectronElastic()) { title << "ve"; }
615  else if (proc.IsGlashowResonance() ) { title << "glres"; }
616  else {
617  LOG("gspl2root", pWARN) << "Process " << proc
618  << " scattering type not recognised: spline not added" ;
619  continue; }
620 
621  if (proc.IsWeakCC()) { title << "_cc"; }
622  else if (proc.IsWeakNC()) { title << "_nc"; }
623  else if (proc.IsWeakMix()) { title << "_ccncmix"; }
624  else if (proc.IsEM() ) { title << "_em"; }
625  else if (proc.IsDarkNeutralCurrent() ) { title << "_dark"; }
626  else {
627  LOG("gspl2root", pWARN) << "Process " << proc
628  << " interaction type has not recongnised: spline not added " ;
629  continue; }
630 
631  if(tgt.HitNucIsSet()) {
632  int hitnuc = tgt.HitNucPdg();
633  if ( pdg::IsProton (hitnuc) ) { title << "_p"; }
634  else if ( pdg::IsNeutron(hitnuc) ) { title << "_n"; }
635  else if ( pdg::Is2NucleonCluster(hitnuc) )
636  {
637  if (hitnuc == kPdgClusterNN) { title << "_nn"; }
638  else if (hitnuc == kPdgClusterNP) { title << "_np"; }
639  else if (hitnuc == kPdgClusterPP) { title << "_pp"; }
640  else {
641  LOG("gspl2root", pWARN) << "Can't handle hit 2-nucleon cluster PDG = " << hitnuc;
642  }
643  }
644  else {
645  LOG("gspl2root", pWARN) << "Can't handle hit nucleon PDG = " << hitnuc;
646  }
647 
648  if(tgt.HitQrkIsSet()) {
649  int qrkpdg = tgt.HitQrkPdg();
650  bool insea = tgt.HitSeaQrk();
651 
652  if ( pdg::IsUQuark(qrkpdg) ) { title << "_u"; }
653  else if ( pdg::IsDQuark(qrkpdg) ) { title << "_d"; }
654  else if ( pdg::IsSQuark(qrkpdg) ) { title << "_s"; }
655  else if ( pdg::IsCQuark(qrkpdg) ) { title << "_c"; }
656  else if ( pdg::IsBQuark(qrkpdg) ) { title << "_b"; }
657  else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
658  else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
659  else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
660  else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
661  else if ( pdg::IsAntiBQuark(qrkpdg) ) { title << "_bbar"; }
662 
663  if(insea) { title << "sea"; }
664  else { title << "val"; }
665  }
666  }
667  if(proc.IsResonant()) {
668  Resonance_t res = xcls.Resonance();
669  string resname = res::AsString(res);
670  resname = str::FilterString(")", resname);
671  resname = str::FilterString("(", resname);
672  title << "_" << resname.substr(3,4) << resname.substr(0,3);
673  }
674 
675  if(xcls.IsStrangeEvent()) {
676  title << "_strange";
677  if(!xcls.IsInclusiveStrange()) { title << xcls.StrangeHadronPdg(); }
678  }
679 
680  if(xcls.IsCharmEvent()) {
681  title << "_charm";
682  if(!xcls.IsInclusiveCharm()) { title << xcls.CharmHadronPdg(); }
683  }
684 
685  if(xcls.IsFinalQuarkEvent()) {
686  int qrkpdg = xcls.FinalQuarkPdg();
687  if ( pdg::IsUQuark(qrkpdg) ) { title << "_u"; }
688  else if ( pdg::IsDQuark(qrkpdg) ) { title << "_d"; }
689  else if ( pdg::IsSQuark(qrkpdg) ) { title << "_s"; }
690  else if ( pdg::IsCQuark(qrkpdg) ) { title << "_c"; }
691  else if ( pdg::IsBQuark(qrkpdg) ) { title << "_b"; }
692  else if ( pdg::IsTQuark(qrkpdg) ) { title << "_t"; }
693  else if ( pdg::IsAntiUQuark(qrkpdg) ) { title << "_ubar"; }
694  else if ( pdg::IsAntiDQuark(qrkpdg) ) { title << "_dbar"; }
695  else if ( pdg::IsAntiSQuark(qrkpdg) ) { title << "_sbar"; }
696  else if ( pdg::IsAntiCQuark(qrkpdg) ) { title << "_cbar"; }
697  else if ( pdg::IsAntiBQuark(qrkpdg) ) { title << "_bbar"; }
698  else if ( pdg::IsAntiTQuark(qrkpdg) ) { title << "_tbar"; }
699  }
700  if(xcls.IsFinalLeptonEvent()) {
701  int leppdg = xcls.FinalLeptonPdg();
702  if ( pdg::IsMuon(leppdg) ) { title << "_mu"; }
703  else if ( pdg::IsElectron(leppdg) ) { title << "_e"; }
704  else if ( pdg::IsTau(leppdg) ) { title << "_tau"; }
705  else if ( pdg::IsPion(leppdg) ) { title << "_had"; }
706  }
707 
708  const Spline * spl = evg_driver.XSecSpline(interaction);
709  for(int i=0; i<kNSplineP; i++) {
710  xs[i] = spl->Evaluate(e[i]) * (1E+38/units::cm2);
711  }
712 
713  TGraph * gr = new TGraph(kNSplineP, e, xs);
714  gr->SetName(title.str().c_str());
715  FormatXSecGraph(gr);
716  gr->SetTitle(spl->GetName());
717 
718  topdir->Add(gr);
719  }
720 
721 
722  //
723  // totals for (anti-)neutrino scattering
724  //
725 
726  bool is_neutrino = pdg::IsNeutralLepton(gOptProbePdgCode);
727 
728  if(is_neutrino) {
729 
730  //
731  // add-up all res channels
732  //
733 
734  double * xsresccp = new double[kNSplineP];
735  double * xsresccn = new double[kNSplineP];
736  double * xsresncp = new double[kNSplineP];
737  double * xsresncn = new double[kNSplineP];
738  for(int i=0; i<kNSplineP; i++) {
739  xsresccp[i] = 0;
740  xsresccn[i] = 0;
741  xsresncp[i] = 0;
742  xsresncn[i] = 0;
743  }
744 
745  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
746  const Interaction * interaction = *ilistiter;
747  const ProcessInfo & proc = interaction->ProcInfo();
748  const InitialState & init = interaction->InitState();
749  const Target & tgt = init.Tgt();
750 
751  const Spline * spl = evg_driver.XSecSpline(interaction);
752 
753  if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
754  for(int i=0; i<kNSplineP; i++) {
755  xsresccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
756  }
757  }
758  if (proc.IsResonant() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
759  for(int i=0; i<kNSplineP; i++) {
760  xsresccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
761  }
762  }
763  if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
764  for(int i=0; i<kNSplineP; i++) {
765  xsresncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
766  }
767  }
768  if (proc.IsResonant() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
769  for(int i=0; i<kNSplineP; i++) {
770  xsresncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
771  }
772  }
773  }
774 
775  TGraph * gr_resccp = new TGraph(kNSplineP, e, xsresccp);
776  gr_resccp->SetName("res_cc_p");
777  FormatXSecGraph(gr_resccp);
778  topdir->Add(gr_resccp);
779  TGraph * gr_resccn = new TGraph(kNSplineP, e, xsresccn);
780  gr_resccn->SetName("res_cc_n");
781  FormatXSecGraph(gr_resccn);
782  topdir->Add(gr_resccn);
783  TGraph * gr_resncp = new TGraph(kNSplineP, e, xsresncp);
784  gr_resncp->SetName("res_nc_p");
785  FormatXSecGraph(gr_resncp);
786  topdir->Add(gr_resncp);
787  TGraph * gr_resncn = new TGraph(kNSplineP, e, xsresncn);
788  gr_resncn->SetName("res_nc_n");
789  FormatXSecGraph(gr_resncn);
790  topdir->Add(gr_resncn);
791 
792  //
793  // add-up all dis channels
794  //
795 
796  double * xsdisccp = new double[kNSplineP];
797  double * xsdisccn = new double[kNSplineP];
798  double * xsdisncp = new double[kNSplineP];
799  double * xsdisncn = new double[kNSplineP];
800  for(int i=0; i<kNSplineP; i++) {
801  xsdisccp[i] = 0;
802  xsdisccn[i] = 0;
803  xsdisncp[i] = 0;
804  xsdisncn[i] = 0;
805  }
806  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
807  const Interaction * interaction = *ilistiter;
808  const ProcessInfo & proc = interaction->ProcInfo();
809  const XclsTag & xcls = interaction->ExclTag();
810  const InitialState & init = interaction->InitState();
811  const Target & tgt = init.Tgt();
812 
813  const Spline * spl = evg_driver.XSecSpline(interaction);
814 
815  if(xcls.IsCharmEvent()) continue;
816 
817  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
818  for(int i=0; i<kNSplineP; i++) {
819  xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
820  }
821  }
822  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
823  for(int i=0; i<kNSplineP; i++) {
824  xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
825  }
826  }
827  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
828  for(int i=0; i<kNSplineP; i++) {
829  xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
830  }
831  }
832  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
833  for(int i=0; i<kNSplineP; i++) {
834  xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
835  }
836  }
837  }
838  TGraph * gr_disccp = new TGraph(kNSplineP, e, xsdisccp);
839  gr_disccp->SetName("dis_cc_p");
840  FormatXSecGraph(gr_disccp);
841  topdir->Add(gr_disccp);
842  TGraph * gr_disccn = new TGraph(kNSplineP, e, xsdisccn);
843  gr_disccn->SetName("dis_cc_n");
844  FormatXSecGraph(gr_disccn);
845  topdir->Add(gr_disccn);
846  TGraph * gr_disncp = new TGraph(kNSplineP, e, xsdisncp);
847  gr_disncp->SetName("dis_nc_p");
848  FormatXSecGraph(gr_disncp);
849  topdir->Add(gr_disncp);
850  TGraph * gr_disncn = new TGraph(kNSplineP, e, xsdisncn);
851  gr_disncn->SetName("dis_nc_n");
852  FormatXSecGraph(gr_disncn);
853  topdir->Add(gr_disncn);
854 
855  //
856  // add-up all charm dis channels
857  //
858 
859  for(int i=0; i<kNSplineP; i++) {
860  xsdisccp[i] = 0;
861  xsdisccn[i] = 0;
862  xsdisncp[i] = 0;
863  xsdisncn[i] = 0;
864  }
865  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
866  const Interaction * interaction = *ilistiter;
867  const ProcessInfo & proc = interaction->ProcInfo();
868  const XclsTag & xcls = interaction->ExclTag();
869  const InitialState & init = interaction->InitState();
870  const Target & tgt = init.Tgt();
871 
872  const Spline * spl = evg_driver.XSecSpline(interaction);
873 
874  if(!xcls.IsCharmEvent()) continue;
875 
876  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) {
877  for(int i=0; i<kNSplineP; i++) {
878  xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
879  }
880  }
881  if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) {
882  for(int i=0; i<kNSplineP; i++) {
883  xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
884  }
885  }
886  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) {
887  for(int i=0; i<kNSplineP; i++) {
888  xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
889  }
890  }
891  if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) {
892  for(int i=0; i<kNSplineP; i++) {
893  xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
894  }
895  }
896  }
897  TGraph * gr_disccp_charm = new TGraph(kNSplineP, e, xsdisccp);
898  gr_disccp_charm->SetName("dis_cc_p_charm");
899  FormatXSecGraph(gr_disccp_charm);
900  topdir->Add(gr_disccp_charm);
901  TGraph * gr_disccn_charm = new TGraph(kNSplineP, e, xsdisccn);
902  gr_disccn_charm->SetName("dis_cc_n_charm");
903  FormatXSecGraph(gr_disccn_charm);
904  topdir->Add(gr_disccn_charm);
905  TGraph * gr_disncp_charm = new TGraph(kNSplineP, e, xsdisncp);
906  gr_disncp_charm->SetName("dis_nc_p_charm");
907  FormatXSecGraph(gr_disncp_charm);
908  topdir->Add(gr_disncp_charm);
909  TGraph * gr_disncn_charm = new TGraph(kNSplineP, e, xsdisncn);
910  gr_disncn_charm->SetName("dis_nc_n_charm");
911  FormatXSecGraph(gr_disncn_charm);
912  topdir->Add(gr_disncn_charm);
913 
914  //
915  // add-up all mec channels
916  //
917 
918  double * xsmeccc = new double[kNSplineP];
919  double * xsmecnc = new double[kNSplineP];
920  for(int i=0; i<kNSplineP; i++) {
921  xsmeccc[i] = 0;
922  xsmecnc[i] = 0;
923  }
924 
925  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
926  const Interaction * interaction = *ilistiter;
927  const ProcessInfo & proc = interaction->ProcInfo();
928 
929  const Spline * spl = evg_driver.XSecSpline(interaction);
930 
931  if (proc.IsMEC() && proc.IsWeakCC()) {
932  for(int i=0; i<kNSplineP; i++) {
933  xsmeccc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
934  }
935  }
936  if (proc.IsMEC() && proc.IsWeakNC()) {
937  for(int i=0; i<kNSplineP; i++) {
938  xsmecnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
939  }
940  }
941  }
942 
943  TGraph * gr_meccc = new TGraph(kNSplineP, e, xsmeccc);
944  gr_meccc->SetName("mec_cc");
945  FormatXSecGraph(gr_meccc);
946  topdir->Add(gr_meccc);
947  TGraph * gr_mecnc = new TGraph(kNSplineP, e, xsmecnc);
948  gr_mecnc->SetName("mec_nc");
949  FormatXSecGraph(gr_mecnc);
950  topdir->Add(gr_mecnc);
951 
952  //
953  // add-up all COH channels
954  //
955 
956  double * xscohcc = new double[kNSplineP];
957  double * xscohnc = new double[kNSplineP];
958  double * xscohtot = new double[kNSplineP];
959  for(int i=0; i<kNSplineP; i++) {
960  xscohcc[i] = 0;
961  xscohnc[i] = 0;
962  xscohtot[i] = 0;
963  }
964 
965  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
966 
967  const Interaction * interaction = *ilistiter;
968  const ProcessInfo & proc = interaction->ProcInfo();
969 
970  const Spline * spl = evg_driver.XSecSpline(interaction);
971 
972  if (proc.IsCoherentProduction() && proc.IsWeakCC()) {
973  for(int i=0; i<kNSplineP; i++) {
974  xscohcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
975  }
976  }
977  if (proc.IsCoherentProduction() && proc.IsWeakNC()) {
978  for(int i=0; i<kNSplineP; i++) {
979  xscohnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
980  }
981  }
982  if ( proc.IsCoherentProduction() ) {
983  for(int i=0; i<kNSplineP; i++) {
984  xscohtot[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
985  }
986  }
987 
988  }
989 
990  TGraph * gr_cohcc = new TGraph(kNSplineP, e, xscohcc);
991  gr_cohcc->SetName("coh_cc");
992  FormatXSecGraph(gr_cohcc);
993  topdir->Add(gr_cohcc);
994 
995  TGraph * gr_cohnc = new TGraph(kNSplineP, e, xscohnc);
996  gr_cohnc->SetName("coh_nc");
997  FormatXSecGraph(gr_cohnc);
998  topdir->Add(gr_cohnc);
999 
1000  TGraph * gr_cohtot = new TGraph(kNSplineP, e, xscohtot);
1001  gr_cohtot->SetName("coh");
1002  FormatXSecGraph(gr_cohtot);
1003  topdir->Add(gr_cohtot);
1004 
1005  //
1006  // total cross sections
1007  //
1008  double * xstotcc = new double[kNSplineP];
1009  double * xstotccp = new double[kNSplineP];
1010  double * xstotccn = new double[kNSplineP];
1011  double * xstotnc = new double[kNSplineP];
1012  double * xstotncp = new double[kNSplineP];
1013  double * xstotncn = new double[kNSplineP];
1014  for(int i=0; i<kNSplineP; i++) {
1015  xstotcc [i] = 0;
1016  xstotccp[i] = 0;
1017  xstotccn[i] = 0;
1018  xstotnc [i] = 0;
1019  xstotncp[i] = 0;
1020  xstotncn[i] = 0;
1021  }
1022  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1023  const Interaction * interaction = *ilistiter;
1024  const ProcessInfo & proc = interaction->ProcInfo();
1025  const InitialState & init = interaction->InitState();
1026  const Target & tgt = init.Tgt();
1027 
1028  const Spline * spl = evg_driver.XSecSpline(interaction);
1029 
1030  bool iscc = proc.IsWeakCC();
1031  bool isnc = proc.IsWeakNC();
1032  bool offp = pdg::IsProton (tgt.HitNucPdg());
1033  bool offn = pdg::IsNeutron(tgt.HitNucPdg());
1034 
1035  if (iscc && offp) {
1036  for(int i=0; i<kNSplineP; i++) {
1037  xstotccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1038  }
1039  }
1040  if (iscc && offn) {
1041  for(int i=0; i<kNSplineP; i++) {
1042  xstotccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1043  }
1044  }
1045  if (isnc && offp) {
1046  for(int i=0; i<kNSplineP; i++) {
1047  xstotncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1048  }
1049  }
1050  if (isnc && offn) {
1051  for(int i=0; i<kNSplineP; i++) {
1052  xstotncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1053  }
1054  }
1055 
1056  if (iscc) {
1057  for(int i=0; i<kNSplineP; i++) {
1058  xstotcc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1059  }
1060  }
1061  if (isnc) {
1062  for(int i=0; i<kNSplineP; i++) {
1063  xstotnc[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1064  }
1065  }
1066  }
1067 
1068  TGraph * gr_totcc = new TGraph(kNSplineP, e, xstotcc);
1069  gr_totcc->SetName("tot_cc");
1070  FormatXSecGraph(gr_totcc);
1071  topdir->Add(gr_totcc);
1072  TGraph * gr_totccp = new TGraph(kNSplineP, e, xstotccp);
1073  gr_totccp->SetName("tot_cc_p");
1074  FormatXSecGraph(gr_totccp);
1075  topdir->Add(gr_totccp);
1076  TGraph * gr_totccn = new TGraph(kNSplineP, e, xstotccn);
1077  gr_totccn->SetName("tot_cc_n");
1078  FormatXSecGraph(gr_totccn);
1079  topdir->Add(gr_totccn);
1080  TGraph * gr_totnc = new TGraph(kNSplineP, e, xstotnc);
1081  gr_totnc->SetName("tot_nc");
1082  FormatXSecGraph(gr_totnc);
1083  topdir->Add(gr_totnc);
1084  TGraph * gr_totncp = new TGraph(kNSplineP, e, xstotncp);
1085  gr_totncp->SetName("tot_nc_p");
1086  FormatXSecGraph(gr_totncp);
1087  topdir->Add(gr_totncp);
1088  TGraph * gr_totncn = new TGraph(kNSplineP, e, xstotncn);
1089  gr_totncn->SetName("tot_nc_n");
1090  FormatXSecGraph(gr_totncn);
1091  topdir->Add(gr_totncn);
1092 
1093  delete [] e;
1094  delete [] xs;
1095  delete [] xsresccp;
1096  delete [] xsresccn;
1097  delete [] xsresncp;
1098  delete [] xsresncn;
1099  delete [] xsdisccp;
1100  delete [] xsdisccn;
1101  delete [] xsdisncp;
1102  delete [] xsdisncn;
1103  delete [] xscohcc;
1104  delete [] xscohnc;
1105  delete [] xscohtot;
1106  delete [] xstotcc;
1107  delete [] xstotccp;
1108  delete [] xstotccn;
1109  delete [] xstotnc;
1110  delete [] xstotncp;
1111  delete [] xstotncn;
1112 
1113  }// neutrinos
1114 
1115 
1116  //
1117  // totals for charged lepton scattering
1118  //
1119 
1120  bool is_charged_lepton = pdg::IsChargedLepton(gOptProbePdgCode);
1121 
1122  if(is_charged_lepton) {
1123 
1124  //
1125  // add-up all res channels
1126  //
1127 
1128  double * xsresemp = new double[kNSplineP];
1129  double * xsresemn = new double[kNSplineP];
1130  for(int i=0; i<kNSplineP; i++) {
1131  xsresemp[i] = 0;
1132  xsresemn[i] = 0;
1133  }
1134 
1135  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1136  const Interaction * interaction = *ilistiter;
1137  const ProcessInfo & proc = interaction->ProcInfo();
1138  const InitialState & init = interaction->InitState();
1139  const Target & tgt = init.Tgt();
1140 
1141  const Spline * spl = evg_driver.XSecSpline(interaction);
1142 
1143  if (proc.IsResonant() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1144  for(int i=0; i<kNSplineP; i++) {
1145  xsresemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1146  }
1147  }
1148  if (proc.IsResonant() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1149  for(int i=0; i<kNSplineP; i++) {
1150  xsresemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1151  }
1152  }
1153  }
1154 
1155  TGraph * gr_resemp = new TGraph(kNSplineP, e, xsresemp);
1156  gr_resemp->SetName("res_em_p");
1157  FormatXSecGraph(gr_resemp);
1158  topdir->Add(gr_resemp);
1159  TGraph * gr_resemn = new TGraph(kNSplineP, e, xsresemn);
1160  gr_resemn->SetName("res_em_n");
1161  FormatXSecGraph(gr_resemn);
1162  topdir->Add(gr_resemn);
1163 
1164  //
1165  // add-up all dis channels
1166  //
1167 
1168  double * xsdisemp = new double[kNSplineP];
1169  double * xsdisemn = new double[kNSplineP];
1170  for(int i=0; i<kNSplineP; i++) {
1171  xsdisemp[i] = 0;
1172  xsdisemn[i] = 0;
1173  }
1174  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1175  const Interaction * interaction = *ilistiter;
1176  const ProcessInfo & proc = interaction->ProcInfo();
1177  const XclsTag & xcls = interaction->ExclTag();
1178  const InitialState & init = interaction->InitState();
1179  const Target & tgt = init.Tgt();
1180 
1181  const Spline * spl = evg_driver.XSecSpline(interaction);
1182 
1183  if(xcls.IsCharmEvent()) continue;
1184 
1185  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1186  for(int i=0; i<kNSplineP; i++) {
1187  xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1188  }
1189  }
1190  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1191  for(int i=0; i<kNSplineP; i++) {
1192  xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1193  }
1194  }
1195  }
1196  TGraph * gr_disemp = new TGraph(kNSplineP, e, xsdisemp);
1197  gr_disemp->SetName("dis_em_p");
1198  FormatXSecGraph(gr_disemp);
1199  topdir->Add(gr_disemp);
1200  TGraph * gr_disemn = new TGraph(kNSplineP, e, xsdisemn);
1201  gr_disemn->SetName("dis_em_n");
1202  FormatXSecGraph(gr_disemn);
1203  topdir->Add(gr_disemn);
1204 
1205  //
1206  // add-up all charm dis channels
1207  //
1208 
1209  for(int i=0; i<kNSplineP; i++) {
1210  xsdisemp[i] = 0;
1211  xsdisemn[i] = 0;
1212  }
1213  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1214  const Interaction * interaction = *ilistiter;
1215  const ProcessInfo & proc = interaction->ProcInfo();
1216  const XclsTag & xcls = interaction->ExclTag();
1217  const InitialState & init = interaction->InitState();
1218  const Target & tgt = init.Tgt();
1219 
1220  const Spline * spl = evg_driver.XSecSpline(interaction);
1221 
1222  if(!xcls.IsCharmEvent()) continue;
1223 
1224  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsProton(tgt.HitNucPdg())) {
1225  for(int i=0; i<kNSplineP; i++) {
1226  xsdisemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1227  }
1228  }
1229  if (proc.IsDeepInelastic() && proc.IsEM() && pdg::IsNeutron(tgt.HitNucPdg())) {
1230  for(int i=0; i<kNSplineP; i++) {
1231  xsdisemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1232  }
1233  }
1234  }
1235  TGraph * gr_disemp_charm = new TGraph(kNSplineP, e, xsdisemp);
1236  gr_disemp_charm->SetName("dis_em_p_charm");
1237  FormatXSecGraph(gr_disemp_charm);
1238  topdir->Add(gr_disemp_charm);
1239  TGraph * gr_disemn_charm = new TGraph(kNSplineP, e, xsdisemn);
1240  gr_disemn_charm->SetName("dis_em_n_charm");
1241  FormatXSecGraph(gr_disemn_charm);
1242  topdir->Add(gr_disemn_charm);
1243 
1244  //
1245  // total cross sections
1246  //
1247  double * xstotem = new double[kNSplineP];
1248  double * xstotemp = new double[kNSplineP];
1249  double * xstotemn = new double[kNSplineP];
1250  for(int i=0; i<kNSplineP; i++) {
1251  xstotem [i] = 0;
1252  xstotemp[i] = 0;
1253  xstotemn[i] = 0;
1254  }
1255  for(ilistiter = ilist->begin(); ilistiter != ilist->end(); ++ilistiter) {
1256  const Interaction * interaction = *ilistiter;
1257  const ProcessInfo & proc = interaction->ProcInfo();
1258  const InitialState & init = interaction->InitState();
1259  const Target & tgt = init.Tgt();
1260 
1261  const Spline * spl = evg_driver.XSecSpline(interaction);
1262 
1263  bool isem = proc.IsEM();
1264  bool offp = pdg::IsProton (tgt.HitNucPdg());
1265  bool offn = pdg::IsNeutron(tgt.HitNucPdg());
1266 
1267  if (isem && offp) {
1268  for(int i=0; i<kNSplineP; i++) {
1269  xstotemp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1270  }
1271  }
1272  if (isem && offn) {
1273  for(int i=0; i<kNSplineP; i++) {
1274  xstotemn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1275  }
1276  }
1277  if (isem) {
1278  for(int i=0; i<kNSplineP; i++) {
1279  xstotem[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2));
1280  }
1281  }
1282  }
1283 
1284  TGraph * gr_totem = new TGraph(kNSplineP, e, xstotem);
1285  gr_totem->SetName("tot_em");
1286  FormatXSecGraph(gr_totem);
1287  topdir->Add(gr_totem);
1288  TGraph * gr_totemp = new TGraph(kNSplineP, e, xstotemp);
1289  gr_totemp->SetName("tot_em_p");
1290  FormatXSecGraph(gr_totemp);
1291  topdir->Add(gr_totemp);
1292  TGraph * gr_totemn = new TGraph(kNSplineP, e, xstotemn);
1293  gr_totemn->SetName("tot_em_n");
1294  FormatXSecGraph(gr_totemn);
1295  topdir->Add(gr_totemn);
1296 
1297  delete [] e;
1298  delete [] xs;
1299  delete [] xsresemp;
1300  delete [] xsresemn;
1301  delete [] xsdisemp;
1302  delete [] xsdisemn;
1303  delete [] xstotem;
1304  delete [] xstotemp;
1305  delete [] xstotemn;
1306 
1307  }// charged leptons
1308 
1309  topdir->Write();
1310 
1311  if(froot) {
1312  froot->Close();
1313  delete froot;
1314  }
1315 }
1316 //____________________________________________________________________________
1318 {
1319 /*
1320  //-- create cross section ntuple
1321  //
1322  double brXSec;
1323  double brEv;
1324  bool brIsQel;
1325  bool brIsRes;
1326  bool brIsDis;
1327  bool brIsCoh;
1328  bool brIsImd;
1329  bool brIsNuEEl;
1330  bool brIsCC;
1331  bool brIsNC;
1332  int brNucleon;
1333  int brQrk;
1334  bool brIsSeaQrk;
1335  int brRes;
1336  bool brCharm;
1337  TTree * xsecnt = new TTree("xsecnt",dtitle.str().c_str());
1338  xsecnt->Branch("xsec", &brXSec, "xsec/D");
1339  xsecnt->Branch("e", &brEv, "e/D");
1340  xsecnt->Branch("qel", &brIsQel, "qel/O");
1341  xsecnt->Branch("res", &brIsRes, "res/O");
1342  xsecnt->Branch("dis", &brIsDis, "dis/O");
1343  xsecnt->Branch("coh", &brIsCoh, "coh/O");
1344  xsecnt->Branch("imd", &brIsImd, "imd/O");
1345  xsecnt->Branch("veel", &brIsNuEEl, "veel/O");
1346  xsecnt->Branch("cc", &brIsCC, "cc/O");
1347  xsecnt->Branch("nc", &brIsNC, "nc/O");
1348  xsecnt->Branch("nuc", &brNucleon, "nuc/I");
1349  xsecnt->Branch("qrk", &brQrk, "qrk/I");
1350  xsecnt->Branch("sea", &brIsSeaQrk, "sea/O");
1351  xsecnt->Branch("res", &brRes, "res/I");
1352  xsecnt->Branch("charm", &brCharm, "charm/O");
1353 */
1354 }
1355 //____________________________________________________________________________
1356 void GetCommandLineArgs(int argc, char ** argv)
1357 {
1358  LOG("gspl2root", pINFO) << "Parsing command line arguments";
1359 
1360  // Common run options. Set defaults and read.
1361  RunOpt::Instance()->ReadFromCommandLine(argc,argv);
1362 
1363  // Parse run options for this app
1364 
1365  CmdLnArgParser parser(argc,argv);
1366 
1367  // input XML file name:
1368  if( parser.OptionExists('f') ) {
1369  LOG("gspl2root", pINFO) << "Reading input XML filename";
1370  gOptXMLFilename = parser.ArgAsString('f');
1371  } else {
1372  LOG("gspl2root", pFATAL) << "Unspecified input XML file!";
1373  PrintSyntax();
1374  exit(1);
1375  }
1376 
1377  // probe PDG code:
1378  if( parser.OptionExists('p') ) {
1379  LOG("gspl2root", pINFO) << "Reading probe PDG code";
1380  gOptProbePdgList = GetPDGCodeListFromString(parser.ArgAsString('p'));
1381  } else {
1382  LOG("gspl2root", pFATAL)
1383  << "Unspecified probe PDG code - Exiting";
1384  PrintSyntax();
1385  exit(1);
1386  }
1387 
1388  // target PDG code:
1389  if( parser.OptionExists('t') ) {
1390  LOG("gspl2root", pINFO) << "Reading target PDG code";
1391  gOptTgtPdgList = GetPDGCodeListFromString(parser.ArgAsString('t'));
1392  } else {
1393  LOG("gspl2root", pFATAL)
1394  << "Unspecified target PDG code - Exiting";
1395  PrintSyntax();
1396  exit(1);
1397  }
1398 
1399  // min,max neutrino energy
1400  if( parser.OptionExists('e') ) {
1401  LOG("gspl2root", pINFO) << "Reading neutrino energy";
1402  string nue = parser.ArgAsString('e');
1403  // is it just a value or a range (comma separated set of values)
1404  if(nue.find(",") != string::npos) {
1405  // split the comma separated list
1406  vector<string> nurange = utils::str::Split(nue, ",");
1407  assert(nurange.size() == 2);
1408  gEmin = atof(nurange[0].c_str());
1409  gEmax = atof(nurange[1].c_str());
1410  } else {
1411  const Registry * val_reg = AlgConfigPool::Instance() -> CommonList( "Param", "Validation" ) ;
1412  gEmin = val_reg -> GetDouble( "GVLD-Emin" ) ;
1413  gEmax = atof(nue.c_str());
1414  LOG("gspl2root", pDEBUG)
1415  << "Unspecified Emin - Setting to " << gEmin << " GeV as per configuration";
1416  }
1417  } else {
1418  const Registry * val_reg = AlgConfigPool::Instance() -> CommonList("Param", "Validation" ) ;
1419  gEmin = val_reg -> GetDouble( "GVLD-Emin" ) ;
1420  gEmax = 100;
1421  LOG("gspl2root", pDEBUG)
1422  << "Unspecified Emin,Emax - Setting to " << gEmin << ",100 GeV ";
1423 
1424  }
1425  assert(gEmin<gEmax);
1426 
1427  // output ROOT file name:
1428  if( parser.OptionExists('o') ) {
1429  LOG("gspl2root", pINFO) << "Reading output ROOT filename";
1430  gOptROOTFilename = parser.ArgAsString('o');
1431  } else {
1432  LOG("gspl2root", pDEBUG)
1433  << "Unspecified output ROOT file. Using default: gxsec.root";
1434  gOptROOTFilename = "gxsec.root";
1435  }
1436 
1437  // write-out a PS file with plots
1438  gWriteOutPlots = parser.OptionExists('w');
1439 
1440  // use same abscissa points as splines
1441  //not yet//gKeepSplineKnots = parser.OptionExists('k');
1442 
1443  gInlogE = parser.OptionExists('l');
1444 
1445  // print the options you got from command line arguments
1446  LOG("gspl2root", pINFO) << "Command line arguments:";
1447  LOG("gspl2root", pINFO) << " Input XML file = " << gOptXMLFilename;
1448  LOG("gspl2root", pINFO) << " Probe PDG code = " << gOptProbePdgCode;
1449  LOG("gspl2root", pINFO) << " Target PDG code = " << gOptTgtPdgCode;
1450  LOG("gspl2root", pINFO) << " Min neutrino E = " << gEmin;
1451  LOG("gspl2root", pINFO) << " Max neutrino E = " << gEmax;
1452  LOG("gspl2root", pINFO) << " In logE = " << gInlogE;
1453  //not yet//LOG("gspl2root", pINFO) << " Keep spline knots = " << (gKeepSplineKnots?"true":"false");
1454 }
1455 //____________________________________________________________________________
1456 void PrintSyntax(void)
1457 {
1458  LOG("gspl2root", pNOTICE)
1459  << "\n\n" << "Syntax:" << "\n"
1460  << " gspl2root -f xml_file -p probe_pdg -t target_pdg"
1461  << " [-e emin,emax] [-o output_root_file] [-w] [-l]\n"
1462  << " [--message-thresholds xml_file]\n";
1463 }
1464 //____________________________________________________________________________
1466 {
1467  vector<string> isvec = utils::str::Split(s,",");
1468 
1469  // fill in the PDG code list
1470  PDGCodeList list;
1472  for(iter = isvec.begin(); iter != isvec.end(); ++iter) {
1473  list.push_back( atoi(iter->c_str()) );
1474  }
1475  return list;
1476 
1477 }
1478 //____________________________________________________________________________
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
bool IsBQuark(int pdgc)
Definition: PDGUtils.cxx:283
int kNP
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:323
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
Definition: GEVGDriver.cxx:440
bool IsFinalLeptonEvent(void) const
Definition: XclsTag.h:73
bool IsWeakMix(void) const
bool HitSeaQrk(void) const
Definition: Target.cxx:299
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:489
bool IsFinalQuarkEvent(void) const
Definition: XclsTag.h:71
int FinalLeptonPdg(void) const
Definition: XclsTag.h:74
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
static constexpr double g
Definition: Units.h:144
bool IsUQuark(int pdgc)
Definition: PDGUtils.cxx:263
int HitNucPdg(void) const
Definition: Target.cxx:304
bool IsDarkNeutralCurrent(void) const
std::string string
Definition: nybbler.cc:12
int HitQrkPdg(void) const
Definition: Target.cxx:242
int NPions(void) const
Definition: XclsTag.h:62
bool IsInverseMuDecay(void) const
const int kPdgClusterNP
Definition: PDGCodes.h:215
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:99
void LoadSplines(void)
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
void FormatXSecGraph(TGraph *g)
const int kPsType
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:46
int CharmHadronPdg(void) const
Definition: XclsTag.h:52
#define pFATAL
Definition: Messenger.h:56
bool IsStrangeEvent(void) const
Definition: XclsTag.h:53
struct vector vector
ChannelGroupService::Name Name
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:98
bool gInlogE
string gOptXMLFilename
bool IsTQuark(int pdgc)
Definition: PDGUtils.cxx:288
void SaveToPsFile(void)
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:273
bool IsDiffractive(void) const
intermediate_table::const_iterator const_iterator
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
bool IsIMDAnnihilation(void) const
const Spline * XSecSumSpline(void) const
Definition: GEVGDriver.h:87
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:303
init
Definition: train.py:42
static XSecSplineList * Instance()
PDGCodeList gOptProbePdgList
double Evaluate(double x) const
Definition: Spline.cxx:361
string filename
Definition: train.py:213
enum genie::EResonance Resonance_t
void PrintSyntax(void)
string AsString(void) const
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:298
double gEmin
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
A list of PDG codes.
Definition: PDGCodeList.h:32
const int kPdgClusterNN
Definition: PDGCodes.h:214
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
string gOptROOTFilename
bool exists(std::string path)
int NSingleGammas(void) const
Definition: XclsTag.h:64
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
Summary information for an interaction.
Definition: Interaction.h:56
int gOptProbePdgCode
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
bool IsWeakNC(void) const
const double e
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
int gOptTgtPdgCode
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsCoherentElastic(void) const
void SetEventGeneratorList(string listname)
Definition: GEVGDriver.cxx:348
bool IsNuElectronElastic(void) const
const InteractionList * Interactions(void) const
Definition: GEVGDriver.cxx:334
static constexpr double cm2
Definition: Units.h:69
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
size_t size
Definition: lodepng.cpp:55
bool IsAntiTQuark(int pdgc)
Definition: PDGUtils.cxx:318
bool IsTau(int pdgc)
Definition: PDGUtils.cxx:205
bool IsAntiBQuark(int pdgc)
Definition: PDGUtils.cxx:313
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:54
#define pINFO
Definition: Messenger.h:62
bool IsMuon(int pdgc)
Definition: PDGUtils.cxx:195
static constexpr double ps
Definition: Units.h:99
GEVGDriver GetEventGenDriver(void)
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:92
void SaveGraphsToRootFile(void)
#define pWARN
Definition: Messenger.h:60
double gEmax
bool IsMEC(void) const
bool IsEM(void) const
int FinalQuarkPdg(void) const
Definition: XclsTag.h:72
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:92
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
bool Is2NucleonCluster(int pdgc)
Definition: PDGUtils.cxx:399
PDGCodeList GetPDGCodeListFromString(std::string s)
string FilterString(string filt, string input)
Definition: StringUtils.cxx:79
bool IsInclusiveStrange(void) const
Definition: XclsTag.cxx:71
int main(int argc, char **argv)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
static RunOpt * Instance(void)
Definition: RunOpt.cxx:54
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
bool IsCQuark(int pdgc)
Definition: PDGUtils.cxx:278
const Spline * XSecSpline(const Interaction *interaction) const
Definition: GEVGDriver.cxx:488
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:32
bool HitNucIsSet(void) const
Definition: Target.cxx:283
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
bool IsAntiCQuark(int pdgc)
Definition: PDGUtils.cxx:308
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:54
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:268
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:137
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:577
list x
Definition: train.py:276
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:99
void GetCommandLineArgs(int argc, char **argv)
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
void SaveNtupleToRootFile(void)
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
const Target & Tgt(void) const
Definition: InitialState.h:66
bool gWriteOutPlots
enum genie::EXmlParseStatus XmlParserStatus_t
PDGCodeList gOptTgtPdgList
bool IsGlashowResonance(void) const
List of cross section vs energy splines.
double GetDouble(xmlDocPtr xml_doc, string node_path)
int NRhos(void) const
Definition: XclsTag.h:63
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?
static QCString * s
Definition: config.cpp:1042
Root of GENIE utility namespaces.
bool IsElectron(int pdgc)
Definition: PDGUtils.cxx:185
int kNSplineP
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
static AlgConfigPool * Instance()
bool IsAntiUQuark(int pdgc)
Definition: PDGUtils.cxx:293
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const int kPdgClusterPP
Definition: PDGCodes.h:216