HNIntranuke2018.cxx
Go to the documentation of this file.
1 
2 //____________________________________________________________________________
3 /*
4  Copyright (c) 2003-2019, The GENIE Collaboration
5  For the full text of the license visit http://copyright.genie-mc.org
6  or see $GENIE/LICENSE
7 
8  Author: Steve Dytman <dytman+@pitt.edu>, Pittsburgh Univ.
9  Aaron Meyer <asm58@pitt.edu>, Pittsburgh Univ.
10  Alex Bell, Pittsburgh Univ.
11  Hugh Gallagher <gallag@minos.phy.tufts.edu>, Tufts Univ.
12  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>, Rutherford Lab.
13  September 20, 2005
14 
15  For the class documentation see the corresponding header file.
16 
17  Important revisions after version 2.0.0 :
18  @ Nov 30, 2007 - SD
19  Changed the hadron tracking algorithm to take into account the radial
20  nuclear density dependence. Using the somewhat empirical approach of
21  increasing the nuclear radius by a const (tunable) number times the tracked
22  particle's de Broglie wavelength as this helps getting the hadron+nucleus
23  cross sections right.
24  @ Mar 08, 2008 - CA
25  Fixed code retrieving the remnant nucleus which stopped working as soon as
26  simulation of nuclear de-excitation started pushing photons in the target
27  nucleus daughter list.
28  @ Jun 20, 2008 - CA
29  Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
30  deleted after it was added at the GHEP event record.
31  @ Jul 15, 2010 - AM
32  The hN mode is now implemented in Intranuke. Similar to hA mode, but particles
33  produced by reactions are stepped through the nucleus like probe particles.
34  Particles react with nucleons instead of the entire nucleus, and final states
35  are determined after reactions are finished, not before.
36  @ Dec 15, 2014 - SD, Nick Geary
37  Update fates to include Compound Nucleus final state correctly.
38  @ Jan 9, 2015 - SD, NG, Tomek Golan
39  Added 2014 version of INTRANUKE codes (new class) for independent development.
40  @ Oct, 2015 - TG
41  Added 2015 version of INTRANUKE codes (new class) for independent development. Include Oset model for medium corrections to piA for Tpi<350 MeV.
42  @ May, 2015 Flor Blaszczyk
43  K+ are now handled.
44  @ July, 2016 Nicholas Suarez, Josh Kleckner, SD
45  fix memory leak, fix fates, improve NNCorr binning
46  & Mar, 2018 Nicholas Suarez, SD
47  add compound nucleus option to populate KE<30 MeV
48 */
49 //____________________________________________________________________________
50 
51 #include <cstdlib>
52 #include <sstream>
53 
54 #include <TMath.h>
55 
82 
83 using std::ostringstream;
84 
85 using namespace genie;
86 using namespace genie::utils;
87 using namespace genie::utils::intranuke2018;
88 using namespace genie::constants;
89 using namespace genie::controls;
90 
91 //___________________________________________________________________________
92 //___________________________________________________________________________
93 // Methods specific to INTRANUKE's HN-mode
94 //___________________________________________________________________________
95 //___________________________________________________________________________
97 Intranuke2018("genie::HNIntranuke2018")
98 {
99 
100 }
101 //___________________________________________________________________________
103 Intranuke2018("genie::HNIntranuke2018",config)
104 {
105 
106 }
107 //___________________________________________________________________________
109 {
110 
111 }
112 //___________________________________________________________________________
114 {
115  LOG("HNIntranuke2018", pNOTICE)
116  << "************ Running hN2018 MODE INTRANUKE ************";
117 
118  /* LOG("HNIntranuke2018", pWARN)
119  << print::PrintFramedMesg(
120  "Experimental code (INTRANUKE/hN model) - Run at your own risk");
121  */
122 
124 
125  LOG("HNIntranuke2018", pINFO) << "Done with this event";
126 }
127 //___________________________________________________________________________
129 {
130 // Simulate a hadron interaction for the input particle p in HN mode
131 //
132  if(!p || !ev)
133  {
134  LOG("HNIntranuke2018", pERROR) << "** Null input!";
135  return;
136  }
137 
138  // check particle id
139  int pdgc = p->Pdg();
140  bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
141  bool is_kaon = (pdgc==kPdgKP);
142  bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
143  bool is_gamma = (pdgc==kPdgGamma);
144  if(!(is_pion || is_baryon || is_gamma || is_kaon))
145  {
146  LOG("HNIntranuke2018", pERROR) << "** Cannot handle particle: " << p->Name();
147  return;
148  }
149  try
150  {
151  // select a fate for the input particle
152  INukeFateHN_t fate = this->HadronFateHN(p);
153 
154  // store the fate
155  ev->Particle(p->FirstMother())->SetRescatterCode((int)fate);
156 
157  if(fate == kIHNFtUndefined)
158  {
159  LOG("HNIntranuke2018", pERROR) << "** Couldn't select a fate";
160  LOG("HNIntranuke2018", pERROR) << "** Num Protons: " << fRemnZ
161  << ", Num Neutrons: "<<(fRemnA-fRemnZ);
162  LOG("HNIntranuke2018", pERROR) << "** Particle: " << "\n" << (*p);
163  //LOG("HNIntranuke2018", pERROR) << "** Event Record: " << "\n" << (*ev);
164  //p->SetStatus(kIStUndefined);
165  p->SetStatus(kIStStableFinalState);
166  ev->AddParticle(*p);
167  return;
168  }
169 
170  LOG("HNIntranuke2018", pNOTICE)
171  << "Selected " << p->Name() << " fate: " << INukeHadroFates::AsString(fate);
172 
173  // handle the reaction
174  if(fate == kIHNFtCEx || fate == kIHNFtElas)
175  {
176  this->ElasHN(ev,p,fate);
177  }
178  else if(fate == kIHNFtAbs) {this-> AbsorbHN(ev,p,fate);}
179  else if(fate == kIHNFtInelas && pdgc != kPdgGamma)
180  {
181 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
182  LOG("HNIntranuke2018", pDEBUG)
183  << "Invoking InelasticHN() for a : " << p->Name()
184  << " whose fate is : " << INukeHadroFates::AsString(fate);
185 #endif
186 
187  this-> InelasticHN(ev,p);
188  }
189  else if(fate == kIHNFtInelas && pdgc == kPdgGamma) {this-> GammaInelasticHN(ev,p,fate);}
191  else if(fate == kIHNFtNoInteraction)
192  {
194  ev->AddParticle(*p);
195  return;
196  }
197  }
199  {
200  this->SimulateHadronicFinalState(ev,p);
201  LOG("HNIntranuke2018", pNOTICE)
202  << "retry call to SimulateHadronicFinalState ";
203  LOG("HNIntranuke2018", pNOTICE) << exception;
204 
205  }
206 }
207 //___________________________________________________________________________
209 {
210 // Select a hadron fate in HN mode
211 //
213 
214  // get pdgc code & kinetic energy in MeV
215  int pdgc = p->Pdg();
216  double ke = p->KinE() / units::MeV;
217 
218  bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
219 
220  if (isPion and fUseOset and ke < 350.0) return HadronFateOset ();
221 
222  LOG("HNIntranuke2018", pNOTICE)
223  << "Selecting hN fate for " << p->Name() << " with KE = " << ke << " MeV";
224 
225  // try to generate a hadron fate
226  unsigned int iter = 0;
227  while(iter++ < kRjMaxIterations) {
228 
229  // handle pions
230  //
231  if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
232 
233  double frac_cex = this->FateWeight(pdgc, kIHNFtCEx)
234  * fHadroData2018->Frac(pdgc, kIHNFtCEx, ke, fRemnA, fRemnZ);
235  double frac_elas = this->FateWeight(pdgc, kIHNFtElas)
236  * fHadroData2018->Frac(pdgc, kIHNFtElas, ke, fRemnA, fRemnZ);
237  double frac_inel = this->FateWeight(pdgc, kIHNFtInelas)
238  * fHadroData2018->Frac(pdgc, kIHNFtInelas, ke, fRemnA, fRemnZ);
239  double frac_abs = this->FateWeight(pdgc, kIHNFtAbs)
240  * fHadroData2018->Frac(pdgc, kIHNFtAbs, ke, fRemnA, fRemnZ);
241 
242  frac_cex *= fNucCEXFac; // scaling factors
243  frac_abs *= fNucAbsFac;
244  frac_elas *= fNucQEFac;
245  if(pdgc==kPdgPi0) frac_abs*= 0.665; //isospin factor
246 
247  LOG("HNIntranuke2018", pNOTICE)
248  << "\n frac{" << INukeHadroFates::AsString(kIHNFtCEx) << "} = " << frac_cex
249  << "\n frac{" << INukeHadroFates::AsString(kIHNFtElas) << "} = " << frac_elas
250  << "\n frac{" << INukeHadroFates::AsString(kIHNFtInelas) << "} = " << frac_inel
251  << "\n frac{" << INukeHadroFates::AsString(kIHNFtAbs) << "} = " << frac_abs;
252 
253  // compute total fraction (can be <1 if fates have been switched off)
254  double tf = frac_cex +
255  frac_elas +
256  frac_inel +
257  frac_abs;
258 
259  double r = tf * rnd->RndFsi().Rndm();
260 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
261  LOG("HNIntranuke2018", pDEBUG) << "r = " << r << " (max = " << tf << ")";
262 #endif
263 
264  double cf=0; // current fraction
265 
266  if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
267  if(r < (cf += frac_elas )) return kIHNFtElas; //elas
268  if(r < (cf += frac_inel )) return kIHNFtInelas; //inelas
269  if(r < (cf += frac_abs )) return kIHNFtAbs; //abs
270 
271  LOG("HNIntranuke2018", pWARN)
272  << "No selection after going through all fates! "
273  << "Total fraction = " << tf << " (r = " << r << ")";
274  ////////////////////////////
275  return kIHNFtUndefined;
276  }
277 
278  // handle nucleons
279  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
280 
281  double frac_elas = this->FateWeight(pdgc, kIHNFtElas)
282  * fHadroData2018->Frac(pdgc, kIHNFtElas, ke, fRemnA, fRemnZ);
283  double frac_inel = this->FateWeight(pdgc, kIHNFtInelas)
284  * fHadroData2018->Frac(pdgc, kIHNFtInelas, ke, fRemnA, fRemnZ);
285  double frac_cmp = this->FateWeight(pdgc, kIHNFtCmp)
286  * fHadroData2018->Frac(pdgc, kIHNFtCmp, ke, fRemnA , fRemnZ);
287 
288  LOG("HNIntranuke2018", pINFO)
289  << "\n frac{" << INukeHadroFates::AsString(kIHNFtElas) << "} = " << frac_elas
290  << "\n frac{" << INukeHadroFates::AsString(kIHNFtInelas) << "} = " << frac_inel;
291 
292  // compute total fraction (can be <1 if fates have been switched off)
293  double tf = frac_elas +
294  frac_inel +
295  frac_cmp;
296 
297  double r = tf * rnd->RndFsi().Rndm();
298 
299 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
300  LOG("HNIntranuke2018", pDEBUG) << "r = " << r << " (max = " << tf << ")";
301 #endif
302 
303  double cf=0; // current fraction
304  if(r < (cf += frac_elas )) return kIHNFtElas; // elas
305  if(r < (cf += frac_inel )) return kIHNFtInelas; // inelas
306  if(r < (cf += frac_cmp )) return kIHNFtCmp; // cmp
307 
308  LOG("HNIntranuke2018", pWARN)
309  << "No selection after going through all fates! "
310  << "Total fraction = " << tf << " (r = " << r << ")";
311  //////////////////////////
312  return kIHNFtUndefined;
313  }
314 
315  // handle gamma -- does not currently consider the elastic case
316  else if (pdgc==kPdgGamma) return kIHNFtInelas;
317  // Handle kaon -- elastic + charge exchange
318  else if (pdgc==kPdgKP){
319  double frac_cex = this->FateWeight(pdgc, kIHNFtCEx)
320  * fHadroData2018->Frac(pdgc, kIHNFtCEx, ke, fRemnA, fRemnZ);
321  double frac_elas = this->FateWeight(pdgc, kIHNFtElas)
322  * fHadroData2018->Frac(pdgc, kIHNFtElas, ke, fRemnA, fRemnZ);
323 
324  // frac_cex *= fNucCEXFac; // scaling factors
325  // frac_elas *= fNucQEFac; // Flor - Correct scaling factors?
326 
327  LOG("HNIntranuke", pINFO)
328  << "\n frac{" << INukeHadroFates::AsString(kIHNFtCEx) << "} = " << frac_cex
329  << "\n frac{" << INukeHadroFates::AsString(kIHNFtElas) << "} = " << frac_elas;
330 
331  // compute total fraction (can be <1 if fates have been switched off)
332  double tf = frac_cex +
333  frac_elas;
334 
335  double r = tf * rnd->RndFsi().Rndm();
336 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
337  LOG("HNIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
338 #endif
339 
340  double cf=0; // current fraction
341 
342  if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
343  if(r < (cf += frac_elas )) return kIHNFtElas; //elas
344 
345  LOG("HNIntranuke", pWARN)
346  << "No selection after going through all fates! "
347  << "Total fraction = " << tf << " (r = " << r << ")";
348  ////////////////////////////
349  return kIHNFtUndefined;
350  }//End K+
351 
352  /*else if (pdgc==kPdgKM){
353 
354  return kIHNFtElas;
355  }//End K-*/
356 
357  }//iterations
358 
359  return kIHNFtUndefined;
360 }
361 
362 //___________________________________________________________________________
363 double HNIntranuke2018::FateWeight(int pdgc, INukeFateHN_t fate) const
364 {
365  // turn fates off if the remnant nucleus does not have the number of p,n
366  // required
367 
368  int np = fRemnZ;
369  int nn = fRemnA - fRemnZ;
370 
371  if (np < 1 && nn < 1)
372  {
373  LOG("HNIntranuke2018", pERROR) << "** Nothing left in nucleus!!! **";
374  return 0;
375  }
376  else
377  {
378  if (fate == kIHNFtCEx && pdgc==kPdgPiP ) { return (nn>=1) ? 1. : 0.; }
379  if (fate == kIHNFtCEx && pdgc==kPdgPiM ) { return (np>=1) ? 1. : 0.; }
380  if (fate == kIHNFtCEx && pdgc==kPdgKP ) { return (nn>=1) ? 1. : 0.; } //Added, changed np to nn
381  if (fate == kIHNFtAbs) { return ((nn>=1) && (np>=1)) ? 1. : 0.; }
382  if (fate == kIHNFtCmp ) { return ((pdgc==kPdgProton||pdgc==kPdgNeutron)&&fDoCompoundNucleus&&fRemnA>5) ? 1. : 0.; }
383 
384  }
385  return 1.;
386 }
387 //___________________________________________________________________________
389  GHepRecord * ev, GHepParticle * p, INukeFateHN_t fate) const
390 {
391  // handles pi+d->2p, pi-d->nn, pi0 d->pn absorbtion, all using pi+d values
392 
393  int pdgc = p->Pdg();
394 
395 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
396  LOG("HNIntranuke2018", pDEBUG)
397  << "AbsorbHN() is invoked for a : " << p->Name()
398  << " whose fate is : " << INukeHadroFates::AsString(fate);
399 #endif
400 
401  // check fate
402  if(fate!=kIHNFtAbs)
403  {
404  LOG("HNIntranuke2018", pWARN)
405  << "AbsorbHN() cannot handle fate: " << INukeHadroFates::AsString(fate);
406  return;
407  }
408 
409  // random number generator
411 
412  // Notes on the kinematics
413  // -- Simple variables are used for efficiency
414  // -- Variables are numbered according to particle
415  // -- -- #1 -> incoming particle
416  // -- -- #2 -> target (here, 2_1 and 2_2 for individual particles)
417  // -- -- #3 -> scattered incoming (Particle tracked in hA mode)
418  // -- -- #4 -> other scattered particle
419  // -- Suffix "L" is for lab frame, suffix "CM" is for center of mass frame
420  // -- Subscript "z" is for parallel component, "t" is for transverse
421 
422  int pcode, t1code, t2code, scode, s2code; // particles
423  double M1, M2_1, M2_2, M3, M4; // rest energies, in GeV
424  double E1L, P1L, E2L, P2L, E3L, P3L, E4L, P4L;
425  double P1zL, P2zL;
426  double beta, gm; // speed and gamma for CM frame in lab
427  double Et, E2CM;
428  double C3CM, S3CM; // cos and sin of scattering angle
429  double Theta1, Theta2, theta5;
430  double PHI3; // transverse scattering angle
431  double E1CM, E3CM, E4CM, P3CM;
432  double P3zL, P3tL, P4zL, P4tL;
433  double E2_1L, E2_2L;
434  TVector3 tP2_1L, tP2_2L, tP1L, tP2L, tPtot, P1zCM, P2zCM;
435  TVector3 tP3L, tP4L;
436  TVector3 bDir, tTrans, tbeta, tVect;
437 
438  // Library instance for reference
439  PDGLibrary * pLib = PDGLibrary::Instance();
440 
441  // Handle fermi target
442  Target target(ev->TargetNucleus()->Pdg());
443 
444  // Target should be a deuteron, but for now
445  // handling it as seperate nucleons
446  if(pdgc==211) // pi-plus
447  {
448  pcode = 211;
449  t1code = 2212; // proton
450  t2code = 2112; // neutron
451  scode = 2212;
452  s2code = 2212;
453  }
454  else if(pdgc==-211) // pi-minus
455  {
456  pcode = -211;
457  t1code = 2212;
458  t2code = 2112;
459  scode = 2112;
460  s2code = 2112;
461  }
462  else if(pdgc==111) // pi-zero
463  {
464  pcode = 111;
465  t1code = 2212;
466  t2code = 2112;
467  scode = 2212;
468  s2code = 2112;
469  }
470  else
471  {
472  LOG("HNIntranuke2018", pWARN)
473  << "AbsorbHN() cannot handle probe: " << pdgc;
474  return;
475  }
476 
477  // assign proper masses
478  M1 = pLib->Find(pcode) ->Mass();
479  M2_1 = pLib->Find(t1code)->Mass();
480  M2_2 = pLib->Find(t2code)->Mass();
481  M3 = pLib->Find(scode) ->Mass();
482  M4 = pLib->Find(s2code)->Mass();
483 
484  // handle fermi momentum
485  if(fDoFermi)
486  {
487  target.SetHitNucPdg(t1code);
489  tP2_1L=fFermiFac * fNuclmodel->Momentum3();
490  E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
491 
492  target.SetHitNucPdg(t2code);
494  tP2_2L=fFermiFac * fNuclmodel->Momentum3();
495  E2_2L = TMath::Sqrt(tP2_2L.Mag2() + M2_2*M2_2);
496  }
497  else
498  {
499  tP2_1L.SetXYZ(0.0, 0.0, 0.0);
500  E2_1L = M2_1;
501 
502  tP2_2L.SetXYZ(0.0, 0.0, 0.0);
503  E2_2L = M2_2;
504  }
505 
506  E2L = E2_1L + E2_2L;
507 
508  // adjust p to reflect scattering
509  // get random scattering angle
510  C3CM = fHadroData2018->IntBounce(p,t1code,scode,fate);
511  if (C3CM<-1.)
512  {
514  ev->AddParticle(*p);
515  return;
516  }
517  S3CM = TMath::Sqrt(1.0 - C3CM*C3CM);
518 
519  // Get lab energy and momenta
520  E1L = p->E();
521  if(E1L<0.001) E1L=0.001;
522  P1L = TMath::Sqrt(E1L*E1L - M1*M1);
523  tP1L = p->P4()->Vect();
524  tP2L = tP2_1L + tP2_2L;
525  P2L = tP2L.Mag();
526  tPtot = tP1L + tP2L;
527 
528  // get unit vectors and angles needed for later
529  bDir = tPtot.Unit();
530  Theta1 = tP1L.Angle(bDir);
531  Theta2 = tP2L.Angle(bDir);
532 
533  // get parallel and transverse components
534  P1zL = P1L*TMath::Cos(Theta1);
535  P2zL = P2L*TMath::Cos(Theta2);
536  tVect.SetXYZ(1,0,0);
537  if(TMath::Abs((tVect - bDir).Mag())<.01) tVect.SetXYZ(0,1,0);
538  theta5 = tVect.Angle(bDir);
539  tTrans = (tVect - TMath::Cos(theta5)*bDir).Unit();
540 
541  // calculate beta and gamma
542  tbeta = tPtot * (1.0 / (E1L + E2L));
543  beta = tbeta.Mag();
544  gm = 1.0 / TMath::Sqrt(1.0 - beta*beta);
545 
546  // boost to CM frame to get scattered particle momenta
547  E1CM = gm*E1L - gm*beta*P1zL;
548  P1zCM = gm*P1zL*bDir - gm*tbeta*E1L;
549  E2CM = gm*E2L - gm*beta*P2zL;
550  P2zCM = gm*P2zL*bDir - gm*tbeta*E2L;
551  Et = E1CM + E2CM;
552  E3CM = (Et*Et + (M3*M3) - (M4*M4)) / (2.0*Et);
553  E4CM = Et - E3CM;
554  P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
555 
556  // boost back to lab
557  P3zL = gm*beta*E3CM + gm*P3CM*C3CM;
558  P3tL = P3CM*S3CM;
559  P4zL = gm*beta*E4CM + gm*P3CM*(-C3CM);
560  P4tL = P3CM*(-S3CM);
561 
562  P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
563  P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
564 
565  // check for too small values
566  // may introduce error, so warn if it occurs
567  if(!(TMath::Finite(P3L))||P3L<.001)
568  {
569  LOG("HNIntranuke2018",pINFO)
570  << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
571  << "\n" << "--> Assigning .001 as new momentum";
572  P3tL = 0;
573  P3zL = .001;
574  P3L = .001;
575  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
576  }
577 
578  if(!(TMath::Finite(P4L))||P4L<.001)
579  {
580  LOG("HNIntranuke2018",pINFO)
581  << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
582  << "\n" << "--> Assigning .001 as new momentum";
583  P4tL = 0;
584  P4zL = .001;
585  P4L = .001;
586  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
587  }
588 
589  // pauli blocking (do not apply PB for Oset)
590  //if(!fUseOset && (P3L < fFermiMomentum || P4L < fFermiMomentum))
591  double ke = p->KinE() / units::MeV;
592  if((!fUseOset || ke > 350.0 ) && (P3L < fFermiMomentum || P4L < fFermiMomentum))
593  {
594 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
595  LOG("HNIntranuke2018",pINFO) << "AbsorbHN failed: Pauli blocking";
596 #endif
597  /*
598  p->SetStatus(kIStHadronInTheNucleus);
599  //disable until needed
600  // utils::intranuke2018::StepParticle(p,fFreeStep,fTrackingRadius);
601  ev->AddParticle(*p);
602  return;
603  */
604  // new attempt at error handling:
605  LOG("HNIntranuke2018", pINFO) << "AbsorbHN failed: Pauli blocking";
607  exception.SetReason("hN absorption failed");
608  throw exception;
609  }
610 
611  // handle remnant nucleus updates
612  fRemnZ--;
613  fRemnA -=2;
614  fRemnP4 -= TLorentzVector(tP2_1L,E2_1L);
615  fRemnP4 -= TLorentzVector(tP2_2L,E2_2L);
616 
617  // get random phi angle, distributed uniformally in 360 deg
618  PHI3 = 2 * kPi * rnd->RndFsi().Rndm();
619 
620  tP3L = P3zL*bDir + P3tL*tTrans;
621  tP4L = P4zL*bDir + P4tL*tTrans;
622 
623  tP3L.Rotate(PHI3,bDir); // randomize transverse components
624  tP4L.Rotate(PHI3,bDir);
625 
626  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
627  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
628 
629  // create t particle w/ appropriate momenta, code, and status
630  // set target's mom to be the mom of the hadron that was cloned
631  GHepParticle * t = new GHepParticle(*p);
632  t->SetFirstMother(p->FirstMother());
633  t->SetLastMother(p->LastMother());
634 
635  TLorentzVector t4P4L(tP4L,E4L);
636  t->SetPdgCode(s2code);
637  t->SetMomentum(t4P4L);
639 
640  // adjust p to reflect scattering
641  TLorentzVector t4P3L(tP3L,E3L);
642  p->SetPdgCode(scode);
643  p->SetMomentum(t4P3L);
645 
646 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
647  LOG("HNIntranuke2018", pDEBUG)
648  << "|p3| = " << (P3L) << ", E3 = " << (E3L);
649  LOG("HNIntranuke2018", pDEBUG)
650  << "|p4| = " << (P4L) << ", E4 = " << (E4L);
651 #endif
652 
653  ev->AddParticle(*p);
654  ev->AddParticle(*t);
655 
656  delete t; // delete particle clone
657 }
658 //___________________________________________________________________________
660  GHepRecord * ev, GHepParticle * p, INukeFateHN_t fate) const
661 {
662  // scatters particle within nucleus
663 
664 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
665  LOG("HNIntranuke2018", pDEBUG)
666  << "ElasHN() is invoked for a : " << p->Name()
667  << " whose fate is : " << INukeHadroFates::AsString(fate);
668 #endif
669 
670  if(fate!=kIHNFtCEx && fate!=kIHNFtElas)
671  {
672  LOG("HNIntranuke2018", pWARN)
673  << "ElasHN() cannot handle fate: " << INukeHadroFates::AsString(fate);
674  return;
675  }
676 
677  // Random number generator
679 
680  // vars for incoming particle, target, and scattered pdg codes
681  int pcode = p->Pdg();
682  int tcode, scode, s2code;
683  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
684 
685  // Select a target randomly, weighted to #
686  // -- Unless, of course, the fate is CEx,
687  // -- in which case the target may be deterministic
688  // Also assign scattered particle code
689  if(fate==kIHNFtCEx)
690  {
691  if(pcode==kPdgPiP) {tcode = kPdgNeutron; scode = kPdgPi0; s2code = kPdgProton;}
692  else if(pcode==kPdgPiM) {tcode = kPdgProton; scode = kPdgPi0; s2code = kPdgNeutron;}
693  else if(pcode==kPdgKP) {tcode = kPdgNeutron; scode = kPdgK0; s2code = kPdgProton;}
694  else
695  {
696  // for pi0
697  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton) :(kPdgNeutron);
698  scode = (tcode == kPdgProton) ?(kPdgPiP) :(kPdgPiM);
699  s2code = (tcode == kPdgProton) ?(kPdgNeutron):(kPdgProton);
700  }
701  }
702  else
703  {
704  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
705  scode = pcode;
706  s2code = tcode;
707  }
708 
709  // get random scattering angle
710  double C3CM = fHadroData2018->IntBounce(p,tcode,scode,fate);
711  if (C3CM<-1.)
712  {
714  ev->AddParticle(*p);
715  return;
716  }
717 
718  // create scattered particle
719  GHepParticle * t = new GHepParticle(*p);
720  t->SetPdgCode(tcode);
721  double Mt = t->Mass();
722  //t->SetMomentum(TLorentzVector(0,0,0,Mt));
723 
724  // handle fermi momentum
725  if(fDoFermi)
726  {
727  // Handle fermi target
728  Target target(ev->TargetNucleus()->Pdg());
729  //LOG("HAIntranuke2018", pNOTICE) << "Nuclmodel= " << fNuclmodel->ModelType(target) ;
730  target.SetHitNucPdg(tcode);
732  TVector3 tP3L = fFermiFac * fNuclmodel->Momentum3();
733  double tE = TMath::Sqrt(tP3L.Mag2() + Mt*Mt);
734  t->SetMomentum(TLorentzVector(tP3L,tE));
735  }
736  else
737  {
738  t->SetMomentum(TLorentzVector(0,0,0,Mt));
739  }
740 
741  bool pass = utils::intranuke2018::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
743 
744 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
745  LOG("HNIntranuke2018",pDEBUG)
746  << "|p3| = " << P3L << ", E3 = " << E3L;
747  LOG("HNIntranuke2018",pDEBUG)
748  << "|p4| = " << P4L << ", E4 = " << E4L;
749 #endif
750 
751  if (pass==true)
752  {
753  ev->AddParticle(*p);
754  ev->AddParticle(*t);
755  } else
756  {
757  delete t; //fixes memory leak
758  LOG("HNIntranuke2018", pINFO) << "Elastic in hN failed calling TwoBodyCollision";
760  exception.SetReason("hN scattering kinematics through TwoBodyCollision failed");
761  throw exception;
762  }
763 
764  delete t;
765 
766 }
767 //___________________________________________________________________________
769 {
770  // Aaron Meyer (Jan 2010)
771  // Updated version of InelasticHN
772 
773  GHepParticle s1(*p);
774  GHepParticle s2(*p);
775  GHepParticle s3(*p);
776 
777 
779  {
780  // set status of particles and return
781 
785 
786  ev->AddParticle(s1);
787  ev->AddParticle(s2);
788  ev->AddParticle(s3);
789  }
790  else
791  {
792  LOG("HNIntranuke2018", pNOTICE) << "Error: could not create pion production final state";
794  exception.SetReason("PionProduction in hN failed");
795  throw exception;
796  }
797  return;
798 
799 }
800 //___________________________________________________________________________
802 {
803  // This function handles pion photoproduction reactions
804 
805 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
806  LOG("HNIntranuke2018", pDEBUG)
807  << "GammaInelasticHN() is invoked for a : " << p->Name()
808  << " whose fate is : " << INukeHadroFates::AsString(fate);
809 #endif
810 
811  if(fate!=kIHNFtInelas && p->Pdg()!=kPdgGamma)
812  {
813  LOG("HNIntranuke2018", pWARN)
814  << "GammaInelasticHN() cannot handle fate: " << INukeHadroFates::AsString(fate);
815  return;
816  }
817 
818  // random number generator
820 
821  // vars for incoming particle, target, and scattered reaction products
822  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
823  int pcode = p->Pdg();
824  int tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
825  int scode, s2code;
826  double ke = p->KinE() / units::MeV;
827 
828  LOG("HNIntranuke2018", pNOTICE)
829  << "Particle code: " << pcode << ", target: " << tcode;
830 
831 
832  if (rnd->RndFsi().Rndm() * (fHadroData2018 -> XSecGamp_fs() -> Evaluate(ke) +
833  fHadroData2018 -> XSecGamn_fs() -> Evaluate(ke) )
834  <= fHadroData2018 -> XSecGamp_fs() -> Evaluate(ke) ) { scode = kPdgProton; }
835  else { scode = kPdgNeutron; }
836 
837  //scode=fHadroData2018->AngleAndProduct(p,tcode,C3CM,fate);
838  //double C3CM = 0.0; // cos of scattering angle
839  double C3CM = fHadroData2018->IntBounce(p,tcode,scode,fate);
840 
841  if ((tcode == kPdgProton ) && (scode==kPdgProton )) {s2code=kPdgPi0;}
842  else if ((tcode == kPdgProton ) && (scode==kPdgNeutron)) {s2code=kPdgPiP;}
843  else if ((tcode == kPdgNeutron) && (scode==kPdgProton )) {s2code=kPdgPiM;}
844  else if ((tcode == kPdgNeutron) && (scode==kPdgNeutron)) {s2code=kPdgPi0;}
845  else {
846  LOG("HNIntranuke2018", pERROR)
847  << "Error: could not determine particle final states";
848  ev->AddParticle(*p);
849  return;
850  }
851 
852  LOG("HNIntranuke2018", pNOTICE)
853  << "GammaInelastic fate: " << INukeHadroFates::AsString(fate);
854  LOG("HNIntranuke2018", pNOTICE)
855  << " final state: " << scode << " and " << s2code;
856  LOG("HNIntranuke2018", pNOTICE)
857  << " scattering angle: " << C3CM;
858 
859  GHepParticle * t = new GHepParticle(*p);
860  t->SetPdgCode(tcode);
861  double Mt = t->Mass();
862 
863  // handle fermi momentum
864  if(fDoFermi)
865  {
866  // Handle fermi target
867  Target target(ev->TargetNucleus()->Pdg());
868 
869  target.SetHitNucPdg(tcode);
871  TVector3 tP3L = fFermiFac * fNuclmodel->Momentum3();
872  double tE = TMath::Sqrt(tP3L.Mag2() + Mt*Mt);
873  t->SetMomentum(TLorentzVector(tP3L,tE));
874  }
875  else
876  {
877  t->SetMomentum(TLorentzVector(0,0,0,Mt));
878  }
879 
880  bool pass = utils::intranuke2018::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
882 
883 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
884  LOG("HNIntranuke2018",pDEBUG)
885  << "|p3| = " << P3L << ", E3 = " << E3L;
886  LOG("HNIntranuke2018",pDEBUG)
887  << "|p4| = " << P4L << ", E4 = " << E4L;
888 #endif
889 
890  if (pass==true)
891  {
892  //p->SetStatus(kIStStableFinalState);
893  //t->SetStatus(kIStStableFinalState);
894  ev->AddParticle(*p);
895  ev->AddParticle(*t);
896  } else
897  {
898  ev->AddParticle(*p);
899  }
900 
901  delete t;
902 
903 }
904 //___________________________________________________________________________
906 {
907 
908  // handle compound nucleus option
909  // -- Call the PreEquilibrium function
911  { // random number generator
912  //unused var - quiet compiler warning//RandomGen * rnd = RandomGen::Instance();
913 
914  if((p->KinE() < fEPreEq) )
915  {
916  if(fRemnA>4) //this needs to be matched to what is in PreEq and Eq
917  {
918  GHepParticle * sp = new GHepParticle(*p);
919  sp->SetFirstMother(mom);
920  // this was PreEquilibrium - now just used for hN
921  //same arguement lists for PreEq and Eq
924 
925  delete sp;
926  return 2;
927  }
928  else
929  {
930  // nothing left to interact with!
931  LOG("HNIntranuke2018", pNOTICE)
932  << "*** Nothing left to interact with, escaping.";
933  GHepParticle * sp = new GHepParticle(*p);
934  sp->SetFirstMother(mom);
936  ev->AddParticle(*sp);
937  delete sp;
938  return 1;
939  }
940  }
941  }
942  return 0;
943 }
944 //___________________________________________________________________________
946 {
947  return (this->HandleCompoundNucleus(ev,p,p->FirstMother())==2);
948 }
949 //___________________________________________________________________________
950 
952 {
953  // load hadronic cross sections
955 
956  // fermi momentum setup
957  // this is specifically set in Intranuke2018::Configure(string)
958  fNuclmodel = dynamic_cast<const NuclearModelI *>( this -> SubAlg("NuclearModel") ) ;
959 
960  // other intranuke config params
961  GetParam( "NUCL-R0", fR0 ); // fm
962  GetParam( "NUCL-NR", fNR );
963 
964  GetParam( "INUKE-NucRemovalE", fNucRmvE ); // GeV
965  GetParam( "INUKE-HadStep", fHadStep ) ;
966  GetParam( "INUKE-NucAbsFac", fNucAbsFac ) ;
967  GetParam( "INUKE-NucQEFac", fNucQEFac ) ;
968  GetParam( "INUKE-NucCEXFac", fNucCEXFac ) ;
969  GetParam( "INUKE-Energy_Pre_Eq", fEPreEq ) ;
970  GetParam( "INUKE-FermiFac", fFermiFac ) ;
971  GetParam( "INUKE-FermiMomentum", fFermiMomentum ) ;
972 
973  GetParam( "INUKE-DoCompoundNucleus", fDoCompoundNucleus ) ;
974  GetParam( "INUKE-DoFermi", fDoFermi ) ;
975  GetParam( "INUKE-XsecNNCorr", fXsecNNCorr ) ;
976  GetParamDef( "AltOset", fAltOset, false ) ;
977 
978  GetParam( "HNINUKE-UseOset", fUseOset ) ;
979  GetParam( "HNINUKE-DelRPion", fDelRPion ) ;
980  GetParam( "HNINUKE-DelRNucleon", fDelRNucleon ) ;
981 
982  // report
983  LOG("HNIntranuke2018", pINFO) << "Settings for Intranuke2018 mode: " << INukeMode::AsString(kIMdHN);
984  LOG("HNIntranuke2018", pWARN) << "R0 = " << fR0 << " fermi";
985  LOG("HNIntranuke2018", pWARN) << "NR = " << fNR;
986  LOG("HNIntranuke2018", pWARN) << "DelRPion = " << fDelRPion;
987  LOG("HNIntranuke2018", pWARN) << "DelRNucleon = " << fDelRNucleon;
988  LOG("HNIntranuke2018", pWARN) << "HadStep = " << fHadStep << " fermi";
989  LOG("HNIntranuke2018", pWARN) << "NucAbsFac = " << fNucAbsFac;
990  LOG("HNIntranuke2018", pWARN) << "NucQEFac = " << fNucQEFac;
991  LOG("HNIntranuke2018", pWARN) << "NucCEXFac = " << fNucCEXFac;
992  LOG("HNIntranuke2018", pWARN) << "EPreEq = " << fEPreEq;
993  LOG("HNIntranuke2018", pWARN) << "FermiFac = " << fFermiFac;
994  LOG("HNIntranuke2018", pWARN) << "FermiMomtm = " << fFermiMomentum;
995  LOG("HNIntranuke2018", pWARN) << "DoFermi? = " << ((fDoFermi)?(true):(false));
996  LOG("HNIntranuke2018", pWARN) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
997  LOG("HNIntranuke2018", pWARN) << "useOset = " << fUseOset;
998  LOG("HNIntranuke2018", pWARN) << "altOset = " << fAltOset;
999  LOG("HNIntranuke2018", pWARN) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
1000 }
1001 //___________________________________________________________________________
1002 
1004 {
1005  //LOG("HNIntranuke2018", pWARN) << "IN HadronFateOset";
1006 
1007  //LOG("HNIntranuke2018", pWARN) << "{ frac abs = " << osetUtils::currentInstance->getAbsorptionFraction();
1008  //LOG("HNIntranuke2018", pWARN) << " frac cex = " << osetUtils::currentInstance->getCexFraction() << " }";
1009 
1010  const double fractionAbsorption = osetUtils::currentInstance->
1011  getAbsorptionFraction();
1012  const double fractionCex = osetUtils::currentInstance->getCexFraction ();
1013 
1014  RandomGen *randomGenerator = RandomGen::Instance();
1015  const double randomNumber = randomGenerator->RndFsi().Rndm();
1016 
1017  //LOG("HNIntranuke2018", pWARN) << "{ frac abs = " << fractionAbsorption;
1018  //LOG("HNIntranuke2018", pWARN) << " frac cex = " << fractionCex;
1019  //LOG("HNIntranuke2018", pWARN) << " frac elas = " << 1-fractionAbsorption-fractionCex << " }";
1020 
1021  if (randomNumber < fractionAbsorption && fRemnA > 1) return kIHNFtAbs;
1022  else if (randomNumber < fractionAbsorption + fractionCex) return kIHNFtCEx;
1023  else return kIHNFtElas;
1024 }
static string AsString(INukeMode_t mode)
Definition: INukeMode.h:41
const double kPi
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
void SetFirstMother(int m)
Definition: GHepParticle.h:133
Basic constants.
void PreEquilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
double fFermiMomentum
whether or not particle collision is pauli blocked
const XML_Char * target
Definition: expat.h:268
INukeOset * currentInstance
Definition: INukeOset.cxx:5
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
double Frac(int hpdgc, INukeFateHN_t fate, double ke, int targA=0, int targZ=0) const
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
double fEPreEq
threshold for pre-equilibrium reaction
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
double fFermiFac
testing parameter to modify fermi momentum
const char * p
Definition: xmltok.h:285
static const double MeV
Definition: Units.h:130
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:99
double FateWeight(int pdgc, INukeFateHN_t fate) const
void ElasHN(GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
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)
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
bool HandleCompoundNucleusHN(GHepRecord *ev, GHepParticle *p) const
bool fXsecNNCorr
use nuclear medium correction for NN cross section
void SetMomentum(const TLorentzVector &p4)
Definition: config.py:1
bool IsInNucleus(const GHepParticle *p) const
void ProcessEventRecord(GHepRecord *event_rec) const
double Mass(void) const
Mass that corresponds to the PDG code.
Double_t beta
Timing fit.
double fNucAbsFac
absorption xsec correction factor (hN Mode)
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2018.h:98
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
const int kPdgK0
Definition: PDGCodes.h:151
int LastMother(void) const
Definition: GHepParticle.h:68
int Pdg(void) const
Definition: GHepParticle.h:64
void InelasticHN(GHepRecord *ev, GHepParticle *p) const
int FirstMother(void) const
Definition: GHepParticle.h:67
string Name(void) const
Name that corresponds to the PDG code.
int HandleCompoundNucleus(GHepRecord *ev, GHepParticle *p, int mom) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
enum genie::EINukeFateHN_t INukeFateHN_t
double getCexFraction() const
return fraction of cex events
Definition: INukeOset.h:53
const int kPdgGamma
Definition: PDGCodes.h:166
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
const int kPdgKP
Definition: PDGCodes.h:149
INukeFateHN_t HadronFateOset() const
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.
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
bool fDoFermi
whether or not to do fermi mom.
INukeFateHN_t HadronFateHN(const GHepParticle *p) const
#define pINFO
Definition: Messenger.h:63
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
TVector3 Unit() const
void Equilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Misc GENIE control constants.
double fR0
effective nuclear size param
#define pWARN
Definition: Messenger.h:61
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
virtual void ProcessEventRecord(GHepRecord *event_rec) const
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static INukeHadroData2018 * Instance(void)
void SetLastMother(int m)
Definition: GHepParticle.h:134
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
static string AsString(INukeFateHN_t fate)
bool fAltOset
NuWro&#39;s table-based implementation (not recommended)
INukeHadroData2018 * fHadroData2018
a collection of h+N,h+A data & calculations
Definition: Intranuke2018.h:96
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:314
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
float Mag() const
double fHadStep
step size for intranuclear hadron transport
TRandom3 r(0)
const int kPdgPiM
Definition: PDGCodes.h:136
virtual bool GenerateNucleon(const Target &) const =0
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const int kPdgProton
Definition: PDGCodes.h:65
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
void GammaInelasticHN(GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
double fNucRmvE
binding energy to subtract from cascade nucleons
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
void AbsorbHN(GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
bool fUseOset
Oset model for low energy pion in hN.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
TGeoManager * gm
Definition: make_fe_box.C:4
Root of GENIE utility namespaces.
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
int fRemnZ
remnant nucleus Z
#define pDEBUG
Definition: Messenger.h:64
void SetPdgCode(int c)
TLorentzVector fRemnP4
P4 of remnant system.
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353