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