GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
INukeHadroData2018.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>, Rutherford Lab.
8  Steve Dytman <dytman+@pitt.edu>, Pittsburgh Univ.
9  Aaron Meyer <asm58@pitt.edu>, Pittsburgh Univ.
10  Alex Bell, Pittsburgh Univ.
11  February 01, 2007
12 
13  For the class documentation see the corresponding header file.
14 
15  Important revisions after version 2.0.0 :
16  @ Dec 06, 2008 - CA
17  Tweak dtor so as not to clutter the output if GENIE exits in err so as to
18  spot the fatal mesg immediately.
19  @ Jul 15, 2010 - AM
20  MeanFreePath now distinguishes between protons and neutrons. To account for
21  this, additional splines were added for each nucleon. Absorption fates
22  condensed into a single fate, splines updated accordingly. Added gamma and kaon
23  splines.Outdated splines were removed. Function IntBounce implemented to calculate
24  a CM scattering angle for given probe, target, product, and fate. AngleAndProduct
25  similar to IntBounce, but also determines the target nucleon.
26  @ May 01, 2012 - CA
27  Pick data from $GENIE/data/evgen/intranuke/
28  @ Jan 9, 2015 - SD, Nick Geary, Tomek Golan
29  Added 2014 version of INTRANUKE codes for independent development.
30  @ Oct, 2015 - TG
31  Added 2015 version of INTRANUKE codes for independent development.
32  Include Oset data files.
33  @ Apr, 2016 - Flor Blasczyk
34  Added K+ cex data files
35 */
36 //____________________________________________________________________________
37 
38 #include <cassert>
39 #include <string>
40 
41 #include <TSystem.h>
42 #include <TNtupleD.h>
43 #include <TGraph2D.h>
44 #include <TTree.h>
45 #include <TMath.h>
46 #include <TFile.h>
47 #include <iostream>
48 
56 
57 using std::ostringstream;
58 using std::ios;
59 
60 using namespace genie;
61 using namespace genie::constants;
62 
63 //____________________________________________________________________________
65 //____________________________________________________________________________
66 double INukeHadroData2018::fMinKinEnergy = 1.0; // MeV
67 double INukeHadroData2018::fMaxKinEnergyHA = 999.0; // MeV
68 double INukeHadroData2018::fMaxKinEnergyHN = 1799.0; // MeV
69 //____________________________________________________________________________
71 {
72  this->LoadCrossSections();
73  fInstance = 0;
74 }
75 //____________________________________________________________________________
77 {
78 
79  // pi+n/p hA x-section splines
80  delete fXSecPipn_Tot;
81  delete fXSecPipn_CEx;
82  delete fXSecPipn_Elas;
83  delete fXSecPipn_Reac;
84  delete fXSecPipp_Tot;
85  delete fXSecPipp_CEx;
86  delete fXSecPipp_Elas;
87  delete fXSecPipp_Reac;
88  delete fXSecPipd_Abs;
89 
90  delete fXSecPp_Cmp; //added to fix memory leak; no noticeable effect, but good convention.
91  delete fXSecPn_Cmp;
92  delete fXSecNn_Cmp;
93  delete fXSecPp_Tot;
94  delete fXSecPp_Elas;
95  delete fXSecPp_Reac;
96  delete fXSecPn_Tot;
97  delete fXSecPn_Elas;
98  delete fXSecPn_Reac;
99  delete fXSecNn_Tot;
100  delete fXSecNn_Elas;
101  delete fXSecNn_Reac;
102 
103  // pi0n/p hA x-section splines
104  delete fXSecPi0n_Tot;
105  delete fXSecPi0n_CEx;
106  delete fXSecPi0n_Elas;
107  delete fXSecPi0n_Reac;
108  delete fXSecPi0p_Tot;
109  delete fXSecPi0p_CEx;
110  delete fXSecPi0p_Elas;
111  delete fXSecPi0p_Reac;
112  delete fXSecPi0d_Abs;
113 
114  // K+N x-section splines
115  delete fXSecKpn_Elas;
116  delete fXSecKpp_Elas;
117  delete fXSecKpn_CEx;
118  delete fXSecKpN_Abs;
119  delete fXSecKpN_Tot;
120 
121  // gamma x-section splines (inelastic only)
122  delete fXSecGamp_fs;
123  delete fXSecGamn_fs;
124  delete fXSecGamN_Tot;
125 
126  // N+A x-section splines
127  delete fFracPA_Tot;
128  // delete fFracPA_Elas;
129  delete fFracPA_Inel;
130  delete fFracPA_CEx;
131  delete fFracPA_Abs;
132  delete fFracPA_PiPro;
133  delete fFracNA_Tot;
134  // delete fFracNA_Elas;
135  delete fFracNA_Inel;
136  delete fFracNA_CEx;
137  delete fFracNA_Abs;
138  delete fFracNA_PiPro;
139 
140  delete fFracPA_Cmp; // cmp - add support later
141  delete fFracNA_Cmp;
142 
143  // hN data
144  delete fhN2dXSecPP_Elas;
145  delete fhN2dXSecNP_Elas;
146  delete fhN2dXSecPipN_Elas;
147  delete fhN2dXSecPi0N_Elas;
148  delete fhN2dXSecPimN_Elas;
149  delete fhN2dXSecKpN_Elas;
150  delete fhN2dXSecKpP_Elas;
151  delete fhN2dXSecPiN_CEx;
152  delete fhN2dXSecPiN_Abs;
153  delete fhN2dXSecGamPi0P_Inelas;
154  delete fhN2dXSecGamPi0N_Inelas;
155  delete fhN2dXSecGamPipN_Inelas;
156  delete fhN2dXSecGamPimP_Inelas;
157  delete fhN2dXSecKpN_CEx;
158 
159 // delete TPipA_Tot;
160  delete TfracPipA_Abs;
161  delete TfracPipA_CEx;
162 // delete TfracPipA_Elas;
163  delete TfracPipA_Inelas;
164  delete TfracPipA_PiPro;
165 
166  // K+A x-section fraction splines
167  delete fFracKA_Tot;
168  delete fFracKA_Elas;
169  delete fFracKA_CEx;
170  delete fFracKA_Inel;
171  delete fFracKA_Abs;
172 
173 }
174 //____________________________________________________________________________
176 {
177  if(fInstance == 0) {
178  LOG("INukeData", pINFO) << "INukeHadroData2018 late initialization";
179  static INukeHadroData2018::Cleaner cleaner;
181  fInstance = new INukeHadroData2018;
182  }
183  return fInstance;
184 }
185 //____________________________________________________________________________
187 {
188 // Loads hadronic x-section data
189 
190  //-- Get the top-level directory with input hadron cross-section data
191  // (search for $GINUKEHADRONDATA or use default location)
192  string data_dir = (gSystem->Getenv("GINUKEHADRONDATA")) ?
193  string(gSystem->Getenv("GINUKEHADRONDATA")) :
194  string(gSystem->Getenv("GENIE")) + string("/data/evgen/intranuke");
195 
196  LOG("INukeData", pINFO)
197  << "Loading INTRANUKE hadron data from: " << data_dir;
198 
199  //-- Build filenames
200 
201  string datafile_NN = data_dir + "/tot_xsec/intranuke-xsections-NN2014.dat";
202  string datafile_pipN = data_dir + "/tot_xsec/intranuke-xsections-pi+N.dat";
203  string datafile_pi0N = data_dir + "/tot_xsec/intranuke-xsections-pi0N.dat";
204  string datafile_NA = data_dir + "/tot_xsec/intranuke-fractions-NA2016.dat";
205  string datafile_KA = data_dir + "/tot_xsec/intranuke-fractions-KA.dat";
206  string datafile_gamN = data_dir + "/tot_xsec/intranuke-xsections-gamN.dat";
207  string datafile_kN = data_dir + "/tot_xsec/intranuke-xsections-kaonN2018.dat";
208 
209  //-- Make sure that all data files are available
210 
211  assert( ! gSystem->AccessPathName(datafile_NN. c_str()) );
212  assert( ! gSystem->AccessPathName(datafile_pipN.c_str()) );
213  assert( ! gSystem->AccessPathName(datafile_pi0N.c_str()) );
214  assert( ! gSystem->AccessPathName(datafile_NA. c_str()) );
215  assert( ! gSystem->AccessPathName(datafile_KA. c_str()) );
216  assert( ! gSystem->AccessPathName(datafile_gamN.c_str()) );
217  assert( ! gSystem->AccessPathName(datafile_kN. c_str()) );
218 
219  LOG("INukeData", pINFO) << "Found all necessary data files...";
220 
221  //-- Load data files
222 
223  TTree data_NN;
224  TTree data_pipN;
225  TTree data_pi0N;
226  TTree data_NA;
227  TTree data_KA;
228  TTree data_gamN;
229  TTree data_kN;
230 
231  data_NN.ReadFile(datafile_NN.c_str(),"ke/D:pp_tot/D:pp_elas/D:pp_reac/D:pn_tot/D:pn_elas/D:pn_reac/D:nn_tot/D:nn_elas/D:nn_reac/D:pp_cmp/D:pn_cmp/D:nn_cmp/D");
232  data_pipN.ReadFile(datafile_pipN.c_str(),
233  "ke/D:pipn_tot/D:pipn_cex/D:pipn_elas/D:pipn_reac/D:pipp_tot/D:pipp_cex/D:pipp_elas/D:pipp_reac/D:pipd_abs");
234  data_pi0N.ReadFile(datafile_pi0N.c_str(),
235  "ke/D:pi0n_tot/D:pi0n_cex/D:pi0n_elas/D:pi0n_reac/D:pi0p_tot/D:pi0p_cex/D:pi0p_elas/D:pi0p_reac/D:pi0d_abs");
236  //data_NA.ReadFile(datafile_NA.c_str(),
237  //"ke/D:pA_tot/D:pA_elas/D:pA_inel/D:pA_cex/D:pA_abs/D:pA_pipro/D"); // add support for cmp here (?)
238  data_NA.ReadFile(datafile_NA.c_str(),
239  "ke/D:pA_tot/D:pA_inel/D:pA_cex/D:pA_abs/D:pA_pipro/D:pA_cmp/D"); // add support for cmp here (?)
240  data_gamN.ReadFile(datafile_gamN.c_str(),
241  "ke/D:pi0p_tot/D:pipn_tot/D:pimp_tot/D:pi0n_tot/D:gamp_fs/D:gamn_fs/D:gamN_tot/D");
242  data_kN.ReadFile(datafile_kN.c_str(),
243  "ke/D:kpp_elas/D:kpn_elas/D:kpn_cex/D:kp_abs/D:kpN_tot/D");
244  data_KA.ReadFile(datafile_KA.c_str(),
245  "ke/D:KA_tot/D:KA_elas/D:KA_inel/D:KA_abs/D");
246 
247  LOG("INukeData", pDEBUG) << "Number of data rows in NN : " << data_NN.GetEntries();
248  LOG("INukeData", pDEBUG) << "Number of data rows in pipN : " << data_pipN.GetEntries();
249  LOG("INukeData", pDEBUG) << "Number of data rows in pi0N : " << data_pi0N.GetEntries();
250  LOG("INukeData", pDEBUG) << "Number of data rows in NA : " << data_NA.GetEntries();
251  LOG("INukeData", pDEBUG) << "Number of data rows in KA : " << data_KA.GetEntries();
252  LOG("INukeData", pDEBUG) << "Number of data rows in gamN : " << data_gamN.GetEntries();
253  LOG("INukeData", pDEBUG) << "Number of data rows in kN : " << data_kN.GetEntries();
254 
255  LOG("INukeData", pINFO) << "Done loading all x-section files...";
256 
257  //-- Build x-section splines
258 
259  // p/n+p/n hA x-section splines
260  fXSecPp_Tot = new Spline(&data_NN, "ke:pp_tot");
261  fXSecPp_Elas = new Spline(&data_NN, "ke:pp_elas");
262  fXSecPp_Reac = new Spline(&data_NN, "ke:pp_reac");
263  fXSecPn_Tot = new Spline(&data_NN, "ke:pn_tot");
264  fXSecPn_Elas = new Spline(&data_NN, "ke:pn_elas");
265  fXSecPn_Reac = new Spline(&data_NN, "ke:pn_reac");
266  fXSecNn_Tot = new Spline(&data_NN, "ke:nn_tot");
267  fXSecNn_Elas = new Spline(&data_NN, "ke:nn_elas");
268  fXSecNn_Reac = new Spline(&data_NN, "ke:nn_reac");
269  fXSecPp_Cmp = new Spline(&data_NN, "ke:pp_cmp"); //for compound nucleus fate
270  fXSecPn_Cmp = new Spline(&data_NN, "ke:pn_cmp");
271  fXSecNn_Cmp = new Spline(&data_NN, "ke:nn_cmp");
272 
273  // pi+n/p hA x-section splines
274  fXSecPipn_Tot = new Spline(&data_pipN, "ke:pipn_tot");
275  fXSecPipn_CEx = new Spline(&data_pipN, "ke:pipn_cex");
276  fXSecPipn_Elas = new Spline(&data_pipN, "ke:pipn_elas");
277  fXSecPipn_Reac = new Spline(&data_pipN, "ke:pipn_reac");
278  fXSecPipp_Tot = new Spline(&data_pipN, "ke:pipp_tot");
279  fXSecPipp_CEx = new Spline(&data_pipN, "ke:pipp_cex");
280  fXSecPipp_Elas = new Spline(&data_pipN, "ke:pipp_elas");
281  fXSecPipp_Reac = new Spline(&data_pipN, "ke:pipp_reac");
282  fXSecPipd_Abs = new Spline(&data_pipN, "ke:pipd_abs");
283 
284  // pi0n/p hA x-section splines
285  fXSecPi0n_Tot = new Spline(&data_pi0N, "ke:pi0n_tot");
286  fXSecPi0n_CEx = new Spline(&data_pi0N, "ke:pi0n_cex");
287  fXSecPi0n_Elas = new Spline(&data_pi0N, "ke:pi0n_elas");
288  fXSecPi0n_Reac = new Spline(&data_pi0N, "ke:pi0n_reac");
289  fXSecPi0p_Tot = new Spline(&data_pi0N, "ke:pi0p_tot");
290  fXSecPi0p_CEx = new Spline(&data_pi0N, "ke:pi0p_cex");
291  fXSecPi0p_Elas = new Spline(&data_pi0N, "ke:pi0p_elas");
292  fXSecPi0p_Reac = new Spline(&data_pi0N, "ke:pi0p_reac");
293  fXSecPi0d_Abs = new Spline(&data_pi0N, "ke:pi0d_abs");
294 
295  // K+N x-section splines
296  fXSecKpn_Elas = new Spline(&data_kN, "ke:kpn_elas");
297  fXSecKpp_Elas = new Spline(&data_kN, "ke:kpp_elas");
298  fXSecKpn_CEx = new Spline(&data_kN, "ke:kpn_cex");
299  fXSecKpN_Abs = 0; // new Spline(&data_kN, "ke:kp_abs"); why not used?
300  fXSecKpN_Tot = new Spline(&data_kN, "ke:kpN_tot");
301 
302  // gamma x-section splines
303  fXSecGamp_fs = new Spline(&data_gamN, "ke:gamp_fs");
304  fXSecGamn_fs = new Spline(&data_gamN, "ke:gamn_fs");
305  fXSecGamN_Tot = new Spline(&data_gamN, "ke:gamN_tot");
306 
307  // N+A x-section fraction splines
308  fFracPA_Tot = new Spline(&data_NA, "ke:pA_tot");
309  // fFracPA_Elas = new Spline(&data_NA, "ke:pA_elas");
310  fFracPA_Inel = new Spline(&data_NA, "ke:pA_inel");
311  fFracPA_CEx = new Spline(&data_NA, "ke:pA_cex");
312  fFracPA_Abs = new Spline(&data_NA, "ke:pA_abs");
313  fFracPA_PiPro = new Spline(&data_NA, "ke:pA_pipro");
314  fFracNA_Tot = new Spline(&data_NA, "ke:pA_tot"); // assuming nA same as pA
315  // fFracNA_Elas = new Spline(&data_NA, "ke:pA_elas");
316  fFracNA_Inel = new Spline(&data_NA, "ke:pA_inel");
317  fFracNA_CEx = new Spline(&data_NA, "ke:pA_cex");
318  fFracNA_Abs = new Spline(&data_NA, "ke:pA_abs");
319  fFracNA_PiPro = new Spline(&data_NA, "ke:pA_pipro");
320 
321  fFracPA_Cmp = new Spline(&data_NA, "ke:pA_cmp"); //cmp - add support later
322  fFracNA_Cmp = new Spline(&data_NA, "ke:pA_cmp");
323 
324  // K+A x-section fraction splines
325  fFracKA_Tot = new Spline(&data_KA, "ke:KA_tot");
326  fFracKA_Elas = new Spline(&data_KA, "ke:KA_elas");
327  fFracKA_CEx = 0; // new Spline(&data_KA, "ke:KA_cex"); //Added, but needs to be computed
328  fFracKA_Inel = new Spline(&data_KA, "ke:KA_inel");
329  fFracKA_Abs = new Spline(&data_KA, "ke:KA_abs");
330  //
331  // hN stuff
332  //
333 
334 
335  // kIHNFtElas
336  // pp, nn --> read from pp/pp%.txt
337  // pn, np --> read from pp/pn%.txt
338  // pi+ N --> read from pip/pip%.txt
339  // pi0 N --> read from pip/pip%.txt
340  // pi- N --> read from pim/pim%.txt
341  // K+ N --> read from kpn/kpn%.txt
342  // K+ P --> read from kpp/kpp%.txt
343  // kIHNFtCEx
344  // pi+, pi0, pi- --> read from pie/pie%.txt (using pip+n->pi0+p data)
345  // kIHNFtAbs
346  // pi+, pi0, pi- --> read from pid2p/pid2p%.txt (using pip+D->2p data)
347  // kIHNFtInelas
348  // gamma p -> p pi0 --> read from gampi0p/%-pi0p.txt
349  // gamma p -> n pi+ --> read from gampi+n/%-pi+n.txt
350  // gamma n -> n pi0 --> read from gampi0n/%-pi0n.txt
351  // gamma n -> p pi- --> read from gampi-p/%-pi-p.txt
352 
353 
354  // kIHNFtElas, pp&nn :
355  {
356  const int hN_ppelas_nfiles = 20;
357  const int hN_ppelas_points_per_file = 21;
358  const int hN_ppelas_npoints = hN_ppelas_points_per_file * hN_ppelas_nfiles;
359 
360  double hN_ppelas_energies[hN_ppelas_nfiles] = {
361  50, 100, 150, 200, 250, 300, 350, 400, 450, 500,
362  550, 600, 650, 700, 750, 800, 850, 900, 950, 1000
363  };
364 
365  double hN_ppelas_costh [hN_ppelas_points_per_file];
366  double hN_ppelas_xsec [hN_ppelas_npoints];
367 
368  int ipoint=0;
369 
370  for(int ifile = 0; ifile < hN_ppelas_nfiles; ifile++) {
371  // build filename
372  ostringstream hN_datafile;
373  double ke = hN_ppelas_energies[ifile];
374  hN_datafile << data_dir << "/diff_ang/pp/pp" << ke << ".txt";
375  // read data
376  ReadhNFile(
377  hN_datafile.str(), ke, hN_ppelas_points_per_file,
378  ipoint, hN_ppelas_costh, hN_ppelas_xsec,2);
379  }//loop over files
380 
381  fhN2dXSecPP_Elas = new BLI2DNonUnifGrid(hN_ppelas_nfiles,hN_ppelas_points_per_file,
382  hN_ppelas_energies,hN_ppelas_costh,hN_ppelas_xsec);
383  }
384 
385  // kIHNFtElas, pn&np :
386  {
387  const int hN_npelas_nfiles = 20;
388  const int hN_npelas_points_per_file = 21;
389  const int hN_npelas_npoints = hN_npelas_points_per_file * hN_npelas_nfiles;
390 
391  double hN_npelas_energies[hN_npelas_nfiles] = {
392  50, 100, 150, 200, 250, 300, 350, 400, 450, 500,
393  550, 600, 650, 700, 750, 800, 850, 900, 950, 1000
394  };
395 
396  double hN_npelas_costh [hN_npelas_points_per_file];
397  double hN_npelas_xsec [hN_npelas_npoints];
398 
399  int ipoint=0;
400 
401  for(int ifile = 0; ifile < hN_npelas_nfiles; ifile++) {
402  // build filename
403  ostringstream hN_datafile;
404  double ke = hN_npelas_energies[ifile];
405  hN_datafile << data_dir << "/diff_ang/pn/pn" << ke << ".txt";
406  // read data
407  ReadhNFile(
408  hN_datafile.str(), ke, hN_npelas_points_per_file,
409  ipoint, hN_npelas_costh, hN_npelas_xsec,2);
410  }//loop over files
411 
412  fhN2dXSecNP_Elas = new BLI2DNonUnifGrid(hN_npelas_nfiles,hN_npelas_points_per_file,
413  hN_npelas_energies,hN_npelas_costh,hN_npelas_xsec);
414  }
415 
416  // kIHNFtElas, pipN :
417  {
418  const int hN_pipNelas_nfiles = 60;
419  const int hN_pipNelas_points_per_file = 21;
420  const int hN_pipNelas_npoints = hN_pipNelas_points_per_file * hN_pipNelas_nfiles;
421 
422  double hN_pipNelas_energies[hN_pipNelas_nfiles] = {
423  10, 20, 30, 40, 50, 60, 70, 80, 90,
424  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
425  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
426  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
427  700, 740, 780, 820, 860, 900, 940, 980,
428  1020, 1060, 1100, 1140, 1180, 1220, 1260,
429  1300, 1340, 1380, 1420, 1460, 1500
430  };
431 
432  double hN_pipNelas_costh [hN_pipNelas_points_per_file];
433  double hN_pipNelas_xsec [hN_pipNelas_npoints];
434 
435  int ipoint=0;
436 
437  for(int ifile = 0; ifile < hN_pipNelas_nfiles; ifile++) {
438  // build filename
439  ostringstream hN_datafile;
440  double ke = hN_pipNelas_energies[ifile];
441  hN_datafile << data_dir << "/diff_ang/pip/pip" << ke << ".txt";
442  // read data
443  ReadhNFile(
444  hN_datafile.str(), ke, hN_pipNelas_points_per_file,
445  ipoint, hN_pipNelas_costh, hN_pipNelas_xsec,2);
446  }//loop over files
447 
448  fhN2dXSecPipN_Elas = new BLI2DNonUnifGrid(hN_pipNelas_nfiles,hN_pipNelas_points_per_file,
449  hN_pipNelas_energies,hN_pipNelas_costh,hN_pipNelas_xsec);
450  }
451 
452  // kIHNFtElas, pi0N :
453  {
454  const int hN_pi0Nelas_nfiles = 60;
455  const int hN_pi0Nelas_points_per_file = 21;
456  const int hN_pi0Nelas_npoints = hN_pi0Nelas_points_per_file * hN_pi0Nelas_nfiles;
457 
458  double hN_pi0Nelas_energies[hN_pi0Nelas_nfiles] = {
459  10, 20, 30, 40, 50, 60, 70, 80, 90,
460  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
461  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
462  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
463  700, 740, 780, 820, 860, 900, 940, 980,
464  1020, 1060, 1100, 1140, 1180, 1220, 1260,
465  1300, 1340, 1380, 1420, 1460, 1500
466  };
467 
468  double hN_pi0Nelas_costh [hN_pi0Nelas_points_per_file];
469  double hN_pi0Nelas_xsec [hN_pi0Nelas_npoints];
470 
471  int ipoint=0;
472 
473  for(int ifile = 0; ifile < hN_pi0Nelas_nfiles; ifile++) {
474  // build filename
475  ostringstream hN_datafile;
476  double ke = hN_pi0Nelas_energies[ifile];
477  hN_datafile << data_dir << "/diff_ang/pip/pip" << ke << ".txt";
478  // read data
479  ReadhNFile(
480  hN_datafile.str(), ke, hN_pi0Nelas_points_per_file,
481  ipoint, hN_pi0Nelas_costh, hN_pi0Nelas_xsec,2);
482  }//loop over files
483 
484  fhN2dXSecPi0N_Elas = new BLI2DNonUnifGrid(hN_pi0Nelas_nfiles,hN_pi0Nelas_points_per_file,
485  hN_pi0Nelas_energies,hN_pi0Nelas_costh,hN_pi0Nelas_xsec);
486  }
487 
488  // kIHNFtElas, pimN :
489  {
490  const int hN_pimNelas_nfiles = 60;
491  const int hN_pimNelas_points_per_file = 21;
492  const int hN_pimNelas_npoints = hN_pimNelas_points_per_file * hN_pimNelas_nfiles;
493 
494  double hN_pimNelas_energies[hN_pimNelas_nfiles] = {
495  10, 20, 30, 40, 50, 60, 70, 80, 90,
496  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
497  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
498  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
499  700, 740, 780, 820, 860, 900, 940, 980,
500  1020, 1060, 1100, 1140, 1180, 1220, 1260,
501  1300, 1340, 1380, 1420, 1460, 1500
502  };
503 
504  double hN_pimNelas_costh [hN_pimNelas_points_per_file];
505  double hN_pimNelas_xsec [hN_pimNelas_npoints];
506 
507  int ipoint=0;
508 
509  for(int ifile = 0; ifile < hN_pimNelas_nfiles; ifile++) {
510  // build filename
511  ostringstream hN_datafile;
512  double ke = hN_pimNelas_energies[ifile];
513  hN_datafile << data_dir << "/diff_ang/pim/pim" << ke << ".txt";
514  // read data
515  ReadhNFile(
516  hN_datafile.str(), ke, hN_pimNelas_points_per_file,
517  ipoint, hN_pimNelas_costh, hN_pimNelas_xsec,2);
518  }//loop over files
519 
520  fhN2dXSecPimN_Elas = new BLI2DNonUnifGrid(hN_pimNelas_nfiles,hN_pimNelas_points_per_file,
521  hN_pimNelas_energies,hN_pimNelas_costh,hN_pimNelas_xsec);
522  }
523 
524  // kIHNFtElas, kpn :
525  {
526  const int hN_kpNelas_nfiles = 18;
527  const int hN_kpNelas_points_per_file = 37;
528  const int hN_kpNelas_npoints = hN_kpNelas_points_per_file * hN_kpNelas_nfiles;
529 
530  double hN_kpNelas_energies[hN_kpNelas_nfiles] = {
531  100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
532  1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800
533  };
534 
535  double hN_kpNelas_costh [hN_kpNelas_points_per_file];
536  double hN_kpNelas_xsec [hN_kpNelas_npoints];
537 
538  int ipoint=0;
539 
540  for(int ifile = 0; ifile < hN_kpNelas_nfiles; ifile++) {
541  // build filename
542  ostringstream hN_datafile;
543  double ke = hN_kpNelas_energies[ifile];
544  hN_datafile << data_dir << "/diff_ang/kpn/kpn" << ke << ".txt";
545  // read data
546  ReadhNFile(
547  hN_datafile.str(), ke, hN_kpNelas_points_per_file,
548  ipoint, hN_kpNelas_costh, hN_kpNelas_xsec,2);
549  }//loop over files
550 
551  fhN2dXSecKpN_Elas = new BLI2DNonUnifGrid(hN_kpNelas_nfiles,hN_kpNelas_points_per_file,
552  hN_kpNelas_energies,hN_kpNelas_costh,hN_kpNelas_xsec);
553  }
554  // kIHNFtCEx, kpn :
555  {
556  const int hN_kpNcex_nfiles = 18;
557  const int hN_kpNcex_points_per_file = 37;
558  const int hN_kpNcex_npoints = hN_kpNcex_points_per_file * hN_kpNcex_nfiles;
559 
560  double hN_kpNcex_energies[hN_kpNcex_nfiles] = {
561  100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
562  1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800
563  };
564 
565  double hN_kpNcex_costh [hN_kpNcex_points_per_file];
566  double hN_kpNcex_xsec [hN_kpNcex_npoints];
567 
568  int ipoint=0;
569 
570  for(int ifile = 0; ifile < hN_kpNcex_nfiles; ifile++) {
571  // build filename
572  ostringstream hN_datafile;
573  double ke = hN_kpNcex_energies[ifile];
574  hN_datafile << data_dir << "/diff_ang/kpncex/kpcex" << ke << ".txt";
575  // read data
576  ReadhNFile(
577  hN_datafile.str(), ke, hN_kpNcex_points_per_file,
578  ipoint, hN_kpNcex_costh, hN_kpNcex_xsec,2);
579  }//loop over files
580 
581  /*double hN_kpNcex_costh_cond [hN_kpNcex_points_per_file];
582  for (int ient = 0; ient < hN_kpNcex_points_per_file; ient++) {
583  hN_kpNcex_costh_cond[ient] = hN_kpNcex_costh[ient];
584  }*/
585 
586  fhN2dXSecKpN_CEx = new BLI2DNonUnifGrid(hN_kpNcex_nfiles,hN_kpNcex_points_per_file,
587  hN_kpNcex_energies,hN_kpNcex_costh,hN_kpNcex_xsec);
588  }
589 //----------------------------------------------------------------------------------------
590 
591 
592  // kIHNFtElas, kpp :
593  {
594  const int hN_kpPelas_nfiles = 18;
595  const int hN_kpPelas_points_per_file = 37;
596  const int hN_kpPelas_npoints = hN_kpPelas_points_per_file * hN_kpPelas_nfiles;
597 
598  double hN_kpPelas_energies[hN_kpPelas_nfiles] = {
599  100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
600  1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800
601  };
602 
603  double hN_kpPelas_costh [hN_kpPelas_points_per_file];
604  double hN_kpPelas_xsec [hN_kpPelas_npoints];
605 
606  int ipoint=0;
607 
608  for(int ifile = 0; ifile < hN_kpPelas_nfiles; ifile++) {
609  // build filename
610  ostringstream hN_datafile;
611  double ke = hN_kpPelas_energies[ifile];
612  hN_datafile << data_dir << "/diff_ang/kpp/kpp" << ke << ".txt";
613  // read data
614  ReadhNFile(
615  hN_datafile.str(), ke, hN_kpPelas_points_per_file,
616  ipoint, hN_kpPelas_costh, hN_kpPelas_xsec,2);
617  }//loop over files
618 
619  fhN2dXSecKpP_Elas = new BLI2DNonUnifGrid(hN_kpPelas_nfiles,hN_kpPelas_points_per_file,
620  hN_kpPelas_energies,hN_kpPelas_costh,hN_kpPelas_xsec);
621  }
622 
623  // kIHNFtCEx, (pi+, pi0, pi-) N
624  {
625  const int hN_piNcex_nfiles = 60;
626  const int hN_piNcex_points_per_file = 21;
627  const int hN_piNcex_npoints = hN_piNcex_points_per_file * hN_piNcex_nfiles;
628 
629  double hN_piNcex_energies[hN_piNcex_nfiles] = {
630  10, 20, 30, 40, 50, 60, 70, 80, 90,
631  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
632  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
633  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
634  700, 740, 780, 820, 860, 900, 940, 980,
635  1020, 1060, 1100, 1140, 1180, 1220, 1260,
636  1300, 1340, 1380, 1420, 1460, 1500
637  };
638 
639  double hN_piNcex_costh [hN_piNcex_points_per_file];
640  double hN_piNcex_xsec [hN_piNcex_npoints];
641 
642  int ipoint=0;
643 
644  for(int ifile = 0; ifile < hN_piNcex_nfiles; ifile++) {
645  // build filename
646  ostringstream hN_datafile;
647  double ke = hN_piNcex_energies[ifile];
648  hN_datafile << data_dir << "/diff_ang/pie/pie" << ke << ".txt";
649  // read data
650  ReadhNFile(
651  hN_datafile.str(), ke, hN_piNcex_points_per_file,
652  ipoint, hN_piNcex_costh, hN_piNcex_xsec,2);
653  }//loop over files
654 
655  fhN2dXSecPiN_CEx = new BLI2DNonUnifGrid(hN_piNcex_nfiles,hN_piNcex_points_per_file,
656  hN_piNcex_energies,hN_piNcex_costh,hN_piNcex_xsec);
657  }
658 
659  // kIHNFtAbs, (pi+, pi0, pi-) N
660  {
661  const int hN_piNabs_nfiles = 19;
662  const int hN_piNabs_points_per_file = 21;
663  const int hN_piNabs_npoints = hN_piNabs_points_per_file * hN_piNabs_nfiles;
664 
665  double hN_piNabs_energies[hN_piNabs_nfiles] = {
666  50, 75, 100, 125, 150, 175, 200, 225, 250, 275,
667  300, 325, 350, 375, 400, 425, 450, 475, 500
668  };
669 
670  double hN_piNabs_costh [hN_piNabs_points_per_file];
671  double hN_piNabs_xsec [hN_piNabs_npoints];
672 
673  int ipoint=0;
674 
675  for(int ifile = 0; ifile < hN_piNabs_nfiles; ifile++) {
676  // build filename
677  ostringstream hN_datafile;
678  double ke = hN_piNabs_energies[ifile];
679  hN_datafile << data_dir << "/diff_ang/pid2p/pid2p" << ke << ".txt";
680  // read data
681  ReadhNFile(
682  hN_datafile.str(), ke, hN_piNabs_points_per_file,
683  ipoint, hN_piNabs_costh, hN_piNabs_xsec,2);
684  }//loop over files
685 
686  fhN2dXSecPiN_Abs = new BLI2DNonUnifGrid(hN_piNabs_nfiles,hN_piNabs_points_per_file,
687  hN_piNabs_energies,hN_piNabs_costh,hN_piNabs_xsec);
688  }
689 
690  // kIHNFtInelas, gamma p -> p pi0
691  {
692  const int hN_gampi0pInelas_nfiles = 29;
693  const int hN_gampi0pInelas_points_per_file = 37;
694  const int hN_gampi0pInelas_npoints = hN_gampi0pInelas_points_per_file * hN_gampi0pInelas_nfiles;
695 
696  double hN_gampi0pInelas_energies[hN_gampi0pInelas_nfiles] = {
697  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
698  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
699  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
700  };
701 
702  double hN_gampi0pInelas_costh [hN_gampi0pInelas_points_per_file];
703  double hN_gampi0pInelas_xsec [hN_gampi0pInelas_npoints];
704 
705  int ipoint=0;
706 
707  for(int ifile = 0; ifile < hN_gampi0pInelas_nfiles; ifile++) {
708  // build filename
709  ostringstream hN_datafile;
710  double ke = hN_gampi0pInelas_energies[ifile];
711  hN_datafile << data_dir << "/diff_ang/gampi0p/" << ke << "-pi0p.txt";
712  // read data
713  ReadhNFile(
714  hN_datafile.str(), ke, hN_gampi0pInelas_points_per_file,
715  ipoint, hN_gampi0pInelas_costh, hN_gampi0pInelas_xsec,3);
716  }//loop over files
717 
718  fhN2dXSecGamPi0P_Inelas = new BLI2DNonUnifGrid(hN_gampi0pInelas_nfiles,hN_gampi0pInelas_points_per_file,
719  hN_gampi0pInelas_energies,hN_gampi0pInelas_costh,hN_gampi0pInelas_xsec);
720  }
721 
722  // kIHNFtInelas, gamma n -> n pi0
723  {
724  const int hN_gampi0nInelas_nfiles = 29;
725  const int hN_gampi0nInelas_points_per_file = 37;
726  const int hN_gampi0nInelas_npoints = hN_gampi0nInelas_points_per_file * hN_gampi0nInelas_nfiles;
727 
728  double hN_gampi0nInelas_energies[hN_gampi0nInelas_nfiles] = {
729  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
730  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
731  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
732  };
733 
734  double hN_gampi0nInelas_costh [hN_gampi0nInelas_points_per_file];
735  double hN_gampi0nInelas_xsec [hN_gampi0nInelas_npoints];
736  int ipoint=0;
737 
738  for(int ifile = 0; ifile < hN_gampi0nInelas_nfiles; ifile++) {
739  // build filename
740  ostringstream hN_datafile;
741  double ke = hN_gampi0nInelas_energies[ifile];
742  hN_datafile << data_dir << "/diff_ang/gampi0n/" << ke << "-pi0n.txt";
743  // read data
744  ReadhNFile(
745  hN_datafile.str(), ke, hN_gampi0nInelas_points_per_file,
746  ipoint, hN_gampi0nInelas_costh, hN_gampi0nInelas_xsec,3);
747  }//loop over files
748 
749  fhN2dXSecGamPi0N_Inelas = new BLI2DNonUnifGrid(hN_gampi0nInelas_nfiles,hN_gampi0nInelas_points_per_file,
750  hN_gampi0nInelas_energies,hN_gampi0nInelas_costh,hN_gampi0nInelas_xsec);
751  }
752 
753  // kIHNFtInelas, gamma p -> n pi+
754  {
755  const int hN_gampipnInelas_nfiles = 29;
756  const int hN_gampipnInelas_points_per_file = 37;
757  const int hN_gampipnInelas_npoints = hN_gampipnInelas_points_per_file * hN_gampipnInelas_nfiles;
758 
759  double hN_gampipnInelas_energies[hN_gampipnInelas_nfiles] = {
760  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
761  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
762  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
763  };
764 
765  double hN_gampipnInelas_costh [hN_gampipnInelas_points_per_file];
766  double hN_gampipnInelas_xsec [hN_gampipnInelas_npoints];
767 
768  int ipoint=0;
769 
770  for(int ifile = 0; ifile < hN_gampipnInelas_nfiles; ifile++) {
771  // build filename
772  ostringstream hN_datafile;
773  double ke = hN_gampipnInelas_energies[ifile];
774  hN_datafile << data_dir << "/diff_ang/gampi+n/" << ke << "-pi+n.txt";
775  // read data
776  ReadhNFile(
777  hN_datafile.str(), ke, hN_gampipnInelas_points_per_file,
778  ipoint, hN_gampipnInelas_costh, hN_gampipnInelas_xsec,3);
779  }//loop over files
780 
781  fhN2dXSecGamPipN_Inelas = new BLI2DNonUnifGrid(hN_gampipnInelas_nfiles,hN_gampipnInelas_points_per_file,
782  hN_gampipnInelas_energies,hN_gampipnInelas_costh,hN_gampipnInelas_xsec);
783  }
784 
785  // kIHNFtInelas, gamma n -> p pi-
786  {
787  const int hN_gampimpInelas_nfiles = 29;
788  const int hN_gampimpInelas_points_per_file = 37;
789  const int hN_gampimpInelas_npoints = hN_gampimpInelas_points_per_file * hN_gampimpInelas_nfiles;
790 
791  double hN_gampimpInelas_energies[hN_gampimpInelas_nfiles] = {
792  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
793  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
794  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
795  };
796 
797  double hN_gampimpInelas_costh [hN_gampimpInelas_points_per_file];
798  double hN_gampimpInelas_xsec [hN_gampimpInelas_npoints];
799 
800  int ipoint=0;
801 
802  for(int ifile = 0; ifile < hN_gampimpInelas_nfiles; ifile++) {
803  // build filename
804  ostringstream hN_datafile;
805  double ke = hN_gampimpInelas_energies[ifile];
806  hN_datafile << data_dir << "/diff_ang/gampi-p/" << ke << "-pi-p.txt";
807  // read data
808  ReadhNFile(
809  hN_datafile.str(), ke, hN_gampimpInelas_points_per_file,
810  ipoint, hN_gampimpInelas_costh, hN_gampimpInelas_xsec,3);
811  }//loop over files
812 
813  fhN2dXSecGamPimP_Inelas = new BLI2DNonUnifGrid(hN_gampimpInelas_nfiles,hN_gampimpInelas_points_per_file,
814  hN_gampimpInelas_energies,hN_gampimpInelas_costh,hN_gampimpInelas_xsec);
815  }
816 
817 
818 
819 
820  TFile TGraphs_file;
821  bool saveTGraphsToFile = false; //true;
822 
823  if (saveTGraphsToFile) {
824  string filename = "TGraphs.root";
825  LOG("INukeHadroData2018", pNOTICE) << "Saving INTRANUKE hadron x-section data to ROOT file: " << filename;
826  TGraphs_file.Open(filename.c_str(), "RECREATE");
827  }
828 
829  /*
830  // kIHNFtTot, pip + A PipA_Tot
831  {
832  const int pipATot_nfiles = 22;
833  const int pipATot_nuclei[pipATot_nfiles] = {1, 2, 3, 4, 6, 7, 9, 12, 16, 27, 28, 32, 40, 48, 56, 58, 63, 93, 120, 165, 181, 209};
834  const int pipATot_npoints = 203;
835 
836  TPipA_Tot = new TGraph2D(pipATot_npoints);
837  TPipA_Tot->SetNameTitle("TPipA_Tot","TPipA_Tot");
838  TPipA_Tot->SetDirectory(0);
839 
840  int ipoint=0;
841  double x, y;
842 
843  for(int ifile=0; ifile < pipATot_nfiles; ifile++) {
844  ostringstream ADep_datafile;
845  int nucleus = pipATot_nuclei[ifile];
846  ADep_datafile << data_dir << "/tot_xsec/pipA_tot/pip" << nucleus << "_tot.txt";
847  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
848  buff->SetNameTitle("buff","buff");
849  for(int i=0; i < buff->GetN(); i++) {
850  buff -> GetPoint(i,x,y);
851  TPipA_Tot -> SetPoint(ipoint,(double)nucleus,x,y);
852  ipoint++;
853  }
854  delete buff;
855  }
856 
857  if (saveTGraphsToFile) {
858  TPipA_Tot -> Write("TPipA_Tot"); // TPipA_Tot will be _key_ name
859  }
860  }
861  */
862 
863  // kIHNFtAbs, pip + A PipA_Abs_frac
864  {
865  const int pipAAbs_f_nfiles = 18;
866  const int pipAAbs_f_nuclei[pipAAbs_f_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 48, 56, 58, 63, 93, 120, 165, 181, 209};
867  const int pipAAbs_f_npoints = 111;
868 
869  TfracPipA_Abs = new TGraph2D(pipAAbs_f_npoints);
870  TfracPipA_Abs->SetNameTitle("TfracPipA_Abs","TfracPipA_Abs");
871  TfracPipA_Abs->SetDirectory(0);
872 
873  int ipoint=0;
874  double x, y;
875  for(int ifile=0; ifile < pipAAbs_f_nfiles; ifile++) {
876  ostringstream ADep_datafile;
877  int nucleus = pipAAbs_f_nuclei[ifile];
878  ADep_datafile << data_dir << "/tot_xsec/pipA_abs_frac/pip" << nucleus << "_abs_frac.txt";
879  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
880  buff->SetNameTitle("buff","buff");
881  for(int i=0; i < buff->GetN(); i++) {
882  buff -> GetPoint(i,x,y);
883  TfracPipA_Abs -> SetPoint(ipoint,(double)nucleus,x,y);
884  ipoint++;
885  }
886  delete buff;
887  }
888  if (saveTGraphsToFile) {
889  TfracPipA_Abs -> Write("TfracPipA_Abs");
890  }
891 
892  }
893 
894 
895  // kIHNFtCEx, pip + A PipA_CEx_frac
896  {
897  const int pipACEx_f_nfiles = 18;
898  const int pipACEx_f_nuclei[pipACEx_f_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 48, 56, 58, 63, 93, 120, 165, 181, 209};
899  const int pipACEx_f_npoints = 129;
900 
901  TfracPipA_CEx = new TGraph2D(pipACEx_f_npoints);
902  TfracPipA_CEx->SetNameTitle("TfracPipA_CEx","TfracPipA_CEx");
903  TfracPipA_CEx->SetDirectory(0);
904 
905  int ipoint=0;
906  double x, y;
907 
908  for(int ifile=0; ifile < pipACEx_f_nfiles; ifile++) {
909  ostringstream ADep_datafile;
910  int nucleus = pipACEx_f_nuclei[ifile];
911  ADep_datafile << data_dir << "/tot_xsec/pipA_cex_frac/pip" << nucleus << "_cex_frac.txt";
912  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
913  buff->SetNameTitle("buff","buff");
914  for(int i=0; i < buff->GetN(); i++) {
915  buff -> GetPoint(i,x,y);
916  TfracPipA_CEx -> SetPoint(ipoint,(double)nucleus,x,y);
917  ipoint++;
918  }
919  delete buff;
920  }
921 
922  if (saveTGraphsToFile) {
923  TfracPipA_CEx -> Write("TfracPipA_CEx");
924  }
925 
926  }
927 
928 
929 
930  // kIHNFtCEx, pip + A PipA_CEx (just for developmental purposes)
931  {
932  TGraph2D * TPipA_CEx;
933 
934  const int pipACEx_nfiles = 18;
935  const int pipACEx_nuclei[pipACEx_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 48, 56, 58, 63, 93, 120, 165, 181, 209};
936  const int pipACEx_npoints = 129;
937 
938  TPipA_CEx = new TGraph2D(pipACEx_npoints);
939  TPipA_CEx->SetNameTitle("TPipA_CEx","TPipA_CEx");
940  TPipA_CEx->SetDirectory(0);
941 
942  int ipoint=0;
943  double x, y;
944 
945  for(int ifile=0; ifile < pipACEx_nfiles; ifile++) {
946  ostringstream ADep_datafile;
947  int nucleus = pipACEx_nuclei[ifile];
948  ADep_datafile << data_dir << "/tot_xsec/pipA_cex/pip" << nucleus << "_cex.txt";
949  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
950  buff->SetNameTitle("buff","buff");
951  for(int i=0; i < buff->GetN(); i++) {
952  buff -> GetPoint(i,x,y);
953  TPipA_CEx -> SetPoint(ipoint,(double)nucleus,x,y);
954  ipoint++;
955  }
956  delete buff;
957  }
958 
959  if (saveTGraphsToFile) {
960  TPipA_CEx -> Write("TPipA_CEx");
961  }
962  delete TPipA_CEx;
963 
964  }
965 
966  // kIHNFtAbs, pip + A PipA_Abs (just for developmental purposes)
967  {
968  TGraph2D * TPipA_Abs;
969 
970  const int pipAAbs_nfiles = 18;
971  const int pipAAbs_nuclei[pipAAbs_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 48, 56, 58, 63, 93, 120, 165, 181, 209};
972  const int pipAAbs_npoints = 111;
973 
974  TPipA_Abs = new TGraph2D(pipAAbs_npoints);
975  TPipA_Abs->SetNameTitle("TPipA_Abs","TPipA_Abs");
976  TPipA_Abs->SetDirectory(0);
977 
978  int ipoint=0;
979  double x, y;
980 
981  for(int ifile=0; ifile < pipAAbs_nfiles; ifile++) {
982  ostringstream ADep_datafile;
983  int nucleus = pipAAbs_nuclei[ifile];
984  ADep_datafile << data_dir << "/tot_xsec/pipA_abs/pip" << nucleus << "_abs.txt";
985  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
986  buff->SetNameTitle("buff","buff");
987  for(int i=0; i < buff->GetN(); i++) {
988  buff -> GetPoint(i,x,y);
989  TPipA_Abs -> SetPoint(ipoint,(double)nucleus,x,y);
990  ipoint++;
991  }
992  delete buff;
993  }
994 
995  if (saveTGraphsToFile) {
996  TPipA_Abs -> Write("TPipA_Abs");
997  }
998  delete TPipA_Abs;
999 
1000  }
1001 
1002  // kIHNFtElas, pip + A PipA_Elas (just for developmental purposes)
1003  {
1004  TGraph2D * TPipA_Elas;
1005 
1006  const int pipAElas_nfiles = 18;
1007  const int pipAElas_nuclei[pipAElas_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 48, 56, 58, 63, 93, 120, 165, 181, 209};
1008  const int pipAElas_npoints = 125;
1009 
1010  TPipA_Elas = new TGraph2D(pipAElas_npoints);
1011  TPipA_Elas->SetNameTitle("TPipA_Elas","TPipA_Elas");
1012  TPipA_Elas->SetDirectory(0);
1013 
1014  int ipoint=0;
1015  double x, y;
1016 
1017  for(int ifile=0; ifile < pipAElas_nfiles; ifile++) {
1018  ostringstream ADep_datafile;
1019  int nucleus = pipAElas_nuclei[ifile];
1020  ADep_datafile << data_dir << "/tot_xsec/pipA_elas/pip" << nucleus << "_elas.txt";
1021  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
1022  buff->SetNameTitle("buff","buff");
1023  for(int i=0; i < buff->GetN(); i++) {
1024  buff -> GetPoint(i,x,y);
1025  TPipA_Elas -> SetPoint(ipoint,(double)nucleus,x,y);
1026  ipoint++;
1027  }
1028  delete buff;
1029  }
1030 
1031  if (saveTGraphsToFile) {
1032  TPipA_Elas -> Write("TPipA_Elas");
1033  }
1034  delete TPipA_Elas;
1035  }
1036 
1037  // kIHNFtInelas, pip + A PipA_Inelas (just for developmental purposes)
1038  {
1039  TGraph2D * TPipA_Inelas;
1040 
1041  const int pipAInelas_nfiles = 20;
1042  const int pipAInelas_nuclei[pipAInelas_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 40, 48, 56, 58, 63, 93, 120, 165, 181, 208, 209};
1043  const int pipAInelas_npoints = 118;
1044 
1045  TPipA_Inelas = new TGraph2D(pipAInelas_npoints);
1046  TPipA_Inelas->SetNameTitle("TPipA_Inelas","TPipA_Inelas");
1047  TPipA_Inelas->SetDirectory(0);
1048 
1049  int ipoint=0;
1050  double x, y;
1051 
1052  for(int ifile=0; ifile < pipAInelas_nfiles; ifile++) {
1053  ostringstream ADep_datafile;
1054  int nucleus = pipAInelas_nuclei[ifile];
1055  ADep_datafile << data_dir << "/tot_xsec/pipA_inelas/pip" << nucleus << "_inelas.txt";
1056  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
1057  buff->SetNameTitle("buff","buff");
1058  for(int i=0; i < buff->GetN(); i++) {
1059  buff -> GetPoint(i,x,y);
1060  TPipA_Inelas -> SetPoint(ipoint,(double)nucleus,x,y);
1061  ipoint++;
1062  }
1063  delete buff;
1064  }
1065 
1066  if (saveTGraphsToFile) {
1067  TPipA_Inelas -> Write("TPipA_Inelas");
1068  }
1069  delete TPipA_Inelas;
1070  }
1071 
1072 
1073  /*
1074  // kIHNFtElas, pip + A PipA_Elas_frac
1075  {
1076  const int pipAElas_f_nfiles = 18;
1077  const int pipAElas_f_nuclei[pipAElas_f_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 48, 56, 58, 63, 93, 120, 165, 181, 209};
1078  const int pipAElas_f_npoints = 125;
1079 
1080  TfracPipA_Elas = new TGraph2D(pipAElas_f_npoints);
1081  TfracPipA_Elas->SetNameTitle("TfracPipA_Elas","TfracPipA_Elas");
1082  TfracPipA_Elas->SetDirectory(0);
1083 
1084  int ipoint=0;
1085  double x, y;
1086 
1087  for(int ifile=0; ifile < pipAElas_f_nfiles; ifile++) {
1088  ostringstream ADep_datafile;
1089  int nucleus = pipAElas_f_nuclei[ifile];
1090  ADep_datafile << data_dir << "/tot_xsec/pipA_elas_frac/pip" << nucleus << "_elas_frac.txt";
1091  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
1092  buff->SetNameTitle("buff","buff");
1093  for(int i=0; i < buff->GetN(); i++) {
1094  buff -> GetPoint(i,x,y);
1095  TfracPipA_Elas -> SetPoint(ipoint,(double)nucleus,x,y);
1096  ipoint++;
1097  }
1098  delete buff;
1099  }
1100 
1101  if (saveTGraphsToFile) {
1102  TfracPipA_Elas -> Write("TfracPipA_Elas");
1103  }
1104 
1105  }
1106  */
1107 
1108  // kIHNFtInelas, pip + A PipA_Inelas_frac
1109  {
1110  const int pipAInelas_f_nfiles = 20;
1111  const int pipAInelas_f_nuclei[pipAInelas_f_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 27, 40, 48, 56, 58, 63, 93, 120, 165, 181, 208, 209};
1112  const int pipAInelas_f_npoints = 118;
1113 
1114  TfracPipA_Inelas = new TGraph2D(pipAInelas_f_npoints);
1115  TfracPipA_Inelas->SetNameTitle("TfracPipA_Inelas","TfracPipA_Inelas");
1116  TfracPipA_Inelas->SetDirectory(0);
1117 
1118  int ipoint=0;
1119  double x, y;
1120 
1121  for(int ifile=0; ifile < pipAInelas_f_nfiles; ifile++) {
1122  ostringstream ADep_datafile;
1123  int nucleus = pipAInelas_f_nuclei[ifile];
1124  ADep_datafile << data_dir << "/tot_xsec/pipA_inelas_frac/pip" << nucleus << "_inelas_frac.txt";
1125  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
1126  buff->SetNameTitle("buff","buff");
1127  for(int i=0; i < buff->GetN(); i++) {
1128  buff -> GetPoint(i,x,y);
1129  TfracPipA_Inelas -> SetPoint(ipoint,(double)nucleus,x,y);
1130  ipoint++;
1131  }
1132  delete buff;
1133  }
1134 
1135  if (saveTGraphsToFile) {
1136  TfracPipA_Inelas -> Write("TfracPipA_Inelas");
1137  }
1138 
1139  }
1140 
1141 
1142  // kIHNFtPiPro, pip + A PipA_PiPro_frac
1143  {
1144  const int pipAPiPro_f_nfiles = 17;
1145  const int pipAPiPro_f_nuclei[pipAPiPro_f_nfiles] = {1, 2, 3, 4, 7, 9, 12, 16, 48, 56, 58, 63, 93, 120, 165, 181, 209};
1146  const int pipAPiPro_f_npoints = 76;
1147 
1148  TfracPipA_PiPro = new TGraph2D(pipAPiPro_f_npoints);
1149  TfracPipA_PiPro->SetNameTitle("TfracPipA_PiPro","TfracPipA_PiPro");
1150  TfracPipA_PiPro->SetDirectory(0);
1151 
1152  int ipoint=0;
1153  double x, y;
1154 
1155  for(int ifile=0; ifile < pipAPiPro_f_nfiles; ifile++) {
1156  ostringstream ADep_datafile;
1157  int nucleus = pipAPiPro_f_nuclei[ifile];
1158  ADep_datafile << data_dir << "/tot_xsec/pipA_pipro_frac/pip" << nucleus << "_pipro_frac.txt";
1159  TGraph * buff = new TGraph(ADep_datafile.str().c_str());
1160  buff->SetNameTitle("buff","buff");
1161  for(int i=0; i < buff->GetN(); i++) {
1162  buff -> GetPoint(i,x,y);
1163  TfracPipA_PiPro -> SetPoint(ipoint,(double)nucleus,x,y);
1164  ipoint++;
1165  }
1166  delete buff;
1167  }
1168 
1169  if (saveTGraphsToFile) {
1170  TfracPipA_PiPro -> Write("TfracPipA_PiPro");
1171  }
1172  }
1173 
1174  TGraphs_file.Close();
1175 
1176  LOG("INukeData", pINFO) << "Done building x-section splines...";
1177 
1178 }
1179 //____________________________________________________________________________
1181  string filename, double ke, int npoints, int & curr_point,
1182  double * costh_array, double * xsec_array, int cols)
1183 {
1184  // open
1185  std::ifstream hN_stream(filename.c_str(), ios::in);
1186  if(!hN_stream.good()) {
1187  LOG("INukeData", pERROR)
1188  << "Error reading INTRANUKE/hN data from: " << filename;
1189  return;
1190  }
1191 
1192  if(cols<2) {
1193  LOG("INukeData", pERROR)
1194  << "Error reading INTRANUKE/hN data from: " << filename;
1195  LOG("INukeData", pERROR)
1196  << "Too few columns: " << cols;
1197  return;
1198  }
1199 
1200  LOG("INukeData", pINFO)
1201  << "Reading INTRANUKE/hN data from: " << filename;
1202 
1203  // skip initial comments
1204  char cbuf[501];
1205  hN_stream.getline(cbuf,400);
1206  hN_stream.getline(cbuf,400);
1207  hN_stream.getline(cbuf,400);
1208 
1209  // read
1210  double angle = 0;
1211  double xsec = 0;
1212  double trash = 0;
1213 
1214  for(int ip = 0; ip < npoints; ip++) {
1215  hN_stream >> angle >> xsec;
1216 
1217  for(int ic = 0; ic < (cols-2); ic++) {
1218  hN_stream >> trash;
1219  }
1220 
1221  LOG("INukeData", pDEBUG)
1222  << "Adding data point: (KE = " << ke << " MeV, angle = "
1223  << angle << ", sigma = " << xsec << " mbarn)";
1224  costh_array[ip] = TMath::Cos(angle*kPi/180.);
1225  xsec_array [curr_point] = xsec;
1226  curr_point++;
1227  }
1228 }
1229 //____________________________________________________________________________
1231  int hpdgc, int tgtpdgc, int nppdgc, INukeFateHN_t fate, double ke, double costh) const
1232 {
1233 // inputs
1234 // fate : h+N fate code
1235 // hpdgc : h PDG code
1236 // tgtpdgc : N PDG code
1237 // nppdgc : product N PDG code
1238 // ke : kinetic energy (MeV)
1239 // costh : cos(scattering angle)
1240 // returns
1241 // xsec : mbarn
1242 
1243  double ke_eval = ke;
1244  double costh_eval = costh;
1245 
1246  costh_eval = TMath::Min(costh, 1.);
1247  costh_eval = TMath::Max(costh_eval, -1.);
1248 
1249  if(fate==kIHNFtElas) {
1250 
1251  if( (hpdgc==kPdgProton && tgtpdgc==kPdgProton) ||
1252  (hpdgc==kPdgNeutron && tgtpdgc==kPdgNeutron) )
1253  {
1254  ke_eval = TMath::Min(ke_eval, 999.);
1255  ke_eval = TMath::Max(ke_eval, 50.);
1256  return fhN2dXSecPP_Elas->Evaluate(ke_eval, costh_eval);
1257  }
1258  else
1259  if( (hpdgc==kPdgProton && tgtpdgc==kPdgNeutron) ||
1260  (hpdgc==kPdgNeutron && tgtpdgc==kPdgProton) )
1261  {
1262  ke_eval = TMath::Min(ke_eval, 999.);
1263  ke_eval = TMath::Max(ke_eval, 50.);
1264  return fhN2dXSecNP_Elas->Evaluate(ke_eval, costh_eval);
1265  }
1266  else
1267  if(hpdgc==kPdgPiP)
1268  {
1269  ke_eval = TMath::Min(ke_eval, 1499.);
1270  ke_eval = TMath::Max(ke_eval, 10.);
1271  return fhN2dXSecPipN_Elas->Evaluate(ke_eval, costh_eval);
1272  }
1273  else
1274  if(hpdgc==kPdgPi0)
1275  {
1276  ke_eval = TMath::Min(ke_eval, 1499.);
1277  ke_eval = TMath::Max(ke_eval, 10.);
1278  return fhN2dXSecPi0N_Elas->Evaluate(ke_eval, costh_eval);
1279  }
1280  else
1281  if(hpdgc==kPdgPiM)
1282  {
1283  ke_eval = TMath::Min(ke_eval, 1499.);
1284  ke_eval = TMath::Max(ke_eval, 10.);
1285  return fhN2dXSecPimN_Elas->Evaluate(ke_eval, costh_eval);
1286  }
1287  else
1288  if(hpdgc==kPdgKP && tgtpdgc==kPdgNeutron)
1289  {
1290  ke_eval = TMath::Min(ke_eval, 1799.);
1291  ke_eval = TMath::Max(ke_eval, 100.);
1292  return fhN2dXSecKpN_Elas->Evaluate(ke_eval, costh_eval);
1293  }
1294  else
1295  if(hpdgc==kPdgKP && tgtpdgc==kPdgProton)
1296  {
1297  ke_eval = TMath::Min(ke_eval, 1799.);
1298  ke_eval = TMath::Max(ke_eval, 100.);
1299  return fhN2dXSecKpP_Elas->Evaluate(ke_eval, costh_eval);
1300  }
1301  }
1302 
1303  else if(fate == kIHNFtCEx) {
1304  if( (hpdgc==kPdgPiP || hpdgc==kPdgPi0 || hpdgc==kPdgPiM) &&
1305  (tgtpdgc==kPdgProton || tgtpdgc==kPdgNeutron) )
1306  {
1307  ke_eval = TMath::Min(ke_eval, 1499.);
1308  ke_eval = TMath::Max(ke_eval, 10.);
1309  return fhN2dXSecPiN_CEx->Evaluate(ke_eval, costh_eval);
1310  }
1311  else if( (hpdgc == kPdgProton && tgtpdgc == kPdgProton) ||
1312  (hpdgc == kPdgNeutron && tgtpdgc == kPdgNeutron) )
1313  {
1314  LOG("INukeData", pWARN) << "Inelastic pp does not exist!";
1315  ke_eval = TMath::Min(ke_eval, 999.);
1316  ke_eval = TMath::Max(ke_eval, 50.);
1317  return fhN2dXSecPP_Elas->Evaluate(ke_eval, costh_eval);
1318  }
1319  else if( (hpdgc == kPdgProton && tgtpdgc == kPdgNeutron) ||
1320  (hpdgc == kPdgNeutron && tgtpdgc == kPdgProton) )
1321  {
1322  ke_eval = TMath::Min(ke_eval, 999.);
1323  ke_eval = TMath::Max(ke_eval, 50.);
1324  return fhN2dXSecNP_Elas->Evaluate(ke_eval, costh_eval);
1325  }
1326  else if(hpdgc == kPdgKP && tgtpdgc == kPdgNeutron) {
1327  ke_eval = TMath::Min(ke_eval, 1799.);
1328  ke_eval = TMath::Max(ke_eval, 100.);
1329  return fhN2dXSecKpN_CEx->Evaluate(ke_eval, costh_eval);
1330  }
1331  }
1332 
1333  else if(fate == kIHNFtAbs) {
1334  if( (hpdgc==kPdgPiP || hpdgc==kPdgPi0 || hpdgc==kPdgPiM) &&
1335  (tgtpdgc==kPdgProton || tgtpdgc==kPdgNeutron) )
1336  {
1337  ke_eval = TMath::Min(ke_eval, 499.);
1338  ke_eval = TMath::Max(ke_eval, 50.);
1339  return fhN2dXSecPiN_Abs->Evaluate(ke_eval, costh_eval);
1340  }
1341  if(hpdgc==kPdgKP) return 1.; //isotropic since no data ???
1342  }
1343 
1344  else if(fate == kIHNFtInelas) {
1345  if( hpdgc==kPdgGamma && tgtpdgc==kPdgProton &&nppdgc==kPdgProton )
1346  {
1347  ke_eval = TMath::Min(ke_eval, 1199.);
1348  ke_eval = TMath::Max(ke_eval, 160.);
1349  return fhN2dXSecGamPi0P_Inelas->Evaluate(ke_eval, costh_eval);
1350  }
1351  else
1352  if( hpdgc==kPdgGamma && tgtpdgc==kPdgProton && nppdgc==kPdgNeutron )
1353  {
1354  ke_eval = TMath::Min(ke_eval, 1199.);
1355  ke_eval = TMath::Max(ke_eval, 160.);
1356  return fhN2dXSecGamPipN_Inelas->Evaluate(ke_eval, costh_eval);
1357  }
1358  else
1359  if( hpdgc==kPdgGamma && tgtpdgc==kPdgNeutron && nppdgc==kPdgProton )
1360  {
1361  ke_eval = TMath::Min(ke_eval, 1199.);
1362  ke_eval = TMath::Max(ke_eval, 160.);
1363  return fhN2dXSecGamPimP_Inelas->Evaluate(ke_eval, costh_eval);
1364  }
1365  else
1366  if( hpdgc==kPdgGamma && tgtpdgc==kPdgNeutron && nppdgc==kPdgNeutron )
1367  {
1368  ke_eval = TMath::Min(ke_eval, 1199.);
1369  ke_eval = TMath::Max(ke_eval, 160.);
1370  return fhN2dXSecGamPi0N_Inelas->Evaluate(ke_eval, costh_eval);
1371  }
1372  }
1373 
1374  return 0;
1375 }
1376 //____________________________________________________________________________
1377 double INukeHadroData2018::FracADep(int hpdgc, INukeFateHA_t fate, double ke, int targA) const
1378 {
1379  // return the x-section fraction for the input fate for the particle with the input pdg
1380  // code and the target with the input mass number at the input kinetic energy
1381 
1382  ke = TMath::Max(fMinKinEnergy, ke); // ke >= 1 MeV
1383  ke = TMath::Min(fMaxKinEnergyHA, ke); // ke <= 999 MeV
1384 
1385  targA = TMath::Min(208, targA); // A <= 208
1386 
1387  LOG("INukeData", pDEBUG) << "Querying hA cross section at ke = " << ke << " and target " << targA;
1388 
1389  // Handle pions (currently the same cross sections are used for pi+, pi-, and pi0)
1390  if ( hpdgc == kPdgPiP || hpdgc == kPdgPiM || hpdgc == kPdgPi0 ) {
1391 
1392  double frac_cex = TfracPipA_CEx->Interpolate(targA, ke);
1393  //double frac_elas = TfracPipA_Elas->Interpolate(targA, ke);
1394  double frac_inelas = TfracPipA_Inelas->Interpolate(targA, ke);
1395  double frac_abs = TfracPipA_Abs->Interpolate(targA, ke);
1396  double frac_pipro = TfracPipA_PiPro->Interpolate(targA, ke);
1397 
1398  // Protect against unitarity violation due to interpolation problems
1399  // by renormalizing all available fate fractions to unity.
1400  double total = frac_cex + frac_inelas + frac_abs + frac_pipro; // + frac_elas
1401 
1402  if ( fate == kIHAFtCEx ) return frac_cex / total;
1403  //else if ( fate == kIHAFtElas ) return frac_elas / total;
1404  else if ( fate == kIHAFtInelas ) return frac_inelas / total;
1405  else if ( fate == kIHAFtAbs ) return frac_abs / total;
1406  else if ( fate == kIHAFtPiProd ) return frac_pipro / total;
1407  else {
1408  std::string sign("+");
1409  if ( hpdgc == kPdgPiM ) sign = "-";
1410  else if ( hpdgc == kPdgPi0 ) sign = "0";
1411  LOG("INukeData", pWARN) << "Pi" << sign << "'s don't have this fate: " << INukeHadroFates::AsString(fate);
1412  return 0.;
1413  }
1414  }
1415 
1416  LOG("INukeData", pWARN) << "Can't handle particles with pdg code = " << hpdgc;
1417  return 0.;
1418 }
1419 //____________________________________________________________________________
1420 double INukeHadroData2018::FracAIndep(int hpdgc, INukeFateHA_t fate, double ke) const
1421 {
1422  // return the x-section fraction for the input fate for the particle with the input pdg
1423  // code at the input kinetic energy
1424  ke = TMath::Max(fMinKinEnergy, ke);
1425  ke = TMath::Min(fMaxKinEnergyHA, ke);
1426 
1427  LOG("INukeData", pDEBUG) << "Querying hA cross section at ke = " << ke;
1428 
1429  // TODO: reduce code duplication here
1430  if (hpdgc == kPdgProton) {
1431  // handle protons
1432  double frac_cex = fFracPA_CEx->Evaluate(ke);
1433  double frac_inelas = fFracPA_Inel->Evaluate(ke);
1434  double frac_abs = fFracPA_Abs->Evaluate(ke);
1435  double frac_pipro = fFracPA_PiPro->Evaluate(ke);
1436  double frac_comp = fFracPA_Cmp->Evaluate(ke);
1437 
1438  // Protect against unitarity violation due to interpolation problems
1439  // by renormalizing all available fate fractions to unity.
1440  double total = frac_cex + frac_inelas + frac_abs + frac_pipro + frac_comp; // + frac_elas
1441 
1442  if ( fate == kIHAFtCEx ) return frac_cex / total;
1443  //else if ( fate == kIHAFtElas ) return frac_elas / total;
1444  else if ( fate == kIHAFtInelas ) return frac_inelas / total;
1445  else if ( fate == kIHAFtAbs ) return frac_abs / total;
1446  else if ( fate == kIHAFtPiProd ) return frac_pipro / total;
1447  else if ( fate == kIHAFtCmp ) return frac_comp / total; // cmp - add support for this later
1448  else {
1449  LOG("INukeData", pWARN)
1450  << "Protons don't have this fate: " << INukeHadroFates::AsString(fate);
1451  return 0;
1452  }
1453  }
1454  else if (hpdgc == kPdgNeutron) {
1455  // handle neutrons
1456  double frac_cex = fFracNA_CEx->Evaluate(ke);
1457  double frac_inelas = fFracNA_Inel->Evaluate(ke);
1458  double frac_abs = fFracNA_Abs->Evaluate(ke);
1459  double frac_pipro = fFracNA_PiPro->Evaluate(ke);
1460  double frac_comp = fFracNA_Cmp->Evaluate(ke);
1461 
1462  // Protect against unitarity violation due to interpolation problems
1463  // by renormalizing all available fate fractions to unity.
1464  double total = frac_cex + frac_inelas + frac_abs + frac_pipro + frac_comp; // + frac_elas
1465 
1466  if ( fate == kIHAFtCEx ) return frac_cex / total;
1467  //else if ( fate == kIHAFtElas ) return frac_elas / total;
1468  else if ( fate == kIHAFtInelas ) return frac_inelas / total;
1469  else if ( fate == kIHAFtAbs ) return frac_abs / total;
1470  else if ( fate == kIHAFtPiProd ) return frac_pipro / total;
1471  else if ( fate == kIHAFtCmp ) return frac_comp / total; // cmp - add support for this later
1472  else {
1473  LOG("INukeData", pWARN)
1474  << "Neutrons don't have this fate: " << INukeHadroFates::AsString(fate);
1475  return 0;
1476  }
1477  }
1478  else if (hpdgc == kPdgKP) {
1479  // handle K+
1480  double frac_inelas = fFracKA_Inel->Evaluate(ke);
1481  //double frac_elas = fFracKA_Elas->Evaluate(ke);
1482  double frac_abs = fFracKA_Abs->Evaluate(ke);
1483 
1484  // Protect against unitarity violation due to interpolation problems
1485  // by renormalizing all available fate fractions to unity.
1486  double total = frac_inelas + frac_abs; // + frac_elas
1487 
1488  if ( fate == kIHAFtInelas ) return frac_inelas / total;
1489  else if ( fate == kIHAFtAbs ) return frac_abs / total;
1490  else {
1491  LOG("INukeData", pWARN)
1492  << "K+'s don't have this fate: " << INukeHadroFates::AsString(fate);
1493  return 0.;
1494  }
1495  }
1496  LOG("INukeData", pWARN) << "Can't handle particles with pdg code = " << hpdgc;
1497  return 0.;
1498 }
1499 //____________________________________________________________________________
1500 double INukeHadroData2018::XSec(int hpdgc, INukeFateHN_t fate, double ke, int targA, int targZ) const
1501 {
1502 // return the x-section for the input fate for the particle with the input pdg
1503 // code at the input kinetic energy
1504 //
1505  ke = TMath::Max(fMinKinEnergy, ke);
1506  ke = TMath::Min(fMaxKinEnergyHN, ke);
1507 
1508  LOG("INukeData", pDEBUG) << "Querying hN cross section at ke = " << ke;
1509 
1510  double xsec=0;
1511 
1512  if (hpdgc == kPdgPiP) {
1513  /* handle pi+ */
1514  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecPipp_CEx -> Evaluate(ke)) * targZ;
1515  xsec+= TMath::Max(0., fXSecPipn_CEx -> Evaluate(ke)) * (targA-targZ);
1516  return xsec;}
1517  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPipp_Elas -> Evaluate(ke)) * targZ;
1518  xsec+= TMath::Max(0., fXSecPipn_Elas -> Evaluate(ke)) * (targA-targZ);
1519  return xsec;}
1520  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPipp_Reac -> Evaluate(ke)) * targZ;
1521  xsec+= TMath::Max(0., fXSecPipn_Reac -> Evaluate(ke)) * (targA-targZ);
1522  return xsec;}
1523  else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecPipd_Abs -> Evaluate(ke)) * targA;
1524  return xsec;}
1525  else {
1526  LOG("INukeData", pWARN)
1527  << "Pi+'s don't have this fate: " << INukeHadroFates::AsString(fate);
1528  return 0;
1529  }
1530 
1531  } else if (hpdgc == kPdgPiM) {
1532  /* handle pi- */
1533  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecPipn_CEx -> Evaluate(ke)) * targZ;
1534  xsec+= TMath::Max(0., fXSecPipp_CEx -> Evaluate(ke)) * (targA-targZ);
1535  return xsec;}
1536  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPipn_Elas -> Evaluate(ke)) * targZ;
1537  xsec+= TMath::Max(0., fXSecPipp_Elas -> Evaluate(ke)) * (targA-targZ);
1538  return xsec;}
1539  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPipn_Reac -> Evaluate(ke)) * targZ;
1540  xsec+= TMath::Max(0., fXSecPipp_Reac -> Evaluate(ke)) * (targA-targZ);
1541  return xsec;}
1542  else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecPipd_Abs -> Evaluate(ke)) * targA;
1543  return xsec;}
1544  else {
1545  LOG("INukeData", pWARN)
1546  << "Pi-'s don't have this fate: " << INukeHadroFates::AsString(fate);
1547  return 0;
1548  }
1549 
1550  } else if (hpdgc == kPdgPi0) {
1551  /* handle pi0 */
1552  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecPi0p_CEx -> Evaluate(ke)) * targZ;
1553  xsec+= TMath::Max(0., fXSecPi0n_CEx -> Evaluate(ke)) * (targA-targZ);
1554  return xsec;}
1555  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPi0p_Elas -> Evaluate(ke)) * targZ;
1556  xsec+= TMath::Max(0., fXSecPi0n_Elas -> Evaluate(ke)) * (targA-targZ);
1557  return xsec;}
1558  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPi0p_Reac -> Evaluate(ke)) * targZ;
1559  xsec+= TMath::Max(0., fXSecPi0n_Reac -> Evaluate(ke)) * (targA-targZ);
1560  return xsec;}
1561  else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecPi0d_Abs -> Evaluate(ke)) * targA;
1562  return xsec;}
1563  else {
1564  LOG("INukeData", pWARN)
1565  << "Pi0's don't have this fate: " << INukeHadroFates::AsString(fate);
1566  return 0;
1567  }
1568 
1569  } else if (hpdgc == kPdgProton) {
1570  /* handle protons */
1571  if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPp_Elas -> Evaluate(ke)) * targZ;
1572  xsec+= TMath::Max(0., fXSecPn_Elas -> Evaluate(ke)) * (targA-targZ);
1573  return xsec;}
1574  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPp_Reac -> Evaluate(ke)) * targZ;
1575  xsec+= TMath::Max(0., fXSecPn_Reac -> Evaluate(ke)) * (targA-targZ);
1576  return xsec;}
1577  else if (fate == kIHNFtCmp) {xsec = TMath::Max(0., fXSecPp_Cmp -> Evaluate(ke)) * targZ;
1578  xsec+= TMath::Max(0., fXSecPn_Cmp -> Evaluate(ke)) * (targA-targZ);
1579  return xsec;}
1580  else {
1581  LOG("INukeData", pWARN)
1582  << "Protons don't have this fate: " << INukeHadroFates::AsString(fate);
1583  return 0;
1584  }
1585 
1586  } else if (hpdgc == kPdgNeutron) {
1587  /* handle protons */
1588  if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPn_Elas -> Evaluate(ke)) * targZ;
1589  xsec+= TMath::Max(0., fXSecNn_Elas -> Evaluate(ke)) * (targA-targZ);
1590  return xsec;}
1591  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPn_Reac -> Evaluate(ke)) * targZ;
1592  xsec+= TMath::Max(0., fXSecNn_Reac -> Evaluate(ke)) * (targA-targZ);
1593  return xsec;}
1594  else if (fate == kIHNFtCmp) {xsec = TMath::Max(0., fXSecPp_Cmp -> Evaluate(ke)) * targZ;
1595  xsec+= TMath::Max(0., fXSecPn_Cmp -> Evaluate(ke)) * (targA-targZ);
1596  return xsec;}
1597  else {
1598  LOG("INukeData", pWARN)
1599  << "Neutrons don't have this fate: " << INukeHadroFates::AsString(fate);
1600  return 0;
1601  }
1602  //Adding here kaons, why elastic only on protons? hA or hN? No _Reac for kaons...
1603  } else if (hpdgc == kPdgKP) {
1604  /* handle K+ */
1605  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecKpn_CEx -> Evaluate(ke)) * targZ;
1606  xsec+= TMath::Max(0., fXSecKpn_CEx -> Evaluate(ke)) * (targA-targZ);
1607  return xsec;}
1608  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecKpn_Elas -> Evaluate(ke)) * targZ;
1609  xsec+= TMath::Max(0., fXSecKpn_Elas -> Evaluate(ke)) * (targA-targZ);
1610  return xsec;}
1611  /*else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecKpd_Abs -> Evaluate(ke)) * targA;
1612  return xsec;}*/
1613  else {
1614  LOG("INukeData", pWARN)
1615  << "K+'s don't have this fate: " << INukeHadroFates::AsString(fate);
1616  return 0;
1617  }
1618  //------------------------------------------------
1619  /* } else if (hpdgc == kPdgGamma) {
1620  / * handle gamma * /
1621  if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecGamp_fs -> Evaluate(ke)) * targZ;
1622  xsec+= TMath::Max(0., fXSecGamn_fs -> Evaluate(ke)) * (targA-targZ);
1623  return xsec;}
1624  else {
1625  LOG("INukeData", pWARN)
1626  << "Gamma's don't have this fate: " << INukeHadroFates::AsString(fate);
1627  return 0;
1628  }*/
1629  }
1630  LOG("INukeData", pWARN)
1631  << "Can't handle particles with pdg code = " << hpdgc;
1632 
1633  return 0;
1634 }
1635 
1636 double INukeHadroData2018::Frac(int hpdgc, INukeFateHN_t fate, double ke, int targA, int targZ) const
1637 {
1638 // return the x-section fraction for the input fate for the particle with the
1639 // input pdg code at the input kinetic energy
1640 
1641  ke = TMath::Max(fMinKinEnergy, ke);
1642  ke = TMath::Min(fMaxKinEnergyHN, ke);
1643 
1644  // get x-section
1645  double xsec = this->XSec(hpdgc,fate,ke,targA,targZ);
1646 
1647  // get max x-section
1648  double xsec_tot = 0;
1649  if (hpdgc == kPdgPiP ){xsec_tot = TMath::Max(0., fXSecPipp_Tot -> Evaluate(ke)) * targZ;
1650  xsec_tot+= TMath::Max(0., fXSecPipn_Tot -> Evaluate(ke)) * (targA-targZ);}
1651  else if (hpdgc == kPdgPiM ){xsec_tot = TMath::Max(0., fXSecPipn_Tot -> Evaluate(ke)) * targZ;
1652  xsec_tot+= TMath::Max(0., fXSecPipp_Tot -> Evaluate(ke)) * (targA-targZ);}
1653  else if (hpdgc == kPdgPi0 ){xsec_tot = TMath::Max(0., fXSecPi0p_Tot -> Evaluate(ke)) * targZ;
1654  xsec_tot+= TMath::Max(0., fXSecPi0n_Tot -> Evaluate(ke)) * (targA-targZ);}
1655  else if (hpdgc == kPdgProton ){xsec_tot = TMath::Max(0., fXSecPp_Tot -> Evaluate(ke)) * targZ;
1656  xsec_tot+= TMath::Max(0., fXSecPn_Tot -> Evaluate(ke)) * (targA-targZ);}
1657  else if (hpdgc == kPdgNeutron){xsec_tot = TMath::Max(0., fXSecPn_Tot -> Evaluate(ke)) * targZ;
1658  xsec_tot+= TMath::Max(0., fXSecNn_Tot -> Evaluate(ke)) * (targA-targZ);}
1659  else if (hpdgc == kPdgGamma ) xsec_tot = TMath::Max(0., fXSecGamN_Tot -> Evaluate(ke));
1660  else if (hpdgc == kPdgKP ) xsec_tot = TMath::Max(0., fXSecKpN_Tot -> Evaluate(ke));
1661 
1662  // compute fraction
1663  double frac = (xsec_tot>0) ? xsec/xsec_tot : 0.;
1664  return frac;
1665 }
1666 //____________________________________________________________________________
1667 double INukeHadroData2018::IntBounce(const GHepParticle* p, int target, int scode, INukeFateHN_t fate)
1668 {
1669  // This method returns a random cos(ang) according to a distribution
1670  // based upon the particle and fate. The sampling uses the
1671  // Accept/Reject method, whereby a distribution is bounded above by
1672  // an envelope, or in this case, a number of envelopes, which can be
1673  // easily sampled (here, we use uniform distributions).
1674  // To get a random value, first the envelope is sampled to
1675  // obtain an x-coordinate (cos(ang)), and then another random value
1676  // is obtained uniformally in the range [0,h(j,0)], where h(j,0)
1677  // is the height of the j-th envelope. If the point is beneath the
1678  // distribution, the x-coordinate is accepted, otherwise, we try
1679  // again.
1680 
1681  RandomGen * rnd = RandomGen::Instance();
1682 
1683  // numEnv is the number of envelopes in the total envelope,
1684  // that is, the number of seperate simple uniform distributions
1685  // that will be fit against the distribution in question in the
1686  // Accept/Reject process of sampling
1687  int numEnv = 4;
1688  int numPoints = 1000; // The number of points to be evaluated
1689  // for the purpose of finding the max
1690  // value of the distribution
1691  assert((numPoints%numEnv)==0); // numPoints/numEnv has to be an integer
1692  double sr = 2.0 / numEnv; // Subrange, i.e., range of an envelope
1693  double cstep = 2.0 / (numPoints); // Magnitude of the step between eval. points
1694 
1695  double ke = (p->E() - p->Mass()) * 1000.0; // ke in MeV
1696  if (TMath::Abs((int)ke-ke)<.01) ke+=.3; // make sure ke isn't an integer,
1697  // otherwise sometimes gives weird results
1698  // due to ROOT's Interpolate() function
1699  double avg = 0.0; // average value in envelop
1700 
1701  // Matrices to hold data; buff holds the distribution
1702  // data per envelope from which the max value is
1703  // obtained. That value is then recorded in dist, where
1704  // the integral of the envelope to that point is
1705  // also recorded
1706 
1707  double * buff = new double[numPoints/numEnv + 1];
1708  double ** dist = new double*[numEnv];
1709  for(int ih=0;ih<numEnv;ih++)
1710  {
1711  dist[ih] = new double[3];
1712  }
1713 
1714  // Acc-Rej Sampling Method
1715  // -- Starting at the beginning of each envelope,
1716  // this loop evaluates (numPoints) amount of points
1717  // on the distribution and holds them in buff;
1718  // then takes the max and places it in the first row
1719  // of h. The second row of h contains the interval of
1720  // the total envelope, up to the current envelope.
1721  // Thus, when properly normalized, the last value
1722  // in the second row of h should be 1.
1723  double totxsec = 0.0;
1724  for(int i=0;i<numEnv;i++)
1725  {
1726  double lbound = -1 + i*sr;
1727 
1728  for(int j=0;j<=numPoints / numEnv; j++)
1729  {
1730  buff[j] = this->XSec(p->Pdg(),target,scode,fate,ke,lbound+j*cstep);
1731  avg += buff[j];
1732  }
1733 
1734  totxsec+=avg;
1735  avg/= (double(numPoints)/double(numEnv));
1736  dist[i][0] = TMath::MaxElement(numPoints/numEnv+1,buff);
1737  dist[i][1] = avg;
1738  dist[i][2] = dist[i][1] + ((i==0)?0.0:dist[i-1][2]);
1739  avg=0.0;
1740  }
1741 
1742 
1743  delete [] buff;
1744 
1745  int iter=1; // keep track of iterations
1746  int env=0; // envelope index
1747  double rval = 0.0; // random value
1748  double val = 0.0; // angle value
1749 
1750  // Get a random point, see if its in the distribution, and if not
1751  // then try again.
1752 
1753  rval = rnd->RndFsi().Rndm()*dist[numEnv-1][2];
1754 
1755  env=0;
1756  // Check which envelope it's in, to
1757  // get proper height
1758  while(env<numEnv)
1759  {
1760  if(rval<=dist[env][2]) break;
1761  else env++;
1762  }
1763  if(env==numEnv) env=numEnv - 1;
1764 
1765 while(iter)
1766  {
1767 
1768  // Obtain the correct x-coordinate from the random sample
1769  val = rnd->RndFsi().Rndm()*sr;
1770  val += sr*env-1;
1771  rval = rnd->RndFsi().Rndm()*dist[env][0];
1772 
1773  // Test to see if point is in distribution, if it is, stop and return
1774  if(rval < this->XSec(p->Pdg(),target,scode,fate,ke,val)) break;
1775 
1776  // Possibly an extremely long loop, don't want to
1777  // hold up the program
1778  if(iter==1000)
1779  {
1780  int NUM_POINTS=2000;
1781  int pvalues=0;
1782  double points[200]={0};
1783  for(int k=0;k<NUM_POINTS;k++)
1784  {
1785  points[int(k/10)]=this->XSec(p->Pdg(),target,scode,fate,ke,-1+(2.0/NUM_POINTS)*k);
1786  if(points[int(k/10)]>0) pvalues++;
1787  }
1788  if(pvalues<(.05*NUM_POINTS))
1789  {
1790  // if it reaches here, one more test...if momenta of particle is
1791  // extremely low, just give it an angle from a uniform distribution
1792  if(p->P4()->P()<.005) // 5 MeV
1793  {
1794  val = 2*rnd->RndFsi().Rndm()-1;
1795  break;
1796  }
1797  else
1798  {
1799  LOG("Intranuke", pWARN) << "Hung-up in IntBounce method - Exiting";
1800  LOG("Intranuke", pWARN) << (*p);
1801  LOG("Intranuke", pWARN) << "Target: " << target << ", Scode: " << scode << ", fate: " << INukeHadroFates::AsString(fate);
1802  for(int ie=0;ie<200;ie+=10) {
1803  LOG("Intranuke", pWARN) << points[ie+0] << ", " << points[ie+1] << ", " << points[ie+2] << ", "
1804  << points[ie+3] << ", " << points[ie+4] << ", " << points[ie+5] << ", " << points[ie+6] << ", "
1805  << points[ie+7] << ", " << points[ie+8] << ", " << points[ie+9];
1806  }
1807  for(int ih=0;ih<numEnv;ih++)
1808  {
1809  delete [] dist[ih];
1810  }
1811  delete [] dist;
1812 
1813  return -2.;
1814  }
1815  }
1816  }
1817  iter++;
1818  }
1819 
1820  for(int ih=0;ih<numEnv;ih++)
1821  {
1822  delete [] dist[ih];
1823  }
1824  delete [] dist;
1825 
1826  return val;
1827 }
1828 //___________________________________________________________________________
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
#define pERROR
Definition: Messenger.h:59
double Frac(int hpdgc, INukeFateHN_t fate, double ke, int targA=0, int targZ=0) const
double E(void) const
Get energy.
Definition: GHepParticle.h:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
double Mass(void) const
Mass that corresponds to the PDG code.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
int Pdg(void) const
Definition: GHepParticle.h:63
double FracADep(int hpdgc, INukeFateHA_t fate, double ke, int targA) const
void ReadhNFile(string filename, double ke, int npoints, int &curr_point, double *costh_array, double *xsec_array, int cols)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
enum genie::EINukeFateHN_t INukeFateHN_t
const int kPdgGamma
Definition: PDGCodes.h:189
const int kPdgKP
Definition: PDGCodes.h:172
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
double FracAIndep(int hpdgc, INukeFateHA_t fate, double ke) const
#define pINFO
Definition: Messenger.h:62
static INukeHadroData2018 * fInstance
#define pWARN
Definition: Messenger.h:60
static INukeHadroData2018 * Instance(void)
static string AsString(INukeFateHN_t fate)
const int kPdgPiM
Definition: PDGCodes.h:159
const int kPdgProton
Definition: PDGCodes.h:81
hadnt Write("hadnt")
#define pNOTICE
Definition: Messenger.h:61
static constexpr double sr
Definition: Units.h:166
const int kPdgNeutron
Definition: PDGCodes.h:83
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
enum genie::EINukeFateHA_t INukeFateHA_t
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
double XSec(int hpdgc, int tgt, int nprod, INukeFateHN_t rxnType, double ke, double costh) const
#define pDEBUG
Definition: Messenger.h:63