GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gXSecComp.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gxscomp
5 
6 \brief A simple utility that plots the pre-calculated cross-sections used as
7  input for event generation. Can also compare against a reference set of
8  such pre-computed cross-sections.
9 
10  Syntax:
11  gxscomp
12  -f xsec_file[,label]
13  [-r reference_xsec_file[,label]]
14  [-o output]
15 
16  Options:
17  [] Denotes an optional argument.
18  -f Specifies a ROOT file with GENIE cross section graphs.
19  -r Specifies a reference file with GENIE cross section graphs.
20  -o Specifies the output filename [default: xsec.ps]
21 
22  Notes:
23  The input ROOT files are the ones generated by GENIE's gspl2root
24  utility. See the GENIE User Manual for more details.
25 
26  Example:
27  To compare cross sections in xsec-v2_4.root and xsec-v2_2.root
28  type:
29  shell$ gxscomp -f xsec-v2_4.root -r xsec-v2_2.root
30 
31 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
32  University of Liverpool
33 
34 \created June 06, 2008
35 
36 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
37  For the full text of the license visit http://copyright.genie-mc.org
38 
39 */
40 //____________________________________________________________________________
41 
42 #include <cassert>
43 #include <sstream>
44 #include <string>
45 #include <vector>
46 
47 #include <TSystem.h>
48 #include <TFile.h>
49 #include <TDirectory.h>
50 #include <TGraph.h>
51 #include <TPostScript.h>
52 #include <TH1D.h>
53 #include <TMath.h>
54 #include <TCanvas.h>
55 #include <TPavesText.h>
56 #include <TText.h>
57 #include <TStyle.h>
58 #include <TLegend.h>
59 #include <TObjString.h>
60 
61 #include "Framework/Conventions/GBuild.h"
62 #include "Framework/Messenger/Messenger.h"
63 #include "Framework/Numerical/RandomGen.h"
64 #include "Framework/ParticleData/PDGUtils.h"
65 #include "Framework/ParticleData/PDGCodes.h"
66 #include "Framework/Utils/CmdLnArgParser.h"
67 #include "Framework/Utils/StringUtils.h"
68 #include "Framework/Utils/Style.h"
69 
70 using std::ostringstream;
71 using std::string;
72 using std::vector;
73 
74 using namespace genie;
75 
76 // function prototypes
77 void Init (void);
78 void End (void);
79 void OpenDir (void);
80 void DirNameToProbe (void);
81 void DirNameToTarget (void);
82 void MakePlots (void);
83 void MakePlotsCurrDir (void);
84 TH1F * DrawFrame (TGraph * gr0, TGraph * gr1, TPad * p, const char * xt, const char * yt, double yminsc, double ymaxsc);
85 TH1F * DrawFrame (double xmin, double xmax, double ymin, double ymax, TPad * p, const char * xt, const char * yt);
86 void Draw (const char * plot, const char * title);
87 void Draw (TGraph* gr, const char * opt);
88 TGraph* TrimGraph (TGraph * gr, int max_np_per_decade);
89 TGraph* DrawRatio (TGraph * gr0, TGraph * gr1);
90 void GetCommandLineArgs (int argc, char ** argv);
91 void PrintSyntax (void);
92 bool CheckRootFilename (string filename);
93 string OutputFileName (string input_file_name);
94 
95 // command-line arguments
96 string gOptXSecFilename_curr = ""; // (-f) input ROOT cross section file
97 string gOptXSecFilename_ref0 = ""; // (-r) input ROOT cross section file (ref.)
98 string gOptOutputFilename = ""; // (-o) output PS file
99 bool gOptHaveRef;
100 
101 // globals
102 TFile * gXSecFile_curr = 0;
103 TFile * gXSecFile_ref0 = 0;
104 TDirectory * gDirCurr = 0;
105 TDirectory * gDirRef0 = 0;
106 string gLabelCurr = "";
107 string gLabelRef0 = "";
108 string gDirName = "";
109 TPostScript * gPS = 0;
110 TCanvas * gC = 0;
111 TPad * gPadTitle = 0;
112 TPad * gPadXSecs = 0;
113 TPad * gPadRatio = 0;
114 TLegend * gLS = 0;
115 string gCurrProbeLbl = "";
116 int gCurrProbePdg = 0;
117 bool gCurrProbeIsNu = false;
118 bool gCurrProbeIsNuBar = false;
119 string gCurrTargetLbl = "";
120 bool gCurrTargetHasP = false;
121 bool gCurrTargetHasN = false;
122 bool gCurrTargetIsFreeNuc = false;
123 
124 //_________________________________________________________________________________
125 int main(int argc, char ** argv)
126 {
127  GetCommandLineArgs (argc,argv);
128  utils::style::SetDefaultStyle();
129  Init();
130  MakePlots();
131  End();
132 
133  LOG("gxscomp", pINFO) << "Done!";
134 
135  return 0;
136 }
137 //_________________________________________________________________________________
138 void Init(void)
139 {
140  gC = new TCanvas("c","",20,20,500,650);
141  gC->SetBorderMode(0);
142  gC->SetFillColor(0);
143 
144  // create output file
145  gPS = new TPostScript(gOptOutputFilename.c_str(), 111);
146 
147  // front page
148  gPS->NewPage();
149  gC->Range(0,0,100,100);
150  TPavesText hdr(10,40,90,70,3,"tr");
151  hdr.AddText("GENIE cross sections");
152  hdr.AddText(" ");
153  hdr.AddText(" ");
154  hdr.AddText("Plotting data from: ");
155  hdr.AddText(gOptXSecFilename_curr.c_str());
156  if(gOptHaveRef) {
157  hdr.AddText(" ");
158  hdr.AddText("Comparing with reference data (red circles) from: ");
159  hdr.AddText(gOptXSecFilename_ref0.c_str());
160  } else {
161  hdr.AddText(" ");
162  }
163  hdr.AddText(" ");
164  hdr.AddText(" ");
165  hdr.Draw();
166  gC->Update();
167 
168  gPS->NewPage();
169 
170  //
171  if(gOptHaveRef) {
172  gPadTitle = new TPad("gPadTitle","",0.05,0.90,0.95,0.97);
173  gPadXSecs = new TPad("gPadXSecs","",0.05,0.35,0.95,0.88);
174  gPadRatio = new TPad("gPadRatio","",0.05,0.03,0.95,0.32);
175  } else {
176  gPadTitle = new TPad("gPadTitle","",0.05,0.90,0.95,0.97);
177  gPadXSecs = new TPad("gPadXSecs","",0.05,0.05,0.95,0.88);
178  gPadRatio = new TPad("gPadRatio","",0.05,0.03,0.95,0.04);
179  }
180 
181  gPadTitle->Range(0,0,100,100);
182 
183  gPadTitle->SetBorderMode(0);
184  gPadTitle->SetFillColor(0);
185  gPadXSecs->SetBorderMode(0);
186  gPadXSecs->SetFillColor(0);
187  gPadRatio->SetBorderMode(0);
188  gPadRatio->SetFillColor(0);
189 
190  gPadXSecs->SetGridx();
191  gPadXSecs->SetGridy();
192  gPadXSecs->SetLogx();
193  gPadXSecs->SetLogy();
194  gPadRatio->SetGridx();
195  gPadRatio->SetGridy();
196  gPadRatio->SetLogx();
197 
198  gPadTitle->Draw();
199  gPadXSecs->Draw();
200  gPadRatio->Draw();
201 
202  gPadXSecs->cd();
203 
204  gLS = new TLegend(0.80,0.25,0.90,0.45);
205  gLS -> SetFillColor(0);
206 //gLS -> SetFillStyle(0);
207  gLS -> SetBorderSize(0);
208 }
209 //_________________________________________________________________________________
210 void End(void)
211 {
212  gPS->Close();
213 
214  delete gC;
215  delete gPS;
216  delete gLS;
217 }
218 //_________________________________________________________________________________
219 void OpenDir(void)
220 {
221  gDirCurr = (TDirectory *) gXSecFile_curr->Get(gDirName.c_str());
222  gDirRef0 = 0;
223  if(gXSecFile_ref0) {
224  gDirRef0 = (TDirectory *) gXSecFile_ref0->Get(gDirName.c_str());
225  }
226 
227  if(!gDirRef0) {
228  LOG("gxscomp", pINFO) << "No reference plots will be shown.";
229  }
230 }
231 //_________________________________________________________________________________
232 void DirNameToProbe(void)
233 {
234 // Figure out the probe type from the input directory name
235 //
236 
237  string dirname = gDirName;
238 
239  int pdg = 0;
240  string label = "";
241 
242  if( dirname.find ("nu_e_bar") != string::npos )
243  {
244  label = "#bar{#nu_{e}}";
245  pdg = kPdgAntiNuE;
246  }
247  else
248  if ( dirname.find ("nu_e") != string::npos )
249  {
250  label = "#nu_{e}";
251  pdg = kPdgNuE;
252  }
253  else
254  if ( dirname.find ("nu_mu_bar") != string::npos )
255  {
256  label = "#bar{#nu_{#mu}}";
257  pdg = kPdgAntiNuMu;
258  }
259  else
260  if ( dirname.find ("nu_mu") != string::npos )
261  {
262  label ="#nu_{#mu}";
263  pdg = kPdgNuMu;
264  }
265  else
266  if ( dirname.find ("nu_tau_bar") != string::npos )
267  {
268  label = "#bar{#nu_{#tau}}";
269  pdg = kPdgAntiNuTau;
270  }
271  else
272  if ( dirname.find ("nu_tau") != string::npos )
273  {
274  label = "#nu_{#tau}";
275  pdg = kPdgNuTau;
276  }
277 
278  gCurrProbeLbl = label;
279  gCurrProbePdg = pdg;
280  gCurrProbeIsNu = pdg::IsNeutrino (pdg);
281  gCurrProbeIsNuBar = pdg::IsAntiNeutrino (pdg);
282 }
283 //_________________________________________________________________________________
284 void DirNameToTarget(void)
285 {
286 // Figure out the target type from the input directory name
287 //
288 
289  string label;
290  string dirname = gDirName;
291 
292  if ( gCurrProbePdg == kPdgAntiNuE )
293  {
294  label = dirname.substr ( dirname.find("nu_e_bar")+9, dirname.length() );
295  }
296  else
297  if ( gCurrProbePdg == kPdgNuE )
298  {
299  label = dirname.substr ( dirname.find("nu_e")+5, dirname.length() );
300  }
301  else
302  if ( gCurrProbePdg == kPdgAntiNuMu )
303  {
304  label = dirname.substr ( dirname.find("nu_mu_bar")+10, dirname.length() );
305  }
306  else
307  if ( gCurrProbePdg == kPdgNuMu )
308  {
309  label = dirname.substr ( dirname.find("nu_mu")+6, dirname.length() );
310  }
311  else
312  if ( gCurrProbePdg == kPdgAntiNuTau )
313  {
314  label = dirname.substr ( dirname.find("nu_tau_bar")+11, dirname.length() );
315  }
316  else
317  if ( gCurrProbePdg == kPdgNuTau )
318  {
319  label = dirname.substr ( dirname.find("nu_tau")+7, dirname.length() );
320  }
321 
322  bool free_n = (label.size()==1 && label.find("n") !=string::npos);
323  bool free_p = (label.size()==2 && label.find("H1")!=string::npos);
324  bool free_nuc = free_p || free_n;
325 
326  gCurrTargetHasP = (free_n) ? false : true;
327  gCurrTargetHasN = (free_p) ? false : true;
328 
329  gCurrTargetLbl = (free_nuc) ? "" : ("(" + label + ")");
330 
331  gCurrTargetIsFreeNuc = free_nuc;
332 }
333 //_________________________________________________________________________________
334 void MakePlots(void)
335 {
336  // Open files
337  gXSecFile_curr = new TFile(gOptXSecFilename_curr.c_str(), "read");
338  gXSecFile_ref0 = 0;
339  if (gOptHaveRef) {
340  gXSecFile_ref0 = new TFile(gOptXSecFilename_ref0.c_str(), "read");
341  }
342 
343  // Get list of directories in the input file
344  TList * directories = gXSecFile_curr->GetListOfKeys();
345 
346  // Loop over directories & generate plots for each one
347  int ndir = directories->GetEntries();
348  for(int idir=0; idir<ndir; idir++) {
349  TObjString * dir = (TObjString*) directories->At(idir);
350  gDirName = dir->GetString().Data();
351  MakePlotsCurrDir();
352  }
353 }
354 //_________________________________________________________________________________
355 void MakePlotsCurrDir(void)
356 {
357  LOG("gxscomp", pINFO) << "Plotting graphs from directory: " << gDirName;
358 
359  OpenDir ();
360  DirNameToProbe ();
361  DirNameToTarget ();
362 
363  const char * probestr = gCurrProbeLbl.c_str();
364  const char * targetstr = gCurrTargetLbl.c_str();
365 
366  LOG("gxscomp", pINFO) << "Probe : " << probestr;
367  LOG("gxscomp", pINFO) << "Target : " << targetstr;
368 
369  //
370  // Start plotting...
371  //
372 
373  if(!gCurrTargetIsFreeNuc) {
374  Draw("tot_cc", Form("%s + %s, TOT CC",probestr,targetstr));
375  Draw("tot_nc", Form("%s + %s, TOT NC",probestr,targetstr));
376  }
377 
378  if(gCurrTargetHasN) {
379  Draw("tot_cc_n", Form("%s + n %s, TOT CC" , probestr, targetstr));
380  Draw("tot_nc_n", Form("%s + n %s, TOT NC" , probestr, targetstr));
381  if(gCurrProbeIsNu) {
382  Draw("qel_cc_n", Form("%s + n %s, QEL CC" , probestr, targetstr));
383  }
384  Draw("qel_nc_n", Form("%s + n %s, NCEL" , probestr, targetstr));
385  Draw("res_cc_n", Form("%s + n %s, RES CC" , probestr, targetstr));
386  Draw("res_nc_n", Form("%s + n %s, RES NC" , probestr, targetstr));
387  Draw("res_cc_n_1232P33", Form("%s + n %s, RES CC, P33(1232)" , probestr, targetstr));
388  Draw("res_cc_n_1535S11", Form("%s + n %s, RES CC, S11(1535)" , probestr, targetstr));
389  Draw("res_cc_n_1520D13", Form("%s + n %s, RES CC, D13(1520)" , probestr, targetstr));
390  Draw("res_cc_n_1650S11", Form("%s + n %s, RES CC, S11(1650)" , probestr, targetstr));
391  Draw("res_cc_n_1700D13", Form("%s + n %s, RES CC, D13(1700)" , probestr, targetstr));
392  Draw("res_cc_n_1675D15", Form("%s + n %s, RES CC, D15(1675)" , probestr, targetstr));
393  Draw("res_cc_n_1620S31", Form("%s + n %s, RES CC, S31(1620)" , probestr, targetstr));
394  Draw("res_cc_n_1700D33", Form("%s + n %s, RES CC, D33(1700)" , probestr, targetstr));
395  Draw("res_cc_n_1440P11", Form("%s + n %s, RES CC, P11(1440)" , probestr, targetstr));
396  Draw("res_cc_n_1720P13", Form("%s + n %s, RES CC, P13(1720)" , probestr, targetstr));
397  Draw("res_cc_n_1680F15", Form("%s + n %s, RES CC, F15(1680)" , probestr, targetstr));
398  Draw("res_cc_n_1910P31", Form("%s + n %s, RES CC, P31(1910)" , probestr, targetstr));
399  Draw("res_cc_n_1920P33", Form("%s + n %s, RES CC, P33(1920)" , probestr, targetstr));
400  Draw("res_cc_n_1905F35", Form("%s + n %s, RES CC, F35(1905)" , probestr, targetstr));
401  Draw("res_cc_n_1950F37", Form("%s + n %s, RES CC, F37(1950)" , probestr, targetstr));
402  Draw("res_cc_n_1710P11", Form("%s + n %s, RES CC, P11(1710)" , probestr, targetstr));
403  Draw("res_nc_n_1232P33", Form("%s + n %s, RES NC, P33(1232)" , probestr, targetstr));
404  Draw("res_nc_n_1535S11", Form("%s + n %s, RES NC, S11(1535)" , probestr, targetstr));
405  Draw("res_nc_n_1520D13", Form("%s + n %s, RES NC, D13(1520)" , probestr, targetstr));
406  Draw("res_nc_n_1650S11", Form("%s + n %s, RES NC, S11(1650)" , probestr, targetstr));
407  Draw("res_nc_n_1700D13", Form("%s + n %s, RES NC, D13(1700)" , probestr, targetstr));
408  Draw("res_nc_n_1675D15", Form("%s + n %s, RES NC, D15(1675)" , probestr, targetstr));
409  Draw("res_nc_n_1620S31", Form("%s + n %s, RES NC, S31(1620)" , probestr, targetstr));
410  Draw("res_nc_n_1700D33", Form("%s + n %s, RES NC, D33(1700)" , probestr, targetstr));
411  Draw("res_nc_n_1440P11", Form("%s + n %s, RES NC, P11(1440)" , probestr, targetstr));
412  Draw("res_nc_n_1720P13", Form("%s + n %s, RES NC, P13(1720)" , probestr, targetstr));
413  Draw("res_nc_n_1680F15", Form("%s + n %s, RES NC, F15(1680)" , probestr, targetstr));
414  Draw("res_nc_n_1910P31", Form("%s + n %s, RES NC, P31(1910)" , probestr, targetstr));
415  Draw("res_nc_n_1920P33", Form("%s + n %s, RES NC, P33(1920)" , probestr, targetstr));
416  Draw("res_nc_n_1905F35", Form("%s + n %s, RES NC, F35(1905)" , probestr, targetstr));
417  Draw("res_nc_n_1950F37", Form("%s + n %s, RES NC, F37(1950)" , probestr, targetstr));
418  Draw("res_nc_n_1710P11", Form("%s + n %s, RES NC, P11(1710)" , probestr, targetstr));
419  Draw("dis_cc_n", Form("%s + n %s, DIS CC" , probestr, targetstr));
420  Draw("dis_nc_n", Form("%s + n %s, DIS NC" , probestr, targetstr));
421  if(gCurrProbeIsNu) {
422  Draw("dis_cc_n_ubarsea", Form("%s + n %s, DIS CC (#bar{u}_{sea})" , probestr, targetstr));
423  Draw("dis_cc_n_dval", Form("%s + n %s, DIS CC (d_{val})" , probestr, targetstr));
424  Draw("dis_cc_n_dsea", Form("%s + n %s, DIS CC (d_{sea})" , probestr, targetstr));
425  Draw("dis_cc_n_ssea", Form("%s + n %s, DIS CC (s_{sea})" , probestr, targetstr));
426  }
427  if(gCurrProbeIsNuBar) {
428  Draw("dis_cc_n_sbarsea", Form("%s + n %s, DIS CC (#bar{s}_{sea})" , probestr, targetstr));
429  Draw("dis_cc_n_dbarsea", Form("%s + n %s, DIS CC (#bar{d}_{val})" , probestr, targetstr));
430  Draw("dis_cc_n_uval", Form("%s + n %s, DIS CC (u_{val})" , probestr, targetstr));
431  Draw("dis_cc_n_usea", Form("%s + n %s, DIS CC (u_{sea})" , probestr, targetstr));
432  }
433  Draw("dis_nc_n_sbarsea", Form("%s + n %s, DIS NC (#bar{s}_{sea})" , probestr, targetstr));
434  Draw("dis_nc_n_ubarsea", Form("%s + n %s, DIS NC (#bar{u}_{sea})" , probestr, targetstr));
435  Draw("dis_nc_n_dbarsea", Form("%s + n %s, DIS NC (#bar{d}_{sea})" , probestr, targetstr));
436  Draw("dis_nc_n_dval", Form("%s + n %s, DIS NC (d_{val})" , probestr, targetstr));
437  Draw("dis_nc_n_dsea", Form("%s + n %s, DIS NC (d_{sea})" , probestr, targetstr));
438  Draw("dis_nc_n_uval", Form("%s + n %s, DIS NC (u_{val})" , probestr, targetstr));
439  Draw("dis_nc_n_usea", Form("%s + n %s, DIS NC (u_{sea})" , probestr, targetstr));
440  Draw("dis_nc_n_ssea", Form("%s + n %s, DIS NC (s_{sea})" , probestr, targetstr));
441  if(gCurrProbeIsNu) {
442  Draw("dis_cc_n_dval_charm", Form("%s + n %s, DIS CC (d_{val} -> c)" , probestr, targetstr));
443  Draw("dis_cc_n_dsea_charm", Form("%s + n %s, DIS CC (d_{sea} -> c)" , probestr, targetstr));
444  Draw("dis_cc_n_ssea_charm", Form("%s + n %s, DIS CC (s_{sea} -> c)" , probestr, targetstr));
445  }
446  if(gCurrProbeIsNuBar) {
447  Draw("dis_cc_n_dbarsea_charm", Form("%s + n %s, DIS CC (#bar{d}_{sea} -> #bar{c})" , probestr, targetstr));
448  Draw("dis_cc_n_sbarsea_charm", Form("%s + n %s, DIS CC (#bar{s}_{sea} -> #bar{c})" , probestr, targetstr));
449  }
450  }//N>0?
451 
452 
453  if(gCurrTargetHasP) {
454  Draw("tot_cc_p", Form("%s + p %s, TOT CC" , probestr, targetstr));
455  Draw("tot_nc_p", Form("%s + p %s, TOT NC" , probestr, targetstr));
456  if(gCurrProbeIsNuBar) {
457  Draw("qel_cc_p", Form("%s + p %s, QEL CC" , probestr, targetstr));
458  }
459  Draw("qel_nc_p", Form("%s + p %s, NCEL" , probestr, targetstr));
460  Draw("res_cc_p", Form("%s + p %s, RES CC" , probestr, targetstr));
461  Draw("res_nc_p", Form("%s + p %s, RES NC" , probestr, targetstr));
462  Draw("res_cc_p_1232P33", Form("%s + p %s, RES CC, P33(1232)" , probestr, targetstr));
463  Draw("res_cc_p_1535S11", Form("%s + p %s, RES CC, S11(1535)" , probestr, targetstr));
464  Draw("res_cc_p_1520D13", Form("%s + p %s, RES CC, D13(1520)" , probestr, targetstr));
465  Draw("res_cc_p_1650S11", Form("%s + p %s, RES CC, S11(1650)" , probestr, targetstr));
466  Draw("res_cc_p_1700D13", Form("%s + p %s, RES CC, D13(1700)" , probestr, targetstr));
467  Draw("res_cc_p_1675D15", Form("%s + p %s, RES CC, D15(1675)" , probestr, targetstr));
468  Draw("res_cc_p_1620S31", Form("%s + p %s, RES CC, S31(1620)" , probestr, targetstr));
469  Draw("res_cc_p_1700D33", Form("%s + p %s, RES CC, D33(1700)" , probestr, targetstr));
470  Draw("res_cc_p_1440P11", Form("%s + p %s, RES CC, P11(1440)" , probestr, targetstr));
471  Draw("res_cc_p_1720P13", Form("%s + p %s, RES CC, P13(1720)" , probestr, targetstr));
472  Draw("res_cc_p_1680F15", Form("%s + p %s, RES CC, F15(1680)" , probestr, targetstr));
473  Draw("res_cc_p_1910P31", Form("%s + p %s, RES CC, P31(1910)" , probestr, targetstr));
474  Draw("res_cc_p_1920P33", Form("%s + p %s, RES CC, P33(1920)" , probestr, targetstr));
475  Draw("res_cc_p_1905F35", Form("%s + p %s, RES CC, F35(1905)" , probestr, targetstr));
476  Draw("res_cc_p_1950F37", Form("%s + p %s, RES CC, F37(1950)" , probestr, targetstr));
477  Draw("res_cc_p_1710P11", Form("%s + p %s, RES CC, P11(1710)" , probestr, targetstr));
478  Draw("res_nc_p_1232P33", Form("%s + p %s, RES NC, P33(1232)" , probestr, targetstr));
479  Draw("res_nc_p_1535S11", Form("%s + p %s, RES NC, S11(1535)" , probestr, targetstr));
480  Draw("res_nc_p_1520D13", Form("%s + p %s, RES NC, D13(1520)" , probestr, targetstr));
481  Draw("res_nc_p_1650S11", Form("%s + p %s, RES NC, S11(1650)" , probestr, targetstr));
482  Draw("res_nc_p_1700D13", Form("%s + p %s, RES NC, D13(1700)" , probestr, targetstr));
483  Draw("res_nc_p_1675D15", Form("%s + p %s, RES NC, D15(1675)" , probestr, targetstr));
484  Draw("res_nc_p_1620S31", Form("%s + p %s, RES NC, S31(1620)" , probestr, targetstr));
485  Draw("res_nc_p_1700D33", Form("%s + p %s, RES NC, D33(1700)" , probestr, targetstr));
486  Draw("res_nc_p_1440P11", Form("%s + p %s, RES NC, P11(1440)" , probestr, targetstr));
487  Draw("res_nc_p_1720P13", Form("%s + p %s, RES NC, P13(1720)" , probestr, targetstr));
488  Draw("res_nc_p_1680F15", Form("%s + p %s, RES NC, F15(1680)" , probestr, targetstr));
489  Draw("res_nc_p_1910P31", Form("%s + p %s, RES NC, P31(1910)" , probestr, targetstr));
490  Draw("res_nc_p_1920P33", Form("%s + p %s, RES NC, P33(1920)" , probestr, targetstr));
491  Draw("res_nc_p_1905F35", Form("%s + p %s, RES NC, F35(1905)" , probestr, targetstr));
492  Draw("res_nc_p_1950F37", Form("%s + p %s, RES NC, F37(1950)" , probestr, targetstr));
493  Draw("res_nc_p_1710P11", Form("%s + p %s, RES NC, P11(1710)" , probestr, targetstr));
494  Draw("dis_cc_p", Form("%s + p %s, DIS CC" , probestr, targetstr));
495  Draw("dis_nc_p", Form("%s + p %s, DIS NC" , probestr, targetstr));
496  if(gCurrProbeIsNu) {
497  Draw("dis_cc_p_ubarsea", Form("%s + p %s, DIS CC (#bar{u}_{sea})" , probestr, targetstr));
498  Draw("dis_cc_p_dval", Form("%s + p %s, DIS CC (d_{val})" , probestr, targetstr));
499  Draw("dis_cc_p_dsea", Form("%s + p %s, DIS CC (d_{sea})" , probestr, targetstr));
500  Draw("dis_cc_p_ssea", Form("%s + p %s, DIS CC (s_{sea})" , probestr, targetstr));
501  }
502  if(gCurrProbeIsNuBar) {
503  Draw("dis_cc_n_sbarsea", Form("%s + n %s, DIS CC (#bar{s}_{sea})" , probestr, targetstr));
504  Draw("dis_cc_n_dbarsea", Form("%s + n %s, DIS CC (#bar{d}_{val})" , probestr, targetstr));
505  Draw("dis_cc_n_uval", Form("%s + n %s, DIS CC (u_{val})" , probestr, targetstr));
506  Draw("dis_cc_n_usea", Form("%s + n %s, DIS CC (u_{sea})" , probestr, targetstr));
507  }
508  Draw("dis_nc_p_sbarsea", Form("%s + p %s, DIS NC (#bar{s}_{sea})" , probestr, targetstr));
509  Draw("dis_nc_p_ubarsea", Form("%s + p %s, DIS NC (#bar{u}_{sea})" , probestr, targetstr));
510  Draw("dis_nc_p_dbarsea", Form("%s + p %s, DIS NC (#bar{d}_{sea})" , probestr, targetstr));
511  Draw("dis_nc_p_dval", Form("%s + p %s, DIS NC (d_{val})" , probestr, targetstr));
512  Draw("dis_nc_p_dsea", Form("%s + p %s, DIS NC (d_{sea})" , probestr, targetstr));
513  Draw("dis_nc_p_uval", Form("%s + p %s, DIS NC (u_{val})" , probestr, targetstr));
514  Draw("dis_nc_p_usea", Form("%s + p %s, DIS NC (u_{sea})" , probestr, targetstr));
515  Draw("dis_nc_p_ssea", Form("%s + p %s, DIS NC (s_{sea})" , probestr, targetstr));
516  if(gCurrProbeIsNu) {
517  Draw("dis_cc_p_dval_charm", Form("%s + p %s, DIS CC (d_{val} -> charm)" , probestr, targetstr));
518  Draw("dis_cc_p_dsea_charm", Form("%s + p %s, DIS CC (d_{sea} -> charm)" , probestr, targetstr));
519  Draw("dis_cc_p_ssea_charm", Form("%s + p %s, DIS CC (s_{sea} -> charm)" , probestr, targetstr));
520  }
521  if(gCurrProbeIsNuBar) {
522  Draw("dis_cc_p_dbarsea_charm", Form("%s + p %s, DIS CC (#bar{d}_{sea} -> #bar{charm})" , probestr, targetstr));
523  Draw("dis_cc_p_sbarsea_charm", Form("%s + p %s, DIS CC (#bar{s}_{sea} -> #bar{charm})" , probestr, targetstr));
524  }
525  }//Z>0
526 
527 }
528 //_________________________________________________________________________________
529 TH1F* DrawFrame(TGraph * gr0, TGraph * gr1, TPad * p, const char * xt, const char * yt, double yminsc, double ymaxsc)
530 {
531  double xmin = 1E-5;
532  double xmax = 1;
533  double ymin = 1E-5;
534  double ymax = 1;
535 
536  if(gr0) {
537  TAxis * x0 = gr0 -> GetXaxis();
538  TAxis * y0 = gr0 -> GetYaxis();
539  xmin = x0 -> GetXmin();
540  xmax = x0 -> GetXmax();
541  ymin = y0 -> GetXmin();
542  ymax = y0 -> GetXmax();
543  }
544  if(gr1) {
545  TAxis * x1 = gr1 -> GetXaxis();
546  TAxis * y1 = gr1 -> GetYaxis();
547  xmin = TMath::Min(xmin, x1 -> GetXmin());
548  xmax = TMath::Max(xmax, x1 -> GetXmax());
549  ymin = TMath::Min(ymin, y1 -> GetXmin());
550  ymax = TMath::Max(ymax, y1 -> GetXmax());
551  }
552  xmin *= 0.5;
553  xmax *= 1.5;
554  ymin *= yminsc;
555  ymax *= ymaxsc;
556  xmin = TMath::Max(0.1, xmin);
557 
558  return DrawFrame(xmin, xmax, ymin, ymax, p, xt, yt);
559 }
560 //_________________________________________________________________________________
561 TH1F* DrawFrame(double xmin, double xmax, double ymin, double ymax, TPad * p, const char * xt, const char * yt)
562 {
563  if(!p) return 0;
564  TH1F * hf = (TH1F*) p->DrawFrame(xmin, ymin, xmax, ymax);
565  hf->GetXaxis()->SetTitle(xt);
566  hf->GetYaxis()->SetTitle(yt);
567  hf->GetYaxis()->SetTitleSize(0.03);
568  hf->GetYaxis()->SetTitleOffset(1.5);
569  hf->GetXaxis()->SetLabelSize(0.03);
570  hf->GetYaxis()->SetLabelSize(0.03);
571  return hf;
572 }
573 //_________________________________________________________________________________
574 void Draw(TGraph* gr, const char * opt)
575 {
576  if(!gr) return;
577  gr->Draw(opt);
578 }
579 //_________________________________________________________________________________
580 void Draw(const char * plot, const char * title)
581 {
582  gPadTitle->cd();
583  TPavesText hdr(10,10,90,90,1,"tr");
584  hdr.AddText(title);
585  hdr.SetFillColor(kWhite);
586  hdr.Draw();
587 
588  gPadXSecs->cd();
589  TGraph * gr_curr = (TGraph*) gDirCurr->Get(plot);
590  TGraph * gr_ref0 = 0;
591  if(gDirRef0) {
592  gr_ref0 = (TGraph*) gDirRef0->Get(plot);
593  }
594 
595  if(gr_curr == 0 && gr_ref0 ==0) return;
596 
597  if(!title) return;
598 
599  // trim points in reference plot (shown with markers) so that the markers
600  // don't hide the current prediction (shown with a line).
601  // Keep, at most, 20 points per decade.
602  TGraph * gr_ref0_trim = TrimGraph(gr_ref0, 20);
603 
604  utils::style::Format(gr_curr, 1,1,1,1,1,1.0);
605  utils::style::Format(gr_ref0_trim, 1,1,1,2,4,0.7);
606 
607  DrawFrame (gr_curr, gr_ref0, gPadXSecs, "E_{#nu} (GeV)", "#sigma (10^{-38} cm^{2})", 0.5, 1.5);
608  Draw (gr_curr, "L");
609  Draw (gr_ref0_trim, "P");
610 
611  gLS->Clear();
612  gLS->AddEntry(gr_curr, gLabelCurr.c_str(), "L");
613  if(gOptHaveRef) {
614  gLS->AddEntry(gr_ref0_trim, gLabelRef0.c_str(), "P");
615  }
616  gLS->Draw();
617 
618  // plot ratio of current and reference models
619  if(gOptHaveRef)
620  {
621  gPadRatio->cd();
622  TGraph * gr_ratio = DrawRatio(gr_curr, gr_ref0);
623  DrawFrame (gr_ratio, 0, gPadRatio, "E_{#nu} (GeV)", Form("%s / %s", gLabelCurr.c_str(), gLabelRef0.c_str()), 0.9, 1.1);
624  Draw (gr_ratio, "L");
625  }
626 
627  gC ->Update();
628 
629  if(gr_ref0_trim) {
630  delete gr_ref0_trim;
631  }
632 
633  gPS->NewPage();
634 }
635 //_________________________________________________________________________________
636 TGraph * TrimGraph(TGraph * gr, int max_np_per_decade)
637 {
638  if(!gr) return 0;
639 
640  const int np = gr->GetN();
641  vector<bool> remv(np);
642 
643  int fp = 0;
644  int lp = 0;
645 
646  while(1) {
647  double xmin = gr->GetX()[fp];
648  double xmax = 10 * xmin;
649  int ndec = 0;
650  for(int i=fp; i<np; i++) {
651  remv[i] = false;
652  lp = i;
653  double x = gr->GetX()[i];
654  if(x > xmin && x <= xmax) ndec++;
655  if(x > xmax) break;
656  }
657  if(ndec > max_np_per_decade) {
658  double keep_prob = (double)max_np_per_decade/ (double) ndec;
659  int keep_rate = TMath::FloorNint(1./keep_prob);
660  int ithrow = 0;
661  for(int i=fp; i<=lp; i++) {
662  if(ithrow % keep_rate) {
663  remv[i] = true;
664  }
665  ithrow++;
666  }
667  }
668  if(lp == np-1) break;
669  fp = lp+1;
670  }
671 
672  int np2 = 0;
673  for(int i=0; i<np; i++) {
674  if(!remv[i]) {
675  np2++;
676  }
677  }
678 
679  double * x = new double[np2];
680  double * y = new double[np2];
681 
682  int j=0;
683  for(int i=0; i<np; i++) {
684  if(!remv[i]) {
685  x[j] = gr->GetX()[i];
686  y[j] = gr->GetY()[i];
687  j++;
688  }
689  }
690  TGraph * gr2 = new TGraph(np2,x,y);
691 
692  delete [] x;
693  delete [] y;
694 
695  return gr2;
696 }
697 //_________________________________________________________________________________
698 TGraph * DrawRatio(TGraph * gr0, TGraph * gr1)
699 {
700 // gr0 / gr1
701 
702  if(!gr0) return 0;
703  if(!gr1) return 0;
704 
705  LOG("gxscomp", pDEBUG) << "Drawing ratio...";
706 
707  const int np = gr0->GetN();
708 
709  double * x = new double[np];
710  double * y = new double[np];
711 
712  for(int i=0; i<np; i++) {
713  x[i] = gr0->GetX()[i];
714  y[i] = 0;
715 
716  if(gr0->Eval(x[i]) != 0. && gr1->Eval(x[i]) != 0.)
717  {
718  y[i] = gr0->Eval(x[i]) / gr1->Eval(x[i]);
719  }
720  else {
721  if(gr0->Eval(x[i]) != 0.)
722  {
723  y[i] = -1;
724  }
725  else
726  {
727  y[i] = 1;
728  }
729  }
730  }
731 
732  TGraph * ratio = new TGraph(np,x,y);
733 
734  delete [] x;
735  delete [] y;
736 
737  return ratio;
738 }
739 //_________________________________________________________________________________
740 void GetCommandLineArgs(int argc, char ** argv)
741 {
742  LOG("gxscomp", pINFO) << "*** Parsing command line arguments";
743 
744  CmdLnArgParser parser(argc,argv);
745 
746  // get input GENIE cross section file
747  if( parser.OptionExists('f') ) {
748  string inp = parser.ArgAsString('f');
749  if(inp.find(",") != string::npos) {
750  vector<string> inpv = utils::str::Split(inp,",");
751  assert(inpv.size()==2);
752  gOptXSecFilename_curr = inpv[0];
753  gLabelCurr = inpv[1];
754  } else {
755  gOptXSecFilename_curr = inp;
756  gLabelCurr = "current";
757  }
758  bool ok = CheckRootFilename(gOptXSecFilename_curr.c_str());
759  if(!ok) {
760  PrintSyntax();
761  exit(1);
762  }
763  } else {
764  PrintSyntax();
765  exit(1);
766  }
767 
768  // get [reference] input GENIE cross section file
769  if( parser.OptionExists('r') ) {
770  string inp = parser.ArgAsString('r');
771  if(inp.find(",") != string::npos) {
772  vector<string> inpv = utils::str::Split(inp,",");
773  assert(inpv.size()==2);
774  gOptXSecFilename_ref0 = inpv[0];
775  gLabelRef0 = inpv[1];
776  } else {
777  gOptXSecFilename_ref0 = inp;
778  gLabelRef0 = "reference";
779  }
780  bool ok = CheckRootFilename(gOptXSecFilename_ref0.c_str());
781  if(!ok) {
782  PrintSyntax();
783  exit(1);
784  }
785  gOptHaveRef = true;
786  } else {
787  LOG("gxscomp", pNOTICE) << "No reference cross section file";
788  gOptHaveRef = false;
789  }
790 
791  // get output filename
792  if( parser.OptionExists('o') ) {
793  gOptOutputFilename = parser.ArgAsString('o');
794  } else {
795  gOptOutputFilename = "xsec.ps";
796  }
797 
798 }
799 //_________________________________________________________________________________
800 void PrintSyntax(void)
801 {
802  LOG("gxscomp", pNOTICE)
803  << "\n\n" << "Syntax:" << "\n"
804  << " gxscomp -f xsec_file [-r reference_xsec_file] [-o output]\n";
805 }
806 //_________________________________________________________________________________
807 bool CheckRootFilename(string filename)
808 {
809  if(filename.size() == 0) return false;
810 
811  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));
812  if (!is_accessible) {
813  LOG("gxscomp", pERROR)
814  << "The input ROOT file [" << filename << "] is not accessible";
815  return false;
816  }
817  return true;
818 }
819 //_________________________________________________________________________________
820 string OutputFileName(string inpname)
821 {
822 // Builds the output filename based on the name of the input filename
823 // Perfors the following conversion: name.root -> name.nuxsec_test.ps
824 
825  unsigned int L = inpname.length();
826 
827  // if the last 4 characters are "root" (ROOT file extension) then
828  // remove them
829  if(inpname.substr(L-4, L).find("root") != string::npos) {
830  inpname.erase(L-4, L);
831  }
832 
833  ostringstream name;
834  name << inpname << "nuxsec_test.ps";
835 
836  return gSystem->BaseName(name.str().c_str());
837 }
838 //_________________________________________________________________________________
839 
gOptOutputFilename
string gOptOutputFilename
Definition: gXSecComp.cxx:98
MakePlotsCurrDir
void MakePlotsCurrDir(void)
Definition: gXSecComp.cxx:355
gC
TCanvas * gC
Definition: gXSecComp.cxx:110
MakePlots
void MakePlots(void)
Definition: gPDFComp.cxx:85
DirNameToTarget
void DirNameToTarget(void)
Definition: gXSecComp.cxx:284
genie::kPdgNuE
const int kPdgNuE
Definition: PDGCodes.h:28
gDirCurr
TDirectory * gDirCurr
Definition: gXSecComp.cxx:104
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
gLabelRef0
string gLabelRef0
Definition: gXSecComp.cxx:107
pERROR
#define pERROR
Definition: Messenger.h:59
gCurrProbePdg
int gCurrProbePdg
Definition: gXSecComp.cxx:116
OpenDir
void OpenDir(void)
Definition: gXSecComp.cxx:219
genie::kPdgAntiNuE
const int kPdgAntiNuE
Definition: PDGCodes.h:29
gPadTitle
TPad * gPadTitle
Definition: gXSecComp.cxx:111
gLabelCurr
string gLabelCurr
Definition: gXSecComp.cxx:106
gCurrProbeLbl
string gCurrProbeLbl
Definition: gXSecComp.cxx:115
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:30
TrimGraph
TGraph * TrimGraph(TGraph *gr, int max_np_per_decade)
Definition: gXSecComp.cxx:636
DirNameToProbe
void DirNameToProbe(void)
Definition: gXSecComp.cxx:232
gCurrTargetHasN
bool gCurrTargetHasN
Definition: gXSecComp.cxx:121
gXSecFile_curr
TFile * gXSecFile_curr
Definition: gXSecComp.cxx:102
dir
string dir
Definition: INukeNucleonCorr.cxx:61
gCurrTargetLbl
string gCurrTargetLbl
Definition: gXSecComp.cxx:119
genie::utils::style::SetDefaultStyle
void SetDefaultStyle(bool black_n_white=false)
Definition: Style.cxx:20
main
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
gLS
TLegend * gLS
Definition: gXSecComp.cxx:114
DrawFrame
TH1F * DrawFrame(TGraph *gr0, TGraph *gr1, TPad *p, const char *xt, const char *yt, double yminsc, double ymaxsc)
Definition: gXSecComp.cxx:529
gCurrProbeIsNuBar
bool gCurrProbeIsNuBar
Definition: gXSecComp.cxx:118
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
gPadXSecs
TPad * gPadXSecs
Definition: gXSecComp.cxx:112
gDirRef0
TDirectory * gDirRef0
Definition: gXSecComp.cxx:105
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
gOptXSecFilename_ref0
string gOptXSecFilename_ref0
Definition: gXSecComp.cxx:97
gOptHaveRef
bool gOptHaveRef
Definition: gXSecComp.cxx:99
pINFO
#define pINFO
Definition: Messenger.h:62
CheckRootFilename
bool CheckRootFilename(string filename)
Definition: gEvComp.cxx:2033
Init
void Init(void)
Definition: gXSecComp.cxx:138
genie::kPdgAntiNuTau
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
genie::kPdgAntiNuMu
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
Draw
void Draw(const char *plot, const char *title)
Definition: gXSecComp.cxx:580
gCurrProbeIsNu
bool gCurrProbeIsNu
Definition: gXSecComp.cxx:117
genie::kPdgNuTau
const int kPdgNuTau
Definition: PDGCodes.h:32
genie::utils::str::Split
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
End
void End(void)
Definition: gXSecComp.cxx:210
gCurrTargetHasP
bool gCurrTargetHasP
Definition: gXSecComp.cxx:120
gCurrTargetIsFreeNuc
bool gCurrTargetIsFreeNuc
Definition: gXSecComp.cxx:122
OutputFileName
string OutputFileName(string input_file_name)
Definition: gEvComp.cxx:1982
DrawRatio
TGraph * DrawRatio(TGraph *gr0, TGraph *gr1)
Definition: gXSecComp.cxx:698
gXSecFile_ref0
TFile * gXSecFile_ref0
Definition: gXSecComp.cxx:103
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:146
genie::CmdLnArgParser
Command line argument parser.
Definition: CmdLnArgParser.h:29
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
gOptXSecFilename_curr
string gOptXSecFilename_curr
Definition: gXSecComp.cxx:96
gPS
TPostScript * gPS
Definition: gXSecComp.cxx:109
gPadRatio
TPad * gPadRatio
Definition: gXSecComp.cxx:113
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
PrintSyntax
void PrintSyntax(void)
Definition: gAtmoEvGen.cxx:1051
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
gDirName
string gDirName
Definition: gXSecComp.cxx:108
Generated on Wed Jul 17 2024 14:46:41 for GENIEGenerator by   doxygen 1.8.5