gPDFComp.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gpdfcomp
5 
6 \brief PDF comparison tool
7 
8 \syntax gpdfcomp --pdf-set pdf_set [-o output]
9 
10  --pdf-set :
11  Specifies a comma separated list of GENIE PDFs.
12  The full algorithm name and configuration should be provided for
13  each GENIE PDF as in `genie::model_name/model_config'.
14  (unique color and line styles are defined for up to 4 sets)
15 
16  -o :
17  Specifies a name to be used in the output files.
18  Default: pdf_comp
19 
20 \example gpdfcomp --pdf-set genie::GRV98LO/Default,genie::BYPDF/Default
21 
22 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
23  University of Liverpool & STFC Rutherford Appleton Lab
24 
25 \created Feb 10, 2016
26 
27 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
28  For the full text of the license visit http://copyright.genie-mc.org
29  or see $GENIE/LICENSE
30 */
31 //____________________________________________________________________________
32 
33 #include <vector>
34 #include <string>
35 
36 #include <TNtuple.h>
37 #include <TFile.h>
38 #include <TMath.h>
39 #include <TPostScript.h>
40 #include <TPavesText.h>
41 #include <TText.h>
42 #include <TStyle.h>
43 #include <TLegend.h>
44 #include <TCanvas.h>
45 #include <TGraph.h>
46 #include <TLatex.h>
47 #include <TPaletteAxis.h>
48 
52 //#include "Physics/PartonDistributions/LHAPDF5.h"
56 #include "Framework/Utils/Style.h"
57 
58 using namespace std;
59 using namespace genie;
60 using namespace genie::utils;
61 
62 // globals
63 string gOptPDFSet = ""; // --pdf-set argument
64 string gOptOutFile = "pdf_comp"; // -o argument
65 vector<const PDFModelI *> gPDFAlgList;
66 
67 // function prototypes
68 void GetCommandLineArgs (int argc, char ** argv);
69 void GetAlgorithms (void);
70 void MakePlots (void);
71 
72 //___________________________________________________________________
73 int main(int argc, char ** argv)
74 {
76 
77  GetCommandLineArgs (argc,argv); // Get command line arguments
78  GetAlgorithms(); // Get requested PDF algorithms
79  MakePlots(); // Produce all output plots and fill output n-tuple
80 
81  return 0;
82 }
83 //_________________________________________________________________________________
84 void MakePlots (void)
85 {
86  const unsigned int nm = 4; // number of models
87  int col [nm] = { kBlack, kRed+1, kBlue-3, kGreen+2 };
88  int sty [nm] = { kSolid, kDashed, kDashed, kDashed };
89  int mrk [nm] = { 20, 20, 20, 20 };
90  double msz [nm] = { 0.7, 0.7, 0.7, 0.7 };
91  const char * opt [nm] = { "ap", "l", "l", "l" };
92  const char * lgopt [nm] = { "P", "L", "L", "L" };
93 
94  // Q2 range and values for 1-D plots
95  const unsigned int nQ2 = 20;
96  const double Q2min = 1E-1; // GeV^2
97  const double Q2max = 1E+3; // GeV^2
98  const double log10Q2min = TMath::Log10(Q2min);
99  const double log10Q2max = TMath::Log10(Q2max);
100  const double dlog10Q2 = (log10Q2max-log10Q2min)/(nQ2-1);
101  double Q2_arr [nQ2];
102  for(unsigned int iq2 = 0; iq2 < nQ2; iq2++) {
103  double Q2 = TMath::Power(10, log10Q2min + iq2*dlog10Q2);
104  Q2_arr[iq2] = Q2;
105  }
106 
107  // x values for 1-D plots
108  const unsigned int nx = 22;
109  double x_arr [nx] = {
110  0.0001, 0.0010, 0.0100, 0.0250, 0.0500,
111  0.0750, 0.1000, 0.1500, 0.2000, 0.2500,
112  0.3500, 0.4000, 0.4500, 0.5000, 0.5500,
113  0.6000, 0.7000, 0.7500, 0.8000, 0.8500,
114  0.9000, 0.9500
115  };
116 
117  // Q2 range and values for 2-D plots
118  const unsigned int nQ2_2d = 50;
119  const double Q2min_2d = 1E-1; // GeV^2
120  const double Q2max_2d = 1E+3; // GeV^2
121  const double log10Q2min_2d = TMath::Log10(Q2min_2d);
122  const double log10Q2max_2d = TMath::Log10(Q2max_2d);
123  const double dlog10Q2_2d = (log10Q2max_2d-log10Q2min_2d)/(nQ2_2d-1);
124  double Q2_bin_edges_2d [nQ2_2d];
125  for(unsigned int iq2 = 0; iq2 < nQ2_2d; iq2++) {
126  double Q2 = TMath::Power(10, log10Q2min_2d + iq2*dlog10Q2_2d);
127  Q2_bin_edges_2d[iq2] = Q2;
128  }
129 
130  // x range and values for 2-D plots
131  const unsigned int nx_2d = 50;
132  const double xmin_2d = 1E-4;
133  const double xmax_2d = 0.95;
134  const double log10xmin_2d = TMath::Log10(xmin_2d);
135  const double log10xmax_2d = TMath::Log10(xmax_2d);
136  const double dlog10x_2d = (log10xmax_2d-log10xmin_2d)/(nx_2d-1);
137  double x_bin_edges_2d [nx_2d];
138  for(unsigned int ix = 0; ix < nx_2d; ix++) {
139  double x = TMath::Power(10, log10xmin_2d + ix*dlog10x_2d);
140  x_bin_edges_2d[ix] = x;
141  }
142 
143 
144  // Output ntuple
145  TNtuple * ntpl = new TNtuple("nt","pdfs","i:uv:dv:us:ds:s:g:x:Q2");
146 
147  // Canvas for output plots
148  TCanvas * cnv = new TCanvas("c","",20,20,500,650);
149  cnv->SetBorderMode(0);
150  cnv->SetFillColor(0);
151 
152  // Create output ps file
153  string ps_filename = gOptOutFile + ".ps";
154  TPostScript * ps = new TPostScript(ps_filename.c_str(), 111);
155 
156  // Front page
157  ps->NewPage();
158  cnv->Range(0,0,100,100);
159  TPavesText hdr(10,40,90,70,3,"tr");
160  hdr.AddText("GENIE parton density function (pdf) comparisons");
161  hdr.AddText(" ");
162  hdr.AddText(" ");
163  hdr.AddText(" ");
164  hdr.AddText(" ");
165  hdr.AddText("Models used:");
166  hdr.AddText(" ");
167  for(unsigned int im=0; im < gPDFAlgList.size(); im++) {
168  const char * label = gPDFAlgList[im]->Id().Key().c_str();
169  hdr.AddText(label);
170  }
171  hdr.AddText(" ");
172  hdr.Draw();
173  cnv->Update();
174 
175  ps->NewPage();
176 
177  TLatex * tex = new TLatex;
178  TLegend * lgnd = 0;
179 
180  //
181  // Plots
182  //
183 
184  // PDFs = f(Q^2) for selected vales of x (1-D plots)
185  TGraph * gr_xuv_Q2 [nm] = { 0 };
186  TGraph * gr_xdv_Q2 [nm] = { 0 };
187  TGraph * gr_xus_Q2 [nm] = { 0 };
188  TGraph * gr_xds_Q2 [nm] = { 0 };
189  TGraph * gr_xstr_Q2 [nm] = { 0 };
190  TGraph * gr_xglu_Q2 [nm] = { 0 };
191 
192  // PDFs = f(x,Q^2) with fine x,Q2 binning (2-D plots)
193  TH2D * h2_xuv [nm] = { 0 };
194  TH2D * h2_xdv [nm] = { 0 };
195  TH2D * h2_xus [nm] = { 0 };
196  TH2D * h2_xds [nm] = { 0 };
197  TH2D * h2_xstr[nm] = { 0 };
198  TH2D * h2_xglu[nm] = { 0 };
199 
200  // PDFs ratios relative to first one specified = f(x,Q^2) with fine x,Q2 binning (2-D plots)
201  TH2D * h2_xuv_r [nm] = { 0 };
202  TH2D * h2_xdv_r [nm] = { 0 };
203  TH2D * h2_xus_r [nm] = { 0 };
204  TH2D * h2_xds_r [nm] = { 0 };
205  TH2D * h2_xstr_r[nm] = { 0 };
206  TH2D * h2_xglu_r[nm] = { 0 };
207 
208  //
209  // Generate PDF plots as f(Q^2) for selected vales of x (1-D plots)
210  //
211 
212  for(unsigned int ix=0; ix < nx; ix++) {
213  double x = x_arr[ix];
214  double xuv_arr [nm][nQ2];
215  double xdv_arr [nm][nQ2];
216  double xus_arr [nm][nQ2];
217  double xds_arr [nm][nQ2];
218  double xstr_arr [nm][nQ2];
219  double xglu_arr [nm][nQ2];
220 
221  double max_gr_xuv_Q2 = -9E9;
222  double max_gr_xdv_Q2 = -9E9;
223  double max_gr_xus_Q2 = -9E9;
224  double max_gr_xds_Q2 = -9E9;
225  double max_gr_xstr_Q2 = -9E9;
226  double max_gr_xglu_Q2 = -9E9;
227 
228  for(unsigned int im=0; im < gPDFAlgList.size(); im++) {
229  PDF pdf;
230  pdf.SetModel(gPDFAlgList[im]);
231  for(unsigned int iq2 = 0; iq2 < nQ2; iq2++) {
232  double Q2 = Q2_arr[iq2];
233  pdf.Calculate(x, Q2);
234  //LOG("gpdfcomp", pINFO) << "PDFs:\n" << pdf;
235  double xuv = pdf.UpValence();
236  double xdv = pdf.DownValence();
237  double xus = pdf.UpSea();
238  double xds = pdf.DownSea();
239  double xstr = pdf.Strange();
240  double xglu = pdf.Gluon();
241  xuv_arr [im][iq2] = x * xuv;
242  xdv_arr [im][iq2] = x * xdv;
243  xus_arr [im][iq2] = x * xus;
244  xds_arr [im][iq2] = x * xds;
245  xstr_arr [im][iq2] = x * xstr;
246  xglu_arr [im][iq2] = x * xglu;
247  ntpl->Fill(im,xuv,xdv,xus,xds,xstr,xglu,x,Q2);
248  }//iq2
249 
250  gr_xuv_Q2 [im] = new TGraph (nQ2, Q2_arr, xuv_arr [im]);
251  gr_xdv_Q2 [im] = new TGraph (nQ2, Q2_arr, xdv_arr [im]);
252  gr_xus_Q2 [im] = new TGraph (nQ2, Q2_arr, xus_arr [im]);
253  gr_xds_Q2 [im] = new TGraph (nQ2, Q2_arr, xds_arr [im]);
254  gr_xstr_Q2 [im] = new TGraph (nQ2, Q2_arr, xstr_arr [im]);
255  gr_xglu_Q2 [im] = new TGraph (nQ2, Q2_arr, xglu_arr [im]);
256 
257  genie::utils::style::Format( gr_xuv_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
258  genie::utils::style::Format( gr_xdv_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
259  genie::utils::style::Format( gr_xus_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
260  genie::utils::style::Format( gr_xds_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
261  genie::utils::style::Format( gr_xstr_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
262  genie::utils::style::Format( gr_xglu_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
263 
264  gr_xuv_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
265  gr_xdv_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
266  gr_xus_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
267  gr_xds_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
268  gr_xstr_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
269  gr_xglu_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
270 
271  gr_xuv_Q2 [im] -> GetYaxis() -> SetTitle("x*u_{val}(x,Q^{2})");
272  gr_xdv_Q2 [im] -> GetYaxis() -> SetTitle("x*d_{val}(x,Q^{2})");
273  gr_xus_Q2 [im] -> GetYaxis() -> SetTitle("x*u_{sea}(x,Q^{2})");
274  gr_xds_Q2 [im] -> GetYaxis() -> SetTitle("x*d_{sea}(x,Q^{2})");
275  gr_xstr_Q2 [im] -> GetYaxis() -> SetTitle("x*s(x,Q^{2})");
276  gr_xglu_Q2 [im] -> GetYaxis() -> SetTitle("x*g(x,Q^{2})");
277 
278  double this_max_gr_xuv_Q2 = TMath::MaxElement(gr_xuv_Q2 [im]->GetN(),gr_xuv_Q2 [im]->GetY());
279  double this_max_gr_xdv_Q2 = TMath::MaxElement(gr_xdv_Q2 [im]->GetN(),gr_xdv_Q2 [im]->GetY());
280  double this_max_gr_xus_Q2 = TMath::MaxElement(gr_xus_Q2 [im]->GetN(),gr_xus_Q2 [im]->GetY());
281  double this_max_gr_xds_Q2 = TMath::MaxElement(gr_xds_Q2 [im]->GetN(),gr_xds_Q2 [im]->GetY());
282  double this_max_gr_xstr_Q2 = TMath::MaxElement(gr_xstr_Q2[im]->GetN(),gr_xstr_Q2[im]->GetY());
283  double this_max_gr_xglu_Q2 = TMath::MaxElement(gr_xglu_Q2[im]->GetN(),gr_xglu_Q2[im]->GetY());
284  max_gr_xuv_Q2 = std::max(max_gr_xuv_Q2 ,this_max_gr_xuv_Q2 );
285  max_gr_xdv_Q2 = std::max(max_gr_xdv_Q2 ,this_max_gr_xdv_Q2 );
286  max_gr_xus_Q2 = std::max(max_gr_xus_Q2 ,this_max_gr_xus_Q2 );
287  max_gr_xds_Q2 = std::max(max_gr_xds_Q2 ,this_max_gr_xds_Q2 );
288  max_gr_xstr_Q2 = std::max(max_gr_xstr_Q2,this_max_gr_xstr_Q2);
289  max_gr_xglu_Q2 = std::max(max_gr_xglu_Q2,this_max_gr_xglu_Q2);
290 
291  }//im
292 
293  // Now loop to set sensible limits
294  for(unsigned int im=0; im < gPDFAlgList.size(); im++) {
295  gr_xuv_Q2 [im] -> SetMinimum(0.);
296  gr_xdv_Q2 [im] -> SetMinimum(0.);
297  gr_xus_Q2 [im] -> SetMinimum(0.);
298  gr_xds_Q2 [im] -> SetMinimum(0.);
299  gr_xstr_Q2 [im] -> SetMinimum(0.);
300  gr_xglu_Q2 [im] -> SetMinimum(0.);
301  gr_xuv_Q2 [im] -> SetMaximum(1.1*max_gr_xuv_Q2 );
302  gr_xdv_Q2 [im] -> SetMaximum(1.1*max_gr_xdv_Q2 );
303  gr_xus_Q2 [im] -> SetMaximum(1.1*max_gr_xus_Q2 );
304  gr_xds_Q2 [im] -> SetMaximum(1.1*max_gr_xds_Q2 );
305  gr_xstr_Q2 [im] -> SetMaximum(1.1*max_gr_xstr_Q2 );
306  gr_xglu_Q2 [im] -> SetMaximum(1.1*max_gr_xglu_Q2 );
307  }//im
308 
309  ps->NewPage();
310 
311  if(x<0.3) {
312  lgnd = new TLegend(0.20, 0.55, 0.50, 0.85);
313  } else {
314  lgnd = new TLegend(0.60, 0.55, 0.80, 0.85);
315  }
316  lgnd -> SetFillColor(0);
317  lgnd -> SetFillStyle(0);
318  lgnd -> SetBorderSize(0);
319 
320  lgnd->Clear();
321  for(unsigned int im=0; im < gPDFAlgList.size(); im++) {
322  std::string label(gPDFAlgList[im]->Id().Key());
323  lgnd->AddEntry(gr_xuv_Q2 [im], label.c_str(), lgopt[im]);
324  }
325 
326  cnv->Clear();
327  cnv->Divide(2,3);
328 
329  cnv->cd(1); gPad->SetLogx();
330  cnv->cd(2); gPad->SetLogx();
331  cnv->cd(3); gPad->SetLogx();
332  cnv->cd(4); gPad->SetLogx();
333  cnv->cd(5); gPad->SetLogx();
334  cnv->cd(6); gPad->SetLogx();
335 
336  for(unsigned int im=0; im < gPDFAlgList.size(); im++) {
337  cnv->cd(1); gr_xuv_Q2 [im]->Draw(opt[im]);
338  cnv->cd(2); gr_xdv_Q2 [im]->Draw(opt[im]);
339  cnv->cd(3); gr_xus_Q2 [im]->Draw(opt[im]);
340  cnv->cd(4); gr_xds_Q2 [im]->Draw(opt[im]);
341  cnv->cd(5); gr_xstr_Q2[im]->Draw(opt[im]);
342  cnv->cd(6); gr_xglu_Q2[im]->Draw(opt[im]);
343  }
344 
345  cnv->cd(1); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, TString::Format("x = %.3e",x).Data());
346  cnv->cd(2); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, TString::Format("x = %.3e",x).Data());
347  cnv->cd(3); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, TString::Format("x = %.3e",x).Data());
348  cnv->cd(4); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, TString::Format("x = %.3e",x).Data());
349  cnv->cd(5); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, TString::Format("x = %.3e",x).Data());
350  cnv->cd(6); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, TString::Format("x = %.3e",x).Data());
351 
352  cnv->Update();
353  }
354 
355  //
356  // Plot PDFs = f(x,Q2)
357  //
358 
359  for(unsigned int im=0; im < gPDFAlgList.size(); im++) {
360  h2_xuv [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
361  h2_xdv [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
362  h2_xus [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
363  h2_xds [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
364  h2_xstr[im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
365  h2_xglu[im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
366  PDF pdf;
367  pdf.SetModel(gPDFAlgList[im]);
368  for(int ibinx = 1;
369  ibinx <= h2_xuv[im]->GetXaxis()->GetNbins(); ibinx++) {
370  double x = h2_xuv[im]->GetXaxis()->GetBinCenter(ibinx);
371  for(int ibinq2 = 1;
372  ibinq2 <= h2_xuv[im]->GetYaxis()->GetNbins(); ibinq2++) {
373  double Q2 = h2_xuv[im]->GetYaxis()->GetBinCenter(ibinq2);
374  pdf.Calculate(x, Q2);
375  //LOG("gpdfcomp", pINFO) << "PDFs:\n" << pdf;
376  double xuv = x * pdf.UpValence();
377  double xdv = x * pdf.DownValence();
378  double xus = x * pdf.UpSea();
379  double xds = x * pdf.DownSea();
380  double xstr = x * pdf.Strange();
381  double xglu = x * pdf.Gluon();
382  h2_xuv [im] -> SetBinContent(ibinx, ibinq2, xuv );
383  h2_xdv [im] -> SetBinContent(ibinx, ibinq2, xdv );
384  h2_xus [im] -> SetBinContent(ibinx, ibinq2, xus );
385  h2_xds [im] -> SetBinContent(ibinx, ibinq2, xds );
386  h2_xstr[im] -> SetBinContent(ibinx, ibinq2, xstr);
387  h2_xglu[im] -> SetBinContent(ibinx, ibinq2, xglu);
388  }
389  }
390 
391  ps->NewPage();
392 
393  h2_xuv [im] -> GetXaxis() -> SetTitle("x");
394  h2_xdv [im] -> GetXaxis() -> SetTitle("x");
395  h2_xus [im] -> GetXaxis() -> SetTitle("x");
396  h2_xds [im] -> GetXaxis() -> SetTitle("x");
397  h2_xstr[im] -> GetXaxis() -> SetTitle("x");
398  h2_xglu[im] -> GetXaxis() -> SetTitle("x");
399 
400  h2_xuv [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
401  h2_xdv [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
402  h2_xus [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
403  h2_xds [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
404  h2_xstr[im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
405  h2_xglu[im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
406 
407  h2_xuv [im] -> SetMinimum(0.);
408  h2_xdv [im] -> SetMinimum(0.);
409  h2_xus [im] -> SetMinimum(0.);
410  h2_xds [im] -> SetMinimum(0.);
411  h2_xstr[im] -> SetMinimum(0.);
412  h2_xglu[im] -> SetMinimum(0.);
413 
414 
415 
416  cnv->Divide(2,3);
417 
418  TPaletteAxis * palette = 0;
419  tex->SetTextSize(0.03);
420 
421  cnv->cd(1); gPad->SetLogx(); gPad->SetLogy();
422  h2_xuv [im] -> Draw("colz");
423  gPad->Update();
424  palette = (TPaletteAxis*)h2_xuv[im]->GetListOfFunctions()->FindObject("palette");
425  if(palette) {
426  palette->SetX1NDC(0.2);
427  palette->SetX2NDC(0.25);
428  palette->SetY1NDC(0.4);
429  palette->SetY2NDC(0.8);
430  }
431  tex->DrawTextNDC(0.30, 0.95, TString::Format("x*uv: %s", gPDFAlgList[im]->Id().Key().c_str()).Data() );
432 
433  cnv->cd(2); gPad->SetLogx(); gPad->SetLogy();
434  h2_xdv [im] -> Draw("colz");
435  gPad->Update();
436  palette = (TPaletteAxis*)h2_xdv[im]->GetListOfFunctions()->FindObject("palette");
437  if(palette) {
438  palette->SetX1NDC(0.2);
439  palette->SetX2NDC(0.25);
440  palette->SetY1NDC(0.4);
441  palette->SetY2NDC(0.8);
442  }
443  tex->DrawTextNDC(0.30, 0.95, TString::Format("x*dv: %s", gPDFAlgList[im]->Id().Key().c_str()).Data() );
444 
445  cnv->cd(3); gPad->SetLogx(); gPad->SetLogy();
446  h2_xus [im] -> Draw("colz");
447  gPad->Update();
448  palette = (TPaletteAxis*)h2_xus[im]->GetListOfFunctions()->FindObject("palette");
449  if(palette) {
450  palette->SetX1NDC(0.2);
451  palette->SetX2NDC(0.25);
452  palette->SetY1NDC(0.4);
453  palette->SetY2NDC(0.8);
454  }
455  tex->DrawTextNDC(0.30, 0.95, TString::Format("x*us: %s", gPDFAlgList[im]->Id().Key().c_str()).Data() );
456 
457  cnv->cd(4); gPad->SetLogx(); gPad->SetLogy();
458  h2_xds [im] -> Draw("colz");
459  gPad->Update();
460  palette = (TPaletteAxis*)h2_xds[im]->GetListOfFunctions()->FindObject("palette");
461  if(palette) {
462  palette->SetX1NDC(0.2);
463  palette->SetX2NDC(0.25);
464  palette->SetY1NDC(0.4);
465  palette->SetY2NDC(0.8);
466  }
467  tex->DrawTextNDC(0.30, 0.95, TString::Format("x*ds: %s", gPDFAlgList[im]->Id().Key().c_str()).Data() );
468 
469  cnv->cd(5); gPad->SetLogx(); gPad->SetLogy();
470  h2_xstr[im] -> Draw("colz");
471  gPad->Update();
472  palette = (TPaletteAxis*)h2_xstr[im]->GetListOfFunctions()->FindObject("palette");
473  if(palette) {
474  palette->SetX1NDC(0.2);
475  palette->SetX2NDC(0.25);
476  palette->SetY1NDC(0.4);
477  palette->SetY2NDC(0.8);
478  }
479  tex->DrawTextNDC(0.30, 0.95, TString::Format("x*str: %s", gPDFAlgList[im]->Id().Key().c_str()).Data() );
480 
481  cnv->cd(6); gPad->SetLogx(); gPad->SetLogy();
482  h2_xglu[im] -> Draw("colz");
483  gPad->Update();
484  palette = (TPaletteAxis*)h2_xglu[im]->GetListOfFunctions()->FindObject("palette");
485  if(palette) {
486  palette->SetX1NDC(0.2);
487  palette->SetX2NDC(0.25);
488  palette->SetY1NDC(0.4);
489  palette->SetY2NDC(0.8);
490  }
491  tex->DrawTextNDC(0.30, 0.95, TString::Format("x*glu: %s", gPDFAlgList[im]->Id().Key().c_str()).Data() );
492 
493  cnv->Update();
494  }
495 
496  //
497  // For multiple PDFs sets, plot the ratio of each PDF = f(x,Q2) to the first one
498  //
499 
500  for(unsigned int im=1; im < gPDFAlgList.size(); im++) {
501  h2_xuv_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
502  h2_xdv_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
503  h2_xus_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
504  h2_xds_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
505  h2_xstr_r[im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
506  h2_xglu_r[im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
507  for(int ibinx = 1;
508  ibinx <= h2_xuv[im]->GetXaxis()->GetNbins(); ibinx++) {
509  for(int ibinq2 = 1;
510  ibinq2 <= h2_xuv[im]->GetYaxis()->GetNbins(); ibinq2++) {
511 
512  double xuv = h2_xuv [im] -> GetBinContent(ibinx, ibinq2);
513  double xdv = h2_xdv [im] -> GetBinContent(ibinx, ibinq2);
514  double xus = h2_xus [im] -> GetBinContent(ibinx, ibinq2);
515  double xds = h2_xds [im] -> GetBinContent(ibinx, ibinq2);
516  double xstr = h2_xstr[im] -> GetBinContent(ibinx, ibinq2);
517  double xglu = h2_xglu[im] -> GetBinContent(ibinx, ibinq2);
518 
519  double xuv0 = h2_xuv [0] -> GetBinContent(ibinx, ibinq2);
520  double xdv0 = h2_xdv [0] -> GetBinContent(ibinx, ibinq2);
521  double xus0 = h2_xus [0] -> GetBinContent(ibinx, ibinq2);
522  double xds0 = h2_xds [0] -> GetBinContent(ibinx, ibinq2);
523  double xstr0 = h2_xstr[0] -> GetBinContent(ibinx, ibinq2);
524  double xglu0 = h2_xglu[0] -> GetBinContent(ibinx, ibinq2);
525 
526  double xuv_r = ((xuv0 > 0.) ? (xuv -xuv0 )/xuv0 : 0.);
527  double xdv_r = ((xdv0 > 0.) ? (xdv -xdv0 )/xdv0 : 0.);
528  double xus_r = ((xus0 > 0.) ? (xus -xus0 )/xus0 : 0.);
529  double xds_r = ((xds0 > 0.) ? (xds -xds0 )/xds0 : 0.);
530  double xstr_r = ((xstr0 > 0.) ? (xstr-xstr0)/xstr0 : 0.);
531  double xglu_r = ((xglu0 > 0.) ? (xglu-xglu0)/xglu0 : 0.);
532 
533  h2_xuv_r [im] -> SetBinContent(ibinx, ibinq2, xuv_r );
534  h2_xdv_r [im] -> SetBinContent(ibinx, ibinq2, xdv_r );
535  h2_xus_r [im] -> SetBinContent(ibinx, ibinq2, xus_r );
536  h2_xds_r [im] -> SetBinContent(ibinx, ibinq2, xds_r );
537  h2_xstr_r[im] -> SetBinContent(ibinx, ibinq2, xstr_r);
538  h2_xglu_r[im] -> SetBinContent(ibinx, ibinq2, xglu_r);
539  }
540  }
541 
542  ps->NewPage();
543 
544  h2_xuv_r [im] -> GetXaxis() -> SetTitle("x");
545  h2_xdv_r [im] -> GetXaxis() -> SetTitle("x");
546  h2_xus_r [im] -> GetXaxis() -> SetTitle("x");
547  h2_xds_r [im] -> GetXaxis() -> SetTitle("x");
548  h2_xstr_r[im] -> GetXaxis() -> SetTitle("x");
549  h2_xglu_r[im] -> GetXaxis() -> SetTitle("x");
550 
551  h2_xuv_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
552  h2_xdv_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
553  h2_xus_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
554  h2_xds_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
555  h2_xstr_r[im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
556  h2_xglu_r[im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
557 
558  cnv->Divide(2,3);
559 
560  TPaletteAxis * palette = 0;
561 
562  cnv->cd(1); gPad->SetLogx(); gPad->SetLogy();
563  h2_xuv_r [im] -> Draw("colz");
564  gPad->Update();
565  palette = (TPaletteAxis*)h2_xuv_r[im]->GetListOfFunctions()->FindObject("palette");
566  if(palette) {
567  palette->SetX1NDC(0.2);
568  palette->SetX2NDC(0.25);
569  palette->SetY1NDC(0.4);
570  palette->SetY2NDC(0.8);
571  }
572  tex->DrawTextNDC(0.1, 0.95, TString::Format("uv: (%s-%s)/(%s)",
573  gPDFAlgList[im]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str()).Data() );
574 
575  cnv->cd(2); gPad->SetLogx(); gPad->SetLogy();
576  h2_xdv_r [im] -> Draw("colz");
577  gPad->Update();
578  palette = (TPaletteAxis*)h2_xdv_r[im]->GetListOfFunctions()->FindObject("palette");
579  if(palette) {
580  palette->SetX1NDC(0.2);
581  palette->SetX2NDC(0.25);
582  palette->SetY1NDC(0.4);
583  palette->SetY2NDC(0.8);
584  }
585  tex->DrawTextNDC(0.1, 0.95, TString::Format("dv: (%s-%s)/(%s)",
586  gPDFAlgList[im]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str()).Data() );
587 
588  cnv->cd(3); gPad->SetLogx(); gPad->SetLogy();
589  h2_xus_r [im] -> Draw("colz");
590  gPad->Update();
591  palette = (TPaletteAxis*)h2_xus_r[im]->GetListOfFunctions()->FindObject("palette");
592  if(palette) {
593  palette->SetX1NDC(0.2);
594  palette->SetX2NDC(0.25);
595  palette->SetY1NDC(0.4);
596  palette->SetY2NDC(0.8);
597  }
598  tex->DrawTextNDC(0.1, 0.95, TString::Format("us: (%s-%s)/(%s)",
599  gPDFAlgList[im]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str()).Data() );
600 
601  cnv->cd(4); gPad->SetLogx(); gPad->SetLogy();
602  h2_xds_r [im] -> Draw("colz");
603  gPad->Update();
604  palette = (TPaletteAxis*)h2_xds_r[im]->GetListOfFunctions()->FindObject("palette");
605  if(palette) {
606  palette->SetX1NDC(0.2);
607  palette->SetX2NDC(0.25);
608  palette->SetY1NDC(0.4);
609  palette->SetY2NDC(0.8);
610  }
611  tex->DrawTextNDC(0.1, 0.95, TString::Format("ds: (%s-%s)/(%s)",
612  gPDFAlgList[im]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str()).Data() );
613 
614  cnv->cd(5); gPad->SetLogx(); gPad->SetLogy();
615  h2_xstr_r[im] -> Draw("colz");
616  gPad->Update();
617  palette = (TPaletteAxis*)h2_xstr_r[im]->GetListOfFunctions()->FindObject("palette");
618  if(palette) {
619  palette->SetX1NDC(0.2);
620  palette->SetX2NDC(0.25);
621  palette->SetY1NDC(0.4);
622  palette->SetY2NDC(0.8);
623  }
624  tex->DrawTextNDC(0.1, 0.95, TString::Format("str: (%s-%s)/(%s)",
625  gPDFAlgList[im]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str()).Data() );
626 
627  cnv->cd(6); gPad->SetLogx(); gPad->SetLogy();
628  h2_xglu_r[im] -> Draw("colz");
629  gPad->Update();
630  palette = (TPaletteAxis*)h2_xglu_r[im]->GetListOfFunctions()->FindObject("palette");
631  if(palette) {
632  palette->SetX1NDC(0.2);
633  palette->SetX2NDC(0.25);
634  palette->SetY1NDC(0.4);
635  palette->SetY2NDC(0.8);
636  }
637  tex->DrawTextNDC(0.1, 0.95, TString::Format("glu: (%s-%s)/(%s)",
638  gPDFAlgList[im]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str(), gPDFAlgList[0]->Id().Key().c_str()).Data() );
639 
640  cnv->Update();
641  }
642 
643  ps->Close();
644 
645  delete cnv;
646  delete ps;
647  delete lgnd;
648 
649  string root_filename = gOptOutFile + ".root";
650  TFile f(root_filename.c_str(),"recreate");
651  ntpl->Write("pdflib");
652 
653  f.Close();
654  delete ntpl;
655 
656  std::cout<<"Done."<<std::endl;
657 }
658 //_________________________________________________________________________________
659 void GetCommandLineArgs(int argc, char ** argv)
660 {
661  CmdLnArgParser parser(argc,argv);
662 
663  if(parser.OptionExists("pdf-set")){
664  gOptPDFSet = parser.Arg("pdf-set");
665  LOG("gpdfcomp", pNOTICE) << "Input PDF set: " << gOptPDFSet;
666  } else {
667  LOG("gpdfcomp", pFATAL)
668  << "Please specify PDF sets using the --pdf-set argument";
669  gAbortingInErr = true;
670  exit(1);
671  }
672 
673  if(parser.OptionExists('o')){
674  gOptOutFile = parser.Arg('o');
675  }
676 
677 }
678 //_________________________________________________________________________________
679 void GetAlgorithms(void)
680 {
681  vector<string> vpdfset = str::Split(gOptPDFSet, ",");
682  LOG("gpdfcomp", pNOTICE)
683  << "Number of input PDF sets: " << vpdfset.size();
684  if(vpdfset.size() == 0) {
685  LOG("gpdfcomp", pFATAL) << "Need at least 1 PDF set!";
686  gAbortingInErr = true;
687  exit(1);
688  }
689  vector<string>::iterator vpdfset_iter = vpdfset.begin();
690  for( ; vpdfset_iter != vpdfset.end(); ++vpdfset_iter) {
691  vector<string> vpdf = str::Split(*vpdfset_iter, "/");
692  if(vpdf.size() != 2) {
693  LOG("gpdfcomp", pFATAL)
694  << "Need to specify both a PDF algorithm name and configuration "
695  << "as in genie::GRV98LO/Default";
696  gAbortingInErr = true;
697  exit(1);
698  }
699  string pdf_alg_name = vpdf[0];
700  string pdf_alg_conf = vpdf[1];
701 
702  AlgFactory * algf = AlgFactory::Instance();
703  const PDFModelI * pdf_alg =
704  dynamic_cast<const PDFModelI *> (
705  algf->GetAlgorithm(pdf_alg_name, pdf_alg_conf));
706 
707  if(!pdf_alg) {
708  LOG("gpdfcomp", pFATAL)
709  << "Couldn't instantiate " << pdf_alg_name << "/" << pdf_alg_conf;
710  gAbortingInErr = true;
711  exit(1);
712  }
713  LOG("gpdfcomp", pNOTICE)
714  << "\n Instantiated: " << pdf_alg->Id()
715  << " with the following configuration: "
716  << pdf_alg->GetConfig();
717 
718  gPDFAlgList.push_back(pdf_alg);
719  }
720 }
721 //_________________________________________________________________________________
722 
T max(const caf::Proxy< T > &a, T b)
tree Draw("slc.nhit")
void MakePlots(void)
Definition: gPDFComp.cxx:84
bin1_2sigma SetFillColor(3)
enum BeamMode kRed
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
h7 SetMaximum(maxY)
double Gluon(void) const
Definition: PDF.h:59
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
correl_xv GetYaxis() -> SetDecimals()
double DownValence(void) const
Definition: PDF.h:52
A class to store PDFs.
Definition: PDF.h:38
#define pFATAL
Definition: Messenger.h:57
static constexpr Double_t nm
Definition: Munits.h:133
double UpSea(void) const
Definition: PDF.h:53
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
Pure abstract base class. Defines the PDFModelI interface to be implemented by wrapper classes to exi...
Definition: PDFModelI.h:29
ntuple SetFillStyle(1001)
void SetDefaultStyle(bool black_n_white=false)
Definition: Style.cxx:28
vector< const PDFModelI * > gPDFAlgList
Definition: gPDFComp.cxx:65
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
const char * label
int main(int argc, char **argv)
Definition: gPDFComp.cxx:73
double DownSea(void) const
Definition: PDF.h:54
string gOptOutFile
Definition: gPDFComp.cxx:64
legend SetBorderSize(0)
double Strange(void) const
Definition: PDF.h:55
void GetAlgorithms(void)
Definition: gPDFComp.cxx:679
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
Int_t col[ntarg]
Definition: Style.C:29
correl_xv GetXaxis() -> SetDecimals()
void SetModel(const PDFModelI *model)
Definition: PDF.cxx:48
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
string gOptPDFSet
Definition: gPDFComp.cxx:63
void GetCommandLineArgs(int argc, char **argv)
Definition: gPDFComp.cxx:659
double GetY(int dcm=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:119
void Calculate(double x, double q2)
Definition: PDF.cxx:55
OStream cout
Definition: OStream.cxx:6
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double UpValence(void) const
Definition: PDF.h:51
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
TLatex * tex
Definition: f2_nu.C:499
exit(0)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
histo SetMinimum(1.0e-9)
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
enum BeamMode kGreen
enum BeamMode kBlue
bool gAbortingInErr
Definition: Messenger.cxx:56
bool OptionExists(char opt)
was option set?
Root of GENIE utility namespaces.
char * Arg(char opt)
return argument following -`opt&#39;
enum BeamMode string