GENIEGenerator
Loading...
Searching...
No Matches
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-2025, 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
70using std::ostringstream;
71using std::string;
72using std::vector;
73
74using namespace genie;
75
76// function prototypes
77void Init (void);
78void End (void);
79void OpenDir (void);
80void DirNameToProbe (void);
81void DirNameToTarget (void);
82void MakePlots (void);
83void MakePlotsCurrDir (void);
84TH1F * DrawFrame (TGraph * gr0, TGraph * gr1, TPad * p, const char * xt, const char * yt, double yminsc, double ymaxsc);
85TH1F * DrawFrame (double xmin, double xmax, double ymin, double ymax, TPad * p, const char * xt, const char * yt);
86void Draw (const char * plot, const char * title);
87void Draw (TGraph* gr, const char * opt);
88TGraph* TrimGraph (TGraph * gr, int max_np_per_decade);
89TGraph* DrawRatio (TGraph * gr0, TGraph * gr1);
90void GetCommandLineArgs (int argc, char ** argv);
91void PrintSyntax (void);
92bool CheckRootFilename (string filename);
93string OutputFileName (string input_file_name);
94
95// command-line arguments
96string gOptXSecFilename_curr = ""; // (-f) input ROOT cross section file
97string gOptXSecFilename_ref0 = ""; // (-r) input ROOT cross section file (ref.)
98string gOptOutputFilename = ""; // (-o) output PS file
99bool gOptHaveRef;
100
101// globals
102TFile * gXSecFile_curr = 0;
103TFile * gXSecFile_ref0 = 0;
104TDirectory * gDirCurr = 0;
105TDirectory * gDirRef0 = 0;
106string gLabelCurr = "";
107string gLabelRef0 = "";
108string gDirName = "";
109TPostScript * gPS = 0;
110TCanvas * gC = 0;
111TPad * gPadTitle = 0;
112TPad * gPadXSecs = 0;
113TPad * gPadRatio = 0;
114TLegend * gLS = 0;
115string gCurrProbeLbl = "";
116int gCurrProbePdg = 0;
117bool gCurrProbeIsNu = false;
118bool gCurrProbeIsNuBar = false;
119string gCurrTargetLbl = "";
120bool gCurrTargetHasP = false;
121bool gCurrTargetHasN = false;
122bool gCurrTargetIsFreeNuc = false;
123
124//_________________________________________________________________________________
125int 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//_________________________________________________________________________________
138void 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//_________________________________________________________________________________
210void End(void)
211{
212 gPS->Close();
213
214 delete gC;
215 delete gPS;
216 delete gLS;
217}
218//_________________________________________________________________________________
219void 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//_________________________________________________________________________________
232void 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//_________________________________________________________________________________
284void 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//_________________________________________________________________________________
334void 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//_________________________________________________________________________________
355void 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//_________________________________________________________________________________
529TH1F* 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//_________________________________________________________________________________
561TH1F* 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//_________________________________________________________________________________
574void Draw(TGraph* gr, const char * opt)
575{
576 if(!gr) return;
577 gr->Draw(opt);
578}
579//_________________________________________________________________________________
580void 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//_________________________________________________________________________________
636TGraph * 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//_________________________________________________________________________________
698TGraph * 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//_________________________________________________________________________________
740void 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//_________________________________________________________________________________
800void 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//_________________________________________________________________________________
807bool 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//_________________________________________________________________________________
820string 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
dir
string dir
Definition INukeNucleonCorr.cxx:61
pNOTICE
#define pNOTICE
Definition Messenger.h:61
pINFO
#define pINFO
Definition Messenger.h:62
pERROR
#define pERROR
Definition Messenger.h:59
pDEBUG
#define pDEBUG
Definition Messenger.h:63
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils.
main
int main()
Definition Validation.cpp:290
genie::CmdLnArgParser
Command line argument parser.
Definition CmdLnArgParser.h:29
genie::CmdLnArgParser::ArgAsString
string ArgAsString(char opt)
Definition CmdLnArgParser.cxx:106
genie::CmdLnArgParser::OptionExists
bool OptionExists(char opt)
was option set?
Definition CmdLnArgParser.cxx:88
DrawFrame
TH1F * DrawFrame(TGraph *gr0, TGraph *gr1, TPad *p, const char *xt, const char *yt, double yminsc, double ymaxsc)
Definition gXSecComp.cxx:529
OutputFileName
string OutputFileName(string input_file_name)
Definition gXSecComp.cxx:820
gC
TCanvas * gC
Definition gXSecComp.cxx:110
gOptXSecFilename_ref0
string gOptXSecFilename_ref0
Definition gXSecComp.cxx:97
DrawRatio
TGraph * DrawRatio(TGraph *gr0, TGraph *gr1)
Definition gXSecComp.cxx:698
Draw
void Draw(const char *plot, const char *title)
Definition gXSecComp.cxx:580
gXSecFile_curr
TFile * gXSecFile_curr
Definition gXSecComp.cxx:102
gLabelRef0
string gLabelRef0
Definition gXSecComp.cxx:107
gPS
TPostScript * gPS
Definition gXSecComp.cxx:109
MakePlotsCurrDir
void MakePlotsCurrDir(void)
Definition gXSecComp.cxx:355
DirNameToProbe
void DirNameToProbe(void)
Definition gXSecComp.cxx:232
gDirRef0
TDirectory * gDirRef0
Definition gXSecComp.cxx:105
gPadRatio
TPad * gPadRatio
Definition gXSecComp.cxx:113
gCurrProbeLbl
string gCurrProbeLbl
Definition gXSecComp.cxx:115
End
void End(void)
Definition gXSecComp.cxx:210
gCurrProbeIsNuBar
bool gCurrProbeIsNuBar
Definition gXSecComp.cxx:118
gPadTitle
TPad * gPadTitle
Definition gXSecComp.cxx:111
MakePlots
void MakePlots(void)
Definition gXSecComp.cxx:334
gXSecFile_ref0
TFile * gXSecFile_ref0
Definition gXSecComp.cxx:103
CheckRootFilename
bool CheckRootFilename(string filename)
Definition gXSecComp.cxx:807
gOptXSecFilename_curr
string gOptXSecFilename_curr
Definition gXSecComp.cxx:96
gCurrTargetIsFreeNuc
bool gCurrTargetIsFreeNuc
Definition gXSecComp.cxx:122
gCurrTargetLbl
string gCurrTargetLbl
Definition gXSecComp.cxx:119
gCurrTargetHasP
bool gCurrTargetHasP
Definition gXSecComp.cxx:120
gLabelCurr
string gLabelCurr
Definition gXSecComp.cxx:106
gCurrTargetHasN
bool gCurrTargetHasN
Definition gXSecComp.cxx:121
Init
void Init(void)
Definition gXSecComp.cxx:138
gLS
TLegend * gLS
Definition gXSecComp.cxx:114
gCurrProbeIsNu
bool gCurrProbeIsNu
Definition gXSecComp.cxx:117
gCurrProbePdg
int gCurrProbePdg
Definition gXSecComp.cxx:116
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition gXSecComp.cxx:740
PrintSyntax
void PrintSyntax(void)
Definition gXSecComp.cxx:800
gPadXSecs
TPad * gPadXSecs
Definition gXSecComp.cxx:112
gDirName
string gDirName
Definition gXSecComp.cxx:108
TrimGraph
TGraph * TrimGraph(TGraph *gr, int max_np_per_decade)
Definition gXSecComp.cxx:636
gDirCurr
TDirectory * gDirCurr
Definition gXSecComp.cxx:104
OpenDir
void OpenDir(void)
Definition gXSecComp.cxx:219
gOptOutputFilename
string gOptOutputFilename
Definition gXSecComp.cxx:98
DirNameToTarget
void DirNameToTarget(void)
Definition gXSecComp.cxx:284
gOptHaveRef
bool gOptHaveRef
Definition gXSecComp.cxx:99
genie::pdg
Utilities for improving the code readability when using PDG codes.
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition PDGUtils.cxx:110
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition PDGUtils.cxx:118
genie::utils::str::Split
vector< string > Split(string input, string delim)
Definition StringUtils.cxx:36
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::utils::style::SetDefaultStyle
void SetDefaultStyle(bool black_n_white=false)
Definition Style.cxx:20
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
genie::kPdgAntiNuE
const int kPdgAntiNuE
Definition PDGCodes.h:29
genie::kPdgAntiNuTau
const int kPdgAntiNuTau
Definition PDGCodes.h:33
genie::kPdgNuE
const int kPdgNuE
Definition PDGCodes.h:28
genie::kPdgNuTau
const int kPdgNuTau
Definition PDGCodes.h:32
genie::kPdgAntiNuMu
const int kPdgAntiNuMu
Definition PDGCodes.h:31
genie::kPdgNuMu
const int kPdgNuMu
Definition PDGCodes.h:30
Generated on for GENIEGenerator by doxygen 1.14.0