GENIEGenerator
Loading...
Searching...
No Matches
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-2025, 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
59using namespace std;
60using namespace genie;
61using namespace genie::utils;
62
63// globals
64string gOptPDFSet = ""; // --pdf-set argument
65string gOptOutFile = "pdf_comp"; // -o argument
66vector<const PDFModelI *> gPDFAlgList;
67
68// function prototypes
69void GetCommandLineArgs (int argc, char ** argv);
70void GetAlgorithms (void);
71void MakePlots (void);
72
73//___________________________________________________________________
74int 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//_________________________________________________________________________________
85void 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//_________________________________________________________________________________
660void 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//_________________________________________________________________________________
683void 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
pNOTICE
#define pNOTICE
Definition Messenger.h:61
pFATAL
#define pFATAL
Definition Messenger.h:56
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
main
int main()
Definition Validation.cpp:290
genie::AlgFactory
The GENIE Algorithm Factory.
Definition AlgFactory.h:39
genie::AlgFactory::GetAlgorithm
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition AlgFactory.cxx:75
genie::AlgFactory::Instance
static AlgFactory * Instance()
Definition AlgFactory.cxx:64
genie::Algorithm::GetConfig
virtual const Registry & GetConfig(void) const
Definition Algorithm.cxx:246
genie::Algorithm::Id
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition Algorithm.h:98
genie::CmdLnArgParser
Command line argument parser.
Definition CmdLnArgParser.h:29
genie::CmdLnArgParser::OptionExists
bool OptionExists(char opt)
was option set?
Definition CmdLnArgParser.cxx:88
genie::CmdLnArgParser::Arg
char * Arg(char opt)
return argument following -‘opt’
Definition CmdLnArgParser.cxx:34
genie::PDFModelI
Pure abstract base class. Defines the PDFModelI interface to be implemented by wrapper classes to exi...
Definition PDFModelI.h:28
genie::PDF
A class to store PDFs.
Definition PDF.h:37
genie::PDF::UpSea
double UpSea(void) const
Definition PDF.h:52
genie::PDF::SetModel
void SetModel(const PDFModelI *model)
Definition PDF.cxx:42
genie::PDF::Gluon
double Gluon(void) const
Definition PDF.h:58
genie::PDF::DownSea
double DownSea(void) const
Definition PDF.h:53
genie::PDF::Calculate
void Calculate(double x, double q2)
Definition PDF.cxx:49
genie::PDF::UpValence
double UpValence(void) const
Definition PDF.h:50
genie::PDF::DownValence
double DownValence(void) const
Definition PDF.h:51
genie::PDF::Strange
double Strange(void) const
Definition PDF.h:54
genie::RunOpt::ReadFromCommandLine
void ReadFromCommandLine(int argc, char **argv)
Definition RunOpt.cxx:99
genie::RunOpt::Instance
static RunOpt * Instance(void)
Definition RunOpt.cxx:54
gPDFAlgList
vector< const PDFModelI * > gPDFAlgList
Definition gPDFComp.cxx:66
GetAlgorithms
void GetAlgorithms(void)
Definition gPDFComp.cxx:683
MakePlots
void MakePlots(void)
Definition gPDFComp.cxx:85
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition gPDFComp.cxx:660
gOptOutFile
string gOptOutFile
Definition gPDFComp.cxx:65
gOptPDFSet
string gOptPDFSet
Definition gPDFComp.cxx:64
Draw
void Draw(const char *plot, const char *title)
Definition gXSecComp.cxx:580
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::utils
Root of GENIE utility namespaces.
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
genie::gAbortingInErr
bool gAbortingInErr
Definition Messenger.cxx:34
Generated on for GENIEGenerator by doxygen 1.14.0