GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
INukeUtils2018.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  Author: Jim Dobson <j.dobson07 \at imperial.ac.uk>
8  Imperial College London
9 
10  Costas Andreopoulos <c.andreopoulos \at cern.ch>
11  University of Liverpool
12 
13  Aaron Meyer <asm58 \at pitt.edu>
14  Pittsburgh University
15 
16  For documentation see the corresponding header file.
17 
18  Important revisions after version 2.0.0 :
19  @ Mar 04, 2009 - JD
20  Was first added in v2.5.1. Adapted from the T2K GENIE reweighting tool.
21  @ Mar 05, 2009 - CA
22  Modified ReconstructHadronFateHA() to work with hadron+A event files in
23  addition to neutrino event files.
24  @ Sep 10, 2009 - CA
25  Added MeanFreePath(), Dist2Exit(), Dist2ExitMFP()
26  @ Sep 30, 2009 - CA
27  Added StepParticle() from Intranuke.cxx
28  @ Oct 02, 2009 - CA
29  Added test MeanFreePath_Delta().
30  @ Jul 15, 2010 - AM
31  Added common utility functions used by both hA and hN mode. Updated
32  MeanFreePath to separate proton and neutron cross sections. Added general
33  utility functions.
34  @ Jan 9, 2015 - SD, NG, TG
35  Added 2014 version of INTRANUKE codes for v2.9.0. Uses INukeHadroData2014,
36  but no changes to mean free path.
37 */
38 //____________________________________________________________________________
39 
40 #include <TLorentzVector.h>
41 #include <TMath.h>
42 #include <TSystem.h>
43 
44 #include "Framework/Algorithm/AlgConfigPool.h"
45 #include "Framework/Conventions/Constants.h"
46 #include "Framework/Conventions/Controls.h"
47 #include "Framework/Conventions/Units.h"
48 #include "Framework/Conventions/GBuild.h"
49 #include "Framework/EventGen/EVGThreadException.h"
50 #include "Framework/GHEP/GHepRecord.h"
51 #include "Framework/GHEP/GHepParticle.h"
52 #include "Physics/HadronTransport/Intranuke2018.h"
53 #include "Physics/HadronTransport/INukeException.h"
54 #include "Physics/HadronTransport/INukeUtils2018.h"
55 #include "Physics/HadronTransport/INukeHadroData2018.h"
56 #include "Framework/Messenger/Messenger.h"
57 #include "Framework/Numerical/RandomGen.h"
58 #include "Framework/Numerical/Spline.h"
59 #include "Framework/ParticleData/PDGLibrary.h"
60 #include "Framework/ParticleData/PDGUtils.h"
61 #include "Framework/ParticleData/PDGCodes.h"
62 #include "Framework/ParticleData/PDGCodeList.h"
63 #include "Framework/ParticleData/PDGUtils.h"
64 #include "Framework/Registry/Registry.h"
65 #include "Physics/NuclearState/NuclearUtils.h"
66 #include "Framework/Utils/PrintUtils.h"
67 #include "Physics/HadronTransport/INukeOset.h"
68 #include "Physics/HadronTransport/INukeOsetTable.h"
69 #include "Physics/HadronTransport/INukeOsetFormula.h"
70 #include "TComplex.h"
71 
72 using std::ostringstream;
73 using namespace genie;
74 using namespace genie::utils;
75 using namespace genie::constants;
76 using namespace genie::controls;
77 
78 //____________________________________________________________________________
79 double genie::utils::intranuke2018::MeanFreePath(
80  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4,
81  double A, double Z, double nRpi, double nRnuc, const bool useOset, const bool altOset, const bool xsecNNCorr, string INukeMode)
82 {
83 // Calculate the mean free path (in fm) for a pions and nucleons in a nucleus
84 //
85 // Inputs
86 // pdgc : Hadron PDG code
87 // x4 : Hadron 4-position in the nucleus coordinate system (units: fm)
88 // p4 : Hadron 4-momentum (units: GeV)
89 // A : Nucleus atomic mass number
90 // nRpi : Controls the pion ring size in terms of de-Broglie wavelengths
91 // nRnuc: Controls the nuclepn ring size in terms of de-Broglie wavelengths
92 //
93  bool is_pion = pdgc == kPdgPiP || pdgc == kPdgPi0 || pdgc == kPdgPiM;
94  bool is_nucleon = pdgc == kPdgProton || pdgc == kPdgNeutron;
95  bool is_kaon = pdgc == kPdgKP;
96  bool is_gamma = pdgc == kPdgGamma;
97 
98  if(!is_pion && !is_nucleon && !is_kaon && !is_gamma) return 0.;
99 
100  // before getting the nuclear density at the current position
101  // check whether the nucleus has to become larger by const times the
102  // de Broglie wavelength -- that is somewhat empirical, but this
103  // is what is needed to get piA total cross sections right.
104  // The ring size is different for light nuclei (using gaus density) /
105  // heavy nuclei (using woods-saxon density).
106  // The ring size is different for pions / nucleons.
107  //
108  double momentum = p4.Vect().Mag(); // hadron momentum in GeV
109  double ring = (momentum>0) ? 1.240/momentum : 0; // de-Broglie wavelength
110 
111  if(A<=20) { ring /= 2.; }
112 
113  /*
114  if (is_pion ) { ring *= nRpi; }
115  else if (is_nucleon ) { ring *= nRnuc; }
116  else if (is_gamma || is_kaon || useOset) { ring = 0.; }
117  */
118  if(INukeMode=="hN2018")
119  {
120  if (is_pion ) { ring *= nRpi; }
121  else if (is_nucleon ) { ring *= nRnuc; }
122  else if (is_gamma || is_kaon || useOset) { ring = 0.;}
123  }
124  else
125  {
126  if (is_pion || is_kaon ) { ring *= nRpi; }
127  else if (is_nucleon ) { ring *= nRnuc; }
128  else if (is_gamma ) { ring = 0.; }
129  }
130 
131  // get the nuclear density at the current position
132  double rnow = x4.Vect().Mag();
133  double rho = A * utils::nuclear::Density(rnow,(int) A,ring);
134 
135  // the hadron+nucleon cross section will be evaluated within the range
136  // of the input spline and assumed to be const outside that range
137  //
138  double ke = (p4.Energy() - p4.M()) / units::MeV; // kinetic energy in MeV
139  ke = TMath::Max(INukeHadroData2018::fMinKinEnergy, ke);
140  ke = TMath::Min(INukeHadroData2018::fMaxKinEnergyHN, ke);
141 
142  // get total xsection for the incident hadron at its current
143  // kinetic energy
144  double sigtot = 0;
145  double ppcnt = (double) Z/ (double) A; // % of protons remaining
146  INukeHadroData2018 * fHadroData2018 = INukeHadroData2018::Instance();
147 
148  if (is_pion and (INukeMode == "hN2018") and useOset and ke < 350.0)
149  sigtot = sigmaTotalOset (ke, rho, pdgc, ppcnt, altOset);
150  else if (pdgc == kPdgPiP)
151  { sigtot = fHadroData2018 -> XSecPipp_Tot() -> Evaluate(ke)*ppcnt;
152  sigtot+= fHadroData2018 -> XSecPipn_Tot() -> Evaluate(ke)*(1-ppcnt);}
153  else if (pdgc == kPdgPi0)
154  { sigtot = fHadroData2018 -> XSecPi0p_Tot() -> Evaluate(ke)*ppcnt;
155  sigtot+= fHadroData2018 -> XSecPi0n_Tot() -> Evaluate(ke)*(1-ppcnt);}
156  else if (pdgc == kPdgPiM)
157  { sigtot = fHadroData2018 -> XSecPipn_Tot() -> Evaluate(ke)*ppcnt;
158  sigtot+= fHadroData2018 -> XSecPipp_Tot() -> Evaluate(ke)*(1-ppcnt);}
159  else if (pdgc == kPdgProton)
160  {
161  sigtot = fHadroData2018 -> XSecPp_Tot() -> Evaluate(ke)*ppcnt;
162  //sigtot+= fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*(1-ppcnt);
163 
164  PDGLibrary * pLib = PDGLibrary::Instance();
165  double hc = 197.327;
166  double R0 = 1.25 * TMath::Power(A,1./3.) + 2.0 * 0.65; // should all be in units of fm
167  double Mp = pLib->Find(2212)->Mass();
168  double M = pLib->Find(pdgc)->Mass();
169  //double E = (p4.Energy() - Mp) * 1000.; // Convert GeV to MeV.
170  double E = ke;
171  if (Z*hc/137./x4.Vect().Mag() > E) // Coulomb correction (Cohen, Concepts of Nuclear Physics, pg. 259-260)
172  {
173  double z = 1.0; // charge for single proton
174  double Bc = z*Z*hc/137./R0;
175  double x = E/Bc;
176  double f = TMath::ACos(TMath::Power(x,0.5)) - TMath::Power(x*(1-x),0.5);
177  double B = 0.63*z*Z*TMath::Power((M/Mp)/E,0.5);
178  double Pc = TMath::Exp(-B*f);
179  sigtot *= Pc;
180  }
181  sigtot+= fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*(1-ppcnt);
182 
183  double E0 = TMath::Power(A,0.2)*12.;
184  if (INukeMode=="hN2018"){if(ke<E0){sigtot=0.0;}} //empirical - needed to cut off large number of low energy nucleons
185  // LOG("INukeUtils",pDEBUG) "sigtot for proton= " << sigtot << "; KE= " << ke << "; E0= " << E0;
186  }
187  else if (pdgc == kPdgNeutron)
188  {
189  sigtot = fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*ppcnt;
190  sigtot+= fHadroData2018 -> XSecNn_Tot() -> Evaluate(ke)*(1-ppcnt);
191  double E0 = TMath::Power(A,0.2)*12.;
192  if (INukeMode=="hN2018"){if(ke<E0){sigtot=0.0;}} //empirical - needed to cut off large number of low energy nucleons
193  // LOG("INukeUtils",pDEBUG) "sigtot for neutron= " << sigtot << "; KE= " << ke;
194  }
195  else if (pdgc == kPdgKP)
196  { sigtot = fHadroData2018 -> XSecKpN_Tot() -> Evaluate(ke);
197  // this factor is used to empirically get agreement with tot xs data, justified historically.
198  sigtot*=1.1;}
199  else if (pdgc == kPdgGamma)
200  { sigtot = fHadroData2018 -> XSecGamp_fs() -> Evaluate(ke)*ppcnt;
201  sigtot+= fHadroData2018 -> XSecGamn_fs() -> Evaluate(ke)*(1-ppcnt);}
202  else {
203  return 0;
204  }
205 
206  // the xsection splines in INukeHadroData return the hadron x-section in
207  // mb -> convert to fm^2
208  sigtot *= (units::mb / units::fm2);
209 
210  // avoid defective error handling
211  //if(sigtot<1E-6){sigtot=1E-6;}
212 
213  if (xsecNNCorr and is_nucleon)
214  sigtot *= INukeNucleonCorr::getInstance()->
215  getAvgCorrection (rho, A, p4.E() - PDGLibrary::Instance()->Find(pdgc)->Mass()); //uses lookup tables
216 
217  // avoid defective error handling
218  if(sigtot<1E-6){sigtot=1E-6;}
219 
220  // compute the mean free path
221  double lamda = 1. / (rho * sigtot);
222 
223  // exits if lamda is InF (if cross section is 0)
224  if( ! TMath::Finite(lamda) ) {
225  return -1;
226  }
227 
228 /*
229  LOG("INukeUtils", pDEBUG)
230  << "sig_total = " << sigtot << " fm^2, rho = " << rho
231  << " fm^-3 => mfp = " << lamda << " fm.";
232 */
233  return lamda;
234 }
235 //____________________________________________________________________________
236 double genie::utils::intranuke2018::MeanFreePath_Delta(
237  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4, double A)
238 {
239 //
240 // **test**
241 //
242 
243 // Calculate the mean free path (in fm) for Delta's in a nucleus
244 //
245 // Inputs
246 // pdgc : Hadron PDG code
247 // x4 : Hadron 4-position in the nucleus coordinate system (units: fm)
248 // p4 : Hadron 4-momentum (units: GeV)
249 // A : Nucleus atomic mass number
250 //
251  bool is_deltapp = (pdgc==kPdgP33m1232_DeltaPP);
252  if(!is_deltapp) return 0.;
253 
254  // get the nuclear density at the current position
255  double rnow = x4.Vect().Mag();
256  double rho = A * utils::nuclear::Density(rnow,(int) A);
257 
258  // the Delta+N->N+N cross section will be evaluated within the range
259  // of the input spline and assumed to be const outside that range
260  double ke = (p4.Energy() - p4.M()) / units::MeV; // kinetic energy in MeV
261  ke = TMath::Max(1500., ke);
262  ke = TMath::Min( 0., ke);
263 
264  // get the Delta+N->N+N cross section
265  double sig = 0;
266  if (ke< 500) sig=20;
267  else if (ke<1000) sig=40;
268  else sig=30;
269 
270  // value is in mb -> convert to fm^2
271  sig *= (units::mb / units::fm2);
272 
273  // compute the mean free path
274  double lamda = 1. / (rho * sig);
275 
276  // exits if lamda is InF (if cross section is 0)
277  if( ! TMath::Finite(lamda) ) {
278  return -1;
279  }
280 
281  return lamda;
282 }
283 //____________________________________________________________________________
284 double genie::utils::intranuke2018::ProbSurvival(
285  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4, double A, double /*Z*/,
286  double mfp_scale_factor, const Intranuke2018& fsi_model )
287 {
288 // Calculate the survival probability for a hadron inside a nucleus
289 //
290 // Inputs
291 // pdgc : Hadron PDG code
292 // x4 : Hadron 4-position in the nucleus coordinate system (units: fm)
293 // p4 : Hadron 4-momentum (units: GeV)
294 // A : Target nucleus atomic mass number
295 // mfp_scale_factor: Tweaks the mean free path (mfp -> mfp*scale). Def: 1.0
296 
297  double prob = 1.0;
298 
299  // Get extra parameters from the FSI model that we need to compute the
300  // mean free path
301 
302  // nRpi: Controls the pion ring size in terms of de-Broglie wavelengths
303  double nRpi = fsi_model.GetDelRPion();
304 
305  // nRnuc: Controls the nuclepn ring size in terms of de-Broglie wavelengths
306  double nRnuc = fsi_model.GetDelRNucleon();
307 
308  // NR: How far away to track the hadron, in terms of the corresponding
309  // nuclear radius. Def: 3
310  double NR = fsi_model.GetNR();
311 
312  // R0: R0 in R=R0*A^1/3 (units:fm). Def. 1.4
313  double R0 = fsi_model.GetR0();
314 
315  // step: Step size (fm) to use when tracking hadrons
316  double step = fsi_model.GetHadStep();
317 
318  // Boolean model configuration options
319  bool useOset = fsi_model.GetUseOset();
320  bool altOset = fsi_model.GetAltOset();
321  bool xsecNNCorr = fsi_model.GetXsecNNCorr();
322 
323  // Intranuke mode setting ("HA2018", etc.)
324  std::string inuke_mode = fsi_model.GetINukeMode();
325 
326  // Maximum radius to use in the stepping loop. Note that Intranuke2018 uses
327  // the *target* mass number to choose this radius, not the value for the
328  // pre- or post-FSI remnant.
329  double R = NR * R0 * TMath::Power(A, 1./3.);
330 
331  TVector3 dr3 = p4.Vect().Unit(); // unit vector along its direction
332  TLorentzVector dr4(dr3,0);
333 
334  int remnA = fsi_model.GetRemnA();
335  int remnZ = fsi_model.GetRemnZ();
336 
337  LOG("INukeUtils", pDEBUG)
338  << "Calculating survival probability for hadron with PDG code = " << pdgc
339  << " and momentum = " << p4.P() << " GeV";
340  LOG("INukeUtils", pDEBUG)
341  << "mfp scale = " << mfp_scale_factor
342  << ", nRpi = " << nRpi << ", nRnuc = " << nRnuc << ", NR = " << NR
343  << ", R0 = " << R0 << " fm";
344 
345  TLorentzVector x4_curr(x4); // current position
346 
347  while(1) {
348  double rnow = x4_curr.Vect().Mag();
349  if (rnow > (R+step)) break;
350 
351  x4_curr += (step*dr4);
352  rnow = x4_curr.Vect().Mag();
353 
354  double mfp = genie::utils::intranuke2018::MeanFreePath( pdgc, x4_curr, p4,
355  remnA, remnZ, nRpi, nRnuc, useOset, altOset, xsecNNCorr, inuke_mode );
356  double mfp_twk = mfp * mfp_scale_factor;
357 
358  double dprob = (mfp_twk>0) ? TMath::Exp(-step/mfp_twk) : 0.;
359  prob*=dprob;
360 /*
361  LOG("INukeUtils", pDEBUG)
362  << "+ step size = " << step << " fm, |r| = " << rnow << " fm, "
363  << "mfp = " << mfp_twk << "fm (nominal mfp = " << mfp << " fm): "
364  << "dPsurv = " << dprob << ", current Psurv = " << prob;
365 */
366  }
367 
368  LOG("INukeUtils", pDEBUG) << "Psurv = " << prob;
369 
370  return prob;
371 }
372 //____________________________________________________________________________
373 double genie::utils::intranuke2018::Dist2Exit(
374  const TLorentzVector & x4, const TLorentzVector & p4,
375  double A, double NR, double R0)
376 {
377 // Calculate distance within a nucleus (units: fm) before we stop tracking
378 // the hadron.
379 // See previous functions for a description of inputs.
380 //
381  double R = NR * R0 * TMath::Power(A, 1./3.);
382  double step = 0.05; // fermi
383 
384  TVector3 dr3 = p4.Vect().Unit(); // unit vector along its direction
385  TLorentzVector dr4(dr3,0);
386 
387  TLorentzVector x4_curr(x4); // current position
388 
389  double d=0;
390  while(1) {
391  double rnow = x4_curr.Vect().Mag();
392  x4_curr += (step*dr4);
393  d+=step;
394  rnow = x4_curr.Vect().Mag();
395  if (rnow > R) break;
396  }
397  return d;
398 }
399 //____________________________________________________________________________
400 double genie::utils::intranuke2018::Dist2ExitMFP(
401  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4,
402  double A, double Z, double NR, double R0)
403 {
404 // Calculate distance within a nucleus (expressed in terms of 'mean free
405 // paths') before we stop tracking the hadron.
406 // See previous functions for a description of inputs.
407 //
408 
409 // distance before exiting in mean free path lengths
410 //
411  double R = NR * R0 * TMath::Power(A, 1./3.);
412  double step = 0.05; // fermi
413 
414  TVector3 dr3 = p4.Vect().Unit(); // unit vector along its direction
415  TLorentzVector dr4(dr3,0);
416 
417  TLorentzVector x4_curr(x4); // current position
418 
419  double d=0;
420  double d_mfp=0;
421  while(1) {
422  double rnow = x4_curr.Vect().Mag();
423  x4_curr += (step*dr4);
424  d+=step;
425  rnow = x4_curr.Vect().Mag();
426 
427  double lambda = genie::utils::intranuke2018::MeanFreePath(pdgc,x4_curr,p4,A,Z);
428  d_mfp += (step/lambda);
429 
430  if (rnow > R) break;
431  }
432  return d_mfp;
433 }
434 //____________________________________________________________________________
435 void genie::utils::intranuke2018::StepParticle(
436  GHepParticle * p, double step, double nuclear_radius)
437 {
438 // Steps a particle starting from its current position (in fm) and moving
439 // along the direction of its current momentum by the input step (in fm).
440 // The particle is stepped in a straight line.
441 // If a nuclear radius is set then the following check is performed:
442 // If the step is too large and takes the the particle far away from the
443 // nucleus then its position is scaled back so that the escaped particles are
444 // always within a ~1fm from the "outer nucleus surface"
445 
446 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
447  LOG("INukeUtils", pDEBUG)
448  << "Stepping particle [" << p->Name() << "] by dr = " << step << " fm";
449 #endif
450 
451  // Step particle
452  TVector3 dr = p->P4()->Vect().Unit(); // unit vector along its direction
453  dr.SetMag(step); // spatial step size
454  double dt = 0; // temporal step:
455  TLorentzVector dx4(dr,dt); // 4-vector step
456  TLorentzVector x4new = *(p->X4()) + dx4; // new position
457 
458  if(nuclear_radius > 0.) {
459  // Check position against nuclear boundary. If the particle was stepped
460  // too far away outside the nuclear boundary bring it back to within
461  // 1fm from that boundary
462  double epsilon = 1; // fm
463  double r = x4new.Vect().Mag(); // fm
464  double rmax = nuclear_radius+epsilon;
465  if(r > rmax) {
466  LOG("INukeUtils", pINFO)
467  << "Particle was stepped too far away (r = " << r << " fm)";
468  LOG("INukeUtils", pINFO)
469  << "Placing it " << epsilon
470  << " fm outside the nucleus (r' = " << rmax << " fm)";
471  double scale = rmax/r;
472  x4new *= scale;
473  }//r>rmax
474  }//nucl radius set
475 
476 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
477  LOG("INukeUtils", pDEBUG)
478  << "\n Init direction = " << print::Vec3AsString(&dr)
479  << "\n Init position (in fm,nsec) = " << print::X4AsString(p->X4())
480  << "\n Fin position (in fm,nsec) = " << print::X4AsString(&x4new);
481 #endif
482 
483  p->SetPosition(x4new);
484 }
485 
486 
487 //___________________________________________________________________________
488 // Method to handle compound nucleus considerations, preequilibrium
489 // and equilibrium
490 // Alex Bell -- 6/17/2008
491 void genie::utils::intranuke2018::PreEquilibrium(
492  GHepRecord * ev, GHepParticle * p,
493  int &RemnA, int &RemnZ, TLorentzVector &RemnP4,
494  bool /* DoFermi */, double /* FermiFac */,
495  const NuclearModelI* /* Nuclmodel */, double NucRmvE, EINukeMode mode)
496 {
497 
498 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
499  LOG("INukeUtils", pDEBUG)
500  << "PreEquilibrium() is invoked for a : " << p->Name()
501  << " whose kinetic energy is : " << p->KinE();
502 #endif
503 
504  // Random number generator
505  RandomGen * rnd = RandomGen::Instance();
506  //unused PDGLibrary * pLib = PDGLibrary::Instance();
507 
508  bool allow_dup = true;
509  PDGCodeList list(allow_dup); // list of final state particles
510 
511  double ppcnt = (double) RemnZ / (double) RemnA; // % of protons left
512 
513  // figure out the final state conditions
514 
515  if(p->Pdg()==kPdgProton) list.push_back(kPdgProton);
516  else if(p->Pdg()==kPdgNeutron) list.push_back(kPdgNeutron);
517 
518  for(int i=0;i<3;i++)
519  {
520  if(rnd->RndFsi().Rndm()<ppcnt)
521  {
522  list.push_back(kPdgProton);
523  RemnZ--;
524  }
525  else list.push_back(kPdgNeutron);
526 
527  RemnA--;
528 
529  ppcnt = (double) RemnZ / (double) RemnA;
530  }
531 
532  // Add the fermi energy of the three nucleons to the phase space
533  /*
534  if(DoFermi)
535  {
536  Target target(ev->TargetNucleus()->Pdg());
537  TVector3 pBuf = p->P4()->Vect();
538  double mBuf = p->Mass();
539  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
540  TLorentzVector tSum(pBuf,eBuf);
541  double mSum = 0.0;
542  vector<int>::const_iterator pdg_iter;
543  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
544  {
545  target.SetHitNucPdg(*pdg_iter);
546  Nuclmodel->GenerateNucleon(target);
547  mBuf = pLib->Find(*pdg_iter)->Mass();
548  mSum += mBuf;
549  pBuf = FermiFac * Nuclmodel->Momentum3();
550  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
551  tSum += TLorentzVector(pBuf,eBuf);
552  RemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
553  }
554  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
555  p->SetMomentum(dP4);
556  }
557  */
558  // do the phase space decay & save all f/s particles to the event record
559  bool success = genie::utils::intranuke2018::PhaseSpaceDecay(ev,p,list,RemnP4,NucRmvE,mode);
560  if(success) LOG("INukeUtils2018",pINFO) << "Successful phase space decay for pre-equilibrium nucleus FSI event";
561  else
562  {
563  exceptions::INukeException exception;
564  exception.SetReason("Phase space generation of pre-equilibrium nucleus final state failed, details above");
565  throw exception;
566  }
567 
568  int p_loc = 0;
569  while(p_loc<ev->GetEntries())
570  {
571  GHepParticle * p_ref = ev->Particle(p_loc);
572  if(!p->ComparePdgCodes(p_ref)) p_loc++;
573  else
574  {
575  if(!p->CompareStatusCodes(p_ref)) p_loc++;
576  else
577  {
578  if(!p->CompareMomentum(p_ref)) p_loc++;
579  else break;
580  }
581  }
582  }
583 
584 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
585  LOG("INukeUtils", pDEBUG)
586  << "Particle at: " << p_loc;
587 #endif
588 
589  // find the appropriate daughter
590  vector<int> * descendants = ev->GetStableDescendants(p_loc);
591 
592  int loc = p_loc + 1;
593  int f_loc = p_loc + 1;
594  double energy = ev->Particle(loc)->E();
595 
596 /* // (1) least energetic
597  double min_en = energy;
598 
599  for(unsigned int j=0;j<descendants->size();j++)
600  {
601  loc = (*descendants)[j];
602  energy = ev->Particle(loc)->E();
603  if(energy<min_en)
604  {
605  f_loc = loc;
606  min_en = energy;
607  }
608  }
609 */
610  // (2) most energetic
611  double max_en = energy;
612 
613  for(unsigned int j=0;j<descendants->size();j++)
614  {
615  loc = (*descendants)[j];
616  energy = ev->Particle(loc)->E();
617  if(energy>max_en)
618  {
619  f_loc = loc;
620  max_en = energy;
621  }
622  }
623 
624  // (3) 1st particle
625  // ...just use the defaulted f_loc
626 
627  delete descendants;
628 
629  // change particle status for decaying particle - take out as test
630  //ev->Particle(f_loc)->SetStatus(kIStIntermediateState);
631  // decay a clone particle
632  GHepParticle * t = new GHepParticle(*(ev->Particle(f_loc)));
633  t->SetFirstMother(f_loc);
634  //next statement was in Alex Bell's original code - PreEq, then Equilibrium using particle with highest energy. Note it gets IST=kIStIntermediateState.
635  //genie::utils::intranuke2018::Equilibrium(ev,t,RemnA,RemnZ,RemnP4,DoFermi,FermiFac,Nuclmodel,NucRmvE,mode);
636 
637  delete t;
638 }
639 //___________________________________________________________________________
640 // Method to handle Equilibrium reaction
641 // Alex Bell -- 6/17/2008
642 void genie::utils::intranuke2018::Equilibrium(
643  GHepRecord * ev, GHepParticle * p,
644  int &RemnA, int &RemnZ, TLorentzVector &RemnP4,
645  bool /* DoFermi */, double /* FermiFac */,
646  const NuclearModelI* /* Nuclmodel */, double NucRmvE, EINukeMode mode)
647 {
648 
649 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
650  LOG("INukeUtils", pDEBUG)
651  << "Equilibrium() is invoked for a : " << p->Name()
652  << " whose kinetic energy is : " << p->KinE();
653 #endif
654 
655  // Random number generator
656  RandomGen * rnd = RandomGen::Instance();
657  //usused PDGLibrary * pLib = PDGLibrary::Instance();
658 
659  bool allow_dup = true;
660  PDGCodeList list(allow_dup); // list of final state particles
661 
662  // % of protons left
663  double ppcnt = (double) RemnZ / (double) RemnA;
664 
665  // figure out the final state conditions
666 
667  if(p->Pdg()==kPdgProton) list.push_back(kPdgProton);
668  else if(p->Pdg()==kPdgNeutron) list.push_back(kPdgNeutron);
669 
670  //add additonal particles to stack
671  for(int i=0;i<4;i++)
672  {
673  if(rnd->RndFsi().Rndm()<ppcnt)
674  {
675  list.push_back(kPdgProton);
676  RemnZ--;
677  }
678  else list.push_back(kPdgNeutron);
679 
680  RemnA--;
681 
682  ppcnt = (double) RemnZ / (double) RemnA;
683  }
684 
685 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
686  LOG("INukeUtils2018", pDEBUG)
687  << "Remnant nucleus (A,Z) = (" << RemnA << ", " << RemnZ << ")";
688 #endif
689 
690  // Add the fermi energy of the three nucleons to the phase space
691  /* if(DoFermi)
692  {
693  Target target(ev->TargetNucleus()->Pdg());
694  TVector3 pBuf = p->P4()->Vect();
695  double mBuf = p->Mass();
696  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
697  TLorentzVector tSum(pBuf,eBuf);
698  double mSum = 0.0;
699  vector<int>::const_iterator pdg_iter;
700  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
701  {
702  target.SetHitNucPdg(*pdg_iter);
703  Nuclmodel->GenerateNucleon(target);
704  mBuf = pLib->Find(*pdg_iter)->Mass();
705  mSum += mBuf;
706  pBuf = FermiFac * Nuclmodel->Momentum3();
707  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
708  tSum += TLorentzVector(pBuf,eBuf);
709  RemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
710  }
711  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
712  p->SetMomentum(dP4);
713  }
714  */
715  // do the phase space decay & save all f/s particles to the record
716  bool success = genie::utils::intranuke2018::PhaseSpaceDecay(ev,p,list,RemnP4,NucRmvE,mode);
717  if (success) LOG("INukeUtils",pINFO) << "successful equilibrium interaction";
718  else
719  {
720  exceptions::INukeException exception;
721  exception.SetReason("Phase space generation of compound nucleus final state failed, details above");
722  throw exception;
723  }
724 
725 }
726 
727 
728 //___________________________________________________________________________
729 bool genie::utils::intranuke2018::TwoBodyCollision(
730  GHepRecord* ev, int pcode, int tcode, int scode, int s2code, double C3CM,
731  GHepParticle* p, GHepParticle* t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode)
732 {
733  // Aaron Meyer (10/29/09)
734  // Adapted from kinematics in other function calls
735  //
736  // C3CM is the cosine of the scattering angle, calculated before calling
737  // p and t are the output particles, must be allocated before calling
738  // pcode,tcode,scode,s2code are initial and final particle PDG codes in scattering
739  // return value used for error checking
740 
741  // Kinematic variables
742 
743  double M1, /* M2, */ M3, M4; // rest energies, in GeV
744  double E3L, P3L, E4L, P4L;
745  TVector3 tP1L, tPtot, tbeta, tbetadir, tTrans, tVect;
746  TVector3 tP1zCM, tP2zCM, tP3L, tP4L;
747 
748  // Library instance for reference
749  PDGLibrary * pLib = PDGLibrary::Instance();
750 
751  // random number generator
752  //RandomGen * rnd = RandomGen::Instance();
753 
754  // handle fermi target
755  Target target(ev->TargetNucleus()->Pdg());
756 
757  // get mass for particles
758  M1 = pLib->Find(pcode)->Mass();
759  // usused // M2 = pLib->Find(tcode)->Mass();
760  M3 = pLib->Find(scode)->Mass();
761  M4 = pLib->Find(s2code)->Mass();
762 
763  // get lab energy and momenta and assign to 4 vectors
764  TLorentzVector t4P1L = *p->P4();
765  TLorentzVector t4P2L = *t->P4();
766 
767  // binding energy
768  double bindE = 0.025; // empirical choice, might need to be improved
769  //double bindE = 0.0;
770 
771  LOG("TwoBodyCollision",pNOTICE) << "M1 = " << t4P1L.M() << " , M2 = " << t4P2L.M();
772  LOG("TwoBodyCollision",pNOTICE) << "E1 = " << t4P1L.E() << " , E2 = " << t4P2L.E();
773 
774  if ( (p->Energy()-p->Mass()) < bindE ){bindE = 0.;} // if the probe's energy is less than the binding energy, set the binding energy to 0.
775 
776  // don't use BE unless kinetic energy >> BE.
777  if((pcode==2112||pcode==2212)&&(t4P1L.E()-M1)<.1) bindE = 0.0;
778  if((pcode==211||pcode==-211||pcode==111)&&(t4P1L.E()-M1)<.08) bindE = 0.0;
779  if((pcode==321)&&(t4P1L.E()-M1)<.1) bindE = 0.0;
780 
781  // carry out scattering
782  TLorentzVector t4P3L, t4P4L;
783  if (!TwoBodyKinematics(M3,M4,t4P1L,t4P2L,t4P3L,t4P4L,C3CM,RemnP4,bindE))
784  {
785  P3L = t4P3L.Vect().Mag();
786  P4L = t4P4L.Vect().Mag();
787  E3L = t4P3L.E();
788  E4L = t4P4L.E();
789 
790  LOG("TwoBodyCollision",pNOTICE)
791  << "TwoBodyKinematics fails: C3CM, P3 = " << C3CM << " "
792  << P3L << " " << E3L << "\n" << " P4 = "
793  << P4L << " " << E4L ;
794  return false; //covers all possiblities for now
795  }
796 
797  // error checking
798  P3L = t4P3L.Vect().Mag();
799  P4L = t4P4L.Vect().Mag();
800  E3L = t4P3L.E();
801  E4L = t4P4L.E();
802  LOG("INukeUtils",pINFO)
803  << "C3CM, P3 = " << C3CM << " "
804  << P3L << " " << E3L << "\n" << " P4 = "
805  << P4L << " " << E4L ;
806 
807  // handle very low momentum particles
808  if(!(TMath::Finite(P3L)) || P3L<.001)
809  {
810  LOG("INukeUtils",pINFO)
811  << "Particle 3 momentum small or non-finite: " << P3L
812  << "\n" << "--> Assigning .001 as new momentum";
813  P3L = .001;
814  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
815  }
816  if(!(TMath::Finite(P4L)) || P4L<.001)
817  {
818  LOG("INukeUtils",pINFO)
819  << "Particle 4 momentum small or non-finite: " << P4L
820  << "\n" << "--> Assigning .001 as new momentum";
821  P4L = .001;
822  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
823  }
824 
825  // if this is going to be on in the future, remember to not apply PB for Oset
826  /*// pauli blocking turn off for now to better match data
827  // if(P3L<fFermiMomentum && pdg::IsNeutronOrProton(scode) ||
828  // P4L<fFermiMomentum && pdg::IsNeutronOrProton(s2code) )
829  if(P3L<.25 && pdg::IsNeutronOrProton(scode) ||
830  P4L<.25 && pdg::IsNeutronOrProton(s2code) )
831  {
832  LOG("INukeUtils",pNOTICE)<< "TwoBodyCollision failed: Pauli blocking";
833  p->SetStatus(kIStHadronInTheNucleus);
834  RemnP4 -= TLorentzVector(0,0,0,bindE);
835  return false;
836  }*/
837 
838  // update remnant nucleus
839  RemnP4 -= t4P2L;
840  LOG("INukeUtils",pINFO)
841  << "t4P2L= " << t4P2L.E() << " " << t4P2L.Z()
842  << " RemnP4= " << RemnP4.E() << " " << RemnP4.Z() ;
843  if (tcode==kPdgProton) {RemnZ--;RemnA--;}
844  else if(tcode==kPdgNeutron) RemnA--;
845 
846  // create t particle w/ appropriate momenta, code, and status
847  // Set target's mom to be the mom of the hadron that was cloned
848  t->SetFirstMother(p->FirstMother());
849  t->SetLastMother(p->LastMother());
850 
851  // adjust p to reflect scattering
852  p->SetPdgCode(scode);
853  p->SetMomentum(t4P3L);
854 
855  t->SetPdgCode(s2code);
856  t->SetMomentum(t4P4L);
857 
858  if (mode==kIMdHN)
859  {
860  p->SetStatus(kIStHadronInTheNucleus);
861  t->SetStatus(kIStHadronInTheNucleus);
862  }
863  else
864  {
865  p->SetStatus(kIStStableFinalState);
866  t->SetStatus(kIStStableFinalState);
867  }
868  LOG("INukeUtils",pINFO) << "Successful 2 body collision";
869  return true;
870 
871 }
872 //___________________________________________________________________________
873 bool genie::utils::intranuke2018::TwoBodyKinematics(
874  double M3, double M4, TLorentzVector t4P1L, TLorentzVector t4P2L,
875  TLorentzVector &t4P3L, TLorentzVector &t4P4L, double C3CM, TLorentzVector &RemnP4, double bindE)
876 {
877  // Aaron Meyer (05/17/10)
878  // Adapted from kinematics in other function calls
879  //
880  // Outgoing particle masses M3,M4
881  // Scatters particles according to normal two body collisions
882  //
883  // bindE is the binding energy (GeV) of a particle removed from the nucleus (default 0)
884  // For nonzero binding energy, remove the binding energy from the total energy,
885  // then put both of the particles back on mass shell by shifting momentum/energy
886  // of target
887  // Momentum only shifted in the direction parallel to the probe's motion
888  //
889  // Rotates final transverse component of momentum by a random angle from 0->2pi
890  // Return value for error checking
891  // Gives outgoing 4-momenta of particles 3 and 4 (t4P3L, t4P4L respectively)
892  //
893  // All 4-momenta should be on mass shell
894 
895  double E1L, E2L, P1L, P2L, E3L, P3L;
896  double beta, gm; // speed and gamma for CM frame in Lab
897  double S3CM; // sin of scattering angle
898  double PHI3;
899  double E1CM, E2CM, E3CM, P3CM;//, E4CM, P4CM;
900  double P3zL, P3tL;//, P4zL, P4tL;
901  double Et;
902  double theta1, theta2, theta5, P1zL, P2zL, P1tL, P2tL;
903  TVector3 tbeta, tbetadir, tTrans, tVect;
904  TVector3 tP1zCM, tP2zCM, vP3L;
905  TLorentzVector t4P1buf, t4P2buf, t4Ptot;
906 
907  // Library instance for reference
908  //PDGLibrary * pLib = PDGLibrary::Instance();
909 
910  // random number generator
911  RandomGen * rnd = RandomGen::Instance();
912 
913  // error checking
914  if (C3CM < -1. || C3CM > 1.) return false;
915 
916  // calculate sine from scattering angle
917  S3CM = TMath::Sqrt(1.0 - C3CM*C3CM);
918 
919  // fill buffers
920  t4P1buf = t4P1L;
921  t4P2buf = t4P2L;
922 
923  // get lab energy and momenta
924  E1L = t4P1buf.E();
925  P1L = t4P1buf.P();
926  E2L = t4P2buf.E();
927  P2L = t4P2buf.P();
928  t4Ptot = t4P1buf + t4P2buf;
929 
930  LOG("INukeUtils",pINFO) <<"M1 "<<t4P1L.M()<< ", M2 "<<t4P2L.M();
931  LOG("INukeUtils",pINFO) <<"E1L "<<E1L<< ", E1CM "<<E1CM;
932  LOG("INukeUtils",pINFO) <<"bindE = " << bindE;
933 
934  // binding energy
935  if (bindE!=0)
936  {
937 
938  E1L -= bindE;
939 
940  if (E1L+E2L < M3+M4)
941  {
942  LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Forbidden by binding energy";
943  LOG("INukeUtils",pNOTICE) <<"E1L, E2L, M3, M4 : "<< E1L <<", "<< E2L <<", "<< M3 <<", "<< M4;
944  t4P3L.SetPxPyPzE(0,0,0,0);
945  t4P4L.SetPxPyPzE(0,0,0,0);
946  return false;
947  }
948  }
949 
950  // calculate beta and gamma
951  tbeta = t4Ptot.Vect() * (1.0 / (E1L + E2L));
952  tbetadir = tbeta.Unit();
953  beta = tbeta.Mag();
954  gm = 1.0 / TMath::Sqrt(1.0 - beta*beta);
955 
956  // get angle and component info
957  theta1 = t4P1buf.Angle(tbeta);
958  theta2 = t4P2buf.Angle(tbeta);
959  P1zL = P1L*TMath::Cos(theta1);
960  P2zL = P2L*TMath::Cos(theta2);
961  P1tL = TMath::Sqrt(P1L*P1L - P1zL*P1zL);
962  P2tL = -TMath::Sqrt(P2L*P2L - P2zL*P2zL);
963  tVect.SetXYZ(1,0,0);
964  if(TMath::Abs((tVect - tbetadir).Mag())<.01) tVect.SetXYZ(0,1,0);
965  theta5 = tVect.Angle(tbetadir);
966  tTrans = (tVect - TMath::Cos(theta5)*tbetadir).Unit();
967 
968  // boost to CM frame to get scattered particle momenta
969  E1CM = gm*E1L - gm*beta*P1zL;
970  tP1zCM = gm*P1zL*tbetadir - gm*tbeta*E1L;
971  E2CM = gm*E2L - gm*beta*P2zL;
972  tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
973  Et = E1CM + E2CM;
974 //-------
975 
976  LOG("INukeUtils",pINFO) <<"M1 "<<t4P1L.M()<< ", M2 "<<t4P2L.M();
977  LOG("INukeUtils",pINFO) <<"E1L "<<E1L<< ", E1CM "<<E1CM;
978  LOG("INukeUtils",pINFO) <<"P1zL "<<P1zL<<", P1zCM "<<tP1zCM.Mag()<<", P1tL "<<P1tL;
979  LOG("INukeUtils",pINFO) <<"E2L "<<E2L<< ", E2CM "<<E2CM;
980  LOG("INukeUtils",pINFO) <<"P2zL "<<P2zL<<", P2zCM "<<tP2zCM.Mag()<<", P2tL "<<P2tL;
981  LOG("INukeUtils",pINFO) <<"C3CM "<<C3CM;
982 
983 //-------
984  E3CM = (Et*Et + M3*M3 - M4*M4) / (2.0 * Et);
985 
986  // check to see if decay is viable
987  if(E3CM*E3CM - M3*M3<0 || E3CM<0 || Et<0)
988  {
989  if (Et<0) LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Total energy is negative";
990  if (E3CM<M3) LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Scattered Particle 3 CM energy is too small";
991  if (E3CM*E3CM - M3*M3<0) LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Scattered Particle 3 CM momentum is nonreal";
992  t4P3L.SetPxPyPzE(0,0,0,0);
993  t4P4L.SetPxPyPzE(0,0,0,0);
994  return false;
995  }
996  P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
997 
998  // boost back to lab
999  P3zL = gm*beta*E3CM + gm*P3CM*C3CM;
1000  P3tL = P3CM*S3CM;
1001 
1002  P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
1003  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1004 
1005  //-------
1006 
1007  double E4CM = Et-E3CM;
1008  double P4zL = gm*beta*E4CM - gm*P3CM*C3CM;
1009  double P4tL = -1.*P3tL;
1010  double P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
1011  double E4L = TMath::Sqrt(P4L*P4L + M4*M4);
1012 
1013  LOG("INukeUtils",pINFO) <<"M3 "<< M3 << ", M4 "<< M4;
1014  LOG("INukeUtils",pINFO) <<"E3L "<<E3L<< ", E3CM "<<E3CM;
1015  LOG("INukeUtils",pINFO) <<"P3zL "<<P3zL<<", P3tL "<<P3tL;
1016  LOG("INukeUtils",pINFO) <<"C3L "<<P3zL/P3L;
1017  LOG("INukeUtils",pINFO) <<"Check:";
1018  LOG("INukeUtils",pINFO) <<"E4L "<<E4L<< ", E4CM "<<E4CM;
1019  LOG("INukeUtils",pINFO) <<"P4zL "<<P4zL<<", P4tL "<<P4tL;
1020  LOG("INukeUtils",pINFO) <<"P4L "<<P4L;
1021  LOG("INukeUtils",pINFO) <<"C4L "<<P4zL/P4L;
1022 
1023  double echeck = E1L + E2L - (E3L + E4L);
1024  double pzcheck = P1zL+ P2zL - (P3zL + P4zL);
1025  double ptcheck = P1tL+ P2tL - (P3tL + P4tL);
1026 
1027  LOG("INukeUtils",pINFO) <<"Check 4-momentum conservation - Energy "<<echeck<<", z momentum "<<pzcheck << ", transverse momentum " << ptcheck ;
1028 
1029  // -------
1030 
1031  // handle very low momentum particles
1032  if(!(TMath::Finite(P3L)) || P3L<.001)
1033  {
1034  LOG("INukeUtils",pINFO)
1035  << "Particle 3 momentum small or non-finite: " << P3L
1036  << "\n" << "--> Assigning .001 as new momentum";
1037  P3tL = 0;
1038  P3zL = .001;
1039  P3L = .001;
1040  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1041  }
1042 
1043  // get random phi angle, distributed uniformally in 360 deg
1044  PHI3 = 2 * kPi * rnd->RndFsi().Rndm();
1045 
1046  vP3L = P3zL*tbetadir + P3tL*tTrans;
1047  vP3L.Rotate(PHI3,tbetadir);
1048 
1049  t4P3L.SetVect(vP3L);
1050  t4P3L.SetE(E3L);
1051 
1052  t4P4L = t4P1buf + t4P2buf - t4P3L;
1053  t4P4L-= TLorentzVector(0,0,0,bindE);
1054  /*LOG("INukeUtils",pINFO) <<"GENIE:";
1055  LOG("INukeUtils",pINFO) <<"E4L "<<t4P4L.E();
1056  LOG("INukeUtils",pINFO) <<"P4zL "<<t4P4L.Vect()*tbetadir<<", P4tL "<<-1.*TMath::Sqrt(t4P4L.Vect().Mag2()-TMath::Power(t4P4L.Vect()*tbetadir,2.));
1057  LOG("INukeUtils",pINFO) <<"P4L "<<t4P4L.Vect().Mag();
1058  LOG("INukeUtils",pINFO) <<"C4L "<<t4P4L.Vect()*tbetadir/t4P4L.Vect().Mag(); */
1059 
1060  if(t4P4L.Mag2()<0 || t4P4L.E()<0)
1061  {
1062  LOG("INukeUtils",pNOTICE)<<"TwoBodyKinematics Failed: Target mass or energy is negative";
1063  t4P3L.SetPxPyPzE(0,0,0,0);
1064  t4P4L.SetPxPyPzE(0,0,0,0);
1065  return false;
1066  }
1067 
1068  if (bindE!=0) RemnP4 += TLorentzVector(0,0,0,bindE);
1069  return true;
1070 }
1071 //___________________________________________________________________________
1072 bool genie::utils::intranuke2018::ThreeBodyKinematics(
1073  GHepRecord* ev, GHepParticle* p, int tcode, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3,
1074  bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
1075 {
1076 
1077  // Aaron Meyer (7/15/10)
1078  //
1079  // Kinematics used in utils::intranuke2018::PionProduction
1080  // Handles the kinematics of three body scattering
1081  //
1082  // s1,s2,and s3 are pointers to particles with the PDG code that needs to be scattered
1083  // the last four variables are for Fermi momentum and pauli blocking, default will use none of them
1084 
1085  // kinematic variables
1086  double M1, M2, M3, M4, M5; // rest energies, in GeV
1087  double P1L, P2L, P3L, P4L, P5L;
1088  double E1L, E2L, E3L, E4L, E5L;
1089  double E1CM, E2CM, P3tL;
1090  double PizL, PitL, PiL, EiL;
1091  double EiCM, P4CM2, E4CM2, E5CM2, P3CM, E3CM;
1092  double beta, gm, beta2, gm2;
1093  double P3zL, P4zL, P4tL, P5zL, P5tL;
1094  double Et, M, theta1, theta2;
1095  double P1zL, P2zL;
1096  double theta3, theta4, phi3, phi4, theta5;
1097  TVector3 tP2L, tP1L, tPtot, tbeta, tbetadir, tTrans, tP4L, tP5L;
1098  TVector3 tP1zCM, tP2zCM, tP3L, tPiL, tbeta2, tbetadir2, tVect, tTrans2;
1099 
1100  // Library instance for reference
1101  PDGLibrary * pLib = PDGLibrary::Instance();
1102 
1103  // random number generator
1104  RandomGen * rnd = RandomGen::Instance();
1105 
1106  M1 = pLib->Find(p->Pdg())->Mass();
1107  M2 = pLib->Find(tcode)->Mass();
1108  M3 = pLib->Find(s1->Pdg())->Mass();
1109  M4 = pLib->Find(s2->Pdg())->Mass();
1110  M5 = pLib->Find(s3->Pdg())->Mass();
1111 
1112  // set up fermi target
1113  Target target(ev->TargetNucleus()->Pdg());
1114 
1115  // handle fermi momentum
1116  if(DoFermi)
1117  {
1118  target.SetHitNucPdg(tcode);
1119  Nuclmodel->GenerateNucleon(target);
1120  tP2L = FermiFac * Nuclmodel->Momentum3();
1121  P2L = tP2L.Mag();
1122  E2L = TMath::Sqrt(tP2L.Mag2() + M2*M2);
1123  }
1124  else
1125  {
1126  tP2L.SetXYZ(0.0, 0.0, 0.0);
1127  P2L = 0;
1128  E2L = M2;
1129  }
1130 
1131  // first sequence, handle 4th and 5th products as composite
1132  E1L = p->E();
1133 
1134  P1L = TMath::Sqrt(E1L*E1L - M1*M1);
1135  tP1L = p->P4()->Vect();
1136  tPtot = tP1L + tP2L;
1137 
1138  tbeta = tPtot * (1.0 / (E1L + E2L));
1139  tbetadir = tbeta.Unit();
1140  beta = tbeta.Mag();
1141  gm = 1.0 / TMath::Sqrt(1.0 - beta*beta);
1142 
1143  theta1 = tP1L.Angle(tbeta);
1144  theta2 = tP2L.Angle(tbeta);
1145  P1zL = P1L*TMath::Cos(theta1);
1146  P2zL = P2L*TMath::Cos(theta2);
1147  tVect.SetXYZ(1,0,0);
1148  if(TMath::Abs((tVect - tbetadir).Mag())<.01) tVect.SetXYZ(0,1,0);
1149  theta5 = tVect.Angle(tbetadir);
1150  tTrans = (tVect - TMath::Cos(theta5)*tbetadir).Unit();
1151 
1152  E1CM = gm*E1L - gm*beta*P1zL;
1153  tP1zCM = gm*P1zL*tbetadir - gm*tbeta*E1L;
1154  E2CM = gm*E2L - gm*beta*P2zL;
1155  tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
1156  Et = E1CM + E2CM;
1157  M = (rnd->RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
1158  E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
1159  EiCM = Et - E3CM;
1160  if(E3CM*E3CM - M3*M3<0)
1161  {
1162  LOG("INukeUtils",pNOTICE)
1163  << "PionProduction P3 has non-real momentum - retry kinematics";
1164  LOG("INukeUtils",pNOTICE) << "Energy, masses of 3 fs particales:"
1165  << E3CM << " " << M3 << " " << " " << M4 << " " << M5;
1166  exceptions::INukeException exception;
1167  exception.SetReason("PionProduction particle 3 has non-real momentum");
1168  throw exception;
1169  return false;
1170  }
1171  P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
1172 
1173  theta3 = kPi * rnd->RndFsi().Rndm();
1174  theta4 = kPi * rnd->RndFsi().Rndm();
1175  phi3 = 2*kPi * rnd->RndFsi().Rndm();
1176  phi4 = 2*kPi * rnd->RndFsi().Rndm();
1177 
1178  P3zL = gm*beta*E3CM + gm*P3CM*TMath::Cos(theta3);
1179  P3tL = P3CM*TMath::Sin(theta3);
1180  PizL = gm*beta*EiCM - gm*P3CM*TMath::Cos(theta3);
1181  PitL = -P3CM*TMath::Sin(theta3);
1182 
1183  P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
1184  PiL = TMath::Sqrt(PizL*PizL + PitL*PitL);
1185  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1186  EiL = TMath::Sqrt(PiL*PiL + M*M);
1187 
1188  // handle very low momentum particles
1189  if(!(TMath::Finite(P3L)) || P3L < .001)
1190  {
1191  LOG("INukeUtils",pINFO)
1192  << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
1193  << "\n" << "--> Assigning .001 as new momentum";
1194  P3tL = 0;
1195  P3zL = .001;
1196  P3L = .001;
1197  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1198  }
1199 
1200  tP3L = P3zL*tbetadir + P3tL*tTrans;
1201  tPiL = PizL*tbetadir + PitL*tTrans;
1202  tP3L.Rotate(phi3,tbetadir);
1203  tPiL.Rotate(phi3,tbetadir);
1204 
1205  // second sequence, handle formally composite particles 4 and 5
1206  tbeta2 = tPiL * (1.0 / EiL);
1207  tbetadir2 = tbeta2.Unit();
1208  beta2 = tbeta2.Mag();
1209  gm2 = 1.0 / TMath::Sqrt(1.0 - beta2*beta2);
1210 
1211  E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
1212  E5CM2 = M - E4CM2;
1213  P4CM2 = TMath::Sqrt(E4CM2*E4CM2 - M4*M4);
1214 
1215  tVect.SetXYZ(1,0,0);
1216  if(TMath::Abs((tVect - tbetadir2).Mag())<.01) tVect.SetXYZ(0,1,0);
1217  theta5 = tVect.Angle(tbetadir2);
1218  tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).Unit();
1219 
1220  P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*TMath::Cos(theta4);
1221  P4tL = P4CM2*TMath::Sin(theta4);
1222  P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*TMath::Cos(theta4);
1223  P5tL = - P4tL;
1224 
1225  P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
1226  P5L = TMath::Sqrt(P5zL*P5zL + P5tL*P5tL);
1227  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
1228  E5L = TMath::Sqrt(P5L*P5L + M5*M5);
1229 
1230  // handle very low momentum particles
1231  if(!(TMath::Finite(P4L)) || P4L < .001)
1232  {
1233  LOG("INukeUtils",pINFO)
1234  << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
1235  << "\n" << "--> Assigning .001 as new momentum";
1236  P4tL = 0;
1237  P4zL = .001;
1238  P4L = .001;
1239  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
1240  }
1241  if(!(TMath::Finite(P5L)) || P5L < .001)
1242  {
1243  LOG("INukeUtils",pINFO)
1244  << "Particle 5 " << M5 << " momentum small or non-finite: " << P5L
1245  << "\n" << "--> Assigning .001 as new momentum";
1246  P5tL = 0;
1247  P5zL = .001;
1248  P5L = .001;
1249  E5L = TMath::Sqrt(P5L*P5L + M5*M5);
1250  }
1251 
1252  tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
1253  tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
1254  tP4L.Rotate(phi4,tbetadir2);
1255  tP5L.Rotate(phi4,tbetadir2);
1256 
1257  // pauli blocking
1258  if(P3L < FermiMomentum || ( pdg::IsNeutronOrProton(s2->Pdg()) && P4L < FermiMomentum ) )
1259  {
1260  LOG("INukeUtils",pNOTICE)
1261  << "PionProduction fails because of Pauli blocking - retry kinematics";
1262  exceptions::INukeException exception;
1263  exception.SetReason("PionProduction final state not determined");
1264  throw exception;
1265  return false;
1266  }
1267 
1268  // create scattered particles w/ appropriate momenta, code, and status
1269  // Set moms to be the moms of the hadron that was cloned
1270  s1->SetFirstMother(p->FirstMother());
1271  s2->SetFirstMother(p->FirstMother());
1272  s3->SetFirstMother(p->FirstMother());
1273  s1->SetLastMother(p->LastMother());
1274  s2->SetLastMother(p->LastMother());
1275  s3->SetLastMother(p->LastMother());
1276 
1277  TLorentzVector(tP3L,E3L);
1278  TLorentzVector(tP4L,E4L);
1279  TLorentzVector(tP5L,E5L);
1280 
1281  s1->SetMomentum(TLorentzVector(tP3L,E3L));
1282  s2->SetMomentum(TLorentzVector(tP4L,E4L));
1283  s3->SetMomentum(TLorentzVector(tP5L,E5L));
1284  int mode = kIMdHA;
1285  LOG ("INukeUtils",pDEBUG) << "in Pi Prod, mode = " << mode;
1286  if (mode==kIMdHN)
1287  {
1288  s1->SetStatus(kIStHadronInTheNucleus);
1289  s2->SetStatus(kIStHadronInTheNucleus);
1290  s3->SetStatus(kIStHadronInTheNucleus);
1291  }
1292  else
1293  {
1294  s1->SetStatus(kIStStableFinalState);
1295  s2->SetStatus(kIStStableFinalState);
1296  s3->SetStatus(kIStStableFinalState);
1297  }
1298  return true;
1299 }
1300 //___________________________________________________________________________
1301 bool genie::utils::intranuke2018::PionProduction(
1302  GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ,
1303  TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
1304 {
1305  // Aaron Meyer (7/15/2010)
1306  //
1307  // Handles pion production reactions in both hA and hN mode
1308  // Calculates fundamental cross sections from fit functions
1309  // Uses isospin relations to determine the rest of cross sections
1310  //
1311  // p is the probe particle
1312  // s1, s2, and s3 are the particles produced in the reaction
1313  // must set the status and add particles to the event record after returning from this method
1314  // return value for error checking
1315 
1316 
1317  // random number generator
1318  RandomGen * rnd = RandomGen::Instance();
1319 
1320  // library reference
1321  PDGLibrary * pLib = PDGLibrary::Instance();
1322 
1323  bool ptarg = false;
1324  int pcode = p->Pdg();
1325 
1326  int p1code = p->Pdg();
1327  // Randomly determine target and 1st product baryons
1328  int p3code = 0, p4code = 0, p5code = 0;
1329 
1330  //
1331  // Uses a fit curve log(sigma) = a - b/(T_pi - c) for pions
1332  // Fit parameters determined by Roman Tacik (4/3/09)
1333  // pi- & p cross sections are assumed to be the same as pi+ & n
1334  //
1335  // Fit curve for nucleons:
1336  // sigma = a*(1+b*e^(-c*(eta-d)^2))*(1-e^(-(f*eta)^g))*(1-e^(-h/eta^2))
1337  // 7 parameters (a,b,c,d,f,g,h)
1338  // eta is maximum kinematically allowed momentum of the pion, normalized by the mass
1339  // Uses isotopic spin decomposition of total cross sections
1340  //
1341 
1342  if ((p1code==kPdgPi0)||(p1code==kPdgPiP)||(p1code==kPdgPiM)) {
1343 
1344  double kine = 1000*p->KinE();
1345 
1346  // Determine cross sections
1347 
1348  // pion
1349  // pi- & p
1350  // -> pi0 & pi0 & n
1351  // a = 8.82; b = 573.2; c = 107.3;
1352  double xsec2pi0n = TMath::Max(0.,TMath::Exp(8.82 - (573.2/(kine-107.3))));
1353  // -> pi- & pi+ & n
1354  // a = 11.06; b = 985.9; c = 88.2;
1355  double xsecpippimn = TMath::Max(0.,TMath::Exp(11.06 - (985.9/(kine-88.2))));
1356  // -> pi- & pi0 & p
1357  // a = 9.58; b = 1229.4; c = 60.5;
1358  double xsecpimpi0p = TMath::Max(0.,TMath::Exp(9.58 - (1229.4/(kine-60.5))));
1359  double totpimp = xsec2pi0n + xsecpippimn + xsecpimpi0p;
1360 
1361 
1362  // pi+ & p
1363  // -> pi+ & pi+ & n
1364  // a = 5.64; b = 222.6; c = 150.0;
1365  double xsec2pipn = TMath::Max(0.,TMath::Exp(5.64 - (222.6/(kine-150.))));
1366  // -> pi+ & pi0 & p
1367  // a = 7.95; b = 852.6; c = 77.8;
1368  double xsecpippi0p = TMath::Max(0.,TMath::Exp(7.95 - (852.6/(kine-77.8))));
1369  double totpipp = xsec2pipn + xsecpippi0p;
1370 
1371  if (totpimp<=0 && totpipp<=0) {
1372  LOG("INukeUtils",pNOTICE) << "InelasticHN called below threshold energy";
1373  p->SetStatus(kIStHadronInTheNucleus);
1374  ev->AddParticle(*p);
1375  return false;
1376  }
1377 
1378  double xsecp, xsecn;
1379  switch (p1code) {
1380  case kPdgPi0: xsecp = 0.5 * (totpimp + totpipp); xsecn = xsecp; break;
1381  case kPdgPiP: xsecp = totpipp; xsecn = totpimp; break;
1382  case kPdgPiM: xsecp = totpimp; xsecn = totpipp; break;
1383  default:
1384  LOG("INukeUtils",pWARN) << "InelasticHN cannot handle probe: "
1385  << PDGLibrary::Instance()->Find(p1code)->GetName();
1386  exceptions::INukeException exception;
1387  exception.SetReason("PionProduction final state not determined");
1388  throw exception;
1389  return false;
1390  break;
1391  }
1392 
1393  // Normalize cross sections by Z or A-Z
1394 
1395  xsecp *= RemnZ;
1396  xsecn *= RemnA-RemnZ;
1397 
1398  // determine target
1399 
1400  double rand = rnd->RndFsi().Rndm() * (xsecp + xsecn);
1401  if (rand < xsecp) // proton target
1402  { rand /= RemnZ; ptarg = true;}
1403  else // neutron target
1404  { rand -= xsecp; rand /= RemnA-RemnZ; ptarg = false;}
1405 
1406 
1407  // determine final state
1408 
1409  if (((ptarg==true)&&(p1code==kPdgPiP))
1410  || ((ptarg==false)&&(p1code==kPdgPiM)))
1411  {
1412  if (rand < xsec2pipn) // pi+ & pi+ & n final state
1413  {
1414  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1415  p4code = p1code;
1416  p5code = p4code;
1417  }
1418  else { // pi+ & pi0 & p final state
1419  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1420  p4code = p1code;
1421  p5code = kPdgPi0;
1422  }
1423  }
1424  else if (((ptarg==false)&&(p1code==kPdgPiP))
1425  || ((ptarg==true)&&(p1code==kPdgPiM)))
1426  {
1427  if (rand < xsec2pi0n) // pi0 & pi0 & n final state
1428  {
1429  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1430  p4code = kPdgPi0;
1431  p5code = p4code;
1432  }
1433  else if (rand < (xsec2pi0n + xsecpippimn)) // pi+ & pi- & n final state
1434  {
1435  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1436  p4code = p1code;
1437  p5code = ((p1code==kPdgPiP) ? kPdgPiM : kPdgPiP);
1438  }
1439  else // pi0 & pi- & p final state
1440  {
1441  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1442  p4code = p1code;
1443  p5code = kPdgPi0;
1444  }
1445  }
1446  else if (p1code==kPdgPi0)
1447  {
1448  rand = rnd->RndFsi().Rndm();
1449  if (rand < 191./270.)
1450  { // pi+ & pi- & p final state
1451  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1452  p4code = kPdgPiP;
1453  p5code = kPdgPiM;
1454  }
1455  else if (rand < 7./135.)
1456  { // pi0 & pi0 & p final state
1457  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1458  p4code = kPdgPi0;
1459  p5code = p4code;
1460  }
1461  else
1462  { // pi+ & pi0 & n final state
1463  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1464  p4code = (ptarg ? kPdgPiP : kPdgPiM);
1465  p5code = kPdgPi0;
1466  }
1467  }
1468  else // unhandled
1469  {
1470  LOG("INukeUtils",pNOTICE) << "Pi production final state unable to be determined, picode, ptarg = " <<PDGLibrary::Instance()->Find(p1code)->GetName() << " " << PDGLibrary::Instance()->Find(ptarg)->GetName();
1471  exceptions::INukeException exception;
1472  exception.SetReason("PionProduction final state not determined");
1473  throw exception;
1474  return false;
1475  }
1476 
1477  } else if(p1code==kPdgProton||p1code==kPdgNeutron) //nucleon probes
1478  {
1479 
1480  double tote = p->Energy();
1481  double pMass = pLib->Find(2212)->Mass();
1482  double nMass = pLib->Find(2112)->Mass();
1483  double etapp2ppPi0 =
1484  utils::intranuke2018::CalculateEta(pMass,tote,pMass,pMass+pMass,pLib->Find(111)->Mass());
1485  double etapp2pnPip =
1486  utils::intranuke2018::CalculateEta(pLib->Find(p1code)->Mass(),tote,((p1code==kPdgProton)?pMass:nMass),
1487  pMass+nMass,pLib->Find(211)->Mass());
1488  double etapn2nnPip =
1489  utils::intranuke2018::CalculateEta(pMass,tote,nMass,nMass+nMass,pLib->Find(211)->Mass());
1490  double etapn2ppPim =
1491  utils::intranuke2018::CalculateEta(pMass,tote,nMass,pMass+pMass,pLib->Find(211)->Mass());
1492 
1493  if ((etapp2ppPi0<=0.)&&(etapp2pnPip<=0.)&&(etapn2nnPip<=0.)&&(etapn2ppPim<=0.)) { // below threshold
1494  LOG("INukeUtils",pNOTICE) << "PionProduction() called below threshold energy";
1495  exceptions::INukeException exception;
1496  exception.SetReason("PionProduction final state not possible - below threshold");
1497  throw exception;
1498  return false;
1499  }
1500 
1501  // calculate cross sections
1502  double xsecppPi0=0,xsecpnPiP=0,xsecnnPiP=0,xsecppPiM=0;
1503  if (etapp2ppPi0>0){
1504  xsecppPi0 = 4511*(1-.91*TMath::Exp(-TMath::Power((etapp2ppPi0-.705),2)));
1505  xsecppPi0 *= (1-TMath::Exp(-TMath::Power((.556*etapp2ppPi0),3.5)));
1506  xsecppPi0 *= (1-TMath::Exp(-56.897/(etapp2ppPi0*etapp2ppPi0)));
1507  xsecppPi0 = TMath::Max(0.,xsecppPi0);}
1508 
1509  if (etapp2pnPip>0){
1510  xsecpnPiP = 18840*(1-.808*TMath::Exp(-TMath::Power((etapp2pnPip-.371),2)));
1511  xsecpnPiP *= (1-TMath::Exp(-TMath::Power((.568*etapp2pnPip),3.2)));
1512  xsecpnPiP *= (1-TMath::Exp(-39.818/(etapp2pnPip*etapp2pnPip)));
1513  xsecpnPiP = TMath::Max(0.,xsecpnPiP);}
1514 
1515  if (etapn2nnPip>0){
1516  xsecnnPiP = 7670*(1-.479*TMath::Exp(-TMath::Power((etapn2nnPip-.947),2)));
1517  xsecnnPiP *= (1-TMath::Exp(-TMath::Power((.35*etapn2nnPip),3.2)));
1518  xsecnnPiP *= (1-TMath::Exp(-71.53/(etapn2nnPip*etapn2nnPip)));
1519  xsecnnPiP = TMath::Max(0.,xsecnnPiP);}
1520 
1521  if (etapn2ppPim>0){
1522  xsecppPiM = 7670*(1-.479*TMath::Exp(-TMath::Power((etapn2ppPim-.947),2)));
1523  xsecppPiM *= (1-TMath::Exp(-TMath::Power((.35*etapn2ppPim),3.2)));
1524  xsecppPiM *= (1-TMath::Exp(-71.53/(etapn2ppPim*etapn2ppPim)));
1525  xsecppPiM = TMath::Max(0.,xsecppPiM);}
1526 
1527  // double sigma11 = xsecppPi0;
1528  double sigma10 = TMath::Max(0.,xsecpnPiP - xsecppPi0); // Fundamental cross sections
1529  double sigma01 = TMath::Max(0.,xsecppPiM + xsecnnPiP - xsecppPi0);
1530 
1531  double xsecpnPi0 = .5*(sigma10 + sigma01);
1532  xsecpnPi0 = TMath::Max(xsecpnPi0,0.);
1533 
1534  LOG("INukeUtils",pDEBUG) << '\n' << "Cross section values: "<<'\n'
1535  << xsecppPi0 << " PP pi0" <<'\n'
1536  << xsecpnPiP << " PN pi+" <<'\n'
1537  << xsecnnPiP << " NN pi+" <<'\n'
1538  << xsecpnPi0 << " PN pi0";
1539 
1540  double xsecp=0,xsecn=0;
1541  switch (p1code) {
1542  case kPdgProton: xsecp=xsecppPi0+xsecpnPiP; xsecn=xsecppPiM+xsecnnPiP+xsecpnPi0; break;
1543  case kPdgNeutron: xsecp=xsecppPiM+xsecnnPiP+xsecpnPi0; xsecn=xsecppPi0+xsecpnPiP; break;
1544  default:
1545  LOG("INukeUtils",pWARN) << "InelasticHN cannot handle probe: "
1546  << PDGLibrary::Instance()->Find(p1code)->GetName();
1547  return false;
1548  break;
1549  }
1550 
1551  // Normalize cross sections by Z or (A-Z)
1552 
1553  xsecp *= RemnZ;
1554  xsecn *= RemnA-RemnZ;
1555 
1556  // determine target
1557 
1558  double rand = rnd->RndFsi().Rndm() * (xsecp + xsecn);
1559  if (rand < xsecp) // proton target
1560  { rand /= RemnZ; ptarg = true;}
1561  else // neutron target
1562  { rand -= xsecp; rand /= RemnA-RemnZ; ptarg = false;}
1563 
1564  if(p1code==kPdgProton) // Cross sections not explicitly given are calculated from isospin relations
1565  {
1566  if(ptarg)
1567  {
1568  if (rand<xsecppPi0) {p3code=kPdgProton; p4code=kPdgProton; p5code=kPdgPi0;}
1569  else {p3code=kPdgProton; p4code=kPdgNeutron; p5code=kPdgPiP;}
1570  }
1571  else
1572  {
1573  if (rand<xsecnnPiP) {p3code=kPdgNeutron; p4code=kPdgNeutron; p5code=kPdgPiP;}
1574  else if (rand<xsecppPiM+xsecnnPiP) {p3code=kPdgProton; p4code=kPdgProton; p5code=kPdgPiM;}
1575  else {p3code=kPdgProton; p4code=kPdgNeutron; p5code=kPdgPi0;}
1576  }
1577  }
1578  else if(p1code==kPdgNeutron)
1579  {
1580  if(ptarg)
1581  {
1582  if (rand<xsecnnPiP) {p3code=kPdgNeutron; p4code=kPdgNeutron; p5code=kPdgPiP;}
1583  else if (rand<xsecppPiM+xsecnnPiP) {p3code=kPdgProton; p4code=kPdgProton; p5code=kPdgPiM;}
1584  else {p3code=kPdgProton; p4code=kPdgNeutron; p5code=kPdgPi0;}
1585  }
1586  else
1587  {
1588  if (rand<xsecpnPiP) {p3code=kPdgNeutron; p4code=kPdgProton; p5code=kPdgPiM;}
1589  else {p3code=kPdgNeutron; p4code=kPdgNeutron; p5code=kPdgPi0;}
1590  }
1591  }
1592  }
1593  else
1594  {
1595  LOG("INukeUtils",pWARN)
1596  << "Unable to handle probe (=" << p1code << ") in InelasticHN()";
1597  return false;
1598  }
1599 
1600  // determine if reaction is allowed
1601  if ( RemnA < 1 )
1602  {
1603  LOG("INukeUtils",pNOTICE) << "PionProduction() failed : no nucleons to produce pions";
1604  return false;
1605  }
1606  else if ( RemnZ + ((pcode==kPdgProton || pcode==kPdgPiP)?1:0) - ((pcode==kPdgPiM)?1:0)
1607  < ((p3code==kPdgProton || p3code==kPdgPiP)?1:0) - ((p3code==kPdgPiM)?1:0)
1608  + ((p4code==kPdgProton || p4code==kPdgPiP)?1:0) - ((p4code==kPdgPiM)?1:0)
1609  + ((p5code==kPdgProton || p5code==kPdgPiP)?1:0) - ((p5code==kPdgPiM)?1:0) )
1610  {
1611  LOG("INukeUtils",pNOTICE) << "PionProduction() failed : too few protons in nucleus";
1612  exceptions::INukeException exception;
1613  exception.SetReason("PionProduction fails - too few protons available");
1614  throw exception;
1615  return false;
1616  }
1617 
1618  s1->SetPdgCode(p3code);
1619  s2->SetPdgCode(p4code);
1620  s3->SetPdgCode(p5code);
1621 
1622  if(genie::utils::intranuke2018::ThreeBodyKinematics(
1623  ev,p,(ptarg?kPdgProton:kPdgNeutron),s1,s2,s3,DoFermi,FermiFac,FermiMomentum,Nuclmodel))
1624  {
1625  // okay, handle remnants and return true
1626  // assumes first particle is always the nucleon,
1627  // second can be either nucleon or pion
1628  // last always pion
1629  if (pcode==kPdgProton || pcode==kPdgPiP) RemnZ++;
1630  if (pcode==kPdgPiM) RemnZ--;
1631  if (pdg::IsPion(pcode)) RemnA--;
1632  if (pdg::IsProton(p3code)) RemnZ--;
1633  if (pdg::IsNeutronOrProton(p4code)) RemnA--;
1634  if (p4code==kPdgPiP || p4code==kPdgProton) RemnZ--;
1635  if (p4code==kPdgPiM) RemnZ++;
1636  if (p5code==kPdgPiP) RemnZ--;
1637  if (p5code==kPdgPiM) RemnZ++;
1638 
1639  LOG("INukeUtils",pDEBUG) << "Remnant (A,Z) = (" <<RemnA<<','<<RemnZ<<')';
1640 
1641  RemnP4 -= *s1->P4() + *s2->P4() + *s3->P4() - *p->P4();
1642  return true;
1643  }
1644  else {
1645  exceptions::INukeException exception;
1646  exception.SetReason("PionProduction final state not determined");
1647  throw exception;
1648  return false;
1649  }
1650 }
1651 //___________________________________________________________________________
1652 double genie::utils::intranuke2018::CalculateEta(double Minc, double nrg, double Mtarg,
1653  double Mtwopart, double Mpi)
1654 {
1655  //Aaron Meyer (1/20/2010)
1656 
1657  //Used to calculate the maximum kinematically allowed CM frame pion momentum
1658  // ke in MeV, eta normalized by pion mass
1659  // approximated by taking two ejected nucleons to be one particle of the same mass
1660  //For pion cross sections, in utils::intranuke2018::PionProduction
1661 
1662  //LOG("INukeUtils",pDEBUG) << "Input values: "<<Minc<<' '<<nrg<<' '<<Mtarg<<' '<<Mtwopart<<' '<<Mpi;
1663  double Scm = Minc*Minc + Mtarg*Mtarg + 2*Mtarg*nrg;
1664  double eta = 0;
1665  eta= TMath::Power(Scm,2) + TMath::Power(Mtwopart,4) + TMath::Power(Mpi,4);
1666  eta-= 2*TMath::Power(Mtwopart*Mpi,2);
1667  eta-= 2*Scm*TMath::Power(Mtwopart,2);
1668  eta-= 2*Scm*TMath::Power(Mpi,2);
1669  eta = TMath::Power(eta/Scm,1./2.);
1670  eta/= (2*Mpi);
1671 
1672  eta=TMath::Max(eta,0.);
1673  return eta;
1674 }
1675 //___________________________________________________________________________
1676 // Generic Phase Space Decay methods
1677 //___________________________________________________________________________
1678 bool genie::utils::intranuke2018::PhaseSpaceDecay(
1679  GHepRecord* ev, GHepParticle* p, const PDGCodeList & pdgv,
1680  TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode)
1681 {
1682 // General method decaying the input particle system 'pdgv' with available 4-p
1683 // given by 'pd'. The decayed system is used to populate the input GHepParticle
1684 // array starting from the slot 'offset'.
1685 //
1686  LOG("INukeUtils", pINFO) << "*** Performing a Phase Space Decay";
1687  assert(pdgv.size() > 1);
1688 
1689  LOG("INukeUtils",pINFO) << "probe mass: M = " << p->Mass();
1690 
1691  // Get the decay product masses & names
1692 
1693  ostringstream state_sstream;
1694  state_sstream << "( ";
1695  vector<int>::const_iterator pdg_iter;
1696  int i = 0;
1697  double * mass = new double[pdgv.size()];
1698  double mass_sum = 0;
1699  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1700  int pdgc = *pdg_iter;
1701  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
1702  string nm = PDGLibrary::Instance()->Find(pdgc)->GetName();
1703  mass[i++] = m;
1704  mass_sum += m;
1705  state_sstream << nm << " ";
1706  }
1707  state_sstream << ")";
1708 
1709  TLorentzVector * pd = p->GetP4(); // incident particle 4p
1710 
1711  bool is_nuc = pdg::IsNeutronOrProton(p->Pdg());
1712  bool is_kaon = p->Pdg()==kPdgKP || p->Pdg()==kPdgKM;
1713  // not used // bool is_pion = p->Pdg()==kPdgPiP || p->Pdg()==kPdgPi0 || p->Pdg()==kPdgPiM;
1714  // update available energy -> init (mass + kinetic) + sum of f/s masses
1715  // for pion only. Probe mass not available for nucleon, kaon
1716  double availE = pd->Energy() + mass_sum;
1717  if(is_nuc||is_kaon) availE -= p->Mass();
1718  pd->SetE(availE);
1719 
1720  LOG("INukeUtils",pNOTICE)
1721  << "size, mass_sum, availE, pd mass, energy = " << pdgv.size() << " "
1722  << mass_sum << " " << availE << " " << p->Mass() << " " << p->Energy() ;
1723 
1724  // compute the 4p transfer to the hadronic blob
1725  double dE = mass_sum;
1726  if(is_nuc||is_kaon) dE -= p->Mass();
1727  TLorentzVector premnsub(0,0,0,dE);
1728  RemnP4 -= premnsub;
1729 
1730  LOG("INukeUtils", pINFO)
1731  << "Final state = " << state_sstream.str() << " has N = " << pdgv.size()
1732  << " particles / total mass = " << mass_sum;
1733  LOG("INukeUtils", pINFO)
1734  << "Composite system p4 = " << utils::print::P4AsString(pd);
1735 
1736  // Set the decay
1737  TGenPhaseSpace GenPhaseSpace;
1738  bool permitted = GenPhaseSpace.SetDecay(*pd, pdgv.size(), mass);
1739  if(!permitted) {
1740  LOG("INukeUtils", pERROR)
1741  << " *** Phase space decay is not permitted \n"
1742  << " Total particle mass = " << mass_sum << "\n"
1743  << " Decaying system p4 = " << utils::print::P4AsString(pd);
1744 
1745  // clean-up and return
1746  RemnP4 += premnsub;
1747  delete [] mass;
1748  delete pd;
1749  return false;
1750  }
1751 
1752  // The decay is permitted - add the incident particle at the event record
1753  // and mark is as 'Nucleon Cluster Target' (used to be confusing 'Decayed State')
1754  p->SetStatus(kIStNucleonClusterTarget); //kIStDecayedState);
1755  p->SetPdgCode(kPdgCompNuclCluster);
1756  ev->AddParticle(*p);
1757  // Get the maximum weight
1758  double wmax = -1;
1759  for(int k=0; k<200; k++) {
1760  double w = GenPhaseSpace.Generate();
1761  wmax = TMath::Max(wmax,w);
1762  }
1763  assert(wmax>0);
1764 
1765  LOG("INukeUtils", pINFO)
1766  << "Max phase space gen. weight @ current hadronic interaction: " << wmax;
1767 
1768  // Generate an unweighted decay
1769 
1770  RandomGen * rnd = RandomGen::Instance();
1771  wmax *= 1.2;
1772 
1773  bool accept_decay=false;
1774  unsigned int itry=0;
1775 
1776  while(!accept_decay)
1777  {
1778  itry++;
1779 
1780  if(itry>kMaxUnweightDecayIterations) {
1781  // report, clean-up and return
1782  LOG("INukeUtils", pNOTICE)
1783  << "Couldn't generate an unweighted phase space decay after "
1784  << itry << " attempts";
1785  delete [] mass;
1786  delete pd;
1787  return false;
1788  }
1789 
1790  double w = GenPhaseSpace.Generate();
1791  double gw = wmax * rnd->RndFsi().Rndm();
1792 
1793  if(w > wmax) {
1794  LOG("INukeUtils", pNOTICE)
1795  << "Decay weight = " << w << " > max decay weight = " << wmax;
1796  }
1797 
1798  LOG("INukeUtils", pNOTICE) << "Decay weight = " << w << " / R = " << gw;
1799  accept_decay = (gw<=w);
1800  }
1801 
1802  // Insert final state products into the event record
1803  // - the particles are added as daughters of the decayed state
1804  // - the particles are marked as final stable state (in hA mode)
1805  i=0;
1806  int mom = ev->ParticlePosition(p);
1807  LOG("INukeUtils", pNOTICE) << "mother index = " << mom;
1808  GHepStatus_t ist = kIStStableFinalState;
1809  GHepStatus_t ist_pi = kIStHadronInTheNucleus;
1810 
1811  TLorentzVector * v4 = p->GetX4();
1812 
1813  double checkpx = p->Px();
1814  double checkpy = p->Py();
1815  double checkpz = p->Pz();
1816  double checkE = p->E();
1817 
1818  //LOG("INukeUtils", PNOTICE)
1819 
1820  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1821 
1822  //-- current PDG code
1823  int pdgc = *pdg_iter;
1824  bool isnuc = pdg::IsNeutronOrProton(pdgc);
1825 
1826  //-- get the 4-momentum of the i-th final state particle
1827  TLorentzVector * p4fin = GenPhaseSpace.GetDecay(i++);
1828 
1829  //-- intranuke no longer throws "bindinos" but adds all the energy
1830  // not going at a simulated f/s particle at a "hadronic blob"
1831  // representing the remnant system: do the binding energy subtraction
1832  // here & update the remnant hadronic system 4p
1833  double M = PDGLibrary::Instance()->Find(pdgc)->Mass();
1834  double En = p4fin->Energy();
1835 
1836  double KE = En-M;
1837 
1838  //double KE = En;
1839  //if(is_pion) KE -= M;
1840 
1841  double dE_leftover = TMath::Min(NucRmvE, KE);
1842  KE -= dE_leftover;
1843  En = KE+M;
1844  double pmag_old = p4fin->P();
1845  double pmag_new = TMath::Sqrt(TMath::Max(0.,En*En-M*M));
1846  double scale = pmag_new / pmag_old;
1847  double pxn = scale * p4fin->Px();
1848  double pyn = scale * p4fin->Py();
1849  double pzn = scale * p4fin->Pz();
1850 
1851  TLorentzVector p4n(pxn,pyn,pzn,En);
1852  // LOG("INukeUtils", pNOTICE) << "Px = " << pxn << " Py = " << pyn
1853  // << " Pz = " << pzn << " E = " << KE;
1854  checkpx -= pxn;
1855  checkpy -= pyn;
1856  checkpz -= pzn;
1857  checkE -= KE;
1858 
1859  if (mode==kIMdHA &&
1860  (pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) )
1861  {
1862  if (p4n.Vect().Mag()>=0.001)
1863  {
1864  GHepParticle new_particle(pdgc, ist_pi, mom,-1,-1,-1, p4n, *v4);
1865  ev->AddParticle(new_particle);
1866  }
1867  else
1868  {
1869  // Momentum too small, assign a non-zero momentum to the particle
1870  // Conserve momentum with the remnant nucleus
1871 
1872  LOG("INukeUtils", pINFO)<<"Momentum too small; assigning 0.001 as new momentum";
1873 
1874  double phi = 2*kPi*rnd->RndFsi().Rndm();
1875  double omega = 2*rnd->RndFsi().Rndm();
1876  // throw number against solid angle for uniform distribution
1877 
1878  double E4n = TMath::Sqrt(0.001*0.001+M*M);
1879  p4n.SetPxPyPzE(0.001,0,0,E4n);
1880  p4n.Rotate(TMath::ACos(1-omega),TVector3(0,0,1));
1881  p4n.Rotate(phi,TVector3(1,0,0));
1882 
1883  RemnP4 -= (p4n - TLorentzVector(0,0,0,M));
1884 
1885  GHepParticle new_particle(pdgc, ist, mom,-1,-1,-1, p4n, *v4);
1886  ev->AddParticle(new_particle);
1887  }
1888  }
1889  else
1890  {
1891  GHepParticle new_particle(pdgc, ist, mom,-1,-1,-1, p4n, *v4);
1892 
1893  if(isnuc) new_particle.SetRemovalEnergy(0.);
1894  ev->AddParticle(new_particle);
1895  }
1896 
1897  double dpx = (1-scale)*p4fin->Px();
1898  double dpy = (1-scale)*p4fin->Py();
1899  double dpz = (1-scale)*p4fin->Pz();
1900  TLorentzVector premnadd(dpx,dpy,dpz,dE_leftover);
1901  RemnP4 += premnadd;
1902  }
1903  //LOG("INukeUtils", pNOTICE) << "TEST: " << p->Mass();
1904  LOG("INukeUtils", pNOTICE) << "check conservation: Px = " << checkpx << " Py = " << checkpy
1905  << " Pz = " << checkpz << " E = " << checkE;
1906 
1907  // Clean-up
1908  delete [] mass;
1909  delete pd;
1910  delete v4;
1911 
1912  return true;
1913 }
1914 
1915 double genie::utils::intranuke2018::sigmaTotalOset (
1916  const double &pionKineticEnergy,
1917  const double &density,
1918  const int &pionPDG,
1919  const double &protonFraction,
1920  const bool &isTableChosen
1921  )
1922 {
1923  // ------ OsetCrossSection init (only first time function is called) ------ //
1924  static INukeOset *iNukeOset = NULL;
1925 
1926  if (iNukeOset == NULL)
1927  {
1928  if (isTableChosen)
1929  {
1930  // set directory with data on first call
1931  static const std::string dataDir = (gSystem->Getenv("GINUKEHADRONDATA")) ?
1932  string(gSystem->Getenv("GINUKEHADRONDATA")) :
1933  string(gSystem->Getenv("GENIE")) +
1934  string("/data/evgen/intranuke/");
1935  // set file with Oset table on first call
1936  static const std::string dataFile = dataDir + "tot_xsec/"
1937  "intranuke-xsection-pi+N-Oset.dat";
1938  // initialize OsetCrossSection on first call
1939  iNukeOset = new INukeOsetTable (dataFile.c_str());
1940  }
1941  else iNukeOset = new INukeOsetFormula();
1942  }
1943  // ------ OsetCrossSection init (only first time function is called) ------ //
1944 
1945  // set up Oset class (assign pion Tk, nuclear density etc)
1946  iNukeOset->setupOset (density, pionKineticEnergy, pionPDG, protonFraction);
1947 
1948  return iNukeOset->getTotalCrossSection();
1949 
1950 }
genie::utils::intranuke2018::MeanFreePath_Delta
double MeanFreePath_Delta(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A)
Mean free path (Delta++ test)
Definition: INukeUtils2018.cxx:236
genie::Intranuke2018::GetNR
double GetNR() const
Definition: Intranuke2018.h:83
genie::kPdgP33m1232_DeltaPP
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::GHepParticle::SetFirstMother
void SetFirstMother(int m)
Definition: GHepParticle.h:132
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
genie::utils::intranuke2018::PreEquilibrium
void PreEquilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Definition: INukeUtils2018.cxx:491
genie::GHepRecord::ParticlePosition
virtual int ParticlePosition(int pdg, GHepStatus_t i, int start=0) const
Definition: GHepRecord.cxx:137
genie::GHepParticle::GetX4
TLorentzVector * GetX4(void) const
Definition: GHepParticle.cxx:212
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
INukeNucleonCorr::getInstance
static INukeNucleonCorr * getInstance()
get single instance of INukeNucleonCorr; create if necessary
Definition: INukeNucleonCorr.h:23
pERROR
#define pERROR
Definition: Messenger.h:59
genie::GHepParticle::E
double E(void) const
Get energy.
Definition: GHepParticle.h:91
genie::Intranuke2018::GetUseOset
bool GetUseOset() const
Definition: Intranuke2018.h:91
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::utils::intranuke2018::PhaseSpaceDecay
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
Definition: INukeUtils2018.cxx:1678
genie::utils::nuclear::Density
double Density(double r, int A, double ring=0.)
Definition: NuclearUtils.cxx:396
genie::Intranuke2018::GetINukeMode
virtual string GetINukeMode() const
Definition: Intranuke2018.h:72
INukeOsetTable
Table-based implementation of Oset model.
Definition: INukeOsetTable.h:25
genie::utils::intranuke2018::StepParticle
void StepParticle(GHepParticle *p, double step, double nuclear_radius=-1.)
Step particle.
Definition: INukeUtils2018.cxx:435
genie::utils::intranuke2018::ThreeBodyKinematics
bool ThreeBodyKinematics(GHepRecord *ev, GHepParticle *p, int tcode, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI *Nuclmodel=(const NuclearModelI *) 0)
Definition: INukeUtils2018.cxx:1072
genie::utils::intranuke2018::PionProduction
bool PionProduction(GHepRecord *ev, GHepParticle *p, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI *Nuclmodel)
Definition: INukeUtils2018.cxx:1301
genie::utils::print::P4AsString
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
genie::controls::kMaxUnweightDecayIterations
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
genie::NuclearModelI
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:436
genie::EINukeMode
EINukeMode
Definition: INukeMode.h:29
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:236
genie::Intranuke2018::GetHadStep
double GetHadStep() const
Definition: Intranuke2018.h:89
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::Intranuke2018::GetRemnA
double GetRemnA() const
Definition: Intranuke2018.h:79
genie::GHepParticle::Pz
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
genie::NuclearModelI::Momentum3
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
genie::GHepParticle::Energy
double Energy(void) const
Get energy.
Definition: GHepParticle.h:92
genie::Intranuke2018::GetDelRPion
double GetDelRPion() const
Definition: Intranuke2018.h:85
genie::PDGCodeList
A list of PDG codes.
Definition: PDGCodeList.h:32
genie::GHepParticle::Px
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
genie::GHepParticle::CompareMomentum
bool CompareMomentum(const GHepParticle *p) const
Definition: GHepParticle.cxx:502
epsilon
const double epsilon
Definition: gCalcHEDISDiffXsec.cxx:76
genie::GHepParticle::LastMother
int LastMother(void) const
Definition: GHepParticle.h:67
genie::utils::intranuke2018::Dist2Exit
double Dist2Exit(const TLorentzVector &x4, const TLorentzVector &p4, double A, double NR=3, double R0=1.4)
Distance to exit.
Definition: INukeUtils2018.cxx:373
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
genie::kPdgCompNuclCluster
const int kPdgCompNuclCluster
Definition: PDGCodes.h:217
genie::GHepParticle::SetPosition
void SetPosition(const TLorentzVector &v4)
Definition: GHepParticle.cxx:252
INukeOset::getTotalCrossSection
double getTotalCrossSection() const
return total = (qel+cex+abs) cross section
Definition: INukeOset.h:35
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::units::mb
static constexpr double mb
Definition: Units.h:79
genie::Intranuke2018::GetR0
double GetR0() const
Definition: Intranuke2018.h:82
genie::kPdgGamma
const int kPdgGamma
Definition: PDGCodes.h:189
genie::GHepParticle::GetP4
TLorentzVector * GetP4(void) const
Definition: GHepParticle.cxx:194
genie::GHepParticle::ComparePdgCodes
bool ComparePdgCodes(const GHepParticle *p) const
Definition: GHepParticle.cxx:481
genie::units::fm2
static constexpr double fm2
Definition: Units.h:76
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::utils::intranuke::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils.cxx:754
genie::utils::intranuke2018::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils2018.cxx:873
genie::utils::intranuke2018::TwoBodyCollision
bool TwoBodyCollision(GHepRecord *ev, int pcode, int tcode, int scode, int s2code, double C3CM, GHepParticle *p, GHepParticle *t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode=kIMdHA)
Intranuke utility functions.
Definition: INukeUtils2018.cxx:729
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
genie::GHepRecord::GetStableDescendants
virtual vector< int > * GetStableDescendants(int position) const
Definition: GHepRecord.cxx:174
genie::utils::intranuke2018::CalculateEta
double CalculateEta(double Minc, double ke, double Mtarg, double Mtwopart, double Mpi)
Definition: INukeUtils2018.cxx:1652
pINFO
#define pINFO
Definition: Messenger.h:62
genie::utils::intranuke2018::MeanFreePath
double MeanFreePath(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double nRpi=0.5, double nRnuc=1.0, const bool useOset=false, const bool altOset=false, const bool xsecNNCorr=false, string INukeMode="XX2018")
Mean free path (pions, nucleons)
Definition: INukeUtils2018.cxx:79
genie::utils::intranuke2018::Equilibrium
void Equilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Definition: INukeUtils2018.cxx:642
genie::Intranuke2018::GetAltOset
bool GetAltOset() const
Definition: Intranuke2018.h:92
pWARN
#define pWARN
Definition: Messenger.h:60
INukeOset::setupOset
virtual void setupOset(const double &density, const double &pionTk, const int &pionPDG, const double &protonFraction)=0
use to set up Oset class (assign pion Tk, nuclear density etc)
genie::Intranuke2018::GetRemnZ
double GetRemnZ() const
Definition: Intranuke2018.h:80
genie::GHepParticle::CompareStatusCodes
bool CompareStatusCodes(const GHepParticle *p) const
Definition: GHepParticle.cxx:486
genie::GHepParticle::SetRemovalEnergy
void SetRemovalEnergy(double Erm)
Definition: GHepParticle.cxx:380
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:159
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::INukeHadroData2018::Instance
static INukeHadroData2018 * Instance(void)
Definition: INukeHadroData2018.cxx:175
genie::GHepParticle::SetLastMother
void SetLastMother(int m)
Definition: GHepParticle.h:133
genie::utils::intranuke2018::ProbSurvival
double ProbSurvival(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double mfp_scale_factor, const Intranuke2018 &fsi_model)
Hadron survival probability.
Definition: INukeUtils2018.cxx:284
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::GHepParticle::X4
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
genie::Target::SetHitNucPdg
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
genie::pdg::IsNeutronOrProton
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:351
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::utils::print::X4AsString
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:57
genie::NuclearModelI::GenerateNucleon
virtual bool GenerateNucleon(const Target &) const =0
genie::utils::intranuke2018::sigmaTotalOset
double sigmaTotalOset(const double &pionKineticEnergy, const double &density, const int &pionPDG, const double &protonFraction, const bool &isTableChosen=true)
Definition: INukeUtils2018.cxx:1915
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::GHepRecord
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
genie::Intranuke2018::GetDelRNucleon
double GetDelRNucleon() const
Definition: Intranuke2018.h:86
genie::units::m
static constexpr double m
Definition: Units.h:71
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::utils::intranuke2018::Dist2ExitMFP
double Dist2ExitMFP(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double NR=3, double R0=1.4)
Distance to exit.
Definition: INukeUtils2018.cxx:400
genie::Intranuke2018::GetXsecNNCorr
bool GetXsecNNCorr() const
Definition: Intranuke2018.h:93
genie::utils::print::Vec3AsString
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::GHepParticle::SetPdgCode
void SetPdgCode(int c)
Definition: GHepParticle.cxx:230
genie::GHepParticle::Py
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
INukeOsetFormula
Formula-based implementation of Oset model.
Definition: INukeOsetFormula.h:17
Generated on Wed Jul 17 2024 14:46:42 for GENIEGenerator by   doxygen 1.8.5