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