GENIEGenerator
Loading...
Searching...
No Matches
INukeUtils2018.cxx
Go to the documentation of this file.
1//____________________________________________________________________________
2/*
3 Copyright (c) 2003-2025, 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
72using std::ostringstream;
73using namespace genie;
74using namespace genie::utils;
75using namespace genie::constants;
76using namespace genie::controls;
77
78//____________________________________________________________________________
79double 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//____________________________________________________________________________
236double 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//____________________________________________________________________________
284double 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//____________________________________________________________________________
373double 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//____________________________________________________________________________
400double 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//____________________________________________________________________________
435void 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
491void 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
642void 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//___________________________________________________________________________
729bool 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//___________________________________________________________________________
873bool 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//___________________________________________________________________________
1072bool 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//___________________________________________________________________________
1301bool 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//___________________________________________________________________________
1652double 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//___________________________________________________________________________
1678bool 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
1915double 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}
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::INukeHadroData2018::Instance
static INukeHadroData2018 * Instance(void)
Definition INukeHadroData2018.cxx:175
genie::Intranuke2018::GetXsecNNCorr
bool GetXsecNNCorr() const
Definition Intranuke2018.h:93
genie::Intranuke2018::GetRemnZ
double GetRemnZ() const
Definition Intranuke2018.h:80
genie::Intranuke2018::GetHadStep
double GetHadStep() const
Definition Intranuke2018.h:89
genie::Intranuke2018::GetR0
double GetR0() const
Definition Intranuke2018.h:82
genie::Intranuke2018::GetAltOset
bool GetAltOset() const
Definition Intranuke2018.h:92
genie::Intranuke2018::GetDelRPion
double GetDelRPion() const
Definition Intranuke2018.h:85
genie::Intranuke2018::GetDelRNucleon
double GetDelRNucleon() const
Definition Intranuke2018.h:86
genie::Intranuke2018::GetUseOset
bool GetUseOset() const
Definition Intranuke2018.h:91
genie::Intranuke2018::GetRemnA
double GetRemnA() const
Definition Intranuke2018.h:79
genie::Intranuke2018::GetINukeMode
virtual string GetINukeMode() const
Definition Intranuke2018.h:72
genie::Intranuke2018::GetNR
double GetNR() const
Definition Intranuke2018.h:83
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::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::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::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::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::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::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::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::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::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::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::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::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::intranuke2018::StepParticle
void StepParticle(GHepParticle *p, double step, double nuclear_radius=-1.)
Step particle.
Definition INukeUtils2018.cxx:435
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::utils::intranuke2018::CalculateEta
double CalculateEta(double Minc, double ke, double Mtarg, double Mtwopart, double Mpi)
Definition INukeUtils2018.cxx:1652
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