GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gEvComp.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gevcomp
5 
6 \brief A GENIE neutrino event sample comparison utility.
7  Reads-in a GENIE GHEP event tree and generates a postscript file
8  with MC truth quantity plots (& comparisons with a reference event
9  sample is specified)
10 
11  Syntax :
12  gevcomp -f sample [-r reference_sample]
13 
14  Options:
15  [] Denotes an optional argument
16  -f Specifies the GENIE/ROOT file with the generated event sample
17  -r Specifies another GENIE/ROOT event sample file for comparison
18  -n Specifies how many events to analyze [default: all]
19 
20  Notes:
21  The input ROOT files are the gst summary ntuples generated by
22  GENIE's gntpc utility
23 
24  Example:
25  gevcomp -f /path/gntp.1.gst.root -r /path/gntp.2.gst.root
26 
27 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
28  University of Liverpool
29 
30 \created June 06, 2008
31 
32 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
33  For the full text of the license visit http://copyright.genie-mc.org
34 
35 */
36 //____________________________________________________________________________
37 
38 #include <cassert>
39 #include <sstream>
40 #include <string>
41 
42 #include <TSystem.h>
43 #include <TFile.h>
44 #include <TDirectory.h>
45 #include <TTree.h>
46 #include <TVector3.h>
47 #include <TLorentzVector.h>
48 #include <TPostScript.h>
49 #include <TH1D.h>
50 #include <TH2D.h>
51 #include <TH3D.h>
52 #include <TMath.h>
53 #include <TCanvas.h>
54 #include <TPavesText.h>
55 #include <TText.h>
56 #include <TStyle.h>
57 #include <TLegend.h>
58 
59 #include "Framework/Messenger/Messenger.h"
60 #include "Framework/ParticleData/PDGUtils.h"
61 #include "Framework/ParticleData/PDGCodes.h"
62 #include "Framework/Utils/CmdLnArgParser.h"
63 #include "Framework/Utils/Style.h"
64 
65 using std::ostringstream;
66 using std::string;
67 
68 using namespace genie;
69 
70 // function prototypes
71 void GetCommandLineArgs (int argc, char ** argv);
72 void PrintSyntax (void);
73 bool CheckRootFilename (string filename);
74 string OutputFileName (string input_file_name);
75 void CreatePlots (string filename, string filename_ref);
76 
77 // command-line arguments
78 string gOptInpFile = ""; // (-f) input GENIE event sample file
79 string gOptInpFileRef = ""; // (-r) input GENIE event sample file (reference)
80 
81 //_________________________________________________________________________________
82 int main(int argc, char ** argv)
83 {
84  GetCommandLineArgs(argc,argv);
85  utils::style::SetDefaultStyle();
86  CreatePlots(gOptInpFile,gOptInpFileRef);
87 
88  LOG("gevcomp", pINFO) << "Done!";
89  return 0;
90 }
91 //_________________________________________________________________________________
92 void CreatePlots(string inp_filename, string inp_filename_ref)
93 {
94  TFile * fin_0 = 0;
95  TFile * fin_1 = 0;
96  TTree * gst_0 = 0;
97  TTree * gst_1 = 0;
98 
99  if(!CheckRootFilename(inp_filename)) {
100  LOG("gevcomp", pERROR) << "Input file: " << inp_filename << " doesn't exist";
101  return;
102  }
103  fin_0 = new TFile(inp_filename.c_str(),"READ");
104  gst_0 = (TTree *) fin_0->Get("gst");
105  assert(gst_0);
106 
107  if(CheckRootFilename(inp_filename_ref)) {
108  fin_1 = new TFile(inp_filename_ref.c_str(),"READ");
109  gst_1 = (TTree *) fin_1->Get("gst");
110  assert(gst_1);
111  }
112 
113  gst_0->SetLineColor(kBlack);
114  gst_0->SetLineWidth(3);
115  if(gst_1) {
116  gst_1->SetLineColor(kRed);
117  gst_1->SetMarkerColor(kRed);
118  gst_1->SetLineWidth(2);
119  gst_1->SetMarkerStyle(20);
120  gst_1->SetMarkerSize(1);
121  }
122 
123  // Set global plot style
124  //
125  gStyle->SetOptTitle(0);
126  gStyle->SetOptStat(0);
127  gStyle->SetHistTopMargin(0.33);
128  gStyle->SetHistMinimumZero(true);
129 
130  // Plotting options
131  //
132  bool monoenergetic_sample = true;
133  bool show_coh_plots = true;
134  bool show_calc_kinematics = true;
135  bool show_mult_per_proc = true;
136  bool show_primary_hadsyst = true;
137 
138  gst_0->Draw("1","tgt>1000010010","GOFF");
139  show_coh_plots = (gst_0->GetSelectedRows() > 0);
140 
141  TCanvas * c = new TCanvas("c","",20,20,500,650);
142  c->SetBorderMode(0);
143  c->SetFillColor(0);
144  c->SetGridx();
145  c->SetGridy();
146 
147  TLegend * ls = new TLegend(0.20,0.94,0.99,0.99);
148  ls->SetFillColor(0);
149  ls->SetBorderSize(0);
150 
151  string ps_filename = OutputFileName(inp_filename);
152  TPostScript * ps = new TPostScript(ps_filename.c_str(), 111);
153 
154  //
155  // SECTION: PS File Header
156  //
157 
158  ps->NewPage();
159  c->Range(0,0,100,100);
160  TPavesText hdr(10,40,90,70,3,"tr");
161  hdr.AddText("GENIE Event Sample Comparisons");
162  hdr.AddText(" ");
163  hdr.AddText(" ");
164  hdr.AddText("Event Sample:");
165  //hdr.AddText(sample);
166  hdr.AddText(" ");
167  hdr.AddText("Notes:");
168  hdr.AddText(" ");
169  hdr.AddText(" ");
170  hdr.Draw();
171  c->Update();
172 
173  //
174  // SECTION: Event Numbers
175  //
176 
177  TH1F * h0num = new TH1F("h0num","", 2, 0., 2.);
178  TH1F * h1num = new TH1F("h1num","", 2, 0., 2.);
179  float n0=0, n1=0;
180  ps->NewPage();
181  c->Range(0,0,100,100);
182  TPavesText evn(10,10,90,90,3,"tr");
183  evn.AddText("Event Numbers:");
184  evn.AddText(" ");
185 
186  gst_0->Draw("1>>h0num","1","goff");
187  if(gst_1) gst_1->Draw("1>>h1num","1","goff");
188  n0 = h0num->GetEntries();
189  n1 = (gst_1) ? h1num->GetEntries() : 0;
190  evn.AddText( Form("ALL : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
191 
192  gst_0->Draw("1>>h0num","qel","goff");
193  if(gst_1) gst_1->Draw("1>>h1num","qel","goff");
194  n0 = h0num->GetEntries();
195  n1 = (gst_1) ? h1num->GetEntries() : 0;
196  evn.AddText( Form("QEL : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
197 
198  gst_0->Draw("1>>h0num","qel&&cc","goff");
199  if(gst_1) gst_1->Draw("1>>h1num","qel&&cc","goff");
200  n0 = h0num->GetEntries();
201  n1 = (gst_1) ? h1num->GetEntries() : 0;
202  evn.AddText( Form("QEL-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
203 
204  gst_0->Draw("1>>h0num","qel&&nc","goff");
205  if(gst_1) gst_1->Draw("1>>h1num","qel&&nc","goff");
206  n0 = h0num->GetEntries();
207  n1 = (gst_1) ? h1num->GetEntries() : 0;
208  evn.AddText( Form("QEL-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
209 
210  gst_0->Draw("1>>h0num","res","goff");
211  if(gst_1) gst_1->Draw("1>>h1num","res","goff");
212  n0 = h0num->GetEntries();
213  n1 = (gst_1) ? h1num->GetEntries() : 0;
214  evn.AddText( Form("RES : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
215 
216  gst_0->Draw("1>>h0num","res&&cc","goff");
217  if(gst_1) gst_1->Draw("1>>h1num","res&&cc","goff");
218  n0 = h0num->GetEntries();
219  n1 = (gst_1) ? h1num->GetEntries() : 0;
220  evn.AddText( Form("RES-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
221 
222  gst_0->Draw("1>>h0num","res&&nc","goff");
223  if(gst_1) gst_1->Draw("1>>h1num","res&&nc","goff");
224  n0 = h0num->GetEntries();
225  n1 = (gst_1) ? h1num->GetEntries() : 0;
226  evn.AddText( Form("RES-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
227 
228  gst_0->Draw("1>>h0num","dis","goff");
229  if(gst_1) gst_1->Draw("1>>h1num","dis","goff");
230  n0 = h0num->GetEntries();
231  n1 = (gst_1) ? h1num->GetEntries() : 0;
232  evn.AddText( Form("DIS : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
233 
234  gst_0->Draw("1>>h0num","dis&&cc","goff");
235  if(gst_1) gst_1->Draw("1>>h1num","dis&&cc","goff");
236  n0 = h0num->GetEntries();
237  n1 = (gst_1) ? h1num->GetEntries() : 0;
238  evn.AddText( Form("DIS-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
239 
240  gst_0->Draw("1>>h0num","dis&&nc","goff");
241  if(gst_1) gst_1->Draw("1>>h1num","dis&&nc","goff");
242  n0 = h0num->GetEntries();
243  n1 = (gst_1) ? h1num->GetEntries() : 0;
244  evn.AddText( Form("DIS-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
245 
246  gst_0->Draw("1>>h0num","coh","goff");
247  if(gst_1) gst_1->Draw("1>>h1num","coh","goff");
248  n0 = h0num->GetEntries();
249  n1 = (gst_1) ? h1num->GetEntries() : 0;
250  evn.AddText( Form("COH : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
251 
252  gst_0->Draw("1>>h0num","coh&&cc","goff");
253  if(gst_1) gst_1->Draw("1>>h1num","coh&&cc","goff");
254  n0 = h0num->GetEntries();
255  n1 = (gst_1) ? h1num->GetEntries() : 0;
256  evn.AddText( Form("COH-CC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
257 
258  gst_0->Draw("1>>h0num","coh&&nc","goff");
259  if(gst_1) gst_1->Draw("1>>h1num","coh&&nc","goff");
260  n0 = h0num->GetEntries();
261  n1 = (gst_1) ? h1num->GetEntries() : 0;
262  evn.AddText( Form("COH-NC : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
263 
264  gst_0->Draw("1>>h0num","imd","goff");
265  if(gst_1) gst_1->Draw("1>>h1num","imd","goff");
266  n0 = h0num->GetEntries();
267  n1 = (gst_1) ? h1num->GetEntries() : 0;
268  evn.AddText( Form("IMD : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
269 
270  gst_0->Draw("1>>h0num","nuel","goff");
271  if(gst_1) gst_1->Draw("1>>h1num","nuel","goff");
272  n0 = h0num->GetEntries();
273  n1 = (gst_1) ? h1num->GetEntries() : 0;
274  evn.AddText( Form("NuE-EL : %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
275 
276  gst_0->Draw("1>>h0num","dis&&cc&&charm","goff");
277  if(gst_1) gst_1->Draw("1>>h1num","dis&&cc&&charm","goff");
278  n0 = h0num->GetEntries();
279  n1 = (gst_1) ? h1num->GetEntries() : 0;
280  evn.AddText( Form("DIS-CHARM: %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
281 
282  gst_0->Draw("1>>h0num","qel&&cc&&charm","goff");
283  if(gst_1) gst_1->Draw("1>>h1num","qel&&cc&&charm","goff");
284  n0 = h0num->GetEntries();
285  n1 = (gst_1) ? h1num->GetEntries() : 0;
286  evn.AddText( Form("QEL-CHARM: %7.0f [test sample], %7.0f [ref sample]", n0, n1) );
287 
288  evn.Draw();
289  c->Update();
290 
291  if(!monoenergetic_sample) {
292  gst_0->Draw("Ev","","");
293  if(gst_1) gst_1->Draw("Ev","","perrsame");
294  ls->Clear();
295  ls->SetHeader("Neutrino Energy Spectrum");
296  ls->Draw();
297  c->Update();
298  }
299 
300  //
301  // SECTION: Kinematics
302  //
303  ps->NewPage();
304  c->Clear();
305  c->Range(0,0,100,100);
306  TPavesText hdrk(10,40,90,70,3,"tr");
307  hdrk.AddText("Selected Kinematical Quantities");
308  hdrk.AddText(" ");
309  hdrk.Draw();
310  c->Update();
311 
312  //------ selected Q2 for all events
313  ps->NewPage();
314  gst_0->Draw("Q2s","Q2s>0","");
315  if(gst_1) gst_1->Draw("Q2s","Q2s>0","perrsame");
316  ls->Clear();
317  ls->SetHeader("selected Q2 for all events");
318  ls->Draw();
319  c->Update();
320 
321  //------ selected Q2 for QEL
322  ps->NewPage();
323  gst_0->Draw("Q2s","qel&&!charm","");
324  if(gst_1) gst_1->Draw("Q2s","qel&&!charm","perrsame");
325  ls->Clear();
326  ls->SetHeader("selected Q2 for QEL events");
327  ls->Draw();
328  c->Update();
329 
330  //------ selected Q2 for QEL CC
331  ps->NewPage();
332  gst_0->Draw("Q2s","qel&&cc&&!charm","");
333  if(gst_1) gst_1->Draw("Q2s","qel&&cc&&!charm","perrsame");
334  ls->Clear();
335  ls->SetHeader("selected Q2 for QEL CC events");
336  ls->Draw();
337  c->Update();
338 
339  //------ selected Q2 for QEL NC
340  ps->NewPage();
341  gst_0->Draw("Q2s","qel&&nc&&!charm","");
342  if(gst_1) gst_1->Draw("Q2s","qel&&nc&&!charm","perrsame");
343  ls->Clear();
344  ls->SetHeader("selected Q2 for QEL NC events");
345  ls->Draw();
346  c->Update();
347 
348  //------ selected Q2 for RES
349  ps->NewPage();
350  gst_0->Draw("Q2s","res","");
351  if(gst_1) gst_1->Draw("Q2s","res","perrsame");
352  ls->Clear();
353  ls->SetHeader("selected Q2 for RES events");
354  ls->Draw();
355  c->Update();
356 
357  //------ selected Q2 for RES CC
358  ps->NewPage();
359  gst_0->Draw("Q2s","res&&cc","");
360  if(gst_1) gst_1->Draw("Q2s","res&&cc","perrsame");
361  ls->Clear();
362  ls->SetHeader("selected Q2 for RES CC events");
363  ls->Draw();
364  c->Update();
365 
366  //------ selected Q2 for RES NC
367  ps->NewPage();
368  gst_0->Draw("Q2s","res&&nc","");
369  if(gst_1) gst_1->Draw("Q2s","res&&nc","perrsame");
370  ls->Clear();
371  ls->SetHeader("selected Q2 for RES NC events");
372  ls->Draw();
373  c->Update();
374 
375  //------ selected Q2 for DIS
376  ps->NewPage();
377  gst_0->Draw("Q2s","dis","");
378  if(gst_1) gst_1->Draw("Q2s","dis","perrsame");
379  ls->Clear();
380  ls->SetHeader("selected Q2 for DIS events");
381  ls->Draw();
382  c->Update();
383 
384  //------ selected Q2 for DIS CC
385  ps->NewPage();
386  gst_0->Draw("Q2s","dis&&cc","");
387  if(gst_1) gst_1->Draw("Q2s","dis&&cc","perrsame");
388  ls->Clear();
389  ls->SetHeader("selected Q2 for DIS CC events");
390  ls->Draw();
391  c->Update();
392 
393  //------ selected Q2 for DIS NC
394  ps->NewPage();
395  gst_0->Draw("Q2s","dis&&nc","");
396  if(gst_1) gst_1->Draw("Q2s","dis&&nc","perrsame");
397  ls->Clear();
398  ls->SetHeader("selected Q2 for DIS NC events");
399  ls->Draw();
400  c->Update();
401 
402  //------ selected Q2 for Charm/DIS
403  ps->NewPage();
404  gst_0->Draw("Q2s","dis&&charm","");
405  if(gst_1) gst_1->Draw("Q2s","dis&&charm","perrsame");
406  ls->Clear();
407  ls->SetHeader("selected Q2 for Charm/DIS events");
408  ls->Draw();
409  c->Update();
410 
411  if(show_coh_plots) {
412  //------ selected Q2 for COH
413  ps->NewPage();
414  gst_0->Draw("Q2s","coh","");
415  if(gst_1) gst_1->Draw("Q2s","coh","perrsame");
416  ls->Clear();
417  ls->SetHeader("selected Q2 for COH events");
418  ls->Draw();
419  c->Update();
420 
421  //------ selected Q2 for COH CC
422  ps->NewPage();
423  gst_0->Draw("Q2s","coh&&cc","");
424  if(gst_1) gst_1->Draw("Q2s","coh&&cc","perrsame");
425  ls->Clear();
426  ls->SetHeader("selected Q2 for COH CC events");
427  ls->Draw();
428  c->Update();
429 
430  //------ selected Q2 for COH NC
431  ps->NewPage();
432  gst_0->Draw("Q2s","coh&&nc","");
433  if(gst_1) gst_1->Draw("Q2s","coh&&nc","perrsame");
434  ls->Clear();
435  ls->SetHeader("selected Q2 for COH NC events");
436  ls->Draw();
437  c->Update();
438  }
439 
440  //------ selected W for all events
441  ps->NewPage();
442  gst_0->Draw("Ws","Ws>0","");
443  if(gst_1) gst_1->Draw("Ws","Ws>0","perrsame");
444  ls->Clear();
445  ls->SetHeader("selected W for all events");
446  ls->Draw();
447  c->Update();
448 
449  //------ selected W for QEL
450  ps->NewPage();
451  gst_0->Draw("Ws","qel&&!charm","");
452  if(gst_1) gst_1->Draw("Ws","qel&&!charm","perrsame");
453  ls->Clear();
454  ls->SetHeader("selected W for QEL events");
455  ls->Draw();
456  c->Update();
457 
458  //------ selected W for RES
459  ps->NewPage();
460  gst_0->Draw("Ws","res","");
461  if(gst_1) gst_1->Draw("Ws","res","perrsame");
462  ls->Clear();
463  ls->SetHeader("selected W for RES events");
464  ls->Draw();
465  c->Update();
466 
467  //------ selected W for DIS
468  ps->NewPage();
469  gst_0->Draw("Ws","dis","");
470  if(gst_1) gst_1->Draw("Ws","dis","perrsame");
471  ls->Clear();
472  ls->SetHeader("selected W for DIS events");
473  ls->Draw();
474  c->Update();
475 
476  //------ selected W for DIS CC
477  ps->NewPage();
478  gst_0->Draw("Ws","dis&&cc","");
479  if(gst_1) gst_1->Draw("Ws","dis&&cc","perrsame");
480  ls->Clear();
481  ls->SetHeader("selected W for DIS CC events");
482  ls->Draw();
483  c->Update();
484 
485  //------ selected W for DIS NC
486  ps->NewPage();
487  gst_0->Draw("Ws","dis&&nc","");
488  if(gst_1) gst_1->Draw("Ws","dis&&nc","perrsame");
489  ls->Clear();
490  ls->SetHeader("selected W for DIS NC events");
491  ls->Draw();
492  c->Update();
493 
494  //------ selected x for all events
495  ps->NewPage();
496  gst_0->Draw("xs","","");
497  if(gst_1) gst_1->Draw("xs","","perrsame");
498  ls->Clear();
499  ls->SetHeader("selected x for all events");
500  ls->Draw();
501  c->Update();
502 
503  //------ selected x for QEL
504  ps->NewPage();
505  gst_0->Draw("xs","qel&&!charm","");
506  if(gst_1) gst_1->Draw("xs","qel&&!charm","perrsame");
507  ls->Clear();
508  ls->SetHeader("selected x for QEL events");
509  ls->Draw();
510  c->Update();
511 
512  //------ selected x for RES
513  ps->NewPage();
514  gst_0->Draw("xs","res","");
515  if(gst_1) gst_1->Draw("xs","res","perrsame");
516  ls->Clear();
517  ls->SetHeader("selected x for RES events");
518  ls->Draw();
519  c->Update();
520 
521  //------ selected x for DIS
522  ps->NewPage();
523  gst_0->Draw("xs","dis","");
524  if(gst_1) gst_1->Draw("xs","dis","perrsame");
525  ls->Clear();
526  ls->SetHeader("selected x for DIS events");
527  ls->Draw();
528  c->Update();
529 
530  //------ selected x for DIS CC
531  ps->NewPage();
532  gst_0->Draw("xs","dis&&cc","");
533  if(gst_1) gst_1->Draw("xs","dis&&cc","perrsame");
534  ls->Clear();
535  ls->SetHeader("selected x for DIS CC events");
536  ls->Draw();
537  c->Update();
538 
539  //------ selected x for DIS NC
540  ps->NewPage();
541  gst_0->Draw("xs","dis&&nc","");
542  if(gst_1) gst_1->Draw("xs","dis&&nc","perrsame");
543  ls->Clear();
544  ls->SetHeader("selected x for DIS NC events");
545  ls->Draw();
546  c->Update();
547 
548  //------ selected x for Charm/DIS
549  ps->NewPage();
550  gst_0->Draw("xs","dis&&charm","");
551  if(gst_1) gst_1->Draw("xs","dis&&charm","perrsame");
552  ls->Clear();
553  ls->SetHeader("selected x for Charm/DIS events");
554  ls->Draw();
555  c->Update();
556 
557  if(show_coh_plots) {
558  //------ selected x for COH
559  ps->NewPage();
560  gst_0->Draw("xs","coh","");
561  if(gst_1) gst_1->Draw("xs","coh","perrsame");
562  ls->Clear();
563  ls->SetHeader("selected x for COH events");
564  ls->Draw();
565  c->Update();
566 
567  //------ selected x for COH CC
568  ps->NewPage();
569  gst_0->Draw("xs","coh&&cc","");
570  if(gst_1) gst_1->Draw("xs","coh&&cc","perrsame");
571  ls->Clear();
572  ls->SetHeader("selected x for COH CC events");
573  ls->Draw();
574  c->Update();
575 
576  //------ selected x for COH NC
577  ps->NewPage();
578  gst_0->Draw("xs","coh&&nc","");
579  if(gst_1) gst_1->Draw("xs","coh&&nc","perrsame");
580  ls->Clear();
581  ls->SetHeader("selected x for COH NC events");
582  ls->Draw();
583  c->Update();
584  }
585 
586  //------ selected y for all events
587  ps->NewPage();
588  gst_0->Draw("ys","","");
589  if(gst_1) gst_1->Draw("ys","","perrsame");
590  ls->Clear();
591  ls->SetHeader("selected y for all events");
592  ls->Draw();
593  c->Update();
594 
595  //------ selected y for QEL
596  ps->NewPage();
597  gst_0->Draw("ys","qel&&!charm","");
598  if(gst_1) gst_1->Draw("ys","qel&&!charm","perrsame");
599  ls->Clear();
600  ls->SetHeader("selected y for QEL events");
601  ls->Draw();
602  c->Update();
603 
604  //------ selected y for RES
605  ps->NewPage();
606  gst_0->Draw("ys","res","");
607  if(gst_1) gst_1->Draw("ys","res","perrsame");
608  ls->Clear();
609  ls->SetHeader("selected y for RES events");
610  ls->Draw();
611  c->Update();
612 
613  //------ selected y for DIS
614  ps->NewPage();
615  gst_0->Draw("ys","dis","");
616  if(gst_1) gst_1->Draw("ys","dis","perrsame");
617  ls->Clear();
618  ls->SetHeader("selected y for DIS events");
619  ls->Draw();
620  c->Update();
621 
622  //------ selected y for DIS CC
623  ps->NewPage();
624  gst_0->Draw("ys","dis&&cc","");
625  if(gst_1) gst_1->Draw("ys","dis&&cc","perrsame");
626  ls->Clear();
627  ls->SetHeader("selected y for DIS CC events");
628  ls->Draw();
629  c->Update();
630 
631  //------ selected y for DIS NC
632  ps->NewPage();
633  gst_0->Draw("ys","dis&&nc","");
634  if(gst_1) gst_1->Draw("ys","dis&&nc","perrsame");
635  ls->Clear();
636  ls->SetHeader("selected y for DIS NC events");
637  ls->Draw();
638  c->Update();
639 
640  //------ selected y for Charm/DIS
641  ps->NewPage();
642  gst_0->Draw("ys","dis&&charm","");
643  if(gst_1) gst_1->Draw("ys","dis&&charm","perrsame");
644  ls->Clear();
645  ls->SetHeader("selected y for Charm/DIS events");
646  ls->Draw();
647  c->Update();
648 
649  if(show_coh_plots) {
650  //------ selected y for COH
651  ps->NewPage();
652  gst_0->Draw("ys","coh","");
653  if(gst_1) gst_1->Draw("ys","coh","perrsame");
654  ls->Clear();
655  ls->SetHeader("selected y for COH events");
656  ls->Draw();
657  c->Update();
658 
659  //------ selected y for COH CC
660  ps->NewPage();
661  gst_0->Draw("ys","coh&&cc","");
662  if(gst_1) gst_1->Draw("ys","coh&&cc","perrsame");
663  ls->Clear();
664  ls->SetHeader("selected y for COH CC events");
665  ls->Draw();
666  c->Update();
667 
668  //------ selected y for COH NC
669  ps->NewPage();
670  gst_0->Draw("ys","coh&&nc","");
671  if(gst_1) gst_1->Draw("ys","coh&&nc","perrsame");
672  ls->Clear();
673  ls->SetHeader("selected y for COH NC events");
674  ls->Draw();
675  c->Update();
676 
677  //------ selected t for COH
678  ps->NewPage();
679  gst_0->Draw("ts","coh","");
680  if(gst_1) gst_1->Draw("ts","coh","perrsame");
681  ls->Clear();
682  ls->SetHeader("selected t for COH events");
683  ls->Draw();
684  c->Update();
685  }
686 
687  if(show_calc_kinematics) {
688 
689  //
690  // SECTION: Computed Kinematics
691  //
692  ps->NewPage();
693  c->Clear();
694  c->Range(0,0,100,100);
695  TPavesText hdrck(10,40,90,70,3,"tr");
696  hdrck.AddText("Kinematical Quantities");
697  hdrck.AddText(" ");
698  hdrck.AddText(" ");
699  hdrck.AddText("Similar to the previous set of plots but");
700  hdrck.AddText("showing 'computed' rather than 'selected' variables");
701  hdrck.Draw();
702  c->Update();
703 
704  //------ Q2 for all events
705  ps->NewPage();
706  gst_0->Draw("Q2","","");
707  if(gst_1) gst_1->Draw("Q2","","perrsame");
708  ls->Clear();
709  ls->SetHeader("computed Q2 for all events");
710  ls->Draw();
711  c->Update();
712 
713  //------ Q2 for QEL
714  ps->NewPage();
715  gst_0->Draw("Q2","qel&&!charm","");
716  if(gst_1) gst_1->Draw("Q2","qel&&!charm","perrsame");
717  ls->Clear();
718  ls->SetHeader("computed Q2 for QEL events");
719  ls->Draw();
720  c->Update();
721 
722  //------ Q2 for RES
723  ps->NewPage();
724  gst_0->Draw("Q2","res","");
725  if(gst_1) gst_1->Draw("Q2","res","perrsame");
726  ls->Clear();
727  ls->SetHeader("computed Q2 for RES events");
728  ls->Draw();
729  c->Update();
730 
731  //------ Q2 for DIS
732  ps->NewPage();
733  gst_0->Draw("Q2","dis","");
734  if(gst_1) gst_1->Draw("Q2","dis","perrsame");
735  ls->Clear();
736  ls->SetHeader("computed Q2 for DIS events");
737  ls->Draw();
738  c->Update();
739 
740  //------ x for all events
741  ps->NewPage();
742  gst_0->Draw("x","","");
743  if(gst_1) gst_1->Draw("x","","perrsame");
744  ls->Clear();
745  ls->SetHeader("computed x for all events");
746  ls->Draw();
747  c->Update();
748 
749  //------ x for QEL
750  ps->NewPage();
751  gst_0->Draw("x","qel&&!charm","");
752  if(gst_1) gst_1->Draw("x","qel&&!charm","perrsame");
753  ls->Clear();
754  ls->SetHeader("computed x for QEL events");
755  ls->Draw();
756  c->Update();
757 
758  //------ x for RES
759  ps->NewPage();
760  gst_0->Draw("x","res","");
761  if(gst_1) gst_1->Draw("x","res","perrsame");
762  ls->Clear();
763  ls->SetHeader("computed x for RES events");
764  ls->Draw();
765  c->Update();
766 
767  //------ x for DIS
768  ps->NewPage();
769  gst_0->Draw("x","dis","");
770  if(gst_1) gst_1->Draw("x","dis","perrsame");
771  ls->Clear();
772  ls->SetHeader("computed x for DIS events");
773  ls->Draw();
774  c->Update();
775 
776  //------ y for all events
777  ps->NewPage();
778  gst_0->Draw("y","","");
779  if(gst_1) gst_1->Draw("y","","perrsame");
780  ls->Clear();
781  ls->SetHeader("computed y for all events");
782  ls->Draw();
783  c->Update();
784 
785  //------ y for QEL
786  ps->NewPage();
787  gst_0->Draw("y","qel&&!charm","");
788  if(gst_1) gst_1->Draw("y","qel&&!charm","perrsame");
789  ls->Clear();
790  ls->SetHeader("computed y for QEL events");
791  ls->Draw();
792  c->Update();
793 
794  //------ y for RES
795  ps->NewPage();
796  gst_0->Draw("y","res","");
797  if(gst_1) gst_1->Draw("y","res","perrsame");
798  ls->Clear();
799  ls->SetHeader("computed y for RES events");
800  ls->Draw();
801  c->Update();
802 
803  //------ y for DIS
804  ps->NewPage();
805  gst_0->Draw("y","dis","");
806  if(gst_1) gst_1->Draw("y","dis","perrsame");
807  ls->Clear();
808  ls->SetHeader("computed y for DIS events");
809  ls->Draw();
810  c->Update();
811 
812  }//show?
813 
814 
815  //
816  // SECTION: Initial State nucleon
817  //
818  ps->NewPage();
819  c->Clear();
820  c->Range(0,0,100,100);
821  TPavesText hdrinuc(10,40,90,70,3,"tr");
822  hdrinuc.AddText("Initial state nucleon 4-Momentum");
823  hdrinuc.Draw();
824  c->Update();
825 
826  //------ selected hit nucleon px
827  ps->NewPage();
828  c->Clear();
829  c->Divide(2,2);
830  c->cd(1);
831  gst_0->Draw("pxn","","");
832  if(gst_1) gst_1->Draw("pxn","","perrsame");
833  c->cd(2);
834  gst_0->Draw("pyn","","");
835  if(gst_1) gst_1->Draw("pyn","","perrsame");
836  c->cd(3);
837  gst_0->Draw("pzn","","");
838  if(gst_1) gst_1->Draw("pzn","","perrsame");
839  c->cd(4);
840  gst_0->Draw("En","En>.2","");
841  if(gst_1) gst_1->Draw("En","En>.2","perrsame");
842  c->Update();
843 
844  //
845  // SECTION: Final State Primary Lepton
846  //
847  ps->NewPage();
848  c->Clear();
849  c->Range(0,0,100,100);
850  TPavesText hdrfsl(10,40,90,70,3,"tr");
851  hdrfsl.AddText("Final State Primary Lepton 4-Momentum");
852  hdrfsl.Draw();
853  c->Update();
854 
855  //------ f/s primary lepton : all events
856  ps->NewPage();
857  c->Divide(2,2);
858  c->cd(1);
859  gst_0->Draw("pxl","","");
860  if(gst_1) gst_1->Draw("pxl","","perrsame");
861  c->cd(2);
862  gst_0->Draw("pyl","","");
863  if(gst_1) gst_1->Draw("pyl","","perrsame");
864  c->cd(3);
865  gst_0->Draw("pzl","","");
866  if(gst_1) gst_1->Draw("pzl","","perrsame");
867  c->cd(4);
868  gst_0->Draw("El","","");
869  if(gst_1) gst_1->Draw("El","","perrsame");
870  c->cd();
871  ls->Clear();
872  ls->SetHeader("Final state primary lepton 4-p: All events");
873  ls->Draw();
874  c->Update();
875 
876  //------ f/s primary lepton : all CC events
877  ps->NewPage();
878  c->Clear();
879  c->Divide(2,2);
880  c->cd(1);
881  gst_0->Draw("pxl","cc","");
882  if(gst_1) gst_1->Draw("pxl","cc","perrsame");
883  c->cd(2);
884  gst_0->Draw("pyl","cc","");
885  if(gst_1) gst_1->Draw("pyl","cc","perrsame");
886  c->cd(3);
887  gst_0->Draw("pzl","cc","");
888  if(gst_1) gst_1->Draw("pzl","cc","perrsame");
889  c->cd(4);
890  gst_0->Draw("El","cc","");
891  if(gst_1) gst_1->Draw("El","cc","perrsame");
892  c->cd();
893  ls->Clear();
894  ls->SetHeader("Final state primary lepton 4-p: All CC events");
895  ls->Draw();
896  c->Update();
897 
898  //------ f/s primary lepton : all NC events
899  ps->NewPage();
900  c->Clear();
901  c->Divide(2,2);
902  c->cd(1);
903  gst_0->Draw("pxl","nc","");
904  if(gst_1) gst_1->Draw("pxl","nc","perrsame");
905  c->cd(2);
906  gst_0->Draw("pyl","nc","");
907  if(gst_1) gst_1->Draw("pyl","nc","perrsame");
908  c->cd(3);
909  gst_0->Draw("pzl","nc","");
910  if(gst_1) gst_1->Draw("pzl","nc","perrsame");
911  c->cd(4);
912  gst_0->Draw("El","nc","");
913  if(gst_1) gst_1->Draw("El","nc","perrsame");
914  c->cd();
915  ls->Clear();
916  ls->SetHeader("Final state primary lepton 4-p: All NC events");
917  ls->Draw();
918  c->Update();
919 
920  //------ f/s primary lepton : QEL events
921  ps->NewPage();
922  c->Clear();
923  c->Divide(2,2);
924  c->cd(1);
925  gst_0->Draw("pxl","qel&&!charm","");
926  if(gst_1) gst_1->Draw("pxl","qel&&!charm","perrsame");
927  c->cd(2);
928  gst_0->Draw("pyl","qel&&!charm","");
929  if(gst_1) gst_1->Draw("pyl","qel&&!charm","perrsame");
930  c->cd(3);
931  gst_0->Draw("pzl","qel&&!charm","");
932  if(gst_1) gst_1->Draw("pzl","qel&&!charm","perrsame");
933  c->cd(4);
934  gst_0->Draw("El","qel&&!charm","");
935  if(gst_1) gst_1->Draw("El","qel&&!charm","perrsame");
936  c->cd();
937  ls->Clear();
938  ls->SetHeader("Final state primary lepton 4-p: QEL events");
939  ls->Draw();
940  c->Update();
941 
942  //------ f/s primary lepton : RES events
943  ps->NewPage();
944  c->Clear();
945  c->Divide(2,2);
946  c->cd(1);
947  gst_0->Draw("pxl","res","");
948  if(gst_1) gst_1->Draw("pxl","res","perrsame");
949  c->cd(2);
950  gst_0->Draw("pyl","res","");
951  if(gst_1) gst_1->Draw("pyl","res","perrsame");
952  c->cd(3);
953  gst_0->Draw("pzl","res","");
954  if(gst_1) gst_1->Draw("pzl","res","perrsame");
955  c->cd(4);
956  gst_0->Draw("El","res","");
957  if(gst_1) gst_1->Draw("El","res","perrsame");
958  c->cd();
959  ls->Clear();
960  ls->SetHeader("Final state primary lepton 4-p: RES events");
961  ls->Draw();
962  c->Update();
963 
964  //------ f/s primary lepton : DIS events
965  ps->NewPage();
966  c->Clear();
967  c->Divide(2,2);
968  c->cd(1);
969  gst_0->Draw("pxl","dis","");
970  if(gst_1) gst_1->Draw("pxl","dis","perrsame");
971  c->cd(2);
972  gst_0->Draw("pyl","dis","");
973  if(gst_1) gst_1->Draw("pyl","dis","perrsame");
974  c->cd(3);
975  gst_0->Draw("pzl","dis","");
976  if(gst_1) gst_1->Draw("pzl","dis","perrsame");
977  c->cd(4);
978  gst_0->Draw("El","dis","");
979  if(gst_1) gst_1->Draw("El","dis","perrsame");
980  c->cd();
981  ls->Clear();
982  ls->SetHeader("Final state primary lepton 4-p: All DIS events");
983  ls->Draw();
984  c->Update();
985 
986  if(show_coh_plots) {
987  //------ f/s primary lepton : COH events
988  ps->NewPage();
989  c->Clear();
990  c->Divide(2,2);
991  c->cd(1);
992  gst_0->Draw("pxl","coh","");
993  if(gst_1) gst_1->Draw("pxl","coh","perrsame");
994  c->cd(2);
995  gst_0->Draw("pyl","coh","");
996  if(gst_1) gst_1->Draw("pyl","coh","perrsame");
997  c->cd(3);
998  gst_0->Draw("pzl","coh","");
999  if(gst_1) gst_1->Draw("pzl","coh","perrsame");
1000  c->cd(4);
1001  gst_0->Draw("El","coh","");
1002  if(gst_1) gst_1->Draw("El","coh","perrsame");
1003  c->cd();
1004  ls->Clear();
1005  ls->SetHeader("Final state primary lepton 4-p: COH events");
1006  ls->Draw();
1007  c->Update();
1008  }
1009 
1010  //
1011  // SECTION: Final State Hadronic System Multiplicities & 4P
1012  //
1013  ps->NewPage();
1014  c->Clear();
1015  c->Range(0,0,100,100);
1016  TPavesText hdrfhad(10,40,90,70,3,"tr");
1017  hdrfhad.AddText("Final State Hadronic System");
1018  hdrfhad.AddText("Multiplicities and 4-Momenta");
1019  hdrfhad.AddText(" ");
1020  hdrfhad.AddText(" ");
1021  hdrfhad.AddText(" ");
1022  hdrfhad.AddText(" ");
1023  hdrfhad.AddText("Note:");
1024  hdrfhad.AddText("For nuclear targets these plots include the effect");
1025  hdrfhad.AddText("of intranuclear hadron transport / rescattering");
1026  hdrfhad.Draw();
1027  c->Update();
1028 
1029  //------ number of final state p
1030  ps->NewPage();
1031  gst_0->Draw("nfp","","");
1032  if(gst_1) gst_1->Draw("nfp","","perrsame");
1033  ls->Clear();
1034  ls->SetHeader("Number of final state protons");
1035  ls->Draw();
1036  c->Update();
1037 
1038  //------ number of final state n
1039  ps->NewPage();
1040  gst_0->Draw("nfn","","");
1041  if(gst_1) gst_1->Draw("nfn","","perrsame");
1042  ls->Clear();
1043  ls->SetHeader("Number of final state neutrons");
1044  ls->Draw();
1045  c->Update();
1046 
1047  //------ number of final state pi+
1048  ps->NewPage();
1049  gst_0->Draw("nfpip","","");
1050  if(gst_1) gst_1->Draw("nfpip","","perrsame");
1051  ls->Clear();
1052  ls->SetHeader("Number of final state pi+");
1053  ls->Draw();
1054  c->Update();
1055 
1056  //------ number of final state pi-
1057  ps->NewPage();
1058  gst_0->Draw("nfpim","","");
1059  if(gst_1) gst_1->Draw("nfpim","","perrsame");
1060  ls->Clear();
1061  ls->SetHeader("Number of final state pi-");
1062  ls->Draw();
1063  c->Update();
1064 
1065  //------ number of final state pi0
1066  ps->NewPage();
1067  gst_0->Draw("nfpi0","","");
1068  if(gst_1) gst_1->Draw("nfpi0","","perrsame");
1069  ls->Clear();
1070  ls->SetHeader("Number of final state pi0");
1071  ls->Draw();
1072  c->Update();
1073 
1074  //------ number of final state K+
1075  ps->NewPage();
1076  gst_0->Draw("nfkp","","");
1077  if(gst_1) gst_1->Draw("nfkp","","perrsame");
1078  ls->Clear();
1079  ls->SetHeader("Number of final state K+");
1080  ls->Draw();
1081  c->Update();
1082 
1083  //------ number of final state K-
1084  ps->NewPage();
1085  gst_0->Draw("nfkm","","");
1086  if(gst_1) gst_1->Draw("nfkm","","perrsame");
1087  ls->Clear();
1088  ls->SetHeader("Number of final state K-");
1089  ls->Draw();
1090  c->Update();
1091 
1092  //------ number of final state K0
1093  ps->NewPage();
1094  gst_0->Draw("nfk0","","");
1095  if(gst_1) gst_1->Draw("nfk0","","perrsame");
1096  ls->Clear();
1097  ls->SetHeader("Number of final state K0");
1098  ls->Draw();
1099  c->Update();
1100 
1101  //------ momentum of final state p
1102  ps->NewPage();
1103  c->Divide(2,2);
1104  c->cd(1);
1105  gst_0->Draw("pxf","pdgf==2212","");
1106  if(gst_1) gst_1->Draw("pxf","pdgf==2212","perrsame");
1107  c->cd(2);
1108  gst_0->Draw("pyf","pdgf==2212","");
1109  if(gst_1) gst_1->Draw("pyf","pdgf==2212","perrsame");
1110  c->cd(3);
1111  gst_0->Draw("pzf","pdgf==2212","");
1112  if(gst_1) gst_1->Draw("pzf","pdgf==2212","perrsame");
1113  c->cd(4);
1114  gst_0->Draw("Ef","pdgf==2212","");
1115  if(gst_1) gst_1->Draw("Ef","pdgf==2212","perrsame");
1116  c->cd();
1117  ls->Clear();
1118  ls->SetHeader("Final state protons 4-momentum");
1119  ls->Draw();
1120  c->Update();
1121 
1122  //------ momentum of final state n
1123  ps->NewPage();
1124  c->Clear();
1125  c->Divide(2,2);
1126  c->cd(1);
1127  gst_0->Draw("pxf","pdgf==2112","");
1128  if(gst_1) gst_1->Draw("pxf","pdgf==2112","perrsame");
1129  c->cd(2);
1130  gst_0->Draw("pyf","pdgf==2112","");
1131  if(gst_1) gst_1->Draw("pyf","pdgf==2112","perrsame");
1132  c->cd(3);
1133  gst_0->Draw("pzf","pdgf==2112","");
1134  if(gst_1) gst_1->Draw("pzf","pdgf==2112","perrsame");
1135  c->cd(4);
1136  gst_0->Draw("Ef","pdgf==2112","");
1137  if(gst_1) gst_1->Draw("Ef","pdgf==2112","perrsame");
1138  c->cd();
1139  ls->Clear();
1140  ls->SetHeader("Final state neutrons 4-momentum");
1141  ls->Draw();
1142  c->Update();
1143 
1144  //------ momentum of final state pi0
1145  ps->NewPage();
1146  c->Clear();
1147  c->Divide(2,2);
1148  c->cd(1);
1149  gst_0->Draw("pxf","pdgf==111","");
1150  if(gst_1) gst_1->Draw("pxf","pdgf==111","perrsame");
1151  c->cd(2);
1152  gst_0->Draw("pyf","pdgf==111","");
1153  if(gst_1) gst_1->Draw("pyf","pdgf==111","perrsame");
1154  c->cd(3);
1155  gst_0->Draw("pzf","pdgf==111","");
1156  if(gst_1) gst_1->Draw("pzf","pdgf==111","perrsame");
1157  c->cd(4);
1158  gst_0->Draw("Ef","pdgf==111","");
1159  if(gst_1) gst_1->Draw("Ef","pdgf==111","perrsame");
1160  c->cd();
1161  ls->Clear();
1162  ls->SetHeader("Final state pi0's 4-momentum");
1163  ls->Draw();
1164  c->Update();
1165 
1166  //------ momentum of final state pi+
1167  ps->NewPage();
1168  c->Clear();
1169  c->Divide(2,2);
1170  c->cd(1);
1171  gst_0->Draw("pxf","pdgf==211","");
1172  if(gst_1) gst_1->Draw("pxf","pdgf==211","perrsame");
1173  c->cd(2);
1174  gst_0->Draw("pyf","pdgf==211","");
1175  if(gst_1) gst_1->Draw("pyf","pdgf==211","perrsame");
1176  c->cd(3);
1177  gst_0->Draw("pzf","pdgf==211","");
1178  if(gst_1) gst_1->Draw("pzf","pdgf==211","perrsame");
1179  c->cd(4);
1180  gst_0->Draw("Ef","pdgf==211","");
1181  if(gst_1) gst_1->Draw("Ef","pdgf==211","perrsame");
1182  c->cd();
1183  ls->Clear();
1184  ls->SetHeader("Final state pi+'s 4-momentum");
1185  ls->Draw();
1186  c->Update();
1187 
1188  //------ momentum of final state pi+
1189  ps->NewPage();
1190  c->Clear();
1191  c->Divide(2,2);
1192  c->cd(1);
1193  gst_0->Draw("pxf","pdgf==-211","");
1194  if(gst_1) gst_1->Draw("pxf","pdgf==-211","perrsame");
1195  c->cd(2);
1196  gst_0->Draw("pyf","pdgf==-211","");
1197  if(gst_1) gst_1->Draw("pyf","pdgf==-211","perrsame");
1198  c->cd(3);
1199  gst_0->Draw("pzf","pdgf==-211","");
1200  if(gst_1) gst_1->Draw("pzf","pdgf==-211","perrsame");
1201  c->cd(4);
1202  gst_0->Draw("Ef","pdgf==-211","");
1203  if(gst_1) gst_1->Draw("Ef","pdgf==-211","perrsame");
1204  c->cd();
1205  ls->Clear();
1206  ls->SetHeader("Final state pi-'s 4-momentum");
1207  ls->Draw();
1208  c->Update();
1209 
1210  if(show_mult_per_proc) {
1211 
1212  //
1213  // similarly but for QEL events only
1214  //
1215 
1216  //------ number of final state p /QEL
1217  ps->NewPage();
1218  if(gst_1) gst_1->Draw("nfp","qel&&!charm","");
1219  gst_0->Draw("nfp","qel&&!charm","perrsame");
1220  ls->Clear();
1221  ls->SetHeader("Number of final state protons / QEL only");
1222  ls->Draw();
1223  c->Update();
1224 
1225  //------ number of final state n /QEL
1226  ps->NewPage();
1227  if(gst_1) gst_1->Draw("nfn","qel&&!charm","");
1228  gst_0->Draw("nfn","qel&&!charm","perrsame");
1229  ls->Clear();
1230  ls->SetHeader("Number of final state neutrons / QEL only");
1231  ls->Draw();
1232  c->Update();
1233 
1234  //------ number of final state pi+ /QEL
1235  ps->NewPage();
1236  gst_0->Draw("nfpip","qel&&!charm","");
1237  if(gst_1) gst_1->Draw("nfpip","qel&&!charm","perrsame");
1238  ls->Clear();
1239  ls->SetHeader("Number of final state pi+ / QEL only");
1240  ls->Draw();
1241  c->Update();
1242 
1243  //------ number of final state pi- /QEL
1244  ps->NewPage();
1245  gst_0->Draw("nfpim","qel&&!charm","");
1246  if(gst_1) gst_1->Draw("nfpim","qel&&!charm","perrsame");
1247  ls->Clear();
1248  ls->SetHeader("Number of final state pi- / QEL only");
1249  ls->Draw();
1250  c->Update();
1251 
1252  //------ number of final state pi0 /QEL
1253  ps->NewPage();
1254  gst_0->Draw("nfpi0","qel&&!charm","");
1255  if(gst_1) gst_1->Draw("nfpi0","qel&&!charm","perrsame");
1256  ls->Clear();
1257  ls->SetHeader("Number of final state pi0 / QEL only");
1258  ls->Draw();
1259  c->Update();
1260 
1261  //------ number of final state K+ /QEL
1262  ps->NewPage();
1263  gst_0->Draw("nfkp","qel&&!charm","");
1264  if(gst_1) gst_1->Draw("nfkp","qel&&!charm","perrsame");
1265  ls->Clear();
1266  ls->SetHeader("Number of final state K+ / QEL only");
1267  ls->Draw();
1268  c->Update();
1269 
1270  //------ number of final state K- /QEL
1271  ps->NewPage();
1272  gst_0->Draw("nfkm","qel&&!charm","");
1273  if(gst_1) gst_1->Draw("nfkm","qel&&!charm","perrsame");
1274  ls->Clear();
1275  ls->SetHeader("Number of final state K- / QEL only");
1276  ls->Draw();
1277  c->Update();
1278 
1279  //------ number of final state K0 /QEL
1280  ps->NewPage();
1281  gst_0->Draw("nfk0","qel&&!charm","");
1282  if(gst_1) gst_1->Draw("nfk0","qel&&!charm","perrsame");
1283  ls->Clear();
1284  ls->SetHeader("Number of final state K0 / QEL only");
1285  ls->Draw();
1286  c->Update();
1287 
1288  //------ momentum of final state p /QEL
1289  ps->NewPage();
1290  c->Divide(2,2);
1291  c->cd(1);
1292  gst_0->Draw("pxf","qel&&!charm&&pdgf==2212","");
1293  if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==2212","perrsame");
1294  c->cd(2);
1295  gst_0->Draw("pyf","qel&&!charm&&pdgf==2212","");
1296  if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==2212","perrsame");
1297  c->cd(3);
1298  gst_0->Draw("pzf","qel&&!charm&&pdgf==2212","");
1299  if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==2212","perrsame");
1300  c->cd(4);
1301  gst_0->Draw("Ef","qel&&!charm&&pdgf==2212","");
1302  if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==2212","perrsame");
1303  c->cd();
1304  ls->Clear();
1305  ls->SetHeader("Final state protons 4-momentum / QEL only");
1306  ls->Draw();
1307  c->Update();
1308 
1309  //------ momentum of final state n /QEL
1310  ps->NewPage();
1311  c->Clear();
1312  c->Divide(2,2);
1313  c->cd(1);
1314  gst_0->Draw("pxf","qel&&!charm&&pdgf==2112","");
1315  if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==2112","perrsame");
1316  c->cd(2);
1317  gst_0->Draw("pyf","qel&&!charm&&pdgf==2112","");
1318  if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==2112","perrsame");
1319  c->cd(3);
1320  gst_0->Draw("pzf","qel&&!charm&&pdgf==2112","");
1321  if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==2112","perrsame");
1322  c->cd(4);
1323  gst_0->Draw("Ef","qel&&!charm&&pdgf==2112","");
1324  if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==2112","perrsame");
1325  c->cd();
1326  ls->Clear();
1327  ls->SetHeader("Final state neutrons 4-momentum / QEL only");
1328  ls->Draw();
1329  c->Update();
1330 
1331  //------ momentum of final state pi0 /QEL
1332  ps->NewPage();
1333  c->Clear();
1334  c->Divide(2,2);
1335  c->cd(1);
1336  gst_0->Draw("pxf","qel&&!charm&&pdgf==111","");
1337  if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==111","perrsame");
1338  c->cd(2);
1339  gst_0->Draw("pyf","qel&&!charm&&pdgf==111","");
1340  if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==111","perrsame");
1341  c->cd(3);
1342  gst_0->Draw("pzf","qel&&!charm&&pdgf==111","");
1343  if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==111","perrsame");
1344  c->cd(4);
1345  gst_0->Draw("Ef","qel&&!charm&&pdgf==111","");
1346  if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==111","perrsame");
1347  c->cd();
1348  ls->Clear();
1349  ls->SetHeader("Final state pi0's 4-momentum / QEL only");
1350  ls->Draw();
1351  c->Update();
1352 
1353  //------ momentum of final state pi+ /QEL
1354  ps->NewPage();
1355  c->Clear();
1356  c->Divide(2,2);
1357  c->cd(1);
1358  gst_0->Draw("pxf","qel&&!charm&&pdgf==211","");
1359  if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==211","perrsame");
1360  c->cd(2);
1361  gst_0->Draw("pyf","qel&&!charm&&pdgf==211","");
1362  if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==211","perrsame");
1363  c->cd(3);
1364  gst_0->Draw("pzf","qel&&!charm&&pdgf==211","");
1365  if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==211","perrsame");
1366  c->cd(4);
1367  gst_0->Draw("Ef","qel&&!charm&&pdgf==211","");
1368  if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==211","perrsame");
1369  c->cd();
1370  ls->Clear();
1371  ls->SetHeader("Final state pi+'s 4-momentum / QEL only");
1372  ls->Draw();
1373  c->Update();
1374 
1375  //------ momentum of final state pi+ /QEL
1376  ps->NewPage();
1377  c->Clear();
1378  c->Divide(2,2);
1379  c->cd(1);
1380  gst_0->Draw("pxf","qel&&!charm&&pdgf==-211","");
1381  if(gst_1) gst_1->Draw("pxf","qel&&!charm&&pdgf==-211","perrsame");
1382  c->cd(2);
1383  gst_0->Draw("pyf","qel&&!charm&&pdgf==-211","");
1384  if(gst_1) gst_1->Draw("pyf","qel&&!charm&&pdgf==-211","perrsame");
1385  c->cd(3);
1386  gst_0->Draw("pzf","qel&&!charm&&pdgf==-211","");
1387  if(gst_1) gst_1->Draw("pzf","qel&&!charm&&pdgf==-211","perrsame");
1388  c->cd(4);
1389  gst_0->Draw("Ef","qel&&!charm&&pdgf==-211","");
1390  if(gst_1) gst_1->Draw("Ef","qel&&!charm&&pdgf==-211","perrsame");
1391  c->cd();
1392  ls->Clear();
1393  ls->SetHeader("Final state pi-'s 4-momentum/ QEL only");
1394  ls->Draw();
1395  c->Update();
1396 
1397  //
1398  // similarly but for RES events only
1399  //
1400 
1401  //------ number of final state p /RES
1402  ps->NewPage();
1403  gst_0->Draw("nfp","res","");
1404  if(gst_1) gst_1->Draw("nfp","res","perrsame");
1405  ls->Clear();
1406  ls->SetHeader("Number of final state protons / RES only");
1407  ls->Draw();
1408  c->Update();
1409 
1410  //------ number of final state n /RES
1411  ps->NewPage();
1412  gst_0->Draw("nfn","res","");
1413  if(gst_1) gst_1->Draw("nfn","res","perrsame");
1414  ls->Clear();
1415  ls->SetHeader("Number of final state neutrons / RES only");
1416  ls->Draw();
1417  c->Update();
1418 
1419  //------ number of final state pi+ /RES
1420  ps->NewPage();
1421  gst_0->Draw("nfpip","res","");
1422  if(gst_1) gst_1->Draw("nfpip","res","perrsame");
1423  ls->Clear();
1424  ls->SetHeader("Number of final state pi+ / RES only");
1425  ls->Draw();
1426  c->Update();
1427 
1428  //------ number of final state pi- /RES
1429  ps->NewPage();
1430  gst_0->Draw("nfpim","res","");
1431  if(gst_1) gst_1->Draw("nfpim","res","perrsame");
1432  ls->Clear();
1433  ls->SetHeader("Number of final state pi- / RES only");
1434  ls->Draw();
1435  c->Update();
1436 
1437  //------ number of final state pi0 /RES
1438  ps->NewPage();
1439  gst_0->Draw("nfpi0","res","");
1440  if(gst_1) gst_1->Draw("nfpi0","res","perrsame");
1441  ls->Clear();
1442  ls->SetHeader("Number of final state pi0 / RES only");
1443  ls->Draw();
1444  c->Update();
1445 
1446  //------ number of final state K+ /RES
1447  ps->NewPage();
1448  gst_0->Draw("nfkp","res","");
1449  if(gst_1) gst_1->Draw("nfkp","res","perrsame");
1450  ls->Clear();
1451  ls->SetHeader("Number of final state K+ / RES only");
1452  ls->Draw();
1453  c->Update();
1454 
1455  //------ number of final state K- /RES
1456  ps->NewPage();
1457  gst_0->Draw("nfkm","res","");
1458  if(gst_1) gst_1->Draw("nfkm","res","perrsame");
1459  ls->Clear();
1460  ls->SetHeader("Number of final state K- / RES only");
1461  ls->Draw();
1462  c->Update();
1463 
1464  //------ number of final state K0 /RES
1465  ps->NewPage();
1466  gst_0->Draw("nfk0","res","");
1467  if(gst_1) gst_1->Draw("nfk0","res","perrsame");
1468  ls->Clear();
1469  ls->SetHeader("Number of final state K0 / RES only");
1470  ls->Draw();
1471  c->Update();
1472 
1473  //------ momentum of final state p /RES
1474  ps->NewPage();
1475  c->Divide(2,2);
1476  c->cd(1);
1477  gst_0->Draw("pxf","res&&pdgf==2212","");
1478  if(gst_1) gst_1->Draw("pxf","res&&pdgf==2212","perrsame");
1479  c->cd(2);
1480  gst_0->Draw("pyf","res&&pdgf==2212","");
1481  if(gst_1) gst_1->Draw("pyf","res&&pdgf==2212","perrsame");
1482  c->cd(3);
1483  gst_0->Draw("pzf","res&&pdgf==2212","");
1484  if(gst_1) gst_1->Draw("pzf","res&&pdgf==2212","perrsame");
1485  c->cd(4);
1486  gst_0->Draw("Ef","res&&pdgf==2212","");
1487  if(gst_1) gst_1->Draw("Ef","res&&pdgf==2212","perrsame");
1488  c->cd();
1489  ls->Clear();
1490  ls->SetHeader("Final state protons 4-momentum / RES only");
1491  ls->Draw();
1492  c->Update();
1493 
1494  //------ momentum of final state n /RES
1495  ps->NewPage();
1496  c->Clear();
1497  c->Divide(2,2);
1498  c->cd(1);
1499  gst_0->Draw("pxf","res&&pdgf==2112","");
1500  if(gst_1) gst_1->Draw("pxf","res&&pdgf==2112","perrsame");
1501  c->cd(2);
1502  gst_0->Draw("pyf","res&&pdgf==2112","");
1503  if(gst_1) gst_1->Draw("pyf","res&&pdgf==2112","perrsame");
1504  c->cd(3);
1505  gst_0->Draw("pzf","res&&pdgf==2112","");
1506  if(gst_1) gst_1->Draw("pzf","res&&pdgf==2112","perrsame");
1507  c->cd(4);
1508  gst_0->Draw("Ef","res&&pdgf==2112","");
1509  if(gst_1) gst_1->Draw("Ef","res&&pdgf==2112","perrsame");
1510  c->cd();
1511  ls->Clear();
1512  ls->SetHeader("Final state neutrons 4-momentum / RES only");
1513  ls->Draw();
1514  c->Update();
1515 
1516  //------ momentum of final state pi0 /RES
1517  ps->NewPage();
1518  c->Clear();
1519  c->Divide(2,2);
1520  c->cd(1);
1521  gst_0->Draw("pxf","res&&pdgf==111","");
1522  if(gst_1) gst_1->Draw("pxf","res&&pdgf==111","perrsame");
1523  c->cd(2);
1524  gst_0->Draw("pyf","res&&pdgf==111","");
1525  if(gst_1) gst_1->Draw("pyf","res&&pdgf==111","perrsame");
1526  c->cd(3);
1527  gst_0->Draw("pzf","res&&pdgf==111","");
1528  if(gst_1) gst_1->Draw("pzf","res&&pdgf==111","perrsame");
1529  c->cd(4);
1530  gst_0->Draw("Ef","res&&pdgf==111","");
1531  if(gst_1) gst_1->Draw("Ef","res&&pdgf==111","perrsame");
1532  c->cd();
1533  ls->Clear();
1534  ls->SetHeader("Final state pi0's 4-momentum / RES only");
1535  ls->Draw();
1536  c->Update();
1537 
1538  //------ momentum of final state pi+ /RES
1539  ps->NewPage();
1540  c->Clear();
1541  c->Divide(2,2);
1542  c->cd(1);
1543  gst_0->Draw("pxf","res&&pdgf==211","");
1544  if(gst_1) gst_1->Draw("pxf","res&&pdgf==211","perrsame");
1545  c->cd(2);
1546  gst_0->Draw("pyf","res&&pdgf==211","");
1547  if(gst_1) gst_1->Draw("pyf","res&&pdgf==211","perrsame");
1548  c->cd(3);
1549  gst_0->Draw("pzf","res&&pdgf==211","");
1550  if(gst_1) gst_1->Draw("pzf","res&&pdgf==211","perrsame");
1551  c->cd(4);
1552  gst_0->Draw("Ef","res&&pdgf==211","");
1553  if(gst_1) gst_1->Draw("Ef","res&&pdgf==211","perrsame");
1554  c->cd();
1555  ls->Clear();
1556  ls->SetHeader("Final state pi+'s 4-momentum / RES only");
1557  ls->Draw();
1558  c->Update();
1559 
1560  //------ momentum of final state pi+ /RES
1561  ps->NewPage();
1562  c->Clear();
1563  c->Divide(2,2);
1564  c->cd(1);
1565  gst_0->Draw("pxf","res&&pdgf==-211","");
1566  if(gst_1) gst_1->Draw("pxf","res&&pdgf==-211","perrsame");
1567  c->cd(2);
1568  gst_0->Draw("pyf","res&&pdgf==-211","");
1569  if(gst_1) gst_1->Draw("pyf","res&&pdgf==-211","perrsame");
1570  c->cd(3);
1571  gst_0->Draw("pzf","res&&pdgf==-211","");
1572  if(gst_1) gst_1->Draw("pzf","res&&pdgf==-211","perrsame");
1573  c->cd(4);
1574  gst_0->Draw("Ef","res&&pdgf==-211","");
1575  if(gst_1) gst_1->Draw("Ef","res&&pdgf==-211","perrsame");
1576  c->cd();
1577  ls->Clear();
1578  ls->SetHeader("Final state pi-'s 4-momentum/ RES only");
1579  ls->Draw();
1580  c->Update();
1581 
1582  //
1583  // similarly but for DIS events only
1584  //
1585 
1586  //------ number of final state p /DIS
1587  ps->NewPage();
1588  gst_0->Draw("nfp","dis","");
1589  if(gst_1) gst_1->Draw("nfp","dis","perrsame");
1590  ls->Clear();
1591  ls->SetHeader("Number of final state protons / DIS only");
1592  ls->Draw();
1593  c->Update();
1594 
1595  //------ number of final state n /DIS
1596  ps->NewPage();
1597  gst_0->Draw("nfn","dis","");
1598  if(gst_1) gst_1->Draw("nfn","dis","perrsame");
1599  ls->Clear();
1600  ls->SetHeader("Number of final state neutrons / DIS only");
1601  ls->Draw();
1602  c->Update();
1603 
1604  //------ number of final state pi+ /DIS
1605  ps->NewPage();
1606  gst_0->Draw("nfpip","dis","");
1607  if(gst_1) gst_1->Draw("nfpip","dis","perrsame");
1608  ls->Clear();
1609  ls->SetHeader("Number of final state pi+ / DIS only");
1610  ls->Draw();
1611  c->Update();
1612 
1613  //------ number of final state pi- /DIS
1614  ps->NewPage();
1615  gst_0->Draw("nfpim","dis","");
1616  if(gst_1) gst_1->Draw("nfpim","dis","perrsame");
1617  ls->Clear();
1618  ls->SetHeader("Number of final state pi- / DIS only");
1619  ls->Draw();
1620  c->Update();
1621 
1622  //------ number of final state pi0 /DIS
1623  ps->NewPage();
1624  gst_0->Draw("nfpi0","dis","");
1625  if(gst_1) gst_1->Draw("nfpi0","dis","perrsame");
1626  ls->Clear();
1627  ls->SetHeader("Number of final state pi0 / DIS only");
1628  ls->Draw();
1629  c->Update();
1630 
1631  //------ number of final state K+ /DIS
1632  ps->NewPage();
1633  gst_0->Draw("nfkp","dis","");
1634  if(gst_1) gst_1->Draw("nfkp","dis","perrsame");
1635  ls->Clear();
1636  ls->SetHeader("Number of final state K+ / DIS only");
1637  ls->Draw();
1638  c->Update();
1639 
1640  //------ number of final state K- /DIS
1641  ps->NewPage();
1642  gst_0->Draw("nfkm","dis","");
1643  if(gst_1) gst_1->Draw("nfkm","dis","perrsame");
1644  ls->Clear();
1645  ls->SetHeader("Number of final state K- / DIS only");
1646  ls->Draw();
1647  c->Update();
1648 
1649  //------ number of final state K0 /DIS
1650  ps->NewPage();
1651  gst_0->Draw("nfk0","dis","");
1652  if(gst_1) gst_1->Draw("nfk0","dis","perrsame");
1653  ls->Clear();
1654  ls->SetHeader("Number of final state K0 / DIS only");
1655  ls->Draw();
1656  c->Update();
1657 
1658  //------ momentum of final state p /DIS
1659  ps->NewPage();
1660  c->Divide(2,2);
1661  c->cd(1);
1662  gst_0->Draw("pxf","dis&&pdgf==2212","");
1663  if(gst_1) gst_1->Draw("pxf","dis&&pdgf==2212","perrsame");
1664  c->cd(2);
1665  gst_0->Draw("pyf","dis&&pdgf==2212","");
1666  if(gst_1) gst_1->Draw("pyf","dis&&pdgf==2212","perrsame");
1667  c->cd(3);
1668  gst_0->Draw("pzf","dis&&pdgf==2212","");
1669  if(gst_1) gst_1->Draw("pzf","dis&&pdgf==2212","perrsame");
1670  c->cd(4);
1671  gst_0->Draw("Ef","dis&&pdgf==2212","");
1672  if(gst_1) gst_1->Draw("Ef","dis&&pdgf==2212","perrsame");
1673  c->cd();
1674  ls->Clear();
1675  ls->SetHeader("Final state protons 4-momentum / DIS only");
1676  ls->Draw();
1677  c->Update();
1678 
1679  //------ momentum of final state n /DIS
1680  ps->NewPage();
1681  c->Clear();
1682  c->Divide(2,2);
1683  c->cd(1);
1684  gst_0->Draw("pxf","dis&&pdgf==2112","");
1685  if(gst_1) gst_1->Draw("pxf","dis&&pdgf==2112","perrsame");
1686  c->cd(2);
1687  gst_0->Draw("pyf","dis&&pdgf==2112","");
1688  if(gst_1) gst_1->Draw("pyf","dis&&pdgf==2112","perrsame");
1689  c->cd(3);
1690  gst_0->Draw("pzf","dis&&pdgf==2112","");
1691  if(gst_1) gst_1->Draw("pzf","dis&&pdgf==2112","perrsame");
1692  c->cd(4);
1693  gst_0->Draw("Ef","dis&&pdgf==2112","");
1694  if(gst_1) gst_1->Draw("Ef","dis&&pdgf==2112","perrsame");
1695  c->cd();
1696  ls->Clear();
1697  ls->SetHeader("Final state neutrons 4-momentum / DIS only");
1698  ls->Draw();
1699  c->Update();
1700 
1701  //------ momentum of final state pi0 /DIS
1702  ps->NewPage();
1703  c->Clear();
1704  c->Divide(2,2);
1705  c->cd(1);
1706  gst_0->Draw("pxf","dis&&pdgf==111","");
1707  if(gst_1) gst_1->Draw("pxf","dis&&pdgf==111","perrsame");
1708  c->cd(2);
1709  gst_0->Draw("pyf","dis&&pdgf==111","");
1710  if(gst_1) gst_1->Draw("pyf","dis&&pdgf==111","perrsame");
1711  c->cd(3);
1712  gst_0->Draw("pzf","dis&&pdgf==111","");
1713  if(gst_1) gst_1->Draw("pzf","dis&&pdgf==111","perrsame");
1714  c->cd(4);
1715  gst_0->Draw("Ef","dis&&pdgf==111","");
1716  if(gst_1) gst_1->Draw("Ef","dis&&pdgf==111","perrsame");
1717  c->cd();
1718  ls->Clear();
1719  ls->SetHeader("Final state pi0's 4-momentum / DIS only");
1720  ls->Draw();
1721  c->Update();
1722 
1723  //------ momentum of final state pi+ /DIS
1724  ps->NewPage();
1725  c->Clear();
1726  c->Divide(2,2);
1727  c->cd(1);
1728  gst_0->Draw("pxf","dis&&pdgf==211","");
1729  if(gst_1) gst_1->Draw("pxf","dis&&pdgf==211","perrsame");
1730  c->cd(2);
1731  gst_0->Draw("pyf","dis&&pdgf==211","");
1732  if(gst_1) gst_1->Draw("pyf","dis&&pdgf==211","perrsame");
1733  c->cd(3);
1734  gst_0->Draw("pzf","dis&&pdgf==211","");
1735  if(gst_1) gst_1->Draw("pzf","dis&&pdgf==211","perrsame");
1736  c->cd(4);
1737  gst_0->Draw("Ef","dis&&pdgf==211","");
1738  if(gst_1) gst_1->Draw("Ef","dis&&pdgf==211","perrsame");
1739  c->cd();
1740  ls->Clear();
1741  ls->SetHeader("Final state pi+'s 4-momentum / DIS only");
1742  ls->Draw();
1743  c->Update();
1744 
1745  //------ momentum of final state pi+ /DIS
1746  ps->NewPage();
1747  c->Clear();
1748  c->Divide(2,2);
1749  c->cd(1);
1750  gst_0->Draw("pxf","dis&&pdgf==-211","");
1751  if(gst_1) gst_1->Draw("pxf","dis&&pdgf==-211","perrsame");
1752  c->cd(2);
1753  gst_0->Draw("pyf","dis&&pdgf==-211","");
1754  if(gst_1) gst_1->Draw("pyf","dis&&pdgf==-211","perrsame");
1755  c->cd(3);
1756  gst_0->Draw("pzf","dis&&pdgf==-211","");
1757  if(gst_1) gst_1->Draw("pzf","dis&&pdgf==-211","perrsame");
1758  c->cd(4);
1759  gst_0->Draw("Ef","dis&&pdgf==-211","");
1760  if(gst_1) gst_1->Draw("Ef","dis&&pdgf==-211","perrsame");
1761  c->cd();
1762  ls->Clear();
1763  ls->SetHeader("Final state pi-'s 4-momentum/ DIS only");
1764  ls->Draw();
1765  c->Update();
1766 
1767  } // per-proc
1768 
1769  //
1770  // SECTION: Primary Hadronic System Multiplicities & 4P
1771  //
1772  if(show_primary_hadsyst) {
1773 
1774  ps->NewPage();
1775  c->Clear();
1776  c->Range(0,0,100,100);
1777  TPavesText hdrihad(10,40,90,70,3,"tr");
1778  hdrihad.AddText("Parimary Hadronic System");
1779  hdrihad.AddText("Multiplicities and 4-Momenta");
1780  hdrihad.AddText(" ");
1781  hdrihad.AddText(" ");
1782  hdrihad.AddText(" ");
1783  hdrihad.AddText(" ");
1784  hdrihad.AddText("Note:");
1785  hdrihad.AddText("For nuclear targets these plots show the hadronic system");
1786  hdrihad.AddText("BEFORE any intranuclear hadron transport / rescattering");
1787  hdrihad.Draw();
1788  c->Update();
1789 
1790  //------ number of prim p
1791  ps->NewPage();
1792  gst_0->Draw("nip","","");
1793  if(gst_1) gst_1->Draw("nip","","perrsame");
1794  ls->Clear();
1795  ls->SetHeader("Primary Hadronic System: Number of protons");
1796  ls->Draw();
1797  c->Update();
1798 
1799  //------ number of prim n
1800  ps->NewPage();
1801  gst_0->Draw("nin","","");
1802  if(gst_1) gst_1->Draw("nin","","perrsame");
1803  ls->Clear();
1804  ls->SetHeader("Primary Hadronic System: Number of neutrons");
1805  ls->Draw();
1806  c->Update();
1807 
1808  //------ number of prim pi+
1809  ps->NewPage();
1810  gst_0->Draw("nipip","","");
1811  if(gst_1) gst_1->Draw("nipip","","perrsame");
1812  ls->Clear();
1813  ls->SetHeader("Primary Hadronic System: Number of pi+");
1814  ls->Draw();
1815  c->Update();
1816 
1817  //------ number of prim pi-
1818  ps->NewPage();
1819  gst_0->Draw("nipim","","");
1820  if(gst_1) gst_1->Draw("nipim","","perrsame");
1821  ls->Clear();
1822  ls->SetHeader("Primary Hadronic System: Number of pi-");
1823  ls->Draw();
1824  c->Update();
1825 
1826  //------ number of prim pi0
1827  ps->NewPage();
1828  gst_0->Draw("nipi0","","");
1829  if(gst_1) gst_1->Draw("nipi0","","perrsame");
1830  ls->Clear();
1831  ls->SetHeader("Primary Hadronic System: Number of pi0");
1832  ls->Draw();
1833  c->Update();
1834 
1835  //------ number of prim K+
1836  ps->NewPage();
1837  gst_0->Draw("nikp","","");
1838  if(gst_1) gst_1->Draw("nikp","","perrsame");
1839  ls->Clear();
1840  ls->SetHeader("Primary Hadronic System: Number of K+");
1841  ls->Draw();
1842  c->Update();
1843 
1844  //------ number of prim K-
1845  ps->NewPage();
1846  gst_0->Draw("nikm","","");
1847  if(gst_1) gst_1->Draw("nikm","","perrsame");
1848  ls->Clear();
1849  ls->SetHeader("Primary Hadronic System: Number of K-");
1850  ls->Draw();
1851  c->Update();
1852 
1853  //------ number of prim K0
1854  ps->NewPage();
1855  gst_0->Draw("nik0","","");
1856  if(gst_1) gst_1->Draw("nik0","","perrsame");
1857  ls->Clear();
1858  ls->SetHeader("Primary Hadronic System: Number of K0");
1859  ls->Draw();
1860  c->Update();
1861 
1862  //------ momentum of prim, p
1863  ps->NewPage();
1864  c->Divide(2,2);
1865  c->cd(1);
1866  gst_0->Draw("pxi","pdgi==2212","");
1867  if(gst_1) gst_1->Draw("pxi","pdgi==2212","perrsame");
1868  c->cd(2);
1869  gst_0->Draw("pyi","pdgi==2212","");
1870  if(gst_1) gst_1->Draw("pyi","pdgi==2212","perrsame");
1871  c->cd(3);
1872  gst_0->Draw("pzi","pdgi==2212","");
1873  if(gst_1) gst_1->Draw("pzi","pdgi==2212","perrsame");
1874  c->cd(4);
1875  gst_0->Draw("Ei","pdgi==2212","");
1876  if(gst_1) gst_1->Draw("Ei","pdgi==2212","perrsame");
1877  c->cd();
1878  ls->Clear();
1879  ls->SetHeader("Primary Hadronic System: proton 4-momentum");
1880  ls->Draw();
1881  c->Update();
1882 
1883  //------ momentum of prim. n
1884  ps->NewPage();
1885  c->Clear();
1886  c->Divide(2,2);
1887  c->cd(1);
1888  gst_0->Draw("pxi","pdgi==2112","");
1889  if(gst_1) gst_1->Draw("pxi","pdgi==2112","perrsame");
1890  c->cd(2);
1891  gst_0->Draw("pyi","pdgi==2112","");
1892  if(gst_1) gst_1->Draw("pyi","pdgi==2112","perrsame");
1893  c->cd(3);
1894  gst_0->Draw("pzi","pdgi==2112","");
1895  if(gst_1) gst_1->Draw("pzi","pdgi==2112","perrsame");
1896  c->cd(4);
1897  gst_0->Draw("Ei","pdgi==2112","");
1898  if(gst_1) gst_1->Draw("Ei","pdgi==2112","perrsame");
1899  c->cd();
1900  ls->Clear();
1901  ls->SetHeader("Primary Hadronic System: neutron 4-momentum");
1902  ls->Draw();
1903  c->Update();
1904 
1905  //------ momentum of prim. pi0
1906  ps->NewPage();
1907  c->Clear();
1908  c->Divide(2,2);
1909  c->cd(1);
1910  gst_0->Draw("pxi","pdgi==111","");
1911  if(gst_1) gst_1->Draw("pxi","pdgi==111","perrsame");
1912  c->cd(2);
1913  gst_0->Draw("pyi","pdgi==111","");
1914  if(gst_1) gst_1->Draw("pyi","pdgi==111","perrsame");
1915  c->cd(3);
1916  gst_0->Draw("pzi","pdgi==111","");
1917  if(gst_1) gst_1->Draw("pzi","pdgi==111","perrsame");
1918  c->cd(4);
1919  gst_0->Draw("Ei","pdgi==111","");
1920  if(gst_1) gst_1->Draw("Ei","pdgi==111","perrsame");
1921  c->cd();
1922  ls->Clear();
1923  ls->SetHeader("Primary Hadronic System: pi0's 4-momentum");
1924  ls->Draw();
1925  c->Update();
1926 
1927  //------ momentum of prim pi+
1928  ps->NewPage();
1929  c->Clear();
1930  c->Divide(2,2);
1931  c->cd(1);
1932  gst_0->Draw("pxi","pdgi==211","");
1933  if(gst_1) gst_1->Draw("pxi","pdgi==211","perrsame");
1934  c->cd(2);
1935  gst_0->Draw("pyi","pdgi==211","");
1936  if(gst_1) gst_1->Draw("pyi","pdgi==211","perrsame");
1937  c->cd(3);
1938  gst_0->Draw("pzi","pdgi==211","");
1939  if(gst_1) gst_1->Draw("pzi","pdgi==211","perrsame");
1940  c->cd(4);
1941  gst_0->Draw("Ei","pdgi==211","");
1942  if(gst_1) gst_1->Draw("Ei","pdgi==211","perrsame");
1943  c->cd();
1944  ls->Clear();
1945  ls->SetHeader("Primary Hadronic System:pi+'s 4-momentum");
1946  ls->Draw();
1947  c->Update();
1948 
1949  //------ momentum of prim. pi+
1950  ps->NewPage();
1951  c->Clear();
1952  c->Divide(2,2);
1953  c->cd(1);
1954  gst_0->Draw("pxi","pdgi==-211","");
1955  if(gst_1) gst_1->Draw("pxi","pdgi==-211","perrsame");
1956  c->cd(2);
1957  gst_0->Draw("pyi","pdgi==-211","");
1958  if(gst_1) gst_1->Draw("pyi","pdgi==-211","perrsame");
1959  c->cd(3);
1960  gst_0->Draw("pzi","pdgi==-211","");
1961  if(gst_1) gst_1->Draw("pzi","pdgi==-211","perrsame");
1962  c->cd(4);
1963  gst_0->Draw("Ei","pdgi==-211","");
1964  if(gst_1) gst_1->Draw("Ei","pdgi==-211","perrsame");
1965  c->cd();
1966  ls->Clear();
1967  ls->SetHeader("Primary Hadronic System: pi-'s 4-momentum");
1968  ls->Draw();
1969  c->Update();
1970  }//show?
1971 
1972  ps->Close();
1973 
1974  fin_0->Close();
1975  delete fin_0;
1976  if(fin_1) {
1977  fin_1->Close();
1978  delete fin_1;
1979  }
1980 }
1981 //_________________________________________________________________________________
1982 string OutputFileName(string inpname)
1983 {
1984 // Builds the output filename based on the name of the input filename
1985 // Perfors the following conversion: name.root -> name.nusample_test.ps
1986 
1987  unsigned int L = inpname.length();
1988 
1989  // if the last 4 characters are "root" (ROOT file extension) then
1990  // remove them
1991  if(inpname.substr(L-4, L).find("root") != string::npos) {
1992  inpname.erase(L-4, L);
1993  }
1994 
1995  ostringstream name;
1996  name << inpname << "sample_test.ps";
1997 
1998  return gSystem->BaseName(name.str().c_str());
1999 }
2000 //_________________________________________________________________________________
2001 void GetCommandLineArgs(int argc, char ** argv)
2002 {
2003  LOG("gevcomp", pNOTICE) << "*** Parsing command line arguments";
2004 
2005  CmdLnArgParser parser(argc,argv);
2006 
2007  // get GENIE summary ntuple
2008  if( parser.OptionExists('f') ) {
2009  LOG("gevcomp", pINFO) << "Reading filename for tested event sample";
2010  gOptInpFile = parser.ArgAsString('f');
2011  } else {
2012  LOG("gevcomp", pFATAL) << "Unspecified input filename - Exiting";
2013  PrintSyntax();
2014  exit(1);
2015  }
2016 
2017  // get another (reference) GENIE summary ntuple
2018  if( parser.OptionExists('r') ) {
2019  LOG("gevcomp", pINFO) << "Reading filename for reference event sample";
2020  gOptInpFileRef = parser.ArgAsString('r');
2021  } else {
2022  LOG("gevcomp", pNOTICE) << "Unspecified 'reference' event sample";
2023  }
2024 }
2025 //_________________________________________________________________________________
2026 void PrintSyntax(void)
2027 {
2028  LOG("gevcomp", pNOTICE)
2029  << "\n\n" << "Syntax:" << "\n"
2030  << " gevcomp -f sample.root [-n nev] [-r reference_sample.root]\n";
2031 }
2032 //_________________________________________________________________________________
2033 bool CheckRootFilename(string filename)
2034 {
2035  if(filename.size() == 0) return false;
2036 
2037  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));
2038  if (!is_accessible) {
2039  LOG("gevcomp", pERROR)
2040  << "The input ROOT file [" << filename << "] is not accessible";
2041  return false;
2042  }
2043  return true;
2044 }
2045 //_________________________________________________________________________________
2046 
pERROR
#define pERROR
Definition: Messenger.h:59
pFATAL
#define pFATAL
Definition: Messenger.h:56
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
gOptInpFile
string gOptInpFile
Definition: gEvComp.cxx:78
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
CreatePlots
void CreatePlots(string filename, string filename_ref)
Definition: gEvComp.cxx:92
pINFO
#define pINFO
Definition: Messenger.h:62
CheckRootFilename
bool CheckRootFilename(string filename)
Definition: gEvComp.cxx:2033
genie::units::ps
static constexpr double ps
Definition: Units.h:99
gOptInpFileRef
string gOptInpFileRef
Definition: gEvComp.cxx:79
OutputFileName
string OutputFileName(string input_file_name)
Definition: gEvComp.cxx:1982
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
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
PrintSyntax
void PrintSyntax(void)
Definition: gAtmoEvGen.cxx:1051
Generated on Wed Jul 17 2024 14:46:40 for GENIEGenerator by   doxygen 1.8.5