GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-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/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 
59 using namespace std;
60 using namespace genie;
61 using namespace genie::utils;
62 
63 // globals
64 string gOptSF = ""; // --structure-func argument
65 string gOptOutFile = "sf_comp"; // -o argument
66 vector<const DISStructureFuncModelI *> gSFAlgList;
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 SF 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  // 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 //_________________________________________________________________________________
618 void 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 //_________________________________________________________________________________
639 void 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 
genie::DISStructureFunc::SetModel
void SetModel(const DISStructureFuncModelI *model)
Attach an algorithm.
Definition: DISStructureFunc.cxx:43
genie::DISStructureFuncModelI
Pure Abstract Base Class. Defines the DISStructureFuncModelI interface to be implemented by any algor...
Definition: DISStructureFuncModelI.h:28
MakePlots
void MakePlots(void)
Definition: gPDFComp.cxx:85
genie::DISStructureFunc::F6
double F6(void) const
Get the computed structure function F6.
Definition: DISStructureFunc.h:68
genie::DISStructureFunc::F2
double F2(void) const
Get the computed structure function F2.
Definition: DISStructureFunc.h:56
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Kinematics::SetQ2
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
genie::Interaction::KinePtr
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
gSFAlgList
vector< const DISStructureFuncModelI * > gSFAlgList
Definition: gSFComp.cxx:66
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:30
genie::utils::style::SetDefaultStyle
void SetDefaultStyle(bool black_n_white=false)
Definition: Style.cxx:20
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::DISStructureFunc::F4
double F4(void) const
Get the computed structure function F4.
Definition: DISStructureFunc.h:62
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
genie::DISStructureFunc::F1
double F1(void) const
Get the computed structure function F1.
Definition: DISStructureFunc.h:53
gOptOutFile
string gOptOutFile
Definition: gPDFComp.cxx:65
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::AlgFactory::GetAlgorithm
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
genie::kPdgTgtFreeP
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
gOptSF
string gOptSF
Definition: gSFComp.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::DISStructureFunc::F5
double F5(void) const
Get the computed structure function F5.
Definition: DISStructureFunc.h:65
genie::Kinematics::Setx
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
genie::DISStructureFunc
A class holding Deep Inelastic Scattering (DIS) Form Factors (invariant structure funstions) ...
Definition: DISStructureFunc.h:39
genie::DISStructureFunc::F3
double F3(void) const
Get the computed structure function F3.
Definition: DISStructureFunc.h:59
genie::Algorithm::Id
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
genie::utils::str::Split
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
genie::DISStructureFunc::Calculate
void Calculate(const Interaction *interaction)
Calculate the S/F&#39;s for the input interaction using the attached algorithm.
Definition: DISStructureFunc.cxx:49
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
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
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Generated on Wed Jul 17 2024 14:46:41 for GENIEGenerator by   doxygen 1.8.5