GENIEGenerator
Loading...
Searching...
No Matches
gSFComp.cxx
Go to the documentation of this file.
1//____________________________________________________________________________
2/*!
3
4\program gsfcomp
5
6\brief Structure function comparison tool
7
8\syntax gsfcomp --structure-func sf_set [-o output]
9
10 --structure-func :
11 Specifies a comma separated list of GENIE structure function models.
12 The full algorithm name and configuration should be provided for
13 each GENIE model 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: sf_comp
19
20\example gsfcomp --structure-func genie::Blah/Default,genie::Blah/Tweaked
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/DeepInelastic/XSection/DISStructureFunc.h"
52#include "Physics/DeepInelastic/XSection/DISStructureFuncModelI.h"
53#include "Framework/Interaction/Interaction.h"
54#include "Framework/ParticleData/PDGCodes.h"
55#include "Framework/Utils/CmdLnArgParser.h"
56#include "Framework/Utils/StringUtils.h"
57#include "Framework/Utils/Style.h"
58
59using namespace std;
60using namespace genie;
61using namespace genie::utils;
62
63// globals
64string gOptSF = ""; // --structure-func argument
65string gOptOutFile = "sf_comp"; // -o argument
66vector<const DISStructureFuncModelI *> gSFAlgList;
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 SF 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 // Output ntuple
145 TNtuple * ntpl = new TNtuple("nt","structure functions","i:F1:F2:F3:F4:F5:F6: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 structure function 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 < gSFAlgList.size(); im++) {
168 const char * label = gSFAlgList[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 // structure functions = f(Q^2) for selected vales of x (1-D plots)
185 TGraph * gr_F1_Q2 [nm] = { 0 };
186 TGraph * gr_F2_Q2 [nm] = { 0 };
187 TGraph * gr_F3_Q2 [nm] = { 0 };
188 TGraph * gr_F4_Q2 [nm] = { 0 };
189 TGraph * gr_F5_Q2 [nm] = { 0 };
190 TGraph * gr_F6_Q2 [nm] = { 0 };
191
192 // structure functions = f(x,Q^2) with fine x,Q2 binning (2-D plots)
193 TH2D * h2_F1 [nm] = { 0 };
194 TH2D * h2_F2 [nm] = { 0 };
195 TH2D * h2_F3 [nm] = { 0 };
196 TH2D * h2_F4 [nm] = { 0 };
197 TH2D * h2_F5 [nm] = { 0 };
198 TH2D * h2_F6 [nm] = { 0 };
199
200 // structure functions ratios relative to first one specified = f(x,Q^2) with fine x,Q2 binning (2-D plots)
201 TH2D * h2_F1_r [nm] = { 0 };
202 TH2D * h2_F2_r [nm] = { 0 };
203 TH2D * h2_F3_r [nm] = { 0 };
204 TH2D * h2_F4_r [nm] = { 0 };
205 TH2D * h2_F5_r [nm] = { 0 };
206 TH2D * h2_F6_r [nm] = { 0 };
207
208 //
209 // Generate tructure function 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 F1_arr [nm][nQ2];
215 double F2_arr [nm][nQ2];
216 double F3_arr [nm][nQ2];
217 double F4_arr [nm][nQ2];
218 double F5_arr [nm][nQ2];
219 double F6_arr [nm][nQ2];
220 for(unsigned int im=0; im < gSFAlgList.size(); im++) {
221 DISStructureFunc sf;
222 sf.SetModel(gSFAlgList[im]);
223 for(unsigned int iq2 = 0; iq2 < nQ2; iq2++) {
224 double Q2 = Q2_arr[iq2];
225 Interaction * interaction = Interaction::DISCC(kPdgTgtFreeP,kPdgProton,kPdgNuMu);
226 //LOG("gsfcomp", pNOTICE) << "Setting x = " << x << ", Q2 = " << Q2;
227 interaction->KinePtr()->Setx(x);
228 interaction->KinePtr()->SetQ2(Q2);
229 sf.Calculate(interaction);
230 double F1 = sf.F1();
231 double F2 = sf.F2();
232 double F3 = sf.F3();
233 double F4 = sf.F4();
234 double F5 = sf.F5();
235 double F6 = sf.F6();
236 F1_arr [im][iq2] = F1;
237 F2_arr [im][iq2] = F2;
238 F3_arr [im][iq2] = F3;
239 F4_arr [im][iq2] = F4;
240 F5_arr [im][iq2] = F5;
241 F6_arr [im][iq2] = F6;
242 ntpl->Fill(im,F1,F2,F3,F4,F5,F6,x,Q2);
243 delete interaction;
244 }//iq2
245
246 gr_F1_Q2 [im] = new TGraph (nQ2, Q2_arr, F1_arr [im]);
247 gr_F2_Q2 [im] = new TGraph (nQ2, Q2_arr, F2_arr [im]);
248 gr_F3_Q2 [im] = new TGraph (nQ2, Q2_arr, F3_arr [im]);
249 gr_F4_Q2 [im] = new TGraph (nQ2, Q2_arr, F4_arr [im]);
250 gr_F5_Q2 [im] = new TGraph (nQ2, Q2_arr, F5_arr [im]);
251 gr_F6_Q2 [im] = new TGraph (nQ2, Q2_arr, F6_arr [im]);
252
253 genie::utils::style::Format( gr_F1_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
254 genie::utils::style::Format( gr_F2_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
255 genie::utils::style::Format( gr_F3_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
256 genie::utils::style::Format( gr_F4_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
257 genie::utils::style::Format( gr_F5_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
258 genie::utils::style::Format( gr_F6_Q2 [im], col[im], sty[im], 2, col[im], mrk[im], msz[im]);
259
260 gr_F1_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
261 gr_F2_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
262 gr_F3_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
263 gr_F4_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
264 gr_F5_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
265 gr_F6_Q2 [im] -> GetXaxis() -> SetTitle("Q^{2} (GeV^{2}/c^{2})");
266
267 gr_F1_Q2 [im] -> GetYaxis() -> SetTitle("F1(x,Q^{2})");
268 gr_F2_Q2 [im] -> GetYaxis() -> SetTitle("F2(x,Q^{2})");
269 gr_F3_Q2 [im] -> GetYaxis() -> SetTitle("F3(x,Q^{2})");
270 gr_F4_Q2 [im] -> GetYaxis() -> SetTitle("F4(x,Q^{2})");
271 gr_F5_Q2 [im] -> GetYaxis() -> SetTitle("F5(x,Q^{2})");
272 gr_F6_Q2 [im] -> GetYaxis() -> SetTitle("F6(x,Q^{2})");
273
274 }//im
275
276 ps->NewPage();
277
278 if(x<0.3) {
279 lgnd = new TLegend(0.20, 0.55, 0.50, 0.85);
280 } else {
281 lgnd = new TLegend(0.60, 0.55, 0.80, 0.85);
282 }
283 lgnd -> SetFillColor(0);
284 lgnd -> SetFillStyle(0);
285 lgnd -> SetBorderSize(0);
286
287 lgnd->Clear();
288 for(unsigned int im=0; im < gSFAlgList.size(); im++) {
289 const char * label = gSFAlgList[im]->Id().Key().c_str();
290 lgnd->AddEntry(gr_F1_Q2 [im], Form("%s",label), lgopt[im]);
291 }
292
293 cnv->Divide(2,3);
294
295 cnv->cd(1); gPad->SetLogx();
296 cnv->cd(2); gPad->SetLogx();
297 cnv->cd(3); gPad->SetLogx();
298 cnv->cd(4); gPad->SetLogx();
299 cnv->cd(5); gPad->SetLogx();
300 cnv->cd(6); gPad->SetLogx();
301
302 for(unsigned int im=0; im < gSFAlgList.size(); im++) {
303 cnv->cd(1); gr_F1_Q2 [im]->Draw(opt[im]);
304 cnv->cd(2); gr_F2_Q2 [im]->Draw(opt[im]);
305 cnv->cd(3); gr_F3_Q2 [im]->Draw(opt[im]);
306 cnv->cd(4); gr_F4_Q2 [im]->Draw(opt[im]);
307 cnv->cd(5); gr_F5_Q2 [im]->Draw(opt[im]);
308 cnv->cd(6); gr_F6_Q2 [im]->Draw(opt[im]);
309 }
310
311 cnv->cd(1); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, Form("x = %.3e",x));
312 cnv->cd(2); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, Form("x = %.3e",x));
313 cnv->cd(3); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, Form("x = %.3e",x));
314 cnv->cd(4); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, Form("x = %.3e",x));
315 cnv->cd(5); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, Form("x = %.3e",x));
316 cnv->cd(6); lgnd->Draw(); tex->DrawTextNDC(0.4, 0.95, Form("x = %.3e",x));
317
318 cnv->Update();
319 }
320
321 //
322 // Plot structure functions = f(x,Q2)
323 //
324
325 for(unsigned int im=0; im < gSFAlgList.size(); im++) {
326 h2_F1 [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
327 h2_F2 [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
328 h2_F3 [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
329 h2_F4 [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
330 h2_F5 [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
331 h2_F6 [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
332 DISStructureFunc sf;
333 sf.SetModel(gSFAlgList[im]);
334 for(int ibinx = 1;
335 ibinx <= h2_F1[im]->GetXaxis()->GetNbins(); ibinx++) {
336 double x = h2_F1[im]->GetXaxis()->GetBinCenter(ibinx);
337 for(int ibinq2 = 1;
338 ibinq2 <= h2_F1[im]->GetYaxis()->GetNbins(); ibinq2++) {
339 double Q2 = h2_F1[im]->GetYaxis()->GetBinCenter(ibinq2);
340 Interaction * interaction = Interaction::DISCC(kPdgTgtFreeP,kPdgProton,kPdgNuMu);
341 //LOG("gsfcomp", pNOTICE) << "Setting x (bin " << ibinx << ") = " << x << ", Q2 (bin " << ibinq2 << ") = " << Q2;
342 interaction->KinePtr()->Setx(x);
343 interaction->KinePtr()->SetQ2(Q2);
344 sf.Calculate(interaction);
345 double F1 = sf.F1();
346 double F2 = sf.F2();
347 double F3 = sf.F3();
348 double F4 = sf.F4();
349 double F5 = sf.F5();
350 double F6 = sf.F6();
351 h2_F1 [im] -> SetBinContent(ibinx, ibinq2, F1);
352 h2_F2 [im] -> SetBinContent(ibinx, ibinq2, F2);
353 h2_F3 [im] -> SetBinContent(ibinx, ibinq2, F3);
354 h2_F4 [im] -> SetBinContent(ibinx, ibinq2, F4);
355 h2_F5 [im] -> SetBinContent(ibinx, ibinq2, F5);
356 h2_F6 [im] -> SetBinContent(ibinx, ibinq2, F6);
357 delete interaction;
358 }
359 }
360
361 ps->NewPage();
362
363 h2_F1 [im] -> GetXaxis() -> SetTitle("x");
364 h2_F2 [im] -> GetXaxis() -> SetTitle("x");
365 h2_F3 [im] -> GetXaxis() -> SetTitle("x");
366 h2_F4 [im] -> GetXaxis() -> SetTitle("x");
367 h2_F5 [im] -> GetXaxis() -> SetTitle("x");
368 h2_F6 [im] -> GetXaxis() -> SetTitle("x");
369
370 h2_F1 [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
371 h2_F2 [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
372 h2_F3 [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
373 h2_F4 [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
374 h2_F5 [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
375 h2_F6 [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
376
377 cnv->Divide(2,3);
378
379 TPaletteAxis * palette = 0;
380 tex->SetTextSize(0.03);
381
382 cnv->cd(1); gPad->SetLogx(); gPad->SetLogy();
383 h2_F1 [im] -> Draw("colz");
384 gPad->Update();
385 palette = (TPaletteAxis*)h2_F1[im]->GetListOfFunctions()->FindObject("palette");
386 if(palette) {
387 palette->SetX1NDC(0.8);
388 palette->SetX2NDC(0.85);
389 palette->SetY1NDC(0.4);
390 palette->SetY2NDC(0.8);
391 }
392 tex->DrawTextNDC(0.30, 0.95, Form("F1: %s", gSFAlgList[im]->Id().Key().c_str()) );
393
394 cnv->cd(2); gPad->SetLogx(); gPad->SetLogy();
395 h2_F2 [im] -> Draw("colz");
396 gPad->Update();
397 palette = (TPaletteAxis*)h2_F2[im]->GetListOfFunctions()->FindObject("palette");
398 if(palette) {
399 palette->SetX1NDC(0.8);
400 palette->SetX2NDC(0.85);
401 palette->SetY1NDC(0.4);
402 palette->SetY2NDC(0.8);
403 }
404 tex->DrawTextNDC(0.30, 0.95, Form("F2: %s", gSFAlgList[im]->Id().Key().c_str()) );
405
406 cnv->cd(3); gPad->SetLogx(); gPad->SetLogy();
407 h2_F3 [im] -> Draw("colz");
408 gPad->Update();
409 palette = (TPaletteAxis*)h2_F3[im]->GetListOfFunctions()->FindObject("palette");
410 if(palette) {
411 palette->SetX1NDC(0.8);
412 palette->SetX2NDC(0.85);
413 palette->SetY1NDC(0.4);
414 palette->SetY2NDC(0.8);
415 }
416 tex->DrawTextNDC(0.30, 0.95, Form("F3: %s", gSFAlgList[im]->Id().Key().c_str()) );
417
418 cnv->cd(4); gPad->SetLogx(); gPad->SetLogy();
419 h2_F4 [im] -> Draw("colz");
420 gPad->Update();
421 palette = (TPaletteAxis*)h2_F4[im]->GetListOfFunctions()->FindObject("palette");
422 if(palette) {
423 palette->SetX1NDC(0.8);
424 palette->SetX2NDC(0.85);
425 palette->SetY1NDC(0.4);
426 palette->SetY2NDC(0.8);
427 }
428 tex->DrawTextNDC(0.30, 0.95, Form("F4: %s", gSFAlgList[im]->Id().Key().c_str()) );
429
430 cnv->cd(5); gPad->SetLogx(); gPad->SetLogy();
431 h2_F5[im] -> Draw("colz");
432 gPad->Update();
433 palette = (TPaletteAxis*)h2_F5[im]->GetListOfFunctions()->FindObject("palette");
434 if(palette) {
435 palette->SetX1NDC(0.8);
436 palette->SetX2NDC(0.85);
437 palette->SetY1NDC(0.4);
438 palette->SetY2NDC(0.8);
439 }
440 tex->DrawTextNDC(0.30, 0.95, Form("F5: %s", gSFAlgList[im]->Id().Key().c_str()) );
441
442 cnv->cd(6); gPad->SetLogx(); gPad->SetLogy();
443 h2_F6[im] -> Draw("colz");
444 gPad->Update();
445 palette = (TPaletteAxis*)h2_F6[im]->GetListOfFunctions()->FindObject("palette");
446 if(palette) {
447 palette->SetX1NDC(0.8);
448 palette->SetX2NDC(0.85);
449 palette->SetY1NDC(0.4);
450 palette->SetY2NDC(0.8);
451 }
452 tex->DrawTextNDC(0.30, 0.95, Form("F6: %s", gSFAlgList[im]->Id().Key().c_str()) );
453
454 cnv->Update();
455 }
456
457 //
458 // In the case of multiple input structure functions,
459 // plot the ratio of each structure function = f(x,Q2) to the first one
460 //
461
462 for(unsigned int im=1; im < gSFAlgList.size(); im++) {
463 h2_F1_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
464 h2_F2_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
465 h2_F3_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
466 h2_F4_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
467 h2_F5_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
468 h2_F6_r [im] = new TH2D("","", nx_2d-1, x_bin_edges_2d, nQ2_2d-1, Q2_bin_edges_2d);
469 for(int ibinx = 1;
470 ibinx <= h2_F1[im]->GetXaxis()->GetNbins(); ibinx++) {
471 for(int ibinq2 = 1;
472 ibinq2 <= h2_F1[im]->GetYaxis()->GetNbins(); ibinq2++) {
473
474 double F1 = h2_F1 [im] -> GetBinContent(ibinx, ibinq2);
475 double F2 = h2_F2 [im] -> GetBinContent(ibinx, ibinq2);
476 double F3 = h2_F3 [im] -> GetBinContent(ibinx, ibinq2);
477 double F4 = h2_F4 [im] -> GetBinContent(ibinx, ibinq2);
478 double F5 = h2_F5 [im] -> GetBinContent(ibinx, ibinq2);
479 double F6 = h2_F6 [im] -> GetBinContent(ibinx, ibinq2);
480
481 double F10 = h2_F1 [0] -> GetBinContent(ibinx, ibinq2);
482 double F20 = h2_F2 [0] -> GetBinContent(ibinx, ibinq2);
483 double F30 = h2_F3 [0] -> GetBinContent(ibinx, ibinq2);
484 double F40 = h2_F4 [0] -> GetBinContent(ibinx, ibinq2);
485 double F50 = h2_F5 [0] -> GetBinContent(ibinx, ibinq2);
486 double F60 = h2_F6 [0] -> GetBinContent(ibinx, ibinq2);
487
488 double F1r = ((F10 > 0.) ? (F1-F10)/F10 : 0.);
489 double F2r = ((F20 > 0.) ? (F2-F20)/F20 : 0.);
490 double F3r = ((F30 > 0.) ? (F3-F30)/F30 : 0.);
491 double F4r = ((F40 > 0.) ? (F4-F40)/F40 : 0.);
492 double F5r = ((F50 > 0.) ? (F5-F50)/F50 : 0.);
493 double F6r = ((F60 > 0.) ? (F6-F60)/F60 : 0.);
494
495 h2_F1_r [im] -> SetBinContent(ibinx, ibinq2, F1r);
496 h2_F2_r [im] -> SetBinContent(ibinx, ibinq2, F2r);
497 h2_F3_r [im] -> SetBinContent(ibinx, ibinq2, F3r);
498 h2_F4_r [im] -> SetBinContent(ibinx, ibinq2, F4r);
499 h2_F5_r [im] -> SetBinContent(ibinx, ibinq2, F5r);
500 h2_F6_r [im] -> SetBinContent(ibinx, ibinq2, F6r);
501 }
502 }
503
504 ps->NewPage();
505
506 h2_F1_r [im] -> GetXaxis() -> SetTitle("x");
507 h2_F2_r [im] -> GetXaxis() -> SetTitle("x");
508 h2_F3_r [im] -> GetXaxis() -> SetTitle("x");
509 h2_F4_r [im] -> GetXaxis() -> SetTitle("x");
510 h2_F5_r [im] -> GetXaxis() -> SetTitle("x");
511 h2_F6_r [im] -> GetXaxis() -> SetTitle("x");
512
513 h2_F1_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
514 h2_F2_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
515 h2_F3_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
516 h2_F4_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
517 h2_F5_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
518 h2_F6_r [im] -> GetXaxis() -> SetTitle("Q^2 (GeV^2/c^2)");
519
520 cnv->Divide(2,3);
521
522 TPaletteAxis * palette = 0;
523
524 cnv->cd(1); gPad->SetLogx(); gPad->SetLogy();
525 h2_F1_r [im] -> Draw("colz");
526 gPad->Update();
527 palette = (TPaletteAxis*)h2_F1_r[im]->GetListOfFunctions()->FindObject("palette");
528 if(palette) {
529 palette->SetX1NDC(0.8);
530 palette->SetX2NDC(0.85);
531 palette->SetY1NDC(0.4);
532 palette->SetY2NDC(0.8);
533 }
534 tex->DrawTextNDC(0.1, 0.95, Form("F1: (%s-%s)/(%s)",
535 gSFAlgList[im]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str()) );
536
537 cnv->cd(2); gPad->SetLogx(); gPad->SetLogy();
538 h2_F2_r [im] -> Draw("colz");
539 gPad->Update();
540 palette = (TPaletteAxis*)h2_F2_r[im]->GetListOfFunctions()->FindObject("palette");
541 if(palette) {
542 palette->SetX1NDC(0.8);
543 palette->SetX2NDC(0.85);
544 palette->SetY1NDC(0.4);
545 palette->SetY2NDC(0.8);
546 }
547 tex->DrawTextNDC(0.1, 0.95, Form("F2: (%s-%s)/(%s)",
548 gSFAlgList[im]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str()) );
549
550 cnv->cd(3); gPad->SetLogx(); gPad->SetLogy();
551 h2_F3_r [im] -> Draw("colz");
552 gPad->Update();
553 palette = (TPaletteAxis*)h2_F3_r[im]->GetListOfFunctions()->FindObject("palette");
554 if(palette) {
555 palette->SetX1NDC(0.8);
556 palette->SetX2NDC(0.85);
557 palette->SetY1NDC(0.4);
558 palette->SetY2NDC(0.8);
559 }
560 tex->DrawTextNDC(0.1, 0.95, Form("F3: (%s-%s)/(%s)",
561 gSFAlgList[im]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str()) );
562
563 cnv->cd(4); gPad->SetLogx(); gPad->SetLogy();
564 h2_F4_r [im] -> Draw("colz");
565 gPad->Update();
566 palette = (TPaletteAxis*)h2_F4_r[im]->GetListOfFunctions()->FindObject("palette");
567 if(palette) {
568 palette->SetX1NDC(0.8);
569 palette->SetX2NDC(0.85);
570 palette->SetY1NDC(0.4);
571 palette->SetY2NDC(0.8);
572 }
573 tex->DrawTextNDC(0.1, 0.95, Form("F4: (%s-%s)/(%s)",
574 gSFAlgList[im]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str()) );
575
576 cnv->cd(5); gPad->SetLogx(); gPad->SetLogy();
577 h2_F5_r[im] -> Draw("colz");
578 gPad->Update();
579 palette = (TPaletteAxis*)h2_F5_r[im]->GetListOfFunctions()->FindObject("palette");
580 if(palette) {
581 palette->SetX1NDC(0.8);
582 palette->SetX2NDC(0.85);
583 palette->SetY1NDC(0.4);
584 palette->SetY2NDC(0.8);
585 }
586 tex->DrawTextNDC(0.1, 0.95, Form("F5: (%s-%s)/(%s)",
587 gSFAlgList[im]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str()) );
588
589 cnv->cd(6); gPad->SetLogx(); gPad->SetLogy();
590 h2_F6_r[im] -> Draw("colz");
591 gPad->Update();
592 palette = (TPaletteAxis*)h2_F6_r[im]->GetListOfFunctions()->FindObject("palette");
593 if(palette) {
594 palette->SetX1NDC(0.8);
595 palette->SetX2NDC(0.85);
596 palette->SetY1NDC(0.4);
597 palette->SetY2NDC(0.8);
598 }
599 tex->DrawTextNDC(0.1, 0.95, Form("F6: (%s-%s)/(%s)",
600 gSFAlgList[im]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str(), gSFAlgList[0]->Id().Key().c_str()) );
601
602 cnv->Update();
603 }
604
605 ps->Close();
606
607 delete cnv;
608 delete ps;
609 delete lgnd;
610
611 string root_filename = gOptOutFile + ".root";
612 TFile f(root_filename.c_str(),"recreate");
613 ntpl->Write("sfntp");
614 f.Close();
615 delete ntpl;
616}
617//_________________________________________________________________________________
618void GetCommandLineArgs(int argc, char ** argv)
619{
620 CmdLnArgParser parser(argc,argv);
621
622 if(parser.OptionExists("structure-function")){
623 gOptSF = parser.Arg("structure-function");
624 LOG("gsfcomp", pNOTICE) << "Input structure functions: " << gOptSF;
625 } else {
626 LOG("gsfcomp", pFATAL)
627 << "Please specify structure function models sets "
628 << "using the --structure-function argument";
629 gAbortingInErr = true;
630 exit(1);
631 }
632
633 if(parser.OptionExists('o')){
634 gOptOutFile = parser.Arg('o');
635 }
636
637}
638//_________________________________________________________________________________
639void GetAlgorithms(void)
640{
641 vector<string> vsfset = str::Split(gOptSF, ",");
642 LOG("gsfcomp", pNOTICE)
643 << "Number of input structure functions: " << vsfset.size();
644 if(vsfset.size() == 0) {
645 LOG("gsfcomp", pFATAL) << "Need at least 1 structure function!";
646 gAbortingInErr = true;
647 exit(1);
648 }
649 vector<string>::iterator vsfset_iter = vsfset.begin();
650 for( ; vsfset_iter != vsfset.end(); ++vsfset_iter) {
651 vector<string> vsf = str::Split(*vsfset_iter, "/");
652 if(vsf.size() != 2) {
653 LOG("gsfcomp", pFATAL)
654 << "Need to specify both a structire function algorithm name and configuration "
655 << "as in genie::BYStrucFunc/Default";
656 gAbortingInErr = true;
657 exit(1);
658 }
659 string sf_alg_name = vsf[0];
660 string sf_alg_conf = vsf[1];
661
662 AlgFactory * algf = AlgFactory::Instance();
663 const DISStructureFuncModelI * sf_alg =
664 dynamic_cast<const DISStructureFuncModelI *> (
665 algf->GetAlgorithm(sf_alg_name, sf_alg_conf));
666
667 if(!sf_alg) {
668 LOG("gsfcomp", pFATAL)
669 << "Couldn't instantiate " << sf_alg_name << "/" << sf_alg_conf;
670 gAbortingInErr = true;
671 exit(1);
672 }
673 LOG("gsfcomp", pNOTICE)
674 << "\n Instantiated: " << sf_alg->Id()
675 << " with the following configuration: "
676 << sf_alg->GetConfig();
677
678 gSFAlgList.push_back(sf_alg);
679 }
680}
681//_________________________________________________________________________________
682
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
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::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::DISStructureFuncModelI
Pure Abstract Base Class. Defines the DISStructureFuncModelI interface to be implemented by any algor...
Definition DISStructureFuncModelI.h:28
genie::DISStructureFunc
A class holding Deep Inelastic Scattering (DIS) Form Factors (invariant structure funstions)
Definition DISStructureFunc.h:39
genie::DISStructureFunc::Calculate
void Calculate(const Interaction *interaction)
Calculate the S/F's for the input interaction using the attached algorithm.
Definition DISStructureFunc.cxx:49
genie::DISStructureFunc::F1
double F1(void) const
Get the computed structure function F1.
Definition DISStructureFunc.h:53
genie::DISStructureFunc::F5
double F5(void) const
Get the computed structure function F5.
Definition DISStructureFunc.h:65
genie::DISStructureFunc::F3
double F3(void) const
Get the computed structure function F3.
Definition DISStructureFunc.h:59
genie::DISStructureFunc::F6
double F6(void) const
Get the computed structure function F6.
Definition DISStructureFunc.h:68
genie::DISStructureFunc::F4
double F4(void) const
Get the computed structure function F4.
Definition DISStructureFunc.h:62
genie::DISStructureFunc::SetModel
void SetModel(const DISStructureFuncModelI *model)
Attach an algorithm.
Definition DISStructureFunc.cxx:43
genie::DISStructureFunc::F2
double F2(void) const
Get the computed structure function F2.
Definition DISStructureFunc.h:56
genie::Interaction
Summary information for an interaction.
Definition Interaction.h:56
genie::Interaction::DISCC
static Interaction * DISCC(int tgt, int nuc, int probe, double E=0)
Definition Interaction.cxx:325
genie::Interaction::KinePtr
Kinematics * KinePtr(void) const
Definition Interaction.h:76
genie::Kinematics::Setx
void Setx(double x, bool selected=false)
Definition Kinematics.cxx:231
genie::Kinematics::SetQ2
void SetQ2(double Q2, bool selected=false)
Definition Kinematics.cxx:255
gOptOutFile
string gOptOutFile
Definition gPDFComp.cxx:65
gSFAlgList
vector< const DISStructureFuncModelI * > gSFAlgList
Definition gSFComp.cxx:66
GetAlgorithms
void GetAlgorithms(void)
Definition gSFComp.cxx:639
MakePlots
void MakePlots(void)
Definition gSFComp.cxx:85
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition gSFComp.cxx:618
gOptSF
string gOptSF
Definition gSFComp.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
genie::kPdgProton
const int kPdgProton
Definition PDGCodes.h:81
genie::kPdgTgtFreeP
const int kPdgTgtFreeP
Definition PDGCodes.h:199
genie::kPdgNuMu
const int kPdgNuMu
Definition PDGCodes.h:30
Generated on for GENIEGenerator by doxygen 1.14.0