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