HAIntranuke.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Steve Dytman <dytman+@pitt.edu>, Pittsburgh Univ.
8  Aaron Meyer <asm58@pitt.edu>, Pittsburgh Univ.
9  Alex Bell, Pittsburgh Univ.
10  Hugh Gallagher <gallag@minos.phy.tufts.edu>, Tufts Univ.
11  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>, Rutherford Lab.
12  September 20, 2005
13 
14  For the class documentation see the corresponding header file.
15 
16  Important revisions after version 2.0.0 :
17  @ Nov 30, 2007 - SD
18  Changed the hadron tracking algorithm to take into account the radial
19  nuclear density dependence. Using the somewhat empirical approach of
20  increasing the nuclear radius by a const (tunable) number times the tracked
21  particle's de Broglie wavelength as this helps getting the hadron+nucleus
22  cross sections right.
23  @ Mar 08, 2008 - CA
24  Fixed code retrieving the remnant nucleus which stopped working as soon as
25  simulation of nuclear de-excitation started pushing photons in the target
26  nucleus daughter list.
27  @ Jun 20, 2008 - CA
28  Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
29  deleted after it was added at the GHEP event record.
30  @ Jul 15, 2010 - AM
31  Major overhaul of the function of each interaction type. Absorption fates
32  changed to allow more than 6 particles at a time (up to 85 now). PiPro fates
33  now allow the pion to rescatter inside the nucleus, will be changed at a
34  later date. HAIntranuke class is now defined as derived from virtual class.
35  Intranuke.
36  @ Oct 10, 2011 - SD
37  Changes to keep reweighting alive. Add exception handling in ElasHA, InelasticHA,
38  and Inelastic.
39  @ Jan 24, 2012 - SD
40  Add option of doing K+.
41  @ Aug 30, 2016 - SD
42  Fix memory leaks found by Igor Kakorin.
43 */
44 //____________________________________________________________________________
45 
46 #include <cstdlib>
47 #include <sstream>
48 #include <exception>
49 
50 #include <TMath.h>
51 
77 
78 using std::ostringstream;
79 
80 using namespace genie;
81 using namespace genie::utils;
82 using namespace genie::utils::intranuke;
83 using namespace genie::constants;
84 using namespace genie::controls;
85 
86 //___________________________________________________________________________
87 //___________________________________________________________________________
88 // Methods specific to INTRANUKE's HA-mode
89 //___________________________________________________________________________
90 //___________________________________________________________________________
92 Intranuke("genie::HAIntranuke")
93 {
94 
95 }
96 //___________________________________________________________________________
98 Intranuke("genie::HAIntranuke",config)
99 {
100 
101 }
102 //___________________________________________________________________________
104 {
105 
106 }
107 //___________________________________________________________________________
109 {
110  LOG("HAIntranuke", pNOTICE)
111  << "************ Running HA MODE INTRANUKE ************";
112 
114 
115  LOG("HAIntranuke", pINFO) << "Done with this event";
116 }
117 //___________________________________________________________________________
119  GHepRecord* ev, GHepParticle* p) const
120 {
121 // Simulate a hadron interaction for the input particle p in HA mode
122 //
123  // check inputs
124  if(!p || !ev) {
125  LOG("HAIntranuke", pERROR) << "** Null input!";
126  return;
127  }
128 
129  // get particle code and check whether this particle can be handled
130  int pdgc = p->Pdg();
131  bool is_gamma = (pdgc==kPdgGamma);
132  bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
133  bool is_kaon = (pdgc==kPdgKP || pdgc==kPdgKM);
134  bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
135  bool is_handled = (is_baryon || is_pion || is_kaon || is_gamma);
136  if(!is_handled) {
137  LOG("HAIntranuke", pERROR) << "** Can not handle particle: " << p->Name();
138  return;
139  }
140 
141  // select a fate for the input particle
142  INukeFateHA_t fate = this->HadronFateHA(p);
143 
144  // store the fate
145  ev->Particle(p->FirstMother())->SetRescatterCode((int)fate);
146 
147  if(fate == kIHAFtUndefined) {
148  LOG("HAIntranuke", pERROR) << "** Couldn't select a fate";
150  ev->AddParticle(*p);
151  return;
152  }
153  LOG("HAIntranuke", pNOTICE)
154  << "Selected "<< p->Name() << " fate: "<< INukeHadroFates::AsString(fate);
155 
156  // try to generate kinematics - repeat till is done (should seldom need >2)
157  fNumIterations = 0;
159 }
160 //___________________________________________________________________________
162  GHepRecord* ev, GHepParticle* p) const
163 {
164  // get stored fate
165  INukeFateHA_t fate = (INukeFateHA_t)
166  ev->Particle(p->FirstMother())->RescatterCode();
167 
168  LOG("HAIntranuke", pINFO)
169  << "Generating kinematics for " << p->Name()
170  << " fate: "<< INukeHadroFates::AsString(fate);
171 
172  // try to generate kinematics for the selected fate
173 
174  try
175  {
176  fNumIterations++;
177  if (fate == kIHAFtElas)
178  {
179  this->ElasHA(ev,p,fate);
180  }
181  else
182  if (fate == kIHAFtInelas || fate == kIHAFtCEx)
183  {
184  this->InelasticHA(ev,p,fate);
185  }
186  else if (fate == kIHAFtAbs || fate == kIHAFtPiProd)
187  {
188  this->Inelastic(ev,p,fate);
189  }
190  }
192  {
193  LOG("HAIntranuke", pNOTICE)
194  << exception;
195  if(fNumIterations <= 100) {
196  LOG("HAIntranuke", pNOTICE)
197  << "Failed attempt to generate kinematics for "
198  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
199  << " - After " << fNumIterations << " tries, still retrying...";
201  } else {
202  LOG("HAIntranuke", pNOTICE)
203  << "Failed attempt to generate kinematics for "
204  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
205  << " after " << fNumIterations-1
206  << " attempts. Trying a new fate...";
207  this->SimulateHadronicFinalState(ev,p);
208  }
209  }
210 }
211 //___________________________________________________________________________
213 {
214 // Select a hadron fate in HA mode
215 //
217 
218  // get pdgc code & kinetic energy in MeV
219  int pdgc = p->Pdg();
220  double ke = p->KinE() / units::MeV;
221 
222  LOG("HAIntranuke", pINFO)
223  << "Selecting hA 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 = fHadroData->Frac(pdgc, kIHAFtCEx, ke);
234  double frac_elas = fHadroData->Frac(pdgc, kIHAFtElas, ke);
235  double frac_inel = fHadroData->Frac(pdgc, kIHAFtInelas, ke);
236  double frac_abs = fHadroData->Frac(pdgc, kIHAFtAbs, ke);
237  double frac_piprod = fHadroData->Frac(pdgc, kIHAFtPiProd, ke);
238 
239  LOG("HAIntranuke", pINFO)
240  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
241  << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
242  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
243  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
244  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
245 
246  // apply external tweaks to fractions
247  frac_cex *= fPionFracCExScale;
248  frac_elas *= fPionFracElasScale;
249  frac_inel *= fPionFracInelScale;
250  frac_abs *= fPionFracAbsScale;
251  frac_piprod *= fPionFracPiProdScale;
252  double frac_rescale = 1./(frac_cex + frac_elas + frac_inel + frac_abs + frac_piprod);
253  frac_cex *= frac_rescale;
254  frac_elas *= frac_rescale;
255  frac_inel *= frac_rescale;
256  frac_abs *= frac_rescale;
257  frac_piprod *= frac_rescale;
258 
259  LOG("HAIntranuke", pINFO)
260  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
261  << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
262  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
263  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
264  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
265 
266  // compute total fraction (can be <1 if fates have been switched off)
267  double tf = frac_cex +
268  frac_elas +
269  frac_inel +
270  frac_abs +
271  frac_piprod;
272 
273  double r = tf * rnd->RndFsi().Rndm();
274 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
275  LOG("HAIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
276 #endif
277  double cf=0; // current fraction
278  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
279  if(r < (cf += frac_elas )) return kIHAFtElas; // elas
280  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
281  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
282  if(r < (cf += frac_piprod )) return kIHAFtPiProd; // pi prod
283 
284  LOG("HAIntranuke", pWARN)
285  << "No selection after going through all fates! "
286  << "Total fraction = " << tf << " (r = " << r << ")";
287  }
288 
289  // handle nucleons
290  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
291 
292  double frac_cex = fHadroData->Frac(pdgc, kIHAFtCEx, ke);
293  double frac_elas = fHadroData->Frac(pdgc, kIHAFtElas, ke);
294  double frac_inel = fHadroData->Frac(pdgc, kIHAFtInelas, ke);
295  double frac_abs = fHadroData->Frac(pdgc, kIHAFtAbs, ke);
296  double frac_pipro = fHadroData->Frac(pdgc, kIHAFtPiProd, ke);
297 
298  // apply external tweaks to fractions
299  frac_cex *= fNucleonFracCExScale;
300  frac_elas *= fNucleonFracElasScale;
301  frac_inel *= fNucleonFracInelScale;
302  frac_abs *= fNucleonFracAbsScale;
303  frac_pipro *= fNucleonFracPiProdScale;
304  double frac_rescale = 1./(frac_cex + frac_elas + frac_inel + frac_abs + frac_pipro);
305  frac_cex *= frac_rescale;
306  frac_elas *= frac_rescale;
307  frac_inel *= frac_rescale;
308  frac_abs *= frac_rescale;
309  frac_pipro *= frac_rescale;
310 
311  LOG("HAIntranuke", pDEBUG)
312  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
313  << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
314  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
315  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
316  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_pipro;
317 
318  // compute total fraction (can be <1 if fates have been switched off)
319  double tf = frac_cex +
320  frac_elas +
321  frac_inel +
322  frac_abs +
323  frac_pipro;
324 
325  double r = tf * rnd->RndFsi().Rndm();
326 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
327  LOG("HAIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
328 #endif
329  double cf=0; // current fraction
330  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
331  if(r < (cf += frac_elas )) return kIHAFtElas; // elas
332  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
333  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
334  if(r < (cf += frac_pipro )) return kIHAFtPiProd; // pi prod
335 
336  LOG("HAIntranuke", pWARN)
337  << "No selection after going through all fates! "
338  << "Total fraction = " << tf << " (r = " << r << ")";
339  }
340  // handle kaons
341  else if (pdgc==kPdgKP || pdgc==kPdgKM) {
342  double frac_inel = fHadroData->Frac(pdgc, kIHAFtInelas, ke);
343  double frac_abs = fHadroData->Frac(pdgc, kIHAFtAbs, ke);
344  LOG("HAIntranuke", pDEBUG)
345  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
346  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs;
347  // compute total fraction (can be <1 if fates have been switched off)
348  double tf = frac_inel +
349  frac_abs;
350  double r = tf * rnd->RndFsi().Rndm();
351 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
352  LOG("HAIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
353 #endif
354  double cf=0; // current fraction
355  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
356  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
357  }
358  }//iterations
359 
360  return kIHAFtUndefined;
361 }
362 //___________________________________________________________________________
363 double HAIntranuke::PiBounce(void) const
364 {
365 // [adapted from neugen3 intranuke_bounce.F]
366 // [is a fortran stub / difficult to understand - needs to be improved]
367 //
368 // Generates theta in radians for elastic pion-nucleus scattering/
369 // Lookup table is based on Fig 17 of Freedman, Miller and Henley, Nucl.Phys.
370 // A389, 457 (1982)
371 //
372  const int nprob = 25;
373  double dintor = 0.0174533;
374  double denom = 47979.453;
375  double rprob[nprob] = {
376  5000., 4200., 3000., 2600., 2100., 1800., 1200., 750., 500., 230., 120.,
377  35., 9., 3., 11., 18., 29., 27., 20., 14., 10., 6., 2., 0.14, 0.19 };
378 
379  double angles[nprob];
380  for(int i=0; i<nprob; i++) angles[i] = 2.5*i;
381 
383  double r = rnd->RndFsi().Rndm();
384 
385  double xsum = 0.;
386  double theta = 0.;
387  double binl = 0.;
388  double binh = 0.;
389  int tj = 0;
390  for(int i=0; i<60; i++) {
391  theta = i+0.5;
392  for(int j=0; j < nprob-1; j++) {
393  binl = angles[j];
394  binh = angles[j+1];
395  tj=j;
396  if(binl<=theta && binh>=theta) break;
397  tj=0;
398  }//j
399  int itj = tj;
400  double tfract = (theta-binl)/2.5;
401  double delp = rprob[itj+1] - rprob[itj];
402  xsum += (rprob[itj] + tfract*delp)/denom;
403  if(xsum>r) break;
404  theta = 0.;
405  }//i
406 
407  theta *= dintor;
408 
409  LOG("HAIntranuke", pNOTICE)
410  << "Generated pi+A elastic scattering angle = " << theta << " radians";
411 
412  return theta;
413 }
414 //___________________________________________________________________________
415 double HAIntranuke::PnBounce(void) const
416 {
417 // [adapted from neugen3 intranuke_pnbounce.F]
418 // [is a fortran stub / difficult to understand - needs to be improved]
419 //
420 // Generates theta in radians for elastic nucleon-nucleus scattering.
421 // Use 800 MeV p+O16 as template in same (highly simplified) spirit as pi+A
422 // from table in Adams et al., PRL 1979. Guess value at 0-2 deg based on Ni
423 // data.
424 //
425  const int nprob = 20;
426  double dintor = 0.0174533;
427  double denom = 11967.0;
428  double rprob[nprob] = {
429  2400., 2350., 2200., 2000., 1728., 1261., 713., 312., 106., 35.,
430  6., 5., 10., 12., 11., 9., 6., 1., 1., 1. };
431 
432  double angles[nprob];
433  for(int i=0; i<nprob; i++) angles[i] = 1.0*i;
434 
436  double r = rnd->RndFsi().Rndm();
437 
438  double xsum = 0.;
439  double theta = 0.;
440  double binl = 0.;
441  double binh = 0.;
442  int tj = 0;
443  for(int i=0; i< nprob; i++) {
444  theta = i+0.5;
445  for(int j=0; j < nprob-1; j++) {
446  binl = angles[j];
447  binh = angles[j+1];
448  tj=j;
449  if(binl<=theta && binh>=theta) break;
450  tj=0;
451  }//j
452  int itj = tj;
453  double tfract = (theta-binl)/2.5;
454  double delp = rprob[itj+1] - rprob[itj];
455  xsum += (rprob[itj] + tfract*delp)/denom;
456  if(xsum>r) break;
457  theta = 0.;
458  }//i
459 
460  theta *= dintor;
461 
462  LOG("HAIntranuke", pNOTICE)
463  << "Generated N+A elastic scattering angle = " << theta << " radians";
464 
465  return theta;
466 }
467 //___________________________________________________________________________
469  GHepRecord* ev, GHepParticle* p, INukeFateHA_t fate) const
470 {
471  // scatters particle within nucleus, copy of hN code meant to run only once
472  // in hA mode
473 
474 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
475  LOG("HAIntranuke", pDEBUG)
476  << "ElasHA() is invoked for a : " << p->Name()
477  << " whose fate is : " << INukeHadroFates::AsString(fate);
478 #endif
479 
480  if(fate!=kIHAFtElas)
481  {
482  LOG("HAIntranuke", pWARN)
483  << "ElasHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
484  return;
485  }
486 
487  // check remnants
488  if(fRemnA<0 || fRemnZ<0) // best to stop it here and not try again.
489  {
490  LOG("HAIntranuke", pWARN) << "Invalid Nucleus! : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
492  ev->AddParticle(*p);
493  return;
494  }
495 
496  // vars for incoming particle, target, and scattered pdg codes
497  int pcode = p->Pdg();
498  double Mp = p->Mass();
499  double Mt = 0.;
500  if (ev->TargetNucleus()->A()==fRemnA)
501  { Mt = PDGLibrary::Instance()->Find(ev->TargetNucleus()->Pdg())->Mass(); }
502  else
503  {
504  Mt = fRemnP4.M();
505  }
506  TLorentzVector t4PpL = *p->P4();
507  TLorentzVector t4PtL = fRemnP4;
508  double C3CM = 0.0;
509 
510  // calculate scattering angle
511  if(pcode==kPdgNeutron||pcode==kPdgProton) C3CM = TMath::Cos(this->PnBounce());
512  else C3CM = TMath::Cos(this->PiBounce());
513 
514  // calculate final 4 momentum of probe
515  TLorentzVector t4P3L, t4P4L;
516 
517  if (!utils::intranuke::TwoBodyKinematics(Mp,Mt,t4PpL,t4PtL,t4P3L,t4P4L,C3CM,fRemnP4))
518  {
519  LOG("HAIntranuke", pNOTICE) << "ElasHA() failed";
521  exception.SetReason("TwoBodyKinematics failed in ElasHA");
522  throw exception;
523  }
524 
525  // Update probe particle
526  p->SetMomentum(t4P3L);
528 
529  // Update Remnant nucleus
530  fRemnP4 = t4P4L;
531  LOG("HAIntranuke",pINFO)
532  << "C3cm = " << C3CM;
533  LOG("HAIntranuke",pINFO)
534  << "|p3| = " << t4P3L.Vect().Mag() << ", E3 = " << t4P3L.E() << ",Mp = " << Mp;
535  LOG("HAIntranuke",pINFO)
536  << "|p4| = " << fRemnP4.Vect().Mag() << ", E4 = " << fRemnP4.E() << ",Mt = " << Mt;
537 
538  ev->AddParticle(*p);
539 
540 }
541 //___________________________________________________________________________
543  GHepRecord* ev, GHepParticle* p, INukeFateHA_t fate) const
544 {
545  // charge exch and inelastic - scatters particle within nucleus, hA version
546  // each are treated as quasielastic, particle scatters off single nucleon
547 
548 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
549  LOG("HAIntranuke", pDEBUG)
550  << "InelasticHA() is invoked for a : " << p->Name()
551  << " whose fate is : " << INukeHadroFates::AsString(fate);
552 #endif
553  if(ev->Probe() ) {
554  LOG("HAIntranuke", pINFO) << " probe KE = " << ev->Probe()->KinE();
555  }
556  if(fate!=kIHAFtCEx && fate!=kIHAFtInelas)
557  {
558  LOG("HAIntranuke", pWARN)
559  << "InelasticHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
560  return;
561  }
562 
563  // Random number generator
565 
566  // vars for incoming particle, target, and scattered pdg codes
567  int pcode = p->Pdg();
568  int tcode, scode, s2code;
569  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
570 
571  // Select a hadron fate in HN mode
572  INukeFateHN_t h_fate;
573  if (fate == kIHAFtCEx) h_fate = kIHNFtCEx;
574  else h_fate = kIHNFtElas;
575 
576  // Select a target randomly, weighted to #
577  // -- Unless, of course, the fate is CEx,
578  // -- in which case the target may be deterministic
579  // Also assign scattered particle code
580  if(fate==kIHAFtCEx)
581  {
582  if(pcode==kPdgPiP) {tcode = kPdgNeutron; scode = kPdgPi0; s2code = kPdgProton;}
583  else if(pcode==kPdgPiM) {tcode = kPdgProton; scode = kPdgPi0; s2code = kPdgNeutron;}
584  else if(pcode==kPdgPi0)
585  {
586  // for pi0
587  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton) :(kPdgNeutron);
588  scode = (tcode == kPdgProton) ?(kPdgPiP) :(kPdgPiM);
589  s2code = (tcode == kPdgProton) ?(kPdgNeutron):(kPdgProton);
590  }
591  else if(pcode==kPdgProton) {tcode = kPdgNeutron; scode = kPdgNeutron; s2code = kPdgProton;}
592  else if(pcode==kPdgNeutron){tcode = kPdgProton; scode = kPdgProton; s2code = kPdgNeutron;}
593  else
594  { LOG("HAIntranuke", pWARN) << "InelasticHA() cannot handle fate: "
596  << " for particle " << p->Name();
597  return;
598  }
599  }
600  else
601  {
602  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
603  // if(pcode == kPdgKP || pcode == kPdgKM) tcode = kPdgProton;
604  scode = pcode;
605  s2code = tcode;
606  }
607 
608  // check remnants
609  if ( fRemnA < 1 ) //we've blown nucleus apart, no need to retry anything - exit
610  {
611  LOG("HAIntranuke",pNOTICE) << "InelasticHA() stops : not enough nucleons";
613  ev->AddParticle(*p);
614  return;
615  }
616  else if ( fRemnZ + (((pcode==kPdgProton)||(pcode==kPdgPiP))?1:0) - (pcode==kPdgPiM?1:0)
617  < ((( scode==kPdgProton)||( scode==kPdgPiP)) ?1:0) - (scode ==kPdgPiM ?1:0)
618  + (((s2code==kPdgProton)||(s2code==kPdgPiP)) ?1:0) - (s2code==kPdgPiM ?1:0) )
619  {
620  LOG("HAIntranuke",pWARN) << "InelasticHA() failed : too few protons in nucleus";
622  ev->AddParticle(*p);
623  return; // another extreme case, best strategy is to exit and go to next event
624  }
625 
626  GHepParticle t(*p);
627  t.SetPdgCode(tcode);
628 
629  // set up fermi target
630  Target target(ev->TargetNucleus()->Pdg());
631  double tM = t.Mass();
632 
633  // handle fermi momentum
634  if(fDoFermi)
635  {
636  target.SetHitNucPdg(tcode);
638  TVector3 tP3 = fFermiFac * fNuclmodel->Momentum3();
639  double tE = TMath::Sqrt(tP3.Mag2()+ tM*tM);
640  t.SetMomentum(TLorentzVector(tP3,tE));
641  }
642  else
643  {
644  t.SetMomentum(0,0,0,tM);
645  }
646 
647  GHepParticle * cl = new GHepParticle(*p); // clone particle, to run IntBounce at proper energy
648  // calculate energy and momentum using invariant mass
649  double pM = p->Mass();
650  double E_p = ((*p->P4() + *t.P4()).Mag2() - tM*tM - pM*pM)/(2.0*tM);
651  double P_p = TMath::Sqrt(E_p*E_p - pM*pM);
652  cl->SetMomentum(TLorentzVector(P_p,0,0,E_p));
653  // momentum doesn't have to be in right direction, only magnitude
654  double C3CM = fHadroData->IntBounce(cl,tcode,scode,h_fate);
655  delete cl;
656  if (C3CM<-1.) // hope this doesn't occur too often - unphysical but we just pass it on
657  {
658  LOG("HAIntranuke", pWARN) << "unphysical angle chosen in InelasicHA - put particle outside nucleus";
660  ev->AddParticle(*p);
661  return;
662  }
663  double KE1L = p->KinE();
664  double KE2L = t.KinE();
665  LOG("HAIntranuke",pINFO)
666  << " KE1L = " << KE1L << " " << KE1L << " KE2L = " << KE2L;
667  GHepParticle cl1(*p);
668  GHepParticle cl2(t);
669  bool success = utils::intranuke::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
670  &cl1,&cl2,fRemnA,fRemnZ,fRemnP4,kIMdHA);
671  if(success)
672  {
673  double P3L = TMath::Sqrt(cl1.Px()*cl1.Px() + cl1.Py()*cl1.Py() + cl1.Pz()*cl1.Pz());
674  double P4L = TMath::Sqrt(cl2.Px()*cl2.Px() + cl2.Py()*cl2.Py() + cl2.Pz()*cl2.Pz());
675  double E3L = cl1.KinE();
676  double E4L = cl2.KinE();
677  LOG ("HAIntranuke",pINFO) << "Successful quasielastic scattering or charge exchange";
678  LOG("HAIntranuke",pINFO)
679  << "C3CM = " << C3CM << "\n P3L, E3L = "
680  << P3L << " " << E3L << " P4L, E4L = "<< P4L << " " << E4L ;
681  if(ev->Probe() ) { LOG("HAIntranuke",pINFO)
682  << "P4L = " << P4L << " ;E4L= " << E4L << "\n probe KE = " << ev->Probe()->KinE() << "\n";
683  if (ev->Probe() && (E3L>ev->Probe()->E()||E4L>ev->Probe()->E())) //is this redundant?
684  {
686  exception.SetReason("TwoBodyCollison gives KE> probe KE in hA simulation");
687  throw exception;
688  }
689  }
690  // always get here, even nucleon decay
691  ev->AddParticle(cl1);
692  ev->AddParticle(cl2);
693  LOG("HAIntranuke", pDEBUG) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
694 
695  } else
696  {
698  exception.SetReason("TwoBodyCollison failed in hA simulation");
699  throw exception;
700  }
701 
702 }
703 //___________________________________________________________________________
705  GHepRecord* ev, GHepParticle* p, INukeFateHA_t fate) const
706 {
707 
708  // Aaron Meyer (05/25/10)
709  //
710  // Called to handle all absorption and pi production reactions
711  //
712  // Nucleons -> Reaction approximated by exponential decay in p+n (sum) space,
713  // gaussian in p-n (difference) space
714  // -fit to hN simulations p C, Fe, Pb at 200, 800 MeV
715  // -get n from isospin, np-nn smaller by 2
716  // Pions -> Reaction approximated with a modified gaussian in p+n space,
717  // normal gaussian in p-n space
718  // -based on fits to multiplicity distributions of hN model
719  // for pi+ C, Fe, Pb at 250, 500 MeV
720  // -fit sum and diff of nn, np to Gaussian
721  // -get pi0 from isospin, np-nn smaller by 2
722  // -get pi- from isospin, np-nn smaller by 4
723  // -add 2-body absorption to better match McKeown data
724  // Kaons -> no guidance, use same code as pions.
725  //
726  // Normally distributed random number generated using Box-Muller transformation
727  //
728  // Pion production reactions rescatter pions in nucleus, otherwise unchanged from
729  // older versions of GENIE
730  //
731 
732 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
733  LOG("HAIntranuke", pDEBUG)
734  << "Inelastic() is invoked for a : " << p->Name()
735  << " whose fate is : " << INukeHadroFates::AsString(fate);
736 #endif
737 
738  bool allow_dup = true;
739  PDGCodeList list(allow_dup); // list of final state particles
740 
741  // only absorption/pipro fates allowed
742  if (fate == kIHAFtPiProd) {
743 
744  GHepParticle s1(*p);
745  GHepParticle s2(*p);
746  GHepParticle s3(*p);
747 
750 
751  if (success){
752  LOG ("HAIntranuke",pINFO) << " successful pion production fate";
753  // set status of particles and conserve charge/baryon number
754  s1.SetStatus(kIStStableFinalState); //should be redundant
755  // if (pdg::IsPion(s2->Pdg())) s2->SetStatus(kIStHadronInTheNucleus);
757  // if (pdg::IsPion(s3->Pdg())) s3->SetStatus(kIStHadronInTheNucleus);
759 
760  ev->AddParticle(s1);
761  ev->AddParticle(s2);
762  ev->AddParticle(s3);
763 
764  return;
765  }
766  else {
767  LOG("HAIntranuke", pNOTICE) << "Error: could not create pion production final state";
769  exception.SetReason("PionProduction kinematics failed - retry kinematics");
770  throw exception;
771  }
772  }
773 
774  else if (fate==kIHAFtAbs)
775 // tuned for pions - mixture of 2-body and many-body
776 // use same for kaons as there is no guidance
777  {
778  // Instances for reference
779  PDGLibrary * pLib = PDGLibrary::Instance();
781 
782  double ke = p->KinE() / units::MeV;
783  int pdgc = p->Pdg();
784 
785  if (fRemnA<2)
786  {
787  LOG("HAIntranuke", pNOTICE) << "stop propagation - could not create absorption final state: too few particles";
789  ev->AddParticle(*p);
790  return;
791  }
792  if (fRemnZ<1 && (pdgc==kPdgPiM || pdgc==kPdgKM))
793  {
794  LOG("HAIntranuke", pNOTICE) << "stop propagation - could not create absorption final state: Pi- or K- cannot be absorbed by only neutrons";
796  ev->AddParticle(*p);
797  return;
798  }
799  if (fRemnA-fRemnZ<1 && (pdgc==kPdgPiP || pdgc==kPdgKP))
800  {
801  LOG("HAIntranuke", pINFO) << "stop propagation - could not create absorption final state: Pi+ or K+ cannot be absorbed by only protons";
803  ev->AddParticle(*p);
804  return;
805  }
806 
807  // for now, empirical split between multi-nucleon absorption and pi d -> N N
808  //
809  // added 03/21/11 - Aaron Meyer
810  //
811  if (pdg::IsPion(pdgc) && rnd->RndFsi().Rndm()<1.14*(.903-0.00189*fRemnA)*(1.35-0.00467*ke))
812  { // pi d -> N N, probability determined empirically with McKeown data
813 
814  INukeFateHN_t fate_hN=kIHNFtAbs;
815  int t1code,t2code,scode,s2code;
816  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
817 
818  // choose target nucleon
819  // -- fates weighted by values from Engel, Mosel...
820  if (pdgc==kPdgPiP) {
821  double Prob_pipd_pp=2.*ppcnt*(1.-ppcnt);
822  double Prob_pipnn_pn=.083*(1.-ppcnt)*(1.-ppcnt);
823  if (rnd->RndFsi().Rndm()*(Prob_pipd_pp+Prob_pipnn_pn)<Prob_pipd_pp){
824  t1code=kPdgNeutron; t2code=kPdgProton;
825  scode=kPdgProton; s2code=kPdgProton;}
826  else{
827  t1code=kPdgNeutron; t2code=kPdgNeutron;
828  scode=kPdgProton; s2code=kPdgNeutron;}
829  }
830  if (pdgc==kPdgPiM) {
831  double Prob_pimd_nn=2.*ppcnt*(1.-ppcnt);
832  double Prob_pimpp_pn=.083*ppcnt*ppcnt;
833  if (rnd->RndFsi().Rndm()*(Prob_pimd_nn+Prob_pimpp_pn)<Prob_pimd_nn){
834  t1code=kPdgProton; t2code=kPdgNeutron;
835  scode=kPdgNeutron; s2code=kPdgNeutron; }
836  else{
837  t1code=kPdgProton; t2code=kPdgProton;
838  scode=kPdgProton; s2code=kPdgNeutron;}
839  }
840  else { // pi0
841  double Prob_pi0d_pn=0.88*ppcnt*(1.-ppcnt); // 2 * .44
842  double Prob_pi0pp_pp=.14*ppcnt*ppcnt;
843  double Prob_pi0nn_nn=.14*(1.-ppcnt)*(1.-ppcnt);
844  if (rnd->RndFsi().Rndm()*(Prob_pi0d_pn+Prob_pi0pp_pp+Prob_pi0nn_nn)<Prob_pi0d_pn){
845  t1code=kPdgNeutron; t2code=kPdgProton;
846  scode=kPdgNeutron; s2code=kPdgProton; }
847  else if (rnd->RndFsi().Rndm()*(Prob_pi0d_pn+Prob_pi0pp_pp+Prob_pi0nn_nn)<(Prob_pi0d_pn+Prob_pi0pp_pp)){
848  t1code=kPdgProton; t2code=kPdgProton;
849  scode=kPdgProton; s2code=kPdgProton; }
850  else {
851  t1code=kPdgNeutron; t2code=kPdgNeutron;
852  scode=kPdgNeutron; s2code=kPdgNeutron; }
853  }
854  LOG("HAIntranuke",pINFO) << "choose 2 body absorption, probe, fs = " << pdgc <<" "<< scode <<" "<<s2code;
855  // assign proper masses
856  //double M1 = pLib->Find(pdgc) ->Mass();
857  double M2_1 = pLib->Find(t1code)->Mass();
858  double M2_2 = pLib->Find(t2code)->Mass();
859  //double M2 = M2_1 + M2_2;
860  double M3 = pLib->Find(scode) ->Mass();
861  double M4 = pLib->Find(s2code)->Mass();
862 
863  // handle fermi momentum
864  double E2_1L, E2_2L;
865  TVector3 tP2_1L, tP2_2L;
866  //TLorentzVector dNucl_P4;
867  Target target(ev->TargetNucleus()->Pdg());
868  if(fDoFermi)
869  {
870  target.SetHitNucPdg(t1code);
872  tP2_1L=fFermiFac * fNuclmodel->Momentum3();
873  E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
874 
875  target.SetHitNucPdg(t2code);
877  tP2_2L=fFermiFac * fNuclmodel->Momentum3();
878  E2_2L = TMath::Sqrt(tP2_2L.Mag2() + M2_2*M2_2);
879 
880  //dNucl_P4=TLorentzVector(tP2_1L+tP2_2L,E2_1L+E2_2L);
881  }
882  else
883  {
884  tP2_1L.SetXYZ(0.0, 0.0, 0.0);
885  E2_1L = M2_1;
886 
887  tP2_2L.SetXYZ(0.0, 0.0, 0.0);
888  E2_2L = M2_2;
889  }
890  TLorentzVector dNucl_P4=TLorentzVector(tP2_1L+tP2_2L,E2_1L+E2_2L);
891 
892  double E2L = E2_1L + E2_2L;
893 
894  // adjust p to reflect scattering
895  // get random scattering angle
896  double C3CM = fHadroData->IntBounce(p,t1code,scode,fate_hN);
897  if (C3CM<-1.)
898  {
899  LOG("HAIntranuke", pWARN) << "Inelastic() failed: IntBounce returned bad angle - try again";
901  exception.SetReason("unphysical angle for hN scattering");
902  throw exception;
903  return;
904  }
905 
906  TLorentzVector t4P1L,t4P2L,t4P3L,t4P4L;
907  t4P1L=*p->P4();
908  t4P2L=TLorentzVector(TVector3(tP2_1L+tP2_2L),E2L);
909  double bindE=0.075; // set to fit McKeown data
910  //double bindE=0.0;
911  if (utils::intranuke::TwoBodyKinematics(M3,M4,t4P1L,t4P2L,t4P3L,t4P4L,C3CM,fRemnP4,bindE))
912  {
913  if (pdgc==kPdgPiP || pdgc==kPdgKP) fRemnZ++;
914  if (pdgc==kPdgPiM || pdgc==kPdgKM) fRemnZ--;
915  if (t1code==kPdgProton) fRemnZ--;
916  if (t2code==kPdgProton) fRemnZ--;
917  fRemnA-=2;
918 
919  fRemnP4-=dNucl_P4;
920 
921  // create t particles w/ appropriate momenta, code, and status
922  // Set target's mom to be the mom of the hadron that was cloned
923  GHepParticle* t1 = new GHepParticle(*p);
924  GHepParticle* t2 = new GHepParticle(*p);
925  t1->SetFirstMother(p->FirstMother());
926  t1->SetLastMother(p->LastMother());
927  t2->SetFirstMother(p->FirstMother());
928  t2->SetLastMother(p->LastMother());
929 
930  // adjust p to reflect scattering
931  t1->SetPdgCode(scode);
932  t1->SetMomentum(t4P3L);
933 
934  t2->SetPdgCode(s2code);
935  t2->SetMomentum(t4P4L);
936 
939 
940  ev->AddParticle(*t1);
941  ev->AddParticle(*t2);
942  delete t1;
943  delete t2;
944 
945  return;
946  }
947  else
948  {
949  LOG("HAIntranuke", pNOTICE) << "Inelastic in hA failed calling TwoBodyKineamtics";
951  exception.SetReason("Pion absorption kinematics through TwoBodyKinematics failed");
952  throw exception;
953 
954  }
955 
956  } // end pi d -> N N
957  else // multi-nucleon
958  {
959 
960  // declare some parameters for double gaussian and determine values chosen
961  // parameters for proton and pi+, others come from isospin transformations
962 
963  double ns0=0; // mean - sum of nucleons
964  double nd0=0; // mean - difference of nucleons
965  double Sig_ns=0; // std dev - sum
966  double Sig_nd=0; // std dev - diff
967  double gam_ns=0; // exponential decay rate (for nucleons)
968 
969  if ( pdg::IsNeutronOrProton (pdgc) ) // nucleon probe
970  {
971  // antisymmetric about Z=N
972  if (fRemnA-fRemnZ > fRemnZ)
973  nd0 = 135.227 * TMath::Exp(-7.124*(fRemnA-fRemnZ)/double(fRemnA)) - 2.762;
974  else
975  nd0 = -135.227 * TMath::Exp(-7.124* fRemnZ /double(fRemnA)) + 4.914;
976 
977  Sig_nd = 2.034 + fRemnA * 0.007846;
978 
979  double c1 = 0.041 + ke * 0.0001525;
980  double c2 = -0.003444 - ke * 0.00002324;
981 //change last factor from 30 to 15 so that gam_ns always larger than 0
982 //add check to be certain
983  double c3 = 0.064 - ke * 0.000015;
984  gam_ns = c1 * TMath::Exp(c2*fRemnA) + c3;
985  if(gam_ns<0.002) gam_ns = 0.002;
986  //gam_ns = 10.;
987  LOG("HAIntranuke", pINFO) << "nucleon absorption";
988  LOG("HAIntranuke", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
989  // LOG("HAIntranuke", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
990  LOG("HAIntranuke", pINFO) << "--> gam_ns = " << gam_ns;
991  }
992  else if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) //pion probe
993  {
994  ns0 = .0001*(1.+ke/250.) * (fRemnA-10)*(fRemnA-10) + 3.5;
995  nd0 = (1.+ke/250.) - ((fRemnA/200.)*(1. + 2.*ke/250.));
996  Sig_ns = (10. + 4. * ke/250.)*TMath::Power(fRemnA/250.,0.9); //(1. - TMath::Exp(-0.02*fRemnA));
997  Sig_nd = 4*(1 - TMath::Exp(-0.03*ke));
998  LOG("HAIntranuke", pINFO) << "pion absorption";
999  LOG("HAIntranuke", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1000  LOG("HAIntranuke", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1001  }
1002  else if (pdgc==kPdgKP || pdgc==kPdgKM) // kaon probe
1003  {
1004  ns0 = (rnd->RndFsi().Rndm()>0.5?3:2);
1005  nd0 = 1.;
1006  Sig_ns = 0.1;
1007  Sig_nd = 0.1;
1008  LOG("HAIntranuke", pINFO) << "kaon absorption - set ns, nd later";
1009  // LOG("HAIntranuke", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1010  // LOG("HAIntranuke", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1011  }
1012  else
1013  {
1014  LOG("HAIntranuke", pWARN) << "Inelastic() cannot handle absorption reaction for " << p->Name();
1015  }
1016 
1017  // account for different isospin
1018  if (pdgc==kPdgPi0 || pdgc==kPdgNeutron) nd0-=2.;
1019  if (pdgc==kPdgPiM) nd0-=4.;
1020 
1021  int iter=0; // counter
1022  int np=0,nn=0; // # of p, # of n
1023  bool not_done=true;
1024  double u1 = 0, u2 = 0;
1025 
1026  while (not_done)
1027  {
1028  // infinite loop check
1029  if (iter>=100) {
1030  LOG("HAIntranuke", pNOTICE) << "Error: could not choose absorption final state";
1031  LOG("HAIntranuke", pNOTICE) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1032  LOG("HAIntranuke", pNOTICE) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1033  LOG("HAIntranuke", pNOTICE) << "--> gam_ns = " << gam_ns;
1034  LOG("HAIntranuke", pNOTICE) << "--> A = " << fRemnA << ", Z = " << fRemnZ << ", Energy = " << ke;
1036  exception.SetReason("Absorption choice of # of p,n failed");
1037  throw exception;
1038  }
1039  //here??
1040 
1041  // Box-Muller transform
1042  // Takes two uniform distribution random variables on (0,1]
1043  // Creates two normally distributed random variables on (0,inf)
1044 
1045  u1 = rnd->RndFsi().Rndm(); // uniform random variable 1
1046  u2 = rnd->RndFsi().Rndm(); // " " 2
1047  if (u1==0) u1 = rnd->RndFsi().Rndm();
1048  if (u2==0) u2 = rnd->RndFsi().Rndm(); // Just in case
1049 
1050  // normally distributed random variable
1051  double x2 = TMath::Sqrt(-2*TMath::Log(u1))*TMath::Sin(2*kPi*u2);
1052 
1053  double ns = 0;
1054 
1055  if ( pdg::IsNeutronOrProton (pdgc) ) //nucleon probe
1056  {
1057  ns = -TMath::Log(rnd->RndFsi().Rndm())/gam_ns; // exponential random variable
1058  }
1059  if ( pdg::IsKaon (pdgc) ) //charged kaon probe - either 2 or 3 nucleons to stay simple
1060  {
1061  ns = (rnd->RndFsi().Rndm()<0.5?2:3);
1062  }
1063  else if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) //pion probe
1064  {
1065  // Pion fit for sum takes for xs*exp((xs-x0)^2 / 2*sig_xs0)
1066  // Find random variable by first finding gaussian random variable
1067  // then throwing the value against a linear P.D.F.
1068  //
1069  // max is the maximum value allowed for the random variable (10 std + mean)
1070  // minimum allowed value is 0
1071 
1072  double max = ns0 + Sig_ns * 10;
1073  if(max>fRemnA) max=fRemnA;
1074  double x1 = 0;
1075  bool not_found = true;
1076  int iter2 = 0;
1077 
1078  while (not_found)
1079  {
1080  // infinite loop check
1081  if (iter2>=100)
1082  {
1083  LOG("HAIntranuke", pNOTICE) << "Error: stuck in random variable loop for ns";
1084  LOG("HAIntranuke", pNOTICE) << "--> mean of sum parent distr = " << ns0 << ", Stand dev = " << Sig_ns;
1085  LOG("HAIntranuke", pNOTICE) << "--> A = " << fRemnA << ", Z = " << fRemnZ << ", Energy = " << ke;
1086 
1088  exception.SetReason("Random number generator for choice of #p,n final state failed - unusual - redo kinematics");
1089  throw exception;
1090  }
1091 
1092  // calculate exponential random variable
1093  u1 = rnd->RndFsi().Rndm();
1094  u2 = rnd->RndFsi().Rndm();
1095  if (u1==0) u1 = rnd->RndFsi().Rndm();
1096  if (u2==0) u2 = rnd->RndFsi().Rndm();
1097  x1 = TMath::Sqrt(-2*TMath::Log(u1))*TMath::Cos(2*kPi*u2);
1098 
1099  ns = ns0 + Sig_ns * x1;
1100  if ( ns>max || ns<0 ) {iter2++; continue;}
1101  else if ( rnd->RndFsi().Rndm() > (ns/max) ) {iter2++; continue;}
1102  else {
1103  // accept this sum value
1104  not_found=false;
1105  }
1106  } //while(not_found)
1107  }//else pion
1108 
1109  double nd = nd0 + Sig_nd * x2; // difference (p-n) for both pion, nucleon probe
1110  if (pdgc==kPdgKP) // special for KP
1111  { if (ns==2) nd=0;
1112  if (ns>2) nd=1; }
1113 
1114  np = int((ns+nd)/2.+.5); // Addition of .5 for rounding correction
1115  nn = int((ns-nd)/2.+.5);
1116 
1117  LOG("HAIntranuke", pINFO) << "ns = "<<ns<<", nd = "<<nd<<", np = "<<np<<", nn = "<<nn;
1118  //LOG("HAIntranuke", pNOTICE) << "RemA = "<<fRemnA<<", RemZ = "<<fRemnZ<<", probe = "<<pdgc;
1119 
1120  /*if ((ns+nd)/2. < 0 || (ns-nd)/2. < 0) {iter++; continue;}
1121  else */
1122  //check for problems befor executing phase space 'decay'
1123  if (np < 0 || nn < 0 ) {iter++; continue;}
1124  else if (np + nn < 2. ) {iter++; continue;}
1125  else if ((np + nn == 2.) && pdg::IsNeutronOrProton (pdgc)) {iter++; continue;}
1126  else if (np > fRemnZ + ((pdg::IsProton(pdgc) || pdgc==kPdgPiP || pdgc==kPdgKP)?1:0)
1127  - ((pdgc==kPdgPiM || pdgc==kPdgKM)?1:0)) {iter++; continue;}
1128  else if (nn > fRemnA-fRemnZ + ((pdg::IsNeutron(pdgc)||pdgc==kPdgPiM||pdgc==kPdgKM)?1:0)
1129  - ((pdgc==kPdgPiP||pdgc==kPdgKP)?1:0)) {iter++; continue;}
1130  else {
1131  not_done=false; //success
1132  LOG("HAIntranuke",pINFO) << "success, iter = " << iter << " np, nn = " << np << " " << nn;
1133  if (np+nn>86) // too many particles, scale down
1134  {
1135  double frac = 85./double(np+nn);
1136  np = int(np*frac);
1137  nn = int(nn*frac);
1138  }
1139 
1140  if ( (np==fRemnZ +((pdg::IsProton (pdgc)||pdgc==kPdgPiP||pdgc==kPdgKP)?1:0)-(pdgc==kPdgPiM||pdgc==kPdgKM?1:0))
1141  &&(nn==fRemnA-fRemnZ+((pdg::IsNeutron(pdgc)||pdgc==kPdgPiM||pdgc==kPdgKM)?1:0)-(pdgc==kPdgPiP||pdgc==kPdgKP?1:0)) )
1142  { // leave at least one nucleon in the nucleus to prevent excess momentum
1143  if (rnd->RndFsi().Rndm()<np/(double)(np+nn)) np--;
1144  else nn--;
1145  }
1146 
1147  LOG("HAIntranuke", pNOTICE) << "Final state chosen; # protons : "
1148  << np << ", # neutrons : " << nn;
1149  }
1150  } //while(not_done)
1151 
1152  // change remnants to reflect probe
1153  if ( pdgc==kPdgProton || pdgc==kPdgPiP || pdgc==kPdgKP) fRemnZ++;
1154  if ( pdgc==kPdgPiM || pdgc==kPdgKM) fRemnZ--;
1155  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA++;
1156 
1157  // PhaseSpaceDecay forbids anything over 18 particles
1158  //
1159  // If over 18 particles, split into 5 groups and perform decay on each group
1160  // Anything over 85 particles reduced to 85 in previous step
1161  //
1162  // 4 of the nucleons are used as "probes" as well as original probe,
1163  // with each one sharing 1/5 of the total incident momentum
1164  //
1165  // Note: choosing 5 groups and distributing momentum evenly was arbitrary
1166  // Needs to be revised later
1167 
1168  if ((np+nn)>18)
1169  {
1170  // code lists
1171  PDGCodeList list0(allow_dup);
1172  PDGCodeList list1(allow_dup);
1173  PDGCodeList list2(allow_dup);
1174  PDGCodeList list3(allow_dup);
1175  PDGCodeList list4(allow_dup);
1176  PDGCodeList* listar[5] = {&list0, &list1, &list2, &list3, &list4};
1177 
1178  //set up HadronClusters
1179  // simple for now, each (of 5) hadron cluster has 1/5 of mom and KE
1180 
1181  double probM = pLib->Find(pdgc) ->Mass();
1182  TVector3 pP3 = p->P4()->Vect() * (1./5.);
1183  double probKE = p->P4()->E() -probM;
1184  double clusKE = probKE * (1./5.);
1185  TLorentzVector clusP4(pP3,clusKE); //no mass
1186 
1187  TLorentzVector X4(*p->X4());
1189 
1190  int mom = p->FirstMother();
1191 
1192  GHepParticle * p0 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1193  GHepParticle * p1 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1194  GHepParticle * p2 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1195  GHepParticle * p3 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1196  GHepParticle * p4 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1197 
1198  // To conserve 4-momenta
1199  // fRemnP4 -= probP4 + protP4*np_p + neutP4*(4-np_p) - *p->P4();
1200  fRemnP4 -= 5.*clusP4 - *p->P4();
1201 
1202  for (int i=0;i<(np+nn);i++)
1203  {
1204  if (i<np)
1205  {
1206  listar[i%5]->push_back(kPdgProton);
1207  fRemnZ--;
1208  }
1209  else listar[i%5]->push_back(kPdgNeutron);
1210  fRemnA--;
1211  }
1212  for (int i=0;i<5;i++)
1213  {
1214  LOG("HAIntranuke", pINFO) << "List" << i << " size: " << listar[i]->size();
1215  if (listar[i]->size() <2)
1216  {
1218  exception.SetReason("too few particles for Phase Space decay - try again");
1219  throw exception;
1220  }
1221  }
1222 
1223  // commented out to better fit with absorption reactions
1224  // Add the fermi energy of the nucleons to the phase space
1225  /*if(fDoFermi)
1226  {
1227  GHepParticle* p_ar[5] = {cl, p1, p2, p3, p4};
1228  for (int i=0;i<5;i++)
1229  {
1230  Target target(ev->TargetNucleus()->Pdg());
1231  TVector3 pBuf = p_ar[i]->P4()->Vect();
1232  double mBuf = p_ar[i]->Mass();
1233  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1234  TLorentzVector tSum(pBuf,eBuf);
1235  double mSum = 0.0;
1236  vector<int>::const_iterator pdg_iter;
1237  for(pdg_iter=++(listar[i]->begin());pdg_iter!=listar[i]->end();++pdg_iter)
1238  {
1239  target.SetHitNucPdg(*pdg_iter);
1240  fNuclmodel->GenerateNucleon(target);
1241  mBuf = pLib->Find(*pdg_iter)->Mass();
1242  mSum += mBuf;
1243  pBuf = fFermiFac * fNuclmodel->Momentum3();
1244  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1245  tSum += TLorentzVector(pBuf,eBuf);
1246  fRemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
1247  }
1248  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
1249  p_ar[i]->SetMomentum(dP4);
1250  }
1251  }*/
1252 
1253  bool success1 = utils::intranuke::PhaseSpaceDecay(ev,p0,*listar[0],fRemnP4,fNucRmvE,kIMdHA);
1254  bool success2 = utils::intranuke::PhaseSpaceDecay(ev,p1,*listar[1],fRemnP4,fNucRmvE,kIMdHA);
1255  bool success3 = utils::intranuke::PhaseSpaceDecay(ev,p2,*listar[2],fRemnP4,fNucRmvE,kIMdHA);
1256  bool success4 = utils::intranuke::PhaseSpaceDecay(ev,p3,*listar[3],fRemnP4,fNucRmvE,kIMdHA);
1257  bool success5 = utils::intranuke::PhaseSpaceDecay(ev,p4,*listar[4],fRemnP4,fNucRmvE,kIMdHA);
1258  if(success1 && success2 && success3 && success4 && success5)
1259  {
1260  LOG("HAIntranuke", pINFO)<<"Successful many-body absorption - n>=18";
1261  }
1262  else
1263  {
1264  // try to recover
1265  LOG("HAIntranuke", pWARN) << "PhaseSpace decay fails for HadrCluster- recovery likely incorrect - rethrow event";
1267  ev->AddParticle(*p);
1268  fRemnA+=np+nn;
1269  fRemnZ+=np;
1270  if ( pdgc==kPdgProton || pdgc==kPdgPiP ) fRemnZ--;
1271  if ( pdgc==kPdgPiM ) fRemnZ++;
1272  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA--;
1273  /* exceptions::INukeException exception;
1274  exception.SetReason("Phase space generation of absorption final state failed");
1275  throw exception;
1276  */
1277  }
1278 
1279  // delete cl;
1280  delete p0;
1281  delete p1;
1282  delete p2;
1283  delete p3;
1284  delete p4;
1285 
1286  }
1287  else // less than 18 particles pion
1288  {
1289  if (pdgc==kPdgKP) list.push_back(kPdgKP); //normally conserve strangeness
1290  if (pdgc==kPdgKM) list.push_back(kPdgKM);
1291  for (int i=0;i<np;i++)
1292  {
1293  list.push_back(kPdgProton);
1294  fRemnA--;
1295  fRemnZ--;
1296  }
1297  for (int i=0;i<nn;i++)
1298  {
1299  list.push_back(kPdgNeutron);
1300  fRemnA--;
1301  }
1302 
1303  // Library instance for reference
1304  //PDGLibrary * pLib = PDGLibrary::Instance();
1305 
1306  // commented out to better fit with absorption reactions
1307  // Add the fermi energy of the nucleons to the phase space
1308  /*if(fDoFermi)
1309  {
1310  Target target(ev->TargetNucleus()->Pdg());
1311  TVector3 pBuf = p->P4()->Vect();
1312  double mBuf = p->Mass();
1313  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1314  TLorentzVector tSum(pBuf,eBuf);
1315  double mSum = 0.0;
1316  vector<int>::const_iterator pdg_iter;
1317  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
1318  {
1319  target.SetHitNucPdg(*pdg_iter);
1320  fNuclmodel->GenerateNucleon(target);
1321  mBuf = pLib->Find(*pdg_iter)->Mass();
1322  mSum += mBuf;
1323  pBuf = fFermiFac * fNuclmodel->Momentum3();
1324  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1325  tSum += TLorentzVector(pBuf,eBuf);
1326  fRemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
1327  }
1328  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
1329  p->SetMomentum(dP4);
1330  }*/
1331 
1332  LOG("HAIntranuke", pDEBUG)
1333  << "Remnant nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")";
1334  LOG("HAIntranuke", pINFO) << " list size: " << np+nn;
1335  if (np+nn <2)
1336  {
1338  exception.SetReason("too few particles for Phase Space decay - try again");
1339  throw exception;
1340  }
1341  // GHepParticle * cl = new GHepParticle(*p);
1342  // cl->SetPdgCode(kPdgDecayNuclCluster);
1344  if (success)
1345  {
1346  LOG ("HAIntranuke",pINFO) << "Successful many-body absorption, n<=18";
1347  }
1348  else {
1349  // recover
1351  ev->AddParticle(*p);
1352  fRemnA+=np+nn;
1353  fRemnZ+=np;
1354  if ( pdgc==kPdgProton || pdgc==kPdgPiP ) fRemnZ--;
1355  if ( pdgc==kPdgPiM ) fRemnZ++;
1356  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA--;
1358  exception.SetReason("Phase space generation of absorption final state failed");
1359  throw exception;
1360  }
1361  }
1362  } // end multi-nucleon FS
1363  }
1364  else // not absorption/pipro
1365  {
1366  LOG("HAIntranuke", pWARN)
1367  << "Inelastic() can not handle fate: " << INukeHadroFates::AsString(fate);
1368  return;
1369  }
1370 }
1371 //___________________________________________________________________________
1373  GHepRecord* /*ev*/, GHepParticle* /*p*/, int /*mom*/) const
1374 {
1375  // only relevant for hN mode
1376  return false;
1377 }
1378 //___________________________________________________________________________
1379 void HAIntranuke::Configure (string param_set) {
1380 
1381  Algorithm::Configure(param_set) ;
1382  this -> LoadConfig() ;
1383 }
1384 //___________________________________________________________________________
1386 {
1387  // load hadronic cross sections
1389 
1390  // fermi momentum setup
1391  fNuclmodel = dynamic_cast<const NuclearModelI *>( this -> SubAlg("NuclearModel") ) ;
1392 
1393  // other intranuke config params
1394  GetParam( "NUCL-R0", fR0 ); // fm
1395  GetParam( "NUCL-NR", fNR );
1396 
1397  GetParam( "INUKE-NucRemovalE", fNucRmvE ); // GeV
1398  GetParam( "INUKE-HadStep", fHadStep ) ;
1399  GetParam( "INUKE-NucAbsFac", fNucAbsFac ) ;
1400  GetParam( "INUKE-NucCEXFac", fNucCEXFac ) ;
1401  GetParam( "INUKE-Energy_Pre_Eq", fEPreEq ) ;
1402  GetParam( "INUKE-FermiFac", fFermiFac ) ;
1403  GetParam( "INUKE-FermiMomentum", fFermiMomentum ) ;
1404 
1405  GetParam( "INUKE-DoCompoundNucleus", fDoCompoundNucleus ) ;
1406  GetParam( "INUKE-DoFermi", fDoFermi ) ;
1407 
1408  GetParam( "HAINUKE-DelRPion", fDelRPion ) ;
1409  GetParam( "HAINUKE-DelRNucleon", fDelRNucleon ) ;
1410 
1411  GetParamDef( "FSI-Pion-MFPScale", fPionMFPScale, 1.0 ) ;
1412  GetParamDef( "FSI-Pion-FracCExScale", fPionFracCExScale, 1.0 ) ;
1413  GetParamDef( "FSI-Pion-FracAbsScale", fPionFracAbsScale, 1.0 ) ;
1414  GetParamDef( "FSI-Pion-FracElasScale", fPionFracElasScale, 1.0 ) ;
1415  GetParamDef( "FSI-Pion-FracInelScale", fPionFracInelScale, 1.0 ) ;
1416  GetParamDef( "FSI-Pion-FracPiProdScale", fPionFracPiProdScale, 1.0 ) ;
1417  GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
1418  GetParamDef( "FSI-Nucleon-FracCExScale", fNucleonFracCExScale , 1.0 ) ;
1419  GetParamDef( "FSI-Nucleon-FracElasScale", fNucleonFracElasScale, 1.0 ) ;
1420  GetParamDef( "FSI-Nucleon-FracInelScale", fNucleonFracInelScale, 1.0 ) ;
1421  GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
1422  GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
1423 
1424  // report
1425  LOG("HAIntranuke", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
1426  LOG("HAIntranuke", pINFO) << "R0 = " << fR0 << " fermi";
1427  LOG("HAIntranuke", pINFO) << "NR = " << fNR;
1428  LOG("HAIntranuke", pINFO) << "DelRPion = " << fDelRPion;
1429  LOG("HAIntranuke", pINFO) << "DelRNucleon = " << fDelRNucleon;
1430  LOG("HAIntranuke", pINFO) << "HadStep = " << fHadStep << " fermi";
1431  LOG("HAIntranuke", pINFO) << "EPreEq = " << fHadStep << " fermi";
1432  LOG("HAIntranuke", pINFO) << "NucAbsFac = " << fNucAbsFac;
1433  LOG("HAIntranuke", pINFO) << "NucCEXFac = " << fNucCEXFac;
1434  LOG("HAIntranuke", pINFO) << "FermiFac = " << fFermiFac;
1435  LOG("HAIntranuke", pINFO) << "FermiMomtm = " << fFermiMomentum;
1436  LOG("HAIntranuke", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1437  LOG("HAIntranuke", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1438 }
1439 //___________________________________________________________________________
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke.h:112
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
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:289
Basic constants.
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
const XML_Char * target
Definition: expat.h:268
double fNucleonFracInelScale
Definition: Intranuke.h:126
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
INukeHadroData * fHadroData
a collection of h+N,h+A data & calculations
Definition: Intranuke.h:93
double E(void) const
Get energy.
Definition: GHepParticle.h:92
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke.h:98
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Float_t x1[n_points_granero]
Definition: compare.C:5
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
static const double MeV
Definition: Units.h:130
double PnBounce(void) const
The INTRANUKE intranuclear hadron transport MC. Is a concrete implementation of the EventRecordVisito...
Definition: Intranuke.h:54
double PiBounce(void) const
double fHadStep
step size for intranuclear hadron transport
Definition: Intranuke.h:107
double fNucleonFracCExScale
Definition: Intranuke.h:124
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
double Mass(Resonance_t res)
resonance mass (GeV)
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition: Intranuke.cxx:114
void SetMomentum(const TLorentzVector &p4)
Definition: config.py:1
bool IsKaon(int pdgc)
Definition: PDGUtils.cxx:294
double Mass(void) const
Mass that corresponds to the PDG code.
double fNucleonFracAbsScale
Definition: Intranuke.h:127
Timing fit.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double Frac(int hpdgc, INukeFateHA_t fate, double ke) const
unsigned int fNumIterations
Definition: HAIntranuke.h:79
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
bool HandleCompoundNucleus(GHepRecord *ev, GHepParticle *p, int mom) const
void InelasticHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
A list of PDG codes.
Definition: PDGCodeList.h:33
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
int LastMother(void) const
Definition: GHepParticle.h:68
int Pdg(void) const
Definition: GHepParticle.h:64
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke.h:106
int FirstMother(void) const
Definition: GHepParticle.h:67
c2
Definition: demo5.py:33
string Name(void) const
Name that corresponds to the PDG code.
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
double fNucleonMFPScale
Definition: Intranuke.h:123
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
const int kPdgCompNuclCluster
Definition: PDGCodes.h:194
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke.h:105
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#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 fPionFracPiProdScale
Definition: Intranuke.h:122
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
const int kPdgKM
Definition: PDGCodes.h:150
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)
const int kPdgGamma
Definition: PDGCodes.h:166
double fPionFracInelScale
Definition: Intranuke.h:120
static INukeHadroData * Instance(void)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const int kPdgKP
Definition: PDGCodes.h:149
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils.cxx:754
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
Definition: Intranuke.h:115
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
double fR0
effective nuclear size param
Definition: Intranuke.h:102
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
const double j
Definition: BetheBloch.cxx:29
def success(message)
Definition: log.py:5
double frac(double x)
Fractional part.
void Inelastic(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
#define pINFO
Definition: Messenger.h:63
int fRemnA
remnant nucleus A
Definition: Intranuke.h:96
double t2
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
Misc GENIE control constants.
INukeFateHA_t HadronFateHA(const GHepParticle *p) const
#define pWARN
Definition: Messenger.h:61
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke.h:104
float Mag2() const
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
double fNucAbsFac
absorption xsec correction factor (hN Mode)
Definition: Intranuke.h:108
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
void SetLastMother(int m)
Definition: GHepParticle.h:134
double fEPreEq
threshold for pre-equilibrium reaction
Definition: Intranuke.h:110
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:324
double fPionFracAbsScale
Definition: Intranuke.h:121
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
static string AsString(INukeFateHN_t fate)
double fPionFracCExScale
Definition: Intranuke.h:118
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:314
double fPionFracElasScale
Definition: Intranuke.h:119
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
Definition: Intranuke.h:109
int fRemnZ
remnant nucleus Z
Definition: Intranuke.h:97
void ElasHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
TRandom3 r(0)
const int kPdgPiM
Definition: PDGCodes.h:136
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke.h:113
virtual bool GenerateNucleon(const Target &) const =0
c1
Definition: demo5.py:24
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke.h:111
const int kPdgProton
Definition: PDGCodes.h:65
int A(void) const
#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
double fNucleonFracPiProdScale
Definition: Intranuke.h:128
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
double fPionMFPScale
tweaking factors for tuning
Definition: Intranuke.h:117
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
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
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: INukeUtils.cxx:625
void ProcessEventRecord(GHepRecord *event_rec) const
enum genie::EINukeFateHA_t INukeFateHA_t
double fNucleonFracElasScale
Definition: Intranuke.h:125
Root of GENIE utility namespaces.
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke.h:95
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
virtual void Configure(string param_set)
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
Definition: Intranuke.h:103
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
#define pDEBUG
Definition: Messenger.h:64
void SetPdgCode(int c)
double Py(void) const
Get Py.
Definition: GHepParticle.h:90
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353