INukeHadroData.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (C) 2003-2015, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>, 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 
29 */
30 //____________________________________________________________________________
31 
32 #include <cassert>
33 #include <string>
34 
35 #include <TSystem.h>
36 #include <TNtupleD.h>
37 #include <TGraph2D.h>
38 #include <TTree.h>
39 #include <TMath.h>
40 
48 
49 using std::ostringstream;
50 using std::ios;
51 
52 using namespace genie;
53 using namespace genie::constants;
54 
55 //____________________________________________________________________________
57 //____________________________________________________________________________
58 double INukeHadroData::fMinKinEnergy = 1.0; // MeV
59 double INukeHadroData::fMaxKinEnergyHA = 999.0; // MeV
60 double INukeHadroData::fMaxKinEnergyHN = 1799.0; // MeV
61 //____________________________________________________________________________
63 {
64  this->LoadCrossSections();
65  fInstance = 0;
66 }
67 //____________________________________________________________________________
69 {
70  // pi+n/p hA x-section splines
71  delete fXSecPipn_Tot;
72  delete fXSecPipn_CEx;
73  delete fXSecPipn_Elas;
74  delete fXSecPipn_Reac;
75  delete fXSecPipp_Tot;
76  delete fXSecPipp_CEx;
77  delete fXSecPipp_Elas;
78  delete fXSecPipp_Reac;
79  delete fXSecPipd_Abs;
80 
81  // pi0n/p hA x-section splines
82  delete fXSecPi0n_Tot;
83  delete fXSecPi0n_CEx;
84  delete fXSecPi0n_Elas;
85  delete fXSecPi0n_Reac;
86  delete fXSecPi0p_Tot;
87  delete fXSecPi0p_CEx;
88  delete fXSecPi0p_Elas;
89  delete fXSecPi0p_Reac;
90  delete fXSecPi0d_Abs;
91 
92  // K+N x-section splines (elastic only)
93  delete fXSecKpn_Elas;
94  delete fXSecKpp_Elas;
95  delete fXSecKpN_Abs;
96  delete fXSecKpN_Tot;
97 
98  // gamma x-section splines (inelastic only)
99  delete fXSecGamp_fs;
100  delete fXSecGamn_fs;
101  delete fXSecGamN_Tot;
102 
103  // N+A x-section splines
104  delete fFracPA_Tot;
105  delete fFracPA_Elas;
106  delete fFracPA_Inel;
107  delete fFracPA_CEx;
108  delete fFracPA_Abs;
109  delete fFracPA_Pipro;
110  delete fFracNA_Tot;
111  delete fFracNA_Elas;
112  delete fFracNA_Inel;
113  delete fFracNA_CEx;
114  delete fFracNA_Abs;
115  delete fFracNA_Pipro;
116 
117  // pi+A x-section splines
118  delete fFracPipA_Tot;
119  delete fFracPipA_Elas;
120  delete fFracPipA_Inel;
121  delete fFracPipA_CEx;
122  delete fFracPipA_Abs;
123  delete fFracPipA_PiProd;
124  delete fFracPimA_Tot;
125  delete fFracPimA_Elas;
126  delete fFracPimA_Inel;
127  delete fFracPimA_CEx;
128  delete fFracPimA_Abs;
129  delete fFracPimA_PiProd;
130  delete fFracPi0A_Tot;
131  delete fFracPi0A_Elas;
132  delete fFracPi0A_Inel;
133  delete fFracPi0A_CEx;
134  delete fFracPi0A_Abs;
135  delete fFracPi0A_PiProd;
136 
137  // hN data
138  delete fhN2dXSecPP_Elas;
139  delete fhN2dXSecNP_Elas;
140  delete fhN2dXSecPipN_Elas;
141  delete fhN2dXSecPi0N_Elas;
142  delete fhN2dXSecPimN_Elas;
143  delete fhN2dXSecKpN_Elas;
144  delete fhN2dXSecKpP_Elas;
145  delete fhN2dXSecPiN_CEx;
146  delete fhN2dXSecPiN_Abs;
147  delete fhN2dXSecGamPi0P_Inelas;
148  delete fhN2dXSecGamPi0N_Inelas;
149  delete fhN2dXSecGamPipN_Inelas;
150  delete fhN2dXSecGamPimP_Inelas;
151 
152  // p/n+p/n hA x-section splines
153  delete fXSecPp_Tot;
154  delete fXSecPp_Elas;
155  delete fXSecPp_Reac;
156  delete fXSecPn_Tot;
157  delete fXSecPn_Elas;
158  delete fXSecPn_Reac;
159  delete fXSecNn_Tot;
160  delete fXSecNn_Elas;
161  delete fXSecNn_Reac;
162 
163 
164  // K+A x-section fraction splines
165  delete fFracKA_Tot;
166  delete fFracKA_Elas;
167  delete fFracKA_Inel;
168  delete fFracKA_Abs;
169 
170 }
171 //____________________________________________________________________________
173 {
174  if(fInstance == 0) {
175  LOG("INukeData", pINFO) << "INukeHadroData late initialization";
176  static INukeHadroData::Cleaner cleaner;
178  fInstance = new INukeHadroData;
179  }
180  return fInstance;
181 }
182 //____________________________________________________________________________
184 {
185 // Loads hadronic x-section data
186 
187  //-- Get the top-level directory with input hadron cross-section data
188  // (search for $GINUKEHADRONDATA or use default location)
189  string data_dir = (gSystem->Getenv("GINUKEHADRONDATA")) ?
190  string(gSystem->Getenv("GINUKEHADRONDATA")) :
191  string(gSystem->Getenv("GENIE")) + string("/data/evgen/intranuke");
192 
193  LOG("INukeData", pINFO)
194  << "Loading INTRANUKE hadron data from: " << data_dir;
195 
196  //-- Build filenames
197 
198  string datafile_NN = data_dir + "/tot_xsec/intranuke-xsections-NN.dat";
199  string datafile_pipN = data_dir + "/tot_xsec/intranuke-xsections-pi+N.dat";
200  string datafile_pi0N = data_dir + "/tot_xsec/intranuke-xsections-pi0N.dat";
201  string datafile_NA = data_dir + "/tot_xsec/intranuke-fractions-NA.dat";
202  string datafile_piA = data_dir + "/tot_xsec/intranuke-fractions-piA.dat";
203  string datafile_KA = data_dir + "/tot_xsec/intranuke-fractions-KA.dat";
204  string datafile_gamN = data_dir + "/tot_xsec/intranuke-xsections-gamN.dat";
205  string datafile_kN = data_dir + "/tot_xsec/intranuke-xsections-kaonN.dat";
206 
207  //-- Make sure that all data files are available
208 
209  assert( ! gSystem->AccessPathName(datafile_NN. c_str()) );
210  assert( ! gSystem->AccessPathName(datafile_pipN.c_str()) );
211  assert( ! gSystem->AccessPathName(datafile_pi0N.c_str()) );
212  assert( ! gSystem->AccessPathName(datafile_NA. c_str()) );
213  assert( ! gSystem->AccessPathName(datafile_piA. c_str()) );
214  assert( ! gSystem->AccessPathName(datafile_KA. c_str()) );
215  assert( ! gSystem->AccessPathName(datafile_gamN.c_str()) );
216  assert( ! gSystem->AccessPathName(datafile_kN. c_str()) );
217 
218  LOG("INukeData", pINFO) << "Found all necessary data files...";
219 
220  //-- Load data files
221 
222  TTree data_NN;
223  TTree data_pipN;
224  TTree data_pi0N;
225  TTree data_NA;
226  TTree data_piA;
227  TTree data_KA;
228  TTree data_gamN;
229  TTree data_kN;
230 
231  data_NN.ReadFile(datafile_NN.c_str(),
232  "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");
233  data_pipN.ReadFile(datafile_pipN.c_str(),
234  "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");
235  data_pi0N.ReadFile(datafile_pi0N.c_str(),
236  "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");
237  data_NA.ReadFile(datafile_NA.c_str(),
238  "ke/D:pA_tot/D:pA_elas/D:pA_inel/D:pA_cex/D:pA_abs/D:pA_pipro/D");
239  data_piA.ReadFile(datafile_piA.c_str(),
240  "ke/D:piA_tot/D:piA_elas/D:piA_inel/D:piA_cex/D:piA_np/D:piA_pp/D:piA_npp/D:piA_nnp/D:piA_2n2p/D:piA_piprod/D");
241  data_gamN.ReadFile(datafile_gamN.c_str(),
242  "ke/D:pi0p_tot/D:pipn_tot/D:pimp_tot/D:pi0n_tot/D:gamp_fs/D:gamn_fs/D:gamN_tot/D");
243  data_kN.ReadFile(datafile_kN.c_str(),
244  "ke/D:kpn_elas/D:kpp_elas/D:kp_abs/D:kpN_tot/D"); //????
245  data_KA.ReadFile(datafile_KA.c_str(),
246  "ke/D:KA_tot/D:KA_elas/D:KA_inel/D:KA_abs/D");
247 
248  LOG("INukeData", pDEBUG) << "Number of data rows in NN : " << data_NN.GetEntries();
249  LOG("INukeData", pDEBUG) << "Number of data rows in pipN : " << data_pipN.GetEntries();
250  LOG("INukeData", pDEBUG) << "Number of data rows in pi0N : " << data_pi0N.GetEntries();
251  LOG("INukeData", pDEBUG) << "Number of data rows in NA : " << data_NA.GetEntries();
252  LOG("INukeData", pDEBUG) << "Number of data rows in piA : " << data_piA.GetEntries();
253  LOG("INukeData", pDEBUG) << "Number of data rows in KA : " << data_KA.GetEntries();
254  LOG("INukeData", pDEBUG) << "Number of data rows in gamN : " << data_gamN.GetEntries();
255  LOG("INukeData", pDEBUG) << "Number of data rows in kN : " << data_kN.GetEntries();
256 
257  LOG("INukeData", pINFO) << "Done loading all x-section files...";
258 
259  //-- Build x-section splines
260 
261  // p/n+p/n hA x-section splines
262  fXSecPp_Tot = new Spline(&data_NN, "ke:pp_tot");
263  fXSecPp_Elas = new Spline(&data_NN, "ke:pp_elas");
264  fXSecPp_Reac = new Spline(&data_NN, "ke:pp_reac");
265  fXSecPn_Tot = new Spline(&data_NN, "ke:pn_tot");
266  fXSecPn_Elas = new Spline(&data_NN, "ke:pn_elas");
267  fXSecPn_Reac = new Spline(&data_NN, "ke:pn_reac");
268  fXSecNn_Tot = new Spline(&data_NN, "ke:nn_tot");
269  fXSecNn_Elas = new Spline(&data_NN, "ke:nn_elas");
270  fXSecNn_Reac = new Spline(&data_NN, "ke:nn_reac");
271 
272  // pi+n/p hA x-section splines
273  fXSecPipn_Tot = new Spline(&data_pipN, "ke:pipn_tot");
274  fXSecPipn_CEx = new Spline(&data_pipN, "ke:pipn_cex");
275  fXSecPipn_Elas = new Spline(&data_pipN, "ke:pipn_elas");
276  fXSecPipn_Reac = new Spline(&data_pipN, "ke:pipn_reac");
277  fXSecPipp_Tot = new Spline(&data_pipN, "ke:pipp_tot");
278  fXSecPipp_CEx = new Spline(&data_pipN, "ke:pipp_cex");
279  fXSecPipp_Elas = new Spline(&data_pipN, "ke:pipp_elas");
280  fXSecPipp_Reac = new Spline(&data_pipN, "ke:pipp_reac");
281  fXSecPipd_Abs = new Spline(&data_pipN, "ke:pipd_abs");
282 
283  // pi0n/p hA x-section splines
284  fXSecPi0n_Tot = new Spline(&data_pi0N, "ke:pi0n_tot");
285  fXSecPi0n_CEx = new Spline(&data_pi0N, "ke:pi0n_cex");
286  fXSecPi0n_Elas = new Spline(&data_pi0N, "ke:pi0n_elas");
287  fXSecPi0n_Reac = new Spline(&data_pi0N, "ke:pi0n_reac");
288  fXSecPi0p_Tot = new Spline(&data_pi0N, "ke:pi0p_tot");
289  fXSecPi0p_CEx = new Spline(&data_pi0N, "ke:pi0p_cex");
290  fXSecPi0p_Elas = new Spline(&data_pi0N, "ke:pi0p_elas");
291  fXSecPi0p_Reac = new Spline(&data_pi0N, "ke:pi0p_reac");
292  fXSecPi0d_Abs = new Spline(&data_pi0N, "ke:pi0d_abs");
293 
294  // K+N x-section splines
295  fXSecKpn_Elas = new Spline(&data_kN, "ke:kpn_elas");
296  fXSecKpp_Elas = new Spline(&data_kN, "ke:kpp_elas");
297  fXSecKpN_Abs = new Spline(&data_kN, "ke:kp_abs");
298  fXSecKpN_Tot = new Spline(&data_kN, "ke:kpN_tot");
299 
300  // gamma x-section splines
301  fXSecGamp_fs = new Spline(&data_gamN, "ke:gamp_fs");
302  fXSecGamn_fs = new Spline(&data_gamN, "ke:gamn_fs");
303  fXSecGamN_Tot = new Spline(&data_gamN, "ke:gamN_tot");
304 
305  // N+A x-section fraction splines
306  fFracPA_Tot = new Spline(&data_NA, "ke:pA_tot");
307  fFracPA_Elas = new Spline(&data_NA, "ke:pA_elas");
308  fFracPA_Inel = new Spline(&data_NA, "ke:pA_inel");
309  fFracPA_CEx = new Spline(&data_NA, "ke:pA_cex");
310  fFracPA_Abs = new Spline(&data_NA, "ke:pA_abs");
311  fFracPA_Pipro = new Spline(&data_NA, "ke:pA_pipro");
312  fFracNA_Tot = new Spline(&data_NA, "ke:pA_tot"); // assuming nA same as pA
313  fFracNA_Elas = new Spline(&data_NA, "ke:pA_elas");
314  fFracNA_Inel = new Spline(&data_NA, "ke:pA_inel");
315  fFracNA_CEx = new Spline(&data_NA, "ke:pA_cex");
316  fFracNA_Abs = new Spline(&data_NA, "ke:pA_abs");
317  fFracNA_Pipro = new Spline(&data_NA, "ke:pA_pipro");
318 
319  // pi+A x-section splines
320  fFracPipA_Tot = new Spline(&data_piA, "ke:piA_tot");
321  fFracPipA_Elas = new Spline(&data_piA, "ke:piA_elas");
322  fFracPipA_Inel = new Spline(&data_piA, "ke:piA_inel");
323  fFracPipA_CEx = new Spline(&data_piA, "ke:piA_cex");
324  fFracPipA_Abs = new Spline(&data_piA, "ke:piA_np+piA_pp+piA_npp+piA_nnp+piA_2n2p");
325  fFracPipA_PiProd = new Spline(&data_piA, "ke:piA_piprod");
326  fFracPimA_Tot = new Spline(&data_piA, "ke:piA_tot");
327  fFracPimA_Elas = new Spline(&data_piA, "ke:piA_elas");
328  fFracPimA_Inel = new Spline(&data_piA, "ke:piA_inel");
329  fFracPimA_CEx = new Spline(&data_piA, "ke:piA_cex");
330  fFracPimA_Abs = new Spline(&data_piA, "ke:piA_np+piA_pp+piA_npp+piA_nnp+piA_2n2p");
331  fFracPimA_PiProd = new Spline(&data_piA, "ke:piA_piprod");
332  fFracPi0A_Tot = new Spline(&data_piA, "ke:piA_tot");
333  fFracPi0A_Elas = new Spline(&data_piA, "ke:piA_elas");
334  fFracPi0A_Inel = new Spline(&data_piA, "ke:piA_inel");
335  fFracPi0A_CEx = new Spline(&data_piA, "ke:piA_cex");
336  fFracPi0A_Abs = new Spline(&data_piA, "ke:piA_np+piA_pp+piA_npp+piA_nnp+piA_2n2p");
337  fFracPi0A_PiProd = new Spline(&data_piA, "ke:piA_piprod");
338  // K+A x-section fraction splines
339  fFracKA_Tot = new Spline(&data_KA, "ke:KA_tot");
340  fFracKA_Elas = new Spline(&data_KA, "ke:KA_elas");
341  fFracKA_Inel = new Spline(&data_KA, "ke:KA_inel");
342  fFracKA_Abs = new Spline(&data_KA, "ke:KA_abs");
343  //
344  // hN stuff
345  //
346 
347 
348  // kIHNFtElas
349  // pp, nn --> read from pp/pp%.txt
350  // pn, np --> read from pp/pn%.txt
351  // pi+ N --> read from pip/pip%.txt
352  // pi0 N --> read from pip/pip%.txt
353  // pi- N --> read from pim/pim%.txt
354  // K+ N --> read from kpn/kpn%.txt
355  // K+ P --> read from kpp/kpp%.txt
356  // kIHNFtCEx
357  // pi+, pi0, pi- --> read from pie/pie%.txt (using pip+n->pi0+p data)
358  // kIHNFtAbs
359  // pi+, pi0, pi- --> read from pid2p/pid2p%.txt (using pip+D->2p data)
360  // kIHNFtInelas
361  // gamma p -> p pi0 --> read from gampi0p/%-pi0p.txt
362  // gamma p -> n pi+ --> read from gampi+n/%-pi+n.txt
363  // gamma n -> n pi0 --> read from gampi0n/%-pi0n.txt
364  // gamma n -> p pi- --> read from gampi-p/%-pi-p.txt
365 
366 
367  // kIHNFtElas, pp&nn :
368  {
369  const int hN_ppelas_nfiles = 20;
370  const int hN_ppelas_points_per_file = 21;
371  const int hN_ppelas_npoints = hN_ppelas_points_per_file * hN_ppelas_nfiles;
372 
373  double hN_ppelas_energies[hN_ppelas_nfiles] = {
374  50, 100, 150, 200, 250, 300, 350, 400, 450, 500,
375  550, 600, 650, 700, 750, 800, 850, 900, 950, 1000
376  };
377 
378  double hN_ppelas_costh [hN_ppelas_points_per_file];
379  double hN_ppelas_xsec [hN_ppelas_npoints];
380 
381  int ipoint=0;
382 
383  for(int ifile = 0; ifile < hN_ppelas_nfiles; ifile++) {
384  // build filename
385  ostringstream hN_datafile;
386  double ke = hN_ppelas_energies[ifile];
387  hN_datafile << data_dir << "/diff_ang/pp/pp" << ke << ".txt";
388  // read data
389  ReadhNFile(
390  hN_datafile.str(), ke, hN_ppelas_points_per_file,
391  ipoint, hN_ppelas_costh, hN_ppelas_xsec,2);
392  }//loop over files
393 
394  /*double hN_ppelas_costh_cond [hN_ppelas_points_per_file];
395  for (int ient = 0; ient < hN_ppelas_points_per_file; ient++) {
396  hN_ppelas_costh_cond[ient] = hN_ppelas_costh[ient];
397  }*/
398 
399  fhN2dXSecPP_Elas = new BLI2DNonUnifGrid(hN_ppelas_nfiles,hN_ppelas_points_per_file,
400  hN_ppelas_energies,hN_ppelas_costh,hN_ppelas_xsec);
401  }
402 
403  // kIHNFtElas, pn&np :
404  {
405  const int hN_npelas_nfiles = 20;
406  const int hN_npelas_points_per_file = 21;
407  const int hN_npelas_npoints = hN_npelas_points_per_file * hN_npelas_nfiles;
408 
409  double hN_npelas_energies[hN_npelas_nfiles] = {
410  50, 100, 150, 200, 250, 300, 350, 400, 450, 500,
411  550, 600, 650, 700, 750, 800, 850, 900, 950, 1000
412  };
413 
414  double hN_npelas_costh [hN_npelas_points_per_file];
415  double hN_npelas_xsec [hN_npelas_npoints];
416 
417  int ipoint=0;
418 
419  for(int ifile = 0; ifile < hN_npelas_nfiles; ifile++) {
420  // build filename
421  ostringstream hN_datafile;
422  double ke = hN_npelas_energies[ifile];
423  hN_datafile << data_dir << "/diff_ang/pn/pn" << ke << ".txt";
424  // read data
425  ReadhNFile(
426  hN_datafile.str(), ke, hN_npelas_points_per_file,
427  ipoint, hN_npelas_costh, hN_npelas_xsec,2);
428  }//loop over files
429 
430  /*double hN_npelas_costh_cond [hN_npelas_points_per_file];
431  for (int ient = 0; ient < hN_npelas_points_per_file; ient++) {
432  hN_npelas_costh_cond[ient] = hN_npelas_costh[ient];
433  }*/
434 
435  fhN2dXSecNP_Elas = new BLI2DNonUnifGrid(hN_npelas_nfiles,hN_npelas_points_per_file,
436  hN_npelas_energies,hN_npelas_costh,hN_npelas_xsec);
437  }
438 
439  // kIHNFtElas, pipN :
440  {
441  const int hN_pipNelas_nfiles = 60;
442  const int hN_pipNelas_points_per_file = 21;
443  const int hN_pipNelas_npoints = hN_pipNelas_points_per_file * hN_pipNelas_nfiles;
444 
445  double hN_pipNelas_energies[hN_pipNelas_nfiles] = {
446  10, 20, 30, 40, 50, 60, 70, 80, 90,
447  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
448  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
449  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
450  700, 740, 780, 820, 860, 900, 940, 980,
451  1020, 1060, 1100, 1140, 1180, 1220, 1260,
452  1300, 1340, 1380, 1420, 1460, 1500
453  };
454 
455  double hN_pipNelas_costh [hN_pipNelas_points_per_file];
456  double hN_pipNelas_xsec [hN_pipNelas_npoints];
457 
458  int ipoint=0;
459 
460  for(int ifile = 0; ifile < hN_pipNelas_nfiles; ifile++) {
461  // build filename
462  ostringstream hN_datafile;
463  double ke = hN_pipNelas_energies[ifile];
464  hN_datafile << data_dir << "/diff_ang/pip/pip" << ke << ".txt";
465  // read data
466  ReadhNFile(
467  hN_datafile.str(), ke, hN_pipNelas_points_per_file,
468  ipoint, hN_pipNelas_costh, hN_pipNelas_xsec,2);
469  }//loop over files
470 
471  /*double hN_pipNelas_costh_cond [hN_pipNelas_points_per_file];
472  for (int ient = 0; ient < hN_pipNelas_points_per_file; ient++) {
473  hN_pipNelas_costh_cond[ient] = hN_pipNelas_costh[ient];
474  }*/
475 
476  fhN2dXSecPipN_Elas = new BLI2DNonUnifGrid(hN_pipNelas_nfiles,hN_pipNelas_points_per_file,
477  hN_pipNelas_energies,hN_pipNelas_costh,hN_pipNelas_xsec);
478  }
479 
480  // kIHNFtElas, pi0N :
481  {
482  const int hN_pi0Nelas_nfiles = 60;
483  const int hN_pi0Nelas_points_per_file = 21;
484  const int hN_pi0Nelas_npoints = hN_pi0Nelas_points_per_file * hN_pi0Nelas_nfiles;
485 
486  double hN_pi0Nelas_energies[hN_pi0Nelas_nfiles] = {
487  10, 20, 30, 40, 50, 60, 70, 80, 90,
488  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
489  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
490  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
491  700, 740, 780, 820, 860, 900, 940, 980,
492  1020, 1060, 1100, 1140, 1180, 1220, 1260,
493  1300, 1340, 1380, 1420, 1460, 1500
494  };
495 
496  double hN_pi0Nelas_costh [hN_pi0Nelas_points_per_file];
497  double hN_pi0Nelas_xsec [hN_pi0Nelas_npoints];
498 
499  int ipoint=0;
500 
501  for(int ifile = 0; ifile < hN_pi0Nelas_nfiles; ifile++) {
502  // build filename
503  ostringstream hN_datafile;
504  double ke = hN_pi0Nelas_energies[ifile];
505  hN_datafile << data_dir << "/diff_ang/pip/pip" << ke << ".txt";
506  // read data
507  ReadhNFile(
508  hN_datafile.str(), ke, hN_pi0Nelas_points_per_file,
509  ipoint, hN_pi0Nelas_costh, hN_pi0Nelas_xsec,2);
510  }//loop over files
511 
512  /*double hN_pi0Nelas_costh_cond [hN_pi0Nelas_points_per_file];
513  for (int ient = 0; ient < hN_pi0Nelas_points_per_file; ient++) {
514  hN_pi0Nelas_costh_cond[ient] = hN_pi0Nelas_costh[ient];
515  }*/
516 
517  fhN2dXSecPi0N_Elas = new BLI2DNonUnifGrid(hN_pi0Nelas_nfiles,hN_pi0Nelas_points_per_file,
518  hN_pi0Nelas_energies,hN_pi0Nelas_costh,hN_pi0Nelas_xsec);
519  }
520 
521  // kIHNFtElas, pimN :
522  {
523  const int hN_pimNelas_nfiles = 60;
524  const int hN_pimNelas_points_per_file = 21;
525  const int hN_pimNelas_npoints = hN_pimNelas_points_per_file * hN_pimNelas_nfiles;
526 
527  double hN_pimNelas_energies[hN_pimNelas_nfiles] = {
528  10, 20, 30, 40, 50, 60, 70, 80, 90,
529  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
530  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
531  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
532  700, 740, 780, 820, 860, 900, 940, 980,
533  1020, 1060, 1100, 1140, 1180, 1220, 1260,
534  1300, 1340, 1380, 1420, 1460, 1500
535  };
536 
537  double hN_pimNelas_costh [hN_pimNelas_points_per_file];
538  double hN_pimNelas_xsec [hN_pimNelas_npoints];
539 
540  int ipoint=0;
541 
542  for(int ifile = 0; ifile < hN_pimNelas_nfiles; ifile++) {
543  // build filename
544  ostringstream hN_datafile;
545  double ke = hN_pimNelas_energies[ifile];
546  hN_datafile << data_dir << "/diff_ang/pim/pim" << ke << ".txt";
547  // read data
548  ReadhNFile(
549  hN_datafile.str(), ke, hN_pimNelas_points_per_file,
550  ipoint, hN_pimNelas_costh, hN_pimNelas_xsec,2);
551  }//loop over files
552 
553  /*double hN_pimNelas_costh_cond [hN_pimNelas_points_per_file];
554  for (int ient = 0; ient < hN_pimNelas_points_per_file; ient++) {
555  hN_pimNelas_costh_cond[ient] = hN_pimNelas_costh[ient];
556  }*/
557 
558  fhN2dXSecPimN_Elas = new BLI2DNonUnifGrid(hN_pimNelas_nfiles,hN_pimNelas_points_per_file,
559  hN_pimNelas_energies,hN_pimNelas_costh,hN_pimNelas_xsec);
560  }
561 
562  // kIHNFtElas, kpn :
563  {
564  const int hN_kpNelas_nfiles = 18;
565  const int hN_kpNelas_points_per_file = 37;
566  const int hN_kpNelas_npoints = hN_kpNelas_points_per_file * hN_kpNelas_nfiles;
567 
568  double hN_kpNelas_energies[hN_kpNelas_nfiles] = {
569  100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
570  1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800
571  };
572 
573  double hN_kpNelas_costh [hN_kpNelas_points_per_file];
574  double hN_kpNelas_xsec [hN_kpNelas_npoints];
575 
576  int ipoint=0;
577 
578  for(int ifile = 0; ifile < hN_kpNelas_nfiles; ifile++) {
579  // build filename
580  ostringstream hN_datafile;
581  double ke = hN_kpNelas_energies[ifile];
582  hN_datafile << data_dir << "/diff_ang/kpn/kpn" << ke << ".txt";
583  // read data
584  ReadhNFile(
585  hN_datafile.str(), ke, hN_kpNelas_points_per_file,
586  ipoint, hN_kpNelas_costh, hN_kpNelas_xsec,2);
587  }//loop over files
588 
589  /*double hN_kpNelas_costh_cond [hN_kpNelas_points_per_file];
590  for (int ient = 0; ient < hN_kpNelas_points_per_file; ient++) {
591  hN_kpNelas_costh_cond[ient] = hN_kpNelas_costh[ient];
592  }*/
593 
594  fhN2dXSecKpN_Elas = new BLI2DNonUnifGrid(hN_kpNelas_nfiles,hN_kpNelas_points_per_file,
595  hN_kpNelas_energies,hN_kpNelas_costh,hN_kpNelas_xsec);
596  }
597 
598  // kIHNFtElas, kpp :
599  {
600  const int hN_kpPelas_nfiles = 18;
601  const int hN_kpPelas_points_per_file = 37;
602  const int hN_kpPelas_npoints = hN_kpPelas_points_per_file * hN_kpPelas_nfiles;
603 
604  double hN_kpPelas_energies[hN_kpPelas_nfiles] = {
605  100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
606  1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800
607  };
608 
609  double hN_kpPelas_costh [hN_kpPelas_points_per_file];
610  double hN_kpPelas_xsec [hN_kpPelas_npoints];
611 
612  int ipoint=0;
613 
614  for(int ifile = 0; ifile < hN_kpPelas_nfiles; ifile++) {
615  // build filename
616  ostringstream hN_datafile;
617  double ke = hN_kpPelas_energies[ifile];
618  hN_datafile << data_dir << "/diff_ang/kpp/kpp" << ke << ".txt";
619  // read data
620  ReadhNFile(
621  hN_datafile.str(), ke, hN_kpPelas_points_per_file,
622  ipoint, hN_kpPelas_costh, hN_kpPelas_xsec,2);
623  }//loop over files
624 
625  /*double hN_kpPelas_costh_cond [hN_kpPelas_points_per_file];
626  for (int ient = 0; ient < hN_kpPelas_points_per_file; ient++) {
627  hN_kpPelas_costh_cond[ient] = hN_kpPelas_costh[ient];
628  }*/
629 
630  fhN2dXSecKpP_Elas = new BLI2DNonUnifGrid(hN_kpPelas_nfiles,hN_kpPelas_points_per_file,
631  hN_kpPelas_energies,hN_kpPelas_costh,hN_kpPelas_xsec);
632  }
633 
634  // kIHNFtCEx, (pi+, pi0, pi-) N
635  {
636  const int hN_piNcex_nfiles = 60;
637  const int hN_piNcex_points_per_file = 21;
638  const int hN_piNcex_npoints = hN_piNcex_points_per_file * hN_piNcex_nfiles;
639 
640  double hN_piNcex_energies[hN_piNcex_nfiles] = {
641  10, 20, 30, 40, 50, 60, 70, 80, 90,
642  100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
643  200, 210, 220, 230, 240, 250, 260, 270, 280, 290,
644  300, 340, 380, 420, 460, 500, 540, 580, 620, 660,
645  700, 740, 780, 820, 860, 900, 940, 980,
646  1020, 1060, 1100, 1140, 1180, 1220, 1260,
647  1300, 1340, 1380, 1420, 1460, 1500
648  };
649 
650  double hN_piNcex_costh [hN_piNcex_points_per_file];
651  double hN_piNcex_xsec [hN_piNcex_npoints];
652 
653  int ipoint=0;
654 
655  for(int ifile = 0; ifile < hN_piNcex_nfiles; ifile++) {
656  // build filename
657  ostringstream hN_datafile;
658  double ke = hN_piNcex_energies[ifile];
659  hN_datafile << data_dir << "/diff_ang/pie/pie" << ke << ".txt";
660  // read data
661  ReadhNFile(
662  hN_datafile.str(), ke, hN_piNcex_points_per_file,
663  ipoint, hN_piNcex_costh, hN_piNcex_xsec,2);
664  }//loop over files
665 
666  /*double hN_piNcex_costh_cond [hN_piNcex_points_per_file];
667  for (int ient = 0; ient < hN_piNcex_points_per_file; ient++) {
668  hN_piNcex_costh_cond[ient] = hN_piNcex_costh[ient];
669  }*/
670 
671  fhN2dXSecPiN_CEx = new BLI2DNonUnifGrid(hN_piNcex_nfiles,hN_piNcex_points_per_file,
672  hN_piNcex_energies,hN_piNcex_costh,hN_piNcex_xsec);
673  }
674 
675  // kIHNFtAbs, (pi+, pi0, pi-) N
676  {
677  const int hN_piNabs_nfiles = 19;
678  const int hN_piNabs_points_per_file = 21;
679  const int hN_piNabs_npoints = hN_piNabs_points_per_file * hN_piNabs_nfiles;
680 
681  double hN_piNabs_energies[hN_piNabs_nfiles] = {
682  50, 75, 100, 125, 150, 175, 200, 225, 250, 275,
683  300, 325, 350, 375, 400, 425, 450, 475, 500
684  };
685 
686  double hN_piNabs_costh [hN_piNabs_points_per_file];
687  double hN_piNabs_xsec [hN_piNabs_npoints];
688 
689  int ipoint=0;
690 
691  for(int ifile = 0; ifile < hN_piNabs_nfiles; ifile++) {
692  // build filename
693  ostringstream hN_datafile;
694  double ke = hN_piNabs_energies[ifile];
695  hN_datafile << data_dir << "/diff_ang/pid2p/pid2p" << ke << ".txt";
696  // read data
697  ReadhNFile(
698  hN_datafile.str(), ke, hN_piNabs_points_per_file,
699  ipoint, hN_piNabs_costh, hN_piNabs_xsec,2);
700  }//loop over files
701 
702  /*double hN_piNabs_costh_cond [hN_piNabs_points_per_file];
703  for (int ient = 0; ient < hN_piNabs_points_per_file; ient++) {
704  hN_piNabs_costh_cond[ient] = hN_piNabs_costh[ient];
705  }*/
706 
707  fhN2dXSecPiN_Abs = new BLI2DNonUnifGrid(hN_piNabs_nfiles,hN_piNabs_points_per_file,
708  hN_piNabs_energies,hN_piNabs_costh,hN_piNabs_xsec);
709  }
710 
711  // kIHNFtInelas, gamma p -> p pi0
712  {
713  const int hN_gampi0pInelas_nfiles = 29;
714  const int hN_gampi0pInelas_points_per_file = 37;
715  const int hN_gampi0pInelas_npoints = hN_gampi0pInelas_points_per_file * hN_gampi0pInelas_nfiles;
716 
717  double hN_gampi0pInelas_energies[hN_gampi0pInelas_nfiles] = {
718  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
719  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
720  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
721  };
722 
723  double hN_gampi0pInelas_costh [hN_gampi0pInelas_points_per_file];
724  double hN_gampi0pInelas_xsec [hN_gampi0pInelas_npoints];
725 
726  int ipoint=0;
727 
728  for(int ifile = 0; ifile < hN_gampi0pInelas_nfiles; ifile++) {
729  // build filename
730  ostringstream hN_datafile;
731  double ke = hN_gampi0pInelas_energies[ifile];
732  hN_datafile << data_dir << "/diff_ang/gampi0p/" << ke << "-pi0p.txt";
733  // read data
734  ReadhNFile(
735  hN_datafile.str(), ke, hN_gampi0pInelas_points_per_file,
736  ipoint, hN_gampi0pInelas_costh, hN_gampi0pInelas_xsec,3);
737  }//loop over files
738 
739  /*double hN_gampi0pInelas_costh_cond [hN_gampi0pInelas_points_per_file];
740  for (int ient = 0; ient < hN_gampi0pInelas_points_per_file; ient++) {
741  hN_gampi0pInelas_costh_cond[ient] = hN_gampi0pInelas_costh[ient];
742  }*/
743 
744  fhN2dXSecGamPi0P_Inelas = new BLI2DNonUnifGrid(hN_gampi0pInelas_nfiles,hN_gampi0pInelas_points_per_file,
745  hN_gampi0pInelas_energies,hN_gampi0pInelas_costh,hN_gampi0pInelas_xsec);
746  }
747 
748  // kIHNFtInelas, gamma n -> n pi0
749  {
750  const int hN_gampi0nInelas_nfiles = 29;
751  const int hN_gampi0nInelas_points_per_file = 37;
752  const int hN_gampi0nInelas_npoints = hN_gampi0nInelas_points_per_file * hN_gampi0nInelas_nfiles;
753 
754  double hN_gampi0nInelas_energies[hN_gampi0nInelas_nfiles] = {
755  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
756  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
757  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
758  };
759 
760  double hN_gampi0nInelas_costh [hN_gampi0nInelas_points_per_file];
761  double hN_gampi0nInelas_xsec [hN_gampi0nInelas_npoints];
762  int ipoint=0;
763 
764  for(int ifile = 0; ifile < hN_gampi0nInelas_nfiles; ifile++) {
765  // build filename
766  ostringstream hN_datafile;
767  double ke = hN_gampi0nInelas_energies[ifile];
768  hN_datafile << data_dir << "/diff_ang/gampi0n/" << ke << "-pi0n.txt";
769  // read data
770  ReadhNFile(
771  hN_datafile.str(), ke, hN_gampi0nInelas_points_per_file,
772  ipoint, hN_gampi0nInelas_costh, hN_gampi0nInelas_xsec,3);
773  }//loop over files
774 
775  /*double hN_gampi0nInelas_costh_cond [hN_gampi0nInelas_points_per_file];
776  for (int ient = 0; ient < hN_gampi0nInelas_points_per_file; ient++) {
777  hN_gampi0nInelas_costh_cond[ient] = hN_gampi0nInelas_costh[ient];
778  }*/
779 
780  fhN2dXSecGamPi0N_Inelas = new BLI2DNonUnifGrid(hN_gampi0nInelas_nfiles,hN_gampi0nInelas_points_per_file,
781  hN_gampi0nInelas_energies,hN_gampi0nInelas_costh,hN_gampi0nInelas_xsec);
782  }
783 
784  // kIHNFtInelas, gamma p -> n pi+
785  {
786  const int hN_gampipnInelas_nfiles = 29;
787  const int hN_gampipnInelas_points_per_file = 37;
788  const int hN_gampipnInelas_npoints = hN_gampipnInelas_points_per_file * hN_gampipnInelas_nfiles;
789 
790  double hN_gampipnInelas_energies[hN_gampipnInelas_nfiles] = {
791  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
792  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
793  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
794  };
795 
796  double hN_gampipnInelas_costh [hN_gampipnInelas_points_per_file];
797  double hN_gampipnInelas_xsec [hN_gampipnInelas_npoints];
798 
799  int ipoint=0;
800 
801  for(int ifile = 0; ifile < hN_gampipnInelas_nfiles; ifile++) {
802  // build filename
803  ostringstream hN_datafile;
804  double ke = hN_gampipnInelas_energies[ifile];
805  hN_datafile << data_dir << "/diff_ang/gampi+n/" << ke << "-pi+n.txt";
806  // read data
807  ReadhNFile(
808  hN_datafile.str(), ke, hN_gampipnInelas_points_per_file,
809  ipoint, hN_gampipnInelas_costh, hN_gampipnInelas_xsec,3);
810  }//loop over files
811 
812  /*double hN_gampipnInelas_costh_cond [hN_gampipnInelas_points_per_file];
813  for (int ient = 0; ient < hN_gampipnInelas_points_per_file; ient++) {
814  hN_gampipnInelas_costh_cond[ient] = hN_gampipnInelas_costh[ient];
815  }*/
816 
817  fhN2dXSecGamPipN_Inelas = new BLI2DNonUnifGrid(hN_gampipnInelas_nfiles,hN_gampipnInelas_points_per_file,
818  hN_gampipnInelas_energies,hN_gampipnInelas_costh,hN_gampipnInelas_xsec);
819  }
820 
821  // kIHNFtInelas, gamma n -> p pi-
822  {
823  const int hN_gampimpInelas_nfiles = 29;
824  const int hN_gampimpInelas_points_per_file = 37;
825  const int hN_gampimpInelas_npoints = hN_gampimpInelas_points_per_file * hN_gampimpInelas_nfiles;
826 
827  double hN_gampimpInelas_energies[hN_gampimpInelas_nfiles] = {
828  160, 180, 200, 220, 240, 260, 280, 300, 320, 340,
829  360, 380, 400, 450, 500, 550, 600, 650, 700, 750,
830  800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200
831  };
832 
833  double hN_gampimpInelas_costh [hN_gampimpInelas_points_per_file];
834  double hN_gampimpInelas_xsec [hN_gampimpInelas_npoints];
835 
836  int ipoint=0;
837 
838  for(int ifile = 0; ifile < hN_gampimpInelas_nfiles; ifile++) {
839  // build filename
840  ostringstream hN_datafile;
841  double ke = hN_gampimpInelas_energies[ifile];
842  hN_datafile << data_dir << "/diff_ang/gampi-p/" << ke << "-pi-p.txt";
843  // read data
844  ReadhNFile(
845  hN_datafile.str(), ke, hN_gampimpInelas_points_per_file,
846  ipoint, hN_gampimpInelas_costh, hN_gampimpInelas_xsec,3);
847  }//loop over files
848 
849  /*double hN_gampimpInelas_costh_cond [hN_gampimpInelas_points_per_file];
850  for (int ient = 0; ient < hN_gampimpInelas_points_per_file; ient++) {
851  hN_gampimpInelas_costh_cond[ient] = hN_gampimpInelas_costh[ient];
852  }*/
853 
854  fhN2dXSecGamPimP_Inelas = new BLI2DNonUnifGrid(hN_gampimpInelas_nfiles,hN_gampimpInelas_points_per_file,
855  hN_gampimpInelas_energies,hN_gampimpInelas_costh,hN_gampimpInelas_xsec);
856  }
857 
858  LOG("INukeData", pINFO) << "Done building x-section splines...";
859 }
860 //____________________________________________________________________________
862  string filename, double ke, int npoints, int & curr_point,
863  double * costh_array, double * xsec_array, int cols)
864 {
865  // open
866  std::ifstream hN_stream(filename.c_str(), ios::in);
867  if(!hN_stream.good()) {
868  LOG("INukeData", pERROR)
869  << "Error reading INTRANUKE/hN data from: " << filename;
870  return;
871  }
872 
873  if(cols<2) {
874  LOG("INukeData", pERROR)
875  << "Error reading INTRANUKE/hN data from: " << filename;
876  LOG("INukeData", pERROR)
877  << "Too few columns: " << cols;
878  return;
879  }
880 
881  LOG("INukeData", pINFO)
882  << "Reading INTRANUKE/hN data from: " << filename;
883 
884  // skip initial comments
885  char cbuf[501];
886  hN_stream.getline(cbuf,400);
887  hN_stream.getline(cbuf,400);
888  hN_stream.getline(cbuf,400);
889 
890  // read
891  double angle = 0;
892  double xsec = 0;
893  double trash = 0;
894 
895  for(int ip = 0; ip < npoints; ip++) {
896  hN_stream >> angle >> xsec;
897 
898  for(int ic = 0; ic < (cols-2); ic++) {
899  hN_stream >> trash;
900  }
901 
902  LOG("INukeData", pDEBUG)
903  << "Adding data point: (KE = " << ke << " MeV, angle = "
904  << angle << ", sigma = " << xsec << " mbarn)";
905  costh_array[ip] = TMath::Cos(angle*kPi/180.);
906  xsec_array [curr_point] = xsec;
907  curr_point++;
908  }
909 }
910 //____________________________________________________________________________
912  int hpdgc, int tgtpdgc, int nppdgc, INukeFateHN_t fate, double ke, double costh) const
913 {
914 // inputs
915 // fate : h+N fate code
916 // hpdgc : h PDG code
917 // tgtpdgc : N PDG code
918 // nppdgc : product N PDG code
919 // ke : kinetic energy (MeV)
920 // costh : cos(scattering angle)
921 // returns
922 // xsec : mbarn
923 
924  double ke_eval = ke;
925  double costh_eval = costh;
926 
927  costh_eval = TMath::Min(costh, 1.);
928  costh_eval = TMath::Max(costh_eval, -1.);
929 
930  if(fate==kIHNFtElas) {
931 
932  if( (hpdgc==kPdgProton && tgtpdgc==kPdgProton) ||
933  (hpdgc==kPdgNeutron && tgtpdgc==kPdgNeutron) )
934  {
935  ke_eval = TMath::Min(ke_eval, 999.);
936  ke_eval = TMath::Max(ke_eval, 50.);
937  return fhN2dXSecPP_Elas->Evaluate(ke_eval, costh_eval);
938  }
939  else
940  if( (hpdgc==kPdgProton && tgtpdgc==kPdgNeutron) ||
941  (hpdgc==kPdgNeutron && tgtpdgc==kPdgProton) )
942  {
943  ke_eval = TMath::Min(ke_eval, 999.);
944  ke_eval = TMath::Max(ke_eval, 50.);
945  return fhN2dXSecNP_Elas->Evaluate(ke_eval, costh_eval);
946  }
947  else
948  if(hpdgc==kPdgPiP)
949  {
950  ke_eval = TMath::Min(ke_eval, 1499.);
951  ke_eval = TMath::Max(ke_eval, 10.);
952  return fhN2dXSecPipN_Elas->Evaluate(ke_eval, costh_eval);
953  }
954  else
955  if(hpdgc==kPdgPi0)
956  {
957  ke_eval = TMath::Min(ke_eval, 1499.);
958  ke_eval = TMath::Max(ke_eval, 10.);
959  return fhN2dXSecPi0N_Elas->Evaluate(ke_eval, costh_eval);
960  }
961  else
962  if(hpdgc==kPdgPiM)
963  {
964  ke_eval = TMath::Min(ke_eval, 1499.);
965  ke_eval = TMath::Max(ke_eval, 10.);
966  return fhN2dXSecPimN_Elas->Evaluate(ke_eval, costh_eval);
967  }
968  else
969  if(hpdgc==kPdgKP && tgtpdgc==kPdgNeutron)
970  {
971  ke_eval = TMath::Min(ke_eval, 1799.);
972  ke_eval = TMath::Max(ke_eval, 100.);
973  return fhN2dXSecKpN_Elas->Evaluate(ke_eval, costh_eval);
974  }
975  else
976  if(hpdgc==kPdgKP && tgtpdgc==kPdgProton)
977  {
978  ke_eval = TMath::Min(ke_eval, 1799.);
979  ke_eval = TMath::Max(ke_eval, 100.);
980  return fhN2dXSecKpP_Elas->Evaluate(ke_eval, costh_eval);
981  }
982  }
983 
984  else if(fate == kIHNFtCEx) {
985  if( (hpdgc==kPdgPiP || hpdgc==kPdgPi0 || hpdgc==kPdgPiM) &&
986  (tgtpdgc==kPdgProton || tgtpdgc==kPdgNeutron) )
987  {
988  ke_eval = TMath::Min(ke_eval, 1499.);
989  ke_eval = TMath::Max(ke_eval, 10.);
990  return fhN2dXSecPiN_CEx->Evaluate(ke_eval, costh_eval);
991  }
992  else if( (hpdgc == kPdgProton && tgtpdgc == kPdgProton) ||
993  (hpdgc == kPdgNeutron && tgtpdgc == kPdgNeutron) )
994  {
995  LOG("INukeData", pWARN) << "Inelastic pp does not exist!";
996  ke_eval = TMath::Min(ke_eval, 999.);
997  ke_eval = TMath::Max(ke_eval, 50.);
998  return fhN2dXSecPP_Elas->Evaluate(ke_eval, costh_eval);
999  }
1000  else if( (hpdgc == kPdgProton && tgtpdgc == kPdgNeutron) ||
1001  (hpdgc == kPdgNeutron && tgtpdgc == kPdgProton) )
1002  {
1003  ke_eval = TMath::Min(ke_eval, 999.);
1004  ke_eval = TMath::Max(ke_eval, 50.);
1005  return fhN2dXSecNP_Elas->Evaluate(ke_eval, costh_eval);
1006  }
1007  }
1008 
1009  else if(fate == kIHNFtAbs) {
1010  if( (hpdgc==kPdgPiP || hpdgc==kPdgPi0 || hpdgc==kPdgPiM) &&
1011  (tgtpdgc==kPdgProton || tgtpdgc==kPdgNeutron) )
1012  {
1013  ke_eval = TMath::Min(ke_eval, 499.);
1014  ke_eval = TMath::Max(ke_eval, 50.);
1015  return fhN2dXSecPiN_Abs->Evaluate(ke_eval, costh_eval);
1016  }
1017  if(hpdgc==kPdgKP) return 1.; //isotropic since no data ???
1018  }
1019 
1020  else if(fate == kIHNFtInelas) {
1021  if( hpdgc==kPdgGamma && tgtpdgc==kPdgProton &&nppdgc==kPdgProton )
1022  {
1023  ke_eval = TMath::Min(ke_eval, 1199.);
1024  ke_eval = TMath::Max(ke_eval, 160.);
1025  return fhN2dXSecGamPi0P_Inelas->Evaluate(ke_eval, costh_eval);
1026  }
1027  else
1028  if( hpdgc==kPdgGamma && tgtpdgc==kPdgProton && nppdgc==kPdgNeutron )
1029  {
1030  ke_eval = TMath::Min(ke_eval, 1199.);
1031  ke_eval = TMath::Max(ke_eval, 160.);
1032  return fhN2dXSecGamPipN_Inelas->Evaluate(ke_eval, costh_eval);
1033  }
1034  else
1035  if( hpdgc==kPdgGamma && tgtpdgc==kPdgNeutron && nppdgc==kPdgProton )
1036  {
1037  ke_eval = TMath::Min(ke_eval, 1199.);
1038  ke_eval = TMath::Max(ke_eval, 160.);
1039  return fhN2dXSecGamPimP_Inelas->Evaluate(ke_eval, costh_eval);
1040  }
1041  else
1042  if( hpdgc==kPdgGamma && tgtpdgc==kPdgNeutron && nppdgc==kPdgNeutron )
1043  {
1044  ke_eval = TMath::Min(ke_eval, 1199.);
1045  ke_eval = TMath::Max(ke_eval, 160.);
1046  return fhN2dXSecGamPi0N_Inelas->Evaluate(ke_eval, costh_eval);
1047  }
1048  }
1049 
1050  return 0;
1051 }
1052 //____________________________________________________________________________
1053 double INukeHadroData::Frac(int hpdgc, INukeFateHA_t fate, double ke) const
1054 {
1055 // return the x-section fraction for the input fate for the particle with the input pdg
1056 // code at the input kinetic energy
1057 //
1058  ke = TMath::Max(fMinKinEnergy, ke);
1059  ke = TMath::Min(fMaxKinEnergyHA, ke);
1060 
1061  LOG("INukeData", pDEBUG) << "Querying hA cross section at ke = " << ke;
1062 
1063  if(hpdgc == kPdgProton) {
1064  /* handle protons */
1065  if (fate == kIHAFtCEx ) return TMath::Max(0., fFracPA_CEx -> Evaluate (ke));
1066  else if (fate == kIHAFtElas ) return TMath::Max(0., fFracPA_Elas -> Evaluate (ke));
1067  else if (fate == kIHAFtInelas ) return TMath::Max(0., fFracPA_Inel -> Evaluate (ke));
1068  else if (fate == kIHAFtAbs ) return TMath::Max(0., fFracPA_Abs -> Evaluate (ke));
1069  else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPA_Pipro -> Evaluate (ke));
1070  else {
1071  LOG("INukeData", pWARN)
1072  << "Protons don't have this fate: " << INukeHadroFates::AsString(fate);
1073  return 0;
1074  }
1075 
1076  } else if (hpdgc == kPdgNeutron) {
1077  /* handle neutrons */
1078  if (fate == kIHAFtCEx ) return TMath::Max(0., fFracNA_CEx -> Evaluate (ke));
1079  else if (fate == kIHAFtElas ) return TMath::Max(0., fFracNA_Elas -> Evaluate (ke));
1080  else if (fate == kIHAFtInelas ) return TMath::Max(0., fFracNA_Inel -> Evaluate (ke));
1081  else if (fate == kIHAFtAbs ) return TMath::Max(0., fFracNA_Abs -> Evaluate (ke));
1082  else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracNA_Pipro -> Evaluate (ke));
1083  else {
1084  LOG("INukeData", pWARN)
1085  << "Neutrons don't have this fate: " << INukeHadroFates::AsString(fate);
1086  return 0;
1087  }
1088 
1089  } else if (hpdgc == kPdgPiP) {
1090  /* handle pi+ */
1091  if (fate == kIHAFtCEx ) return TMath::Max(0., fFracPipA_CEx -> Evaluate (ke));
1092  else if (fate == kIHAFtElas ) return TMath::Max(0., fFracPipA_Elas -> Evaluate (ke));
1093  else if (fate == kIHAFtInelas ) return TMath::Max(0., fFracPipA_Inel -> Evaluate (ke));
1094  else if (fate == kIHAFtAbs ) return TMath::Max(0., fFracPipA_Abs -> Evaluate (ke));
1095  // else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPipA_Pipro -> Evaluate (ke));
1096  else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPipA_PiProd -> Evaluate (ke));
1097  else if (fate == kIHAFtPiProd) return TMath::Max(0., fFracPipA_PiProd -> Evaluate (ke));
1098  else {
1099  LOG("INukeData", pWARN)
1100  << "Pi+'s don't have this fate: " << INukeHadroFates::AsString(fate);
1101  return 0;
1102  }
1103 
1104  } else if (hpdgc == kPdgPiM) {
1105  /* handle pi- */
1106  if (fate == kIHAFtCEx ) return TMath::Max(0., fFracPimA_CEx -> Evaluate (ke));
1107  else if (fate == kIHAFtElas ) return TMath::Max(0., fFracPimA_Elas -> Evaluate (ke));
1108  else if (fate == kIHAFtInelas ) return TMath::Max(0., fFracPimA_Inel -> Evaluate (ke));
1109  else if (fate == kIHAFtAbs ) return TMath::Max(0., fFracPimA_Abs -> Evaluate (ke));
1110  // else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPimA_Pipro -> Evaluate (ke));
1111  else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPimA_PiProd -> Evaluate (ke));
1112  else if (fate == kIHAFtPiProd) return TMath::Max(0., fFracPimA_PiProd -> Evaluate (ke));
1113  else {
1114  LOG("INukeData", pWARN)
1115  << "Pi-'s don't have this fate: " << INukeHadroFates::AsString(fate);
1116  return 0;
1117  }
1118 
1119  } else if (hpdgc == kPdgPi0) {
1120  /* handle pi0 */
1121  if (fate == kIHAFtCEx ) return TMath::Max(0., fFracPi0A_CEx -> Evaluate (ke));
1122  else if (fate == kIHAFtElas ) return TMath::Max(0., fFracPi0A_Elas -> Evaluate (ke));
1123  else if (fate == kIHAFtInelas ) return TMath::Max(0., fFracPi0A_Inel -> Evaluate (ke));
1124  else if (fate == kIHAFtAbs ) return TMath::Max(0., fFracPi0A_Abs -> Evaluate (ke));
1125  // else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPi0A_Pipro -> Evaluate (ke));
1126  else if (fate == kIHAFtPiProd ) return TMath::Max(0., fFracPi0A_PiProd -> Evaluate (ke));
1127  else if (fate == kIHAFtPiProd) return TMath::Max(0., fFracPi0A_PiProd -> Evaluate (ke));
1128  else {
1129  LOG("INukeData", pWARN)
1130  << "Pi0's don't have this fate: " << INukeHadroFates::AsString(fate);
1131  return 0;
1132  }
1133  } else if (hpdgc == kPdgKP) {
1134  /* handle K+ */
1135  if (fate == kIHAFtInelas ) return TMath::Max(0., fFracKA_Inel -> Evaluate (ke));
1136  else if (fate == kIHAFtAbs ) return TMath::Max(0., fFracKA_Abs -> Evaluate (ke));
1137  // else if (fate == kIHAFtElas ) return TMath::Max(0., fFracKA_Elas -> Evaluate (ke));
1138  else {
1139  LOG("INukeData", pWARN)
1140  << "K+'s don't have this fate: " << INukeHadroFates::AsString(fate);
1141  return 0;
1142  }
1143  }
1144  LOG("INukeData", pWARN)
1145  << "Can't handle particles with pdg code = " << hpdgc;
1146 
1147  return 0;
1148 }
1149 //____________________________________________________________________________
1150 double INukeHadroData::XSec(int hpdgc, INukeFateHN_t fate, double ke, int targA, int targZ) const
1151 {
1152 // return the x-section for the input fate for the particle with the input pdg
1153 // code at the input kinetic energy
1154 //
1155  ke = TMath::Max(fMinKinEnergy, ke);
1156  ke = TMath::Min(fMaxKinEnergyHN, ke);
1157 
1158  LOG("INukeData", pDEBUG) << "Querying hN cross section at ke = " << ke;
1159 
1160  double xsec=0;
1161 
1162  if (hpdgc == kPdgPiP) {
1163  /* handle pi+ */
1164  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecPipp_CEx -> Evaluate(ke)) * targZ;
1165  xsec+= TMath::Max(0., fXSecPipn_CEx -> Evaluate(ke)) * (targA-targZ);
1166  return xsec;}
1167  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPipp_Elas -> Evaluate(ke)) * targZ;
1168  xsec+= TMath::Max(0., fXSecPipn_Elas -> Evaluate(ke)) * (targA-targZ);
1169  return xsec;}
1170  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPipp_Reac -> Evaluate(ke)) * targZ;
1171  xsec+= TMath::Max(0., fXSecPipn_Reac -> Evaluate(ke)) * (targA-targZ);
1172  return xsec;}
1173  else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecPipd_Abs -> Evaluate(ke)) * targA;
1174  return xsec;}
1175  else {
1176  LOG("INukeData", pWARN)
1177  << "Pi+'s don't have this fate: " << INukeHadroFates::AsString(fate);
1178  return 0;
1179  }
1180 
1181  } else if (hpdgc == kPdgPiM) {
1182  /* handle pi- */
1183  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecPipn_CEx -> Evaluate(ke)) * targZ;
1184  xsec+= TMath::Max(0., fXSecPipp_CEx -> Evaluate(ke)) * (targA-targZ);
1185  return xsec;}
1186  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPipn_Elas -> Evaluate(ke)) * targZ;
1187  xsec+= TMath::Max(0., fXSecPipp_Elas -> Evaluate(ke)) * (targA-targZ);
1188  return xsec;}
1189  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPipn_Reac -> Evaluate(ke)) * targZ;
1190  xsec+= TMath::Max(0., fXSecPipp_Reac -> Evaluate(ke)) * (targA-targZ);
1191  return xsec;}
1192  else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecPipd_Abs -> Evaluate(ke)) * targA;
1193  return xsec;}
1194  else {
1195  LOG("INukeData", pWARN)
1196  << "Pi-'s don't have this fate: " << INukeHadroFates::AsString(fate);
1197  return 0;
1198  }
1199 
1200  } else if (hpdgc == kPdgPi0) {
1201  /* handle pi0 */
1202  if (fate == kIHNFtCEx ) {xsec = TMath::Max(0., fXSecPi0p_CEx -> Evaluate(ke)) * targZ;
1203  xsec+= TMath::Max(0., fXSecPi0n_CEx -> Evaluate(ke)) * (targA-targZ);
1204  return xsec;}
1205  else if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPi0p_Elas -> Evaluate(ke)) * targZ;
1206  xsec+= TMath::Max(0., fXSecPi0n_Elas -> Evaluate(ke)) * (targA-targZ);
1207  return xsec;}
1208  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPi0p_Reac -> Evaluate(ke)) * targZ;
1209  xsec+= TMath::Max(0., fXSecPi0n_Reac -> Evaluate(ke)) * (targA-targZ);
1210  return xsec;}
1211  else if (fate == kIHNFtAbs ) {xsec = TMath::Max(0., fXSecPi0d_Abs -> Evaluate(ke)) * targA;
1212  return xsec;}
1213  else {
1214  LOG("INukeData", pWARN)
1215  << "Pi0's don't have this fate: " << INukeHadroFates::AsString(fate);
1216  return 0;
1217  }
1218 
1219  } else if (hpdgc == kPdgProton) {
1220  /* handle protons */
1221  if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPp_Elas -> Evaluate(ke)) * targZ;
1222  xsec+= TMath::Max(0., fXSecPn_Elas -> Evaluate(ke)) * (targA-targZ);
1223  return xsec;}
1224  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPp_Reac -> Evaluate(ke)) * targZ;
1225  xsec+= TMath::Max(0., fXSecPn_Reac -> Evaluate(ke)) * (targA-targZ);
1226  return xsec;}
1227  else {
1228  LOG("INukeData", pWARN)
1229  << "Protons don't have this fate: " << INukeHadroFates::AsString(fate);
1230  return 0;
1231  }
1232 
1233  } else if (hpdgc == kPdgNeutron) {
1234  /* handle protons */
1235  if (fate == kIHNFtElas ) {xsec = TMath::Max(0., fXSecPn_Elas -> Evaluate(ke)) * targZ;
1236  xsec+= TMath::Max(0., fXSecNn_Elas -> Evaluate(ke)) * (targA-targZ);
1237  return xsec;}
1238  else if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecPn_Reac -> Evaluate(ke)) * targZ;
1239  xsec+= TMath::Max(0., fXSecNn_Reac -> Evaluate(ke)) * (targA-targZ);
1240  return xsec;}
1241  else {
1242  LOG("INukeData", pWARN)
1243  << "Neutrons don't have this fate: " << INukeHadroFates::AsString(fate);
1244  return 0;
1245  }
1246  /* } else if (hpdgc == kPdgGamma) {
1247  / * handle gamma * /
1248  if (fate == kIHNFtInelas) {xsec = TMath::Max(0., fXSecGamp_fs -> Evaluate(ke)) * targZ;
1249  xsec+= TMath::Max(0., fXSecGamn_fs -> Evaluate(ke)) * (targA-targZ);
1250  return xsec;}
1251  else {
1252  LOG("INukeData", pWARN)
1253  << "Gamma's don't have this fate: " << INukeHadroFates::AsString(fate);
1254  return 0;
1255  }*/
1256  }
1257  LOG("INukeData", pWARN)
1258  << "Can't handle particles with pdg code = " << hpdgc;
1259 
1260  return 0;
1261 }
1262 
1263 double INukeHadroData::Frac(int hpdgc, INukeFateHN_t fate, double ke, int targA, int targZ) const
1264 {
1265 // return the x-section fraction for the input fate for the particle with the
1266 // input pdg code at the input kinetic energy
1267 
1268  ke = TMath::Max(fMinKinEnergy, ke);
1269  ke = TMath::Min(fMaxKinEnergyHN, ke);
1270 
1271  // get x-section
1272  double xsec = this->XSec(hpdgc,fate,ke,targA,targZ);
1273 
1274  // get max x-section
1275  double xsec_tot = 0;
1276  if (hpdgc == kPdgPiP ){xsec_tot = TMath::Max(0., fXSecPipp_Tot -> Evaluate(ke)) * targZ;
1277  xsec_tot+= TMath::Max(0., fXSecPipn_Tot -> Evaluate(ke)) * (targA-targZ);}
1278  else if (hpdgc == kPdgPiM ){xsec_tot = TMath::Max(0., fXSecPipn_Tot -> Evaluate(ke)) * targZ;
1279  xsec_tot+= TMath::Max(0., fXSecPipp_Tot -> Evaluate(ke)) * (targA-targZ);}
1280  else if (hpdgc == kPdgPi0 ){xsec_tot = TMath::Max(0., fXSecPi0p_Tot -> Evaluate(ke)) * targZ;
1281  xsec_tot+= TMath::Max(0., fXSecPi0n_Tot -> Evaluate(ke)) * (targA-targZ);}
1282  else if (hpdgc == kPdgProton ){xsec_tot = TMath::Max(0., fXSecPp_Tot -> Evaluate(ke)) * targZ;
1283  xsec_tot+= TMath::Max(0., fXSecPn_Tot -> Evaluate(ke)) * (targA-targZ);}
1284  else if (hpdgc == kPdgNeutron){xsec_tot = TMath::Max(0., fXSecPn_Tot -> Evaluate(ke)) * targZ;
1285  xsec_tot+= TMath::Max(0., fXSecNn_Tot -> Evaluate(ke)) * (targA-targZ);}
1286  else if (hpdgc == kPdgGamma ) xsec_tot = TMath::Max(0., fXSecGamN_Tot -> Evaluate(ke));
1287  else if (hpdgc == kPdgKP ) xsec_tot = TMath::Max(0., fXSecKpN_Tot -> Evaluate(ke));
1288 
1289  // compute fraction
1290  double frac = (xsec_tot>0) ? xsec/xsec_tot : 0.;
1291  return frac;
1292 }
1293 //____________________________________________________________________________
1294 double INukeHadroData::IntBounce(const GHepParticle* p, int target, int scode, INukeFateHN_t fate)
1295 {
1296  // This method returns a random cos(ang) according to a distribution
1297  // based upon the particle and fate. The sampling uses the
1298  // Accept/Reject method, whereby a distribution is bounded above by
1299  // an envelope, or in this case, a number of envelopes, which can be
1300  // easily sampled (here, we use uniform distributions).
1301  // To get a random value, first the envelope is sampled to
1302  // obtain an x-coordinate (cos(ang)), and then another random value
1303  // is obtained uniformally in the range [0,h(j,0)], where h(j,0)
1304  // is the height of the j-th envelope. If the point is beneath the
1305  // distribution, the x-coordinate is accepted, otherwise, we try
1306  // again.
1307 
1309 
1310  // numEnv is the number of envelopes in the total envelope,
1311  // that is, the number of seperate simple uniform distributions
1312  // that will be fit against the distribution in question in the
1313  // Accept/Reject process of sampling
1314  int numEnv = 4;
1315  int numPoints = 1000; // The number of points to be evaluated
1316  // for the purpose of finding the max
1317  // value of the distribution
1318  assert((numPoints%numEnv)==0); // numPoints/numEnv has to be an integer
1319  double sr = 2.0 / numEnv; // Subrange, i.e., range of an envelope
1320  double cstep = 2.0 / (numPoints); // Magnitude of the step between eval. points
1321 
1322  double ke = (p->E() - p->Mass()) * 1000.0; // ke in MeV
1323  if (TMath::Abs((int)ke-ke)<.01) ke+=.3; // make sure ke isn't an integer,
1324  // otherwise sometimes gives weird results
1325  // due to ROOT's Interpolate() function
1326  double avg = 0.0; // average value in envelop
1327 
1328  // Matrices to hold data; buff holds the distribution
1329  // data per envelope from which the max value is
1330  // obtained. That value is then recorded in dist, where
1331  // the integral of the envelope to that point is
1332  // also recorded
1333 
1334  double * buff = new double[numPoints/numEnv + 1];
1335  double ** dist = new double*[numEnv];
1336  for(int ih=0;ih<numEnv;ih++)
1337  {
1338  dist[ih] = new double[3];
1339  }
1340 
1341  // Acc-Rej Sampling Method
1342  // -- Starting at the beginning of each envelope,
1343  // this loop evaluates (numPoints) amount of points
1344  // on the distribution and holds them in buff;
1345  // then takes the max and places it in the first row
1346  // of h. The second row of h contains the interval of
1347  // the total envelope, up to the current envelope.
1348  // Thus, when properly normalized, the last value
1349  // in the second row of h should be 1.
1350  double totxsec = 0.0;
1351  for(int i=0;i<numEnv;i++)
1352  {
1353  double lbound = -1 + i*sr;
1354 
1355  for(int j=0;j<=numPoints / numEnv; j++)
1356  {
1357  buff[j] = this->XSec(p->Pdg(),target,scode,fate,ke,lbound+j*cstep);
1358  avg += buff[j];
1359  }
1360 
1361  totxsec+=avg;
1362  avg/= (double(numPoints)/double(numEnv));
1363  dist[i][0] = TMath::MaxElement(numPoints/numEnv+1,buff);
1364  dist[i][1] = avg;
1365  dist[i][2] = dist[i][1] + ((i==0)?0.0:dist[i-1][2]);
1366  avg=0.0;
1367  }
1368 
1369 
1370  delete [] buff;
1371 
1372  int iter=1; // keep track of iterations
1373  int env=0; // envelope index
1374  double rval = 0.0; // random value
1375  double val = 0.0; // angle value
1376 
1377  // Get a random point, see if its in the distribution, and if not
1378  // then try again.
1379 
1380  rval = rnd->RndFsi().Rndm()*dist[numEnv-1][2];
1381 
1382  env=0;
1383  // Check which envelope it's in, to
1384  // get proper height
1385  while(env<numEnv)
1386  {
1387  if(rval<=dist[env][2]) break;
1388  else env++;
1389  }
1390  if(env==numEnv) env=numEnv - 1;
1391 
1392 while(iter)
1393  {
1394 
1395  // Obtain the correct x-coordinate from the random sample
1396  val = rnd->RndFsi().Rndm()*sr;
1397  val += sr*env-1;
1398  rval = rnd->RndFsi().Rndm()*dist[env][0];
1399 
1400  // Test to see if point is in distribution, if it is, stop and return
1401  if(rval < this->XSec(p->Pdg(),target,scode,fate,ke,val)) break;
1402 
1403  // Possibly an extremely long loop, don't want to
1404  // hold up the program
1405  if(iter==1000)
1406  {
1407  int NUM_POINTS=2000;
1408  int pvalues=0;
1409  double points[200]={0};
1410  for(int k=0;k<NUM_POINTS;k++)
1411  {
1412  points[int(k/10)]=this->XSec(p->Pdg(),target,scode,fate,ke,-1+(2.0/NUM_POINTS)*k);
1413  if(points[int(k/10)]>0) pvalues++;
1414  }
1415  if(pvalues<(.05*NUM_POINTS))
1416  {
1417  // if it reaches here, one more test...if momenta of particle is
1418  // extremely low, just give it an angle from a uniform distribution
1419  if(p->P4()->P()<.005) // 5 MeV
1420  {
1421  val = 2*rnd->RndFsi().Rndm()-1;
1422  break;
1423  }
1424  else
1425  {
1426  LOG("Intranuke", pWARN) << "Hung-up in IntBounce method - Exiting";
1427  LOG("Intranuke", pWARN) << (*p);
1428  LOG("Intranuke", pWARN) << "Target: " << target << ", Scode: " << scode << ", fate: " << INukeHadroFates::AsString(fate);
1429  for(int ie=0;ie<200;ie+=10) {
1430  LOG("Intranuke", pWARN) << points[ie+0] << ", " << points[ie+1] << ", " << points[ie+2] << ", "
1431  << points[ie+3] << ", " << points[ie+4] << ", " << points[ie+5] << ", " << points[ie+6] << ", "
1432  << points[ie+7] << ", " << points[ie+8] << ", " << points[ie+9];
1433  }
1434  for(int ih=0;ih<numEnv;ih++)
1435  {
1436  delete [] dist[ih];
1437  }
1438  delete [] dist;
1439 
1440  return -2.;
1441  }
1442  }
1443  }
1444  iter++;
1445  }
1446 
1447  for(int ih=0;ih<numEnv;ih++)
1448  {
1449  delete [] dist[ih];
1450  }
1451  delete [] dist;
1452 
1453  return val;
1454 }
1455 //___________________________________________________________________________
1456 /*int INukeHadroData::AngleAndProduct(const GHepParticle* p, int target, double &costh, INukeFateHN_t fate)
1457 {
1458  //
1459  // Adapted from genie::utils::intranuke::IntBounce()
1460  // by Aaron Meyer (7/10/09)
1461  //
1462  // This method returns a random cos(ang) and a final product
1463  // PDG number according to a distribution
1464  // based upon the probe particle and fate. The sampling uses the
1465  // Accept/Reject method, described in the notes for IntBounce.
1466  // The final state particle is chosen by comparing a random value
1467  // that was obtained and accepted within the range to one of the
1468  // particle cross section envelops.
1469  // If the value is less than that of the first particle
1470  // cross section, the first particle is chosen; if not, the second
1471  // particle is chosen.
1472  //
1473  // The argument costh should be declared before calling. The cosine
1474  // value is determined in the function and returned via argument by
1475  // reference.
1476 
1477  RandomGen * rnd = RandomGen::Instance();
1478 
1479  // numEnv is the number of envelopes in the total envelope,
1480  // that is, the number of seperate simple uniform distributions
1481  // that will be fit against the distribution in question in the
1482  // Accept/Reject process of sampling
1483 
1484  int numEnv = 4;
1485  int numPoints = 1000; // The number of points to be evaluated
1486  // for the purpose of finding the max
1487  // value of the distribution
1488  assert((numPoints%numEnv)==0); // numPoints/numEnv has to be an integer
1489  double sr = 2.0 / numEnv; // Subrange, i.e., range of an envelope
1490  double cstep = 2.0 / (numPoints); // Magnitude of the step between eval. points
1491  int fsPart1 = 0; // First possible particle final state
1492  int fsPart2 = 0; // Second possible particle final state
1493  int endPart = 0; // Particle to be returned at completion
1494 
1495  if ((p->Pdg()==kPdgGamma) && (target==kPdgProton)) {
1496  fsPart1 = kPdgProton;
1497  fsPart2 = kPdgNeutron;
1498  } else if ((p->Pdg()==kPdgGamma) && (target==kPdgNeutron)) {
1499  fsPart1 = kPdgProton;
1500  fsPart2 = kPdgNeutron;
1501  } else {
1502  LOG("INukeHadroData", pERROR) << "Cannot handle particle: " << p->Pdg()
1503  << " and target: " << target;
1504  return endPart;
1505  }
1506 
1507  LOG("INukeHadroData", pNOTICE) << "Possible final state particles :" << fsPart1 << " or " <<fsPart2;
1508 
1509 
1510  double ke = (p->E() - p->Mass()) * 1000.0; // ke in MeV
1511  double avg = 0.0; // Temp value used to find average of envelope
1512 
1513  double * buff = new double[numPoints/numEnv + 1];
1514  double ** dist = new double*[numEnv];
1515  for(int ih=0;ih<numEnv;ih++)
1516  {
1517  dist[ih] = new double[3];
1518  }
1519 
1520  // Acc-Rej Sampling Method
1521  // see IntBounce above for description
1522  for(int i=0;i<numEnv;i++)
1523  {
1524  double lbound = -1 + i*sr;
1525 
1526  for(int j=0;j<=numPoints / numEnv; j++)
1527  {
1528  buff[j] = this->XSec(p->Pdg(),target,fsPart1,fate,ke,lbound+j*cstep);
1529  buff[j]+= this->XSec(p->Pdg(),target,fsPart2,fate,ke,lbound+j*cstep);
1530  avg += buff[j];
1531  }
1532 
1533  avg /= (double(numPoints)/double(numEnv));
1534  dist[i][0] = TMath::MaxElement(numPoints/numEnv+1,buff);
1535  dist[i][1] = avg;
1536  dist[i][2] = dist[i][1] + ((i==0)?0.0:dist[i-1][2]);
1537  avg = 0.0;
1538  LOG("INukeHadroData", pNOTICE) << "max xsec value for env "<<i<<": " << dist[i][0];
1539  LOG("INukeHadroData", pNOTICE) << "xsec avg for env "<<i<<": " << dist[i][1];
1540  LOG("INukeHadroData", pNOTICE) << "xsec avg integral at env "<<i<<": " << dist[i][2];
1541 
1542  }
1543 
1544  delete [] buff;
1545 
1546  int iter=1; // keep track of iterations
1547  int env=0; // envelope index
1548  double rval = 0.0; // random value
1549  double val = 0.0; // costh value
1550 
1551 
1552  // Get a random point, see if its in the distribution, and if not
1553  // then try again.
1554 
1555  rval = rnd->RndFsi().Rndm()*dist[numEnv-1][2];
1556  LOG("INukeHadroData", pNOTICE) << "random variable for envelope decision: " << rval;
1557 
1558  env=0;
1559  // Check which envelope it's in, to
1560  // get proper height
1561  while(env<numEnv)
1562  {
1563  if(rval<=dist[env][2]) break;
1564  else env++;
1565  }
1566  if(env==numEnv) env=numEnv - 1;
1567  LOG("INukeHadroData", pNOTICE) << "value in envelope: " << env;
1568  LOG("INukeHadroData", pNOTICE) << "weighted envelope bound: " << dist[env][2];
1569 
1570  while(iter)
1571  {
1572  // Obtain the correct x-coordinate from the random sample
1573  val = rnd->RndFsi().Rndm()*sr;
1574  val+= (env*sr - 1);
1575  rval = rnd->RndFsi().Rndm()*dist[env][0];
1576 
1577  // Test to see if point is in distribution, if it is, stop and return
1578  if(rval < (this->XSec(p->Pdg(),target,fsPart1,fate,ke,val) +
1579  this->XSec(p->Pdg(),target,fsPart2,fate,ke,val)))
1580  {
1581  // Determine final state particle
1582  LOG("INukeHadroData", pNOTICE) << "particle 1 bound: " << this->XSec(p->Pdg(),target,fsPart1,fate,ke,val);
1583  LOG("INukeHadroData", pNOTICE) << "xsec bound: " << (this->XSec(p->Pdg(),target,fsPart1,fate,ke,val) +\
1584  this->XSec(p->Pdg(),target,fsPart2,fate,ke,val));
1585 
1586  if (rval < this->XSec(p->Pdg(),target,fsPart1,fate,ke,val))
1587  {endPart = fsPart1;
1588  LOG("INukeHadroData", pNOTICE) << "particle 1 taken: " << endPart;
1589  }
1590  else
1591  {endPart = fsPart2;
1592  LOG("INukeHadroData", pNOTICE) << "particle 2 taken: " << endPart;
1593  }
1594  costh = val;
1595  break;
1596  }
1597 
1598  // Possibly an extremely long loop, don't want to
1599  // hold up the program
1600  if(iter==1000)
1601  {
1602  int NUM_POINTS=2000;
1603  int pvalues=0;
1604  double points[200]={0};
1605  for(int k=0;k<NUM_POINTS;k++)
1606  {
1607  points[int(k/10)]=this->XSec(p->Pdg(),target,fsPart1,fate,ke,-1+(2.0/NUM_POINTS)*k) +
1608  this->XSec(p->Pdg(),target,fsPart2,fate,ke,-1+(2.0/NUM_POINTS)*k);
1609  if(points[int(k/10)]>0) pvalues++;
1610  }
1611  if(pvalues<(.05*NUM_POINTS))
1612  {
1613  // if it reaches here, one more test...if momenta of particle is
1614  // extremely low, just give it an angle from a uniform distribution
1615  if(p->P4()->P()<.005) // 5 MeV
1616  {
1617  val = 2*rnd->RndFsi().Rndm()-1;
1618  break;
1619  }
1620  else
1621  {
1622  LOG("INukeHadroData", pWARN) << "Hung-up in AngleAndProduct method - Exiting";
1623  LOG("INukeHadroData", pWARN) << (*p);
1624  LOG("INukeHadroData", pWARN) << "Target: " << target << ", Potential final state particles: "
1625  << fsPart1 << " " << fsPart2 <<", fate: " << INukeHadroFates::AsString(fate);
1626  for(int ie=0;ie<200;ie+=10) {
1627  LOG("INukeHadroData", pWARN) << points[ie+0] << ", " << points[ie+1] << ", " << points[ie+2] << ", "
1628  << points[ie+3] << ", " << points[ie+4] << ", " << points[ie+5] << ", " << points[ie+6] << ", "
1629  << points[ie+7] << ", " << points[ie+8] << ", " << points[ie+9];
1630  }
1631 
1632  for(int ih=0;ih<numEnv;ih++)
1633  {
1634  delete [] dist[ih];
1635  }
1636  delete [] dist;
1637 
1638  return 0;
1639  }
1640  }
1641  }
1642  iter++;
1643  }
1644  for(int ih=0;ih<numEnv;ih++)
1645  {
1646  delete [] dist[ih];
1647  }
1648  delete [] dist;
1649 
1650  // LOG("INukeHadroData", pNOTICE) << "return value: " << endPart;
1651 
1652  return endPart;
1653 }*/
1654 //___________________________________________________________________________
static INukeHadroData * fInstance
Basic constants.
Double_t angle
Definition: plot.C:86
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
const char * p
Definition: xmltok.h:285
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
TString ip
Definition: loadincs.C:5
string filename
Definition: shutoffs.py:106
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:30
double Frac(int hpdgc, INukeFateHA_t fate, double ke) const
double dist
Definition: runWimpSim.h:113
int Pdg(void) const
Definition: GHepParticle.h:64
cstep
Definition: test2.py:21
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
enum genie::EINukeFateHN_t INukeFateHN_t
const int cols[3]
static double fMinKinEnergy
const int kPdgGamma
Definition: PDGCodes.h:166
static INukeHadroData * Instance(void)
const int kPdgKP
Definition: PDGCodes.h:149
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
const double j
Definition: BetheBloch.cxx:29
double frac(double x)
Fractional part.
#define pINFO
Definition: Messenger.h:63
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
void ReadhNFile(string filename, double ke, int npoints, int &curr_point, double *costh_array, double *xsec_array, int cols)
ifstream in
Definition: comparison.C:7
static string AsString(INukeFateHN_t fate)
assert(nhit_max >=nhit_nbins)
const int kPdgPiM
Definition: PDGCodes.h:136
static double fMaxKinEnergyHN
double XSec(int hpdgc, int tgt, int nprod, INukeFateHN_t rxnType, double ke, double costh) const
const int kPdgProton
Definition: PDGCodes.h:65
Singleton class to load & serve hadron x-section splines used by GENIE&#39;s version of the INTRANUKE cascade...
static double fMaxKinEnergyHA
const int kPdgNeutron
Definition: PDGCodes.h:67
static const double kPi
Definition: Constants.h:38
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:40
enum genie::EINukeFateHA_t INukeFateHA_t
static constexpr Double_t sr
Definition: Munits.h:164
#define pDEBUG
Definition: Messenger.h:64