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