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