GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HAIntranuke2018.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, 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 <c.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 
54 #include "Framework/Algorithm/AlgConfigPool.h"
55 #include "Framework/Algorithm/AlgFactory.h"
56 #include "Framework/Conventions/GBuild.h"
57 #include "Framework/Conventions/Constants.h"
58 #include "Framework/Conventions/Controls.h"
59 #include "Framework/EventGen/EVGThreadException.h"
60 #include "Framework/GHEP/GHepFlags.h"
61 #include "Framework/GHEP/GHepStatus.h"
62 #include "Framework/GHEP/GHepRecord.h"
63 #include "Framework/GHEP/GHepParticle.h"
64 #include "Physics/HadronTransport/INukeException.h"
65 #include "Physics/HadronTransport/Intranuke2018.h"
66 #include "Physics/HadronTransport/HAIntranuke2018.h"
67 #include "Physics/HadronTransport/INukeHadroData2018.h"
68 #include "Physics/HadronTransport/INukeUtils2018.h"
69 #include "Framework/Interaction/Interaction.h"
70 #include "Framework/Messenger/Messenger.h"
71 #include "Framework/Numerical/RandomGen.h"
72 #include "Framework/Numerical/Spline.h"
73 #include "Framework/ParticleData/PDGLibrary.h"
74 #include "Framework/ParticleData/PDGCodes.h"
75 #include "Framework/ParticleData/PDGCodeList.h"
76 #include "Framework/ParticleData/PDGUtils.h"
77 #include "Framework/Utils/PrintUtils.h"
78 #include "Physics/NuclearState/NuclearUtils.h"
79 #include "Physics/NuclearState/NuclearModelI.h"
80 #include "Physics/NuclearState/NuclearModelMap.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 //___________________________________________________________________________
96 HAIntranuke2018::HAIntranuke2018() :
97 Intranuke2018("genie::HAIntranuke2018")
98 {
99 
100 }
101 //___________________________________________________________________________
102 HAIntranuke2018::HAIntranuke2018(string config) :
103 Intranuke2018("genie::HAIntranuke2018",config)
104 {
105 
106 }
107 //___________________________________________________________________________
108 HAIntranuke2018::~HAIntranuke2018()
109 {
110 
111 }
112 //___________________________________________________________________________
113 void HAIntranuke2018::ProcessEventRecord(GHepRecord * evrec) const
114 {
115  LOG("HAIntranuke2018", pNOTICE)
116  << "************ Running hA2018 MODE INTRANUKE ************";
117  GHepParticle * nuclearTarget = evrec -> TargetNucleus();
118  nuclA = nuclearTarget -> A();
119 
120  Intranuke2018::ProcessEventRecord(evrec);
121 
122  LOG("HAIntranuke2018", pINFO) << "Done with this event";
123 }
124 //___________________________________________________________________________
125 void HAIntranuke2018::SimulateHadronicFinalState(
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";
156  p->SetStatus(kIStStableFinalState);
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;
165  this->SimulateHadronicFinalStateKinematics(ev,p);
166 }
167 //___________________________________________________________________________
168 void HAIntranuke2018::SimulateHadronicFinalStateKinematics(
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";
200  utils::intranuke2018::PreEquilibrium(ev,p,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fNuclmodel,fNucRmvE,kIMdHA); //should be kiMdHA or HN?
201  }
202  }
203  catch(exceptions::INukeException exception)
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...";
212  this->SimulateHadronicFinalStateKinematics(ev,p);
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 //___________________________________________________________________________
224 INukeFateHA_t HAIntranuke2018::HadronFateHA(const GHepParticle * p) const
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 //___________________________________________________________________________
481 void HAIntranuke2018::ElasHA(GHepRecord* ev, GHepParticle* p,
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<<')';
502  p->SetStatus(kIStStableFinalState);
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";
531  exceptions::INukeException exception;
532  exception.SetReason("TwoBodyKinematics failed in ElasHA");
533  throw exception;
534  }
535 
536  // Update probe particle
537  p->SetMomentum(t4P3L);
538  p->SetStatus(kIStStableFinalState);
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 //___________________________________________________________________________
553 void HAIntranuke2018::InelasticHA(
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: "
606  << INukeHadroFates::AsString(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";
623  p->SetStatus(kIStStableFinalState);
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";
632  p->SetStatus(kIStStableFinalState);
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);
648  fNuclmodel->GenerateNucleon(target);
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";
670  p->SetStatus(kIStStableFinalState);
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";
722  exceptions::INukeException exception;
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  {
732  exceptions::INukeException exception;
733  exception.SetReason("TwoBodyCollison failed in hA simulation");
734  throw exception;
735  }
736 
737 }
738 //___________________________________________________________________________
739 void HAIntranuke2018::Inelastic(
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 
783  bool success = utils::intranuke2018::PionProduction(
784  ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel);
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);
791  s2.SetStatus(kIStStableFinalState);
792  // if (pdg::IsPion(s3->Pdg())) s3->SetStatus(kIStHadronInTheNucleus);
793  s3.SetStatus(kIStStableFinalState);
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";
803  exceptions::INukeException exception;
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";
823  p->SetStatus(kIStStableFinalState);
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";
830  p->SetStatus(kIStStableFinalState);
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";
837  p->SetStatus(kIStStableFinalState);
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  else 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  double random_number = rnd->RndFsi().Rndm();
880  if (random_number*(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 (random_number*(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);
907  fNuclmodel->GenerateNucleon(target);
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);
913  fNuclmodel->GenerateNucleon(target);
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";
937  exceptions::INukeException exception;
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);
1000  t2->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";
1012  exceptions::INukeException exception;
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;
1097  exceptions::INukeException exception;
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 
1149  exceptions::INukeException exception;
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());
1251  GHepStatus_t ist = kIStNucleonClusterTarget;
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  {
1280  exceptions::INukeException exception;
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";
1353  p->SetStatus(kIStStableFinalState);
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());
1399  GHepStatus_t ist = kIStNucleonClusterTarget;
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  {
1453  exceptions::INukeException exception;
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);
1460  bool success = utils::intranuke2018::PhaseSpaceDecay(ev,p0,list,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
1490  p->SetStatus(kIStStableFinalState);
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--;
1497  exceptions::INukeException exception;
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 //___________________________________________________________________________
1521 void HAIntranuke2018::LoadConfig(void)
1522 {
1523 
1524  // load hadronic cross sections
1525  fHadroData2018 = INukeHadroData2018::Instance();
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-ChargedPion-MFPScale", fChPionMFPScale, 1.0 ) ;
1553  GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
1554  GetParamDef( "FSI-Pion-FracCExScale", fPionFracCExScale, 1.0 ) ;
1555  GetParamDef( "FSI-ChargedPion-FracAbsScale", fChPionFracAbsScale, 1.0 ) ;
1556  GetParamDef( "FSI-NeutralPion-FracAbsScale", fNeutralPionFracAbsScale,1.0 ) ;
1557  GetParamDef( "FSI-Pion-FracInelScale", fPionFracInelScale, 1.0 ) ;
1558  GetParamDef( "FSI-Pion-FracPiProdScale", fPionFracPiProdScale, 1.0 ) ;
1559  GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
1560  GetParamDef( "FSI-Nucleon-FracCExScale", fNucleonFracCExScale , 1.0 ) ;
1561  GetParamDef( "FSI-Nucleon-FracInelScale", fNucleonFracInelScale, 1.0 ) ;
1562  GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
1563  GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
1564 
1565  // report
1566  LOG("HAIntranuke2018", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
1567  LOG("HAIntranuke2018", pINFO) << "R0 = " << fR0 << " fermi";
1568  LOG("HAIntranuke2018", pINFO) << "NR = " << fNR;
1569  LOG("HAIntranuke2018", pINFO) << "DelRPion = " << fDelRPion;
1570  LOG("HAIntranuke2018", pINFO) << "DelRNucleon = " << fDelRNucleon;
1571  LOG("HAIntranuke2018", pINFO) << "HadStep = " << fHadStep << " fermi";
1572  LOG("HAIntranuke2018", pINFO) << "EPreEq = " << fHadStep << " fermi";
1573  LOG("HAIntranuke2018", pINFO) << "NucAbsFac = " << fNucAbsFac;
1574  LOG("HAIntranuke2018", pINFO) << "NucCEXFac = " << fNucCEXFac;
1575  LOG("HAIntranuke2018", pINFO) << "FermiFac = " << fFermiFac;
1576  LOG("HAIntranuke2018", pINFO) << "FermiMomtm = " << fFermiMomentum;
1577  LOG("HAIntranuke2018", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1578  LOG("HAIntranuke2018", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1579  LOG("HAIntranuke2018", pINFO) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
1580 }
1581 //___________________________________________________________________________
1582 /*
1583 INukeFateHA_t HAIntranuke2018::HadronFateOset () const
1584 {
1585  const double fractionAbsorption = osetUtils::currentInstance->
1586  getAbsorptionFraction();
1587  const double fractionCex = osetUtils::currentInstance->getCexFraction ();
1588 
1589  RandomGen *randomGenerator = RandomGen::Instance();
1590  const double randomNumber = randomGenerator->RndFsi().Rndm();
1591 
1592  LOG("HAIntranuke2018", pINFO)
1593  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << fractionCex
1594  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << 1-fractionCex-fractionAbsorption
1595  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << fractionAbsorption;
1596  if (randomNumber < fractionAbsorption && fRemnA > 1) return kIHAFtAbs;
1597  else if (randomNumber < fractionAbsorption + fractionCex) return kIHAFtCEx;
1598  else return kIHAFtInelas;
1599 }
1600 */
genie::INukeMode::AsString
static string AsString(INukeMode_t mode)
Definition: INukeMode.h:41
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::GHepParticle::SetFirstMother
void SetFirstMother(int m)
Definition: GHepParticle.h:132
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
genie::utils::intranuke2018::PreEquilibrium
void PreEquilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Definition: INukeUtils2018.cxx:491
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
genie::Intranuke2018::fFermiMomentum
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke2018.h:135
pERROR
#define pERROR
Definition: Messenger.h:59
genie::HAIntranuke2018::SimulateHadronicFinalState
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2018.cxx:125
genie::HAIntranuke2018::HandleCompoundNucleus
int HandleCompoundNucleus(GHepRecord *ev, GHepParticle *p, int mom) const
Definition: HAIntranuke2018.cxx:1513
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::Intranuke2018::fDelRPion
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke2018.h:128
genie::utils::intranuke2018::PhaseSpaceDecay
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
Definition: INukeUtils2018.cxx:1678
genie::Intranuke2018::fEPreEq
double fEPreEq
threshold for pre-equilibrium reaction
Definition: Intranuke2018.h:133
genie::Intranuke2018::fNR
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the &quot;nuclear bounda...
Definition: Intranuke2018.h:126
genie::Intranuke2018::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2018.h:134
genie::HAIntranuke2018::fNumIterations
unsigned int fNumIterations
Definition: HAIntranuke2018.h:82
genie::Intranuke2018::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:119
genie::utils::intranuke2018::PionProduction
bool PionProduction(GHepRecord *ev, GHepParticle *p, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI *Nuclmodel)
Definition: INukeUtils2018.cxx:1301
genie::Intranuke2018::fChPionMFPScale
double fChPionMFPScale
tweaking factors for tuning
Definition: Intranuke2018.h:143
genie::HAIntranuke2018::PiBounce
double PiBounce(void) const
Definition: HAIntranuke2018.cxx:376
genie::HAIntranuke2018::Inelastic
void Inelastic(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2018.cxx:739
genie::NuclearModelI
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:436
genie::Intranuke2018::fXsecNNCorr
bool fXsecNNCorr
use nuclear medium correction for NN cross section
Definition: Intranuke2018.h:141
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:236
genie::units::ns
static constexpr double ns
Definition: Units.h:98
genie::pdg::IsKaon
bool IsKaon(int pdgc)
Definition: PDGUtils.cxx:331
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
genie::HAIntranuke2018::ProcessEventRecord
void ProcessEventRecord(GHepRecord *event_rec) const
Definition: HAIntranuke2018.cxx:113
genie::HAIntranuke2018::InelasticHA
void InelasticHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2018.cxx:553
genie::Intranuke2018::fNucAbsFac
double fNucAbsFac
absorption xsec correction factor (hN Mode)
Definition: Intranuke2018.h:131
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::GHepParticle::Pz
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
genie::Intranuke2018::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2018.h:118
genie::NuclearModelI::Momentum3
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
genie::PDGCodeList
A list of PDG codes.
Definition: PDGCodeList.h:32
genie::GHepParticle::Px
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
genie::GHepRecord::Probe
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
genie::GHepParticle::LastMother
int LastMother(void) const
Definition: GHepParticle.h:67
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::INukeHadroData2018::FracADep
double FracADep(int hpdgc, INukeFateHA_t fate, double ke, int targA) const
Definition: INukeHadroData2018.cxx:1377
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::kPdgCompNuclCluster
const int kPdgCompNuclCluster
Definition: PDGCodes.h:217
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::INukeFateHN_t
enum genie::EINukeFateHN_t INukeFateHN_t
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::kPdgGamma
const int kPdgGamma
Definition: PDGCodes.h:189
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::utils::intranuke2018::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils2018.cxx:873
genie::utils::intranuke2018::TwoBodyCollision
bool TwoBodyCollision(GHepRecord *ev, int pcode, int tcode, int scode, int s2code, double C3CM, GHepParticle *p, GHepParticle *t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode=kIMdHA)
Intranuke utility functions.
Definition: INukeUtils2018.cxx:729
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::INukeHadroData2018::FracAIndep
double FracAIndep(int hpdgc, INukeFateHA_t fate, double ke) const
Definition: INukeHadroData2018.cxx:1420
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
genie::Intranuke2018::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2018.h:136
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke2018::fNucCEXFac
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
Definition: Intranuke2018.h:132
genie::Intranuke2018::fR0
double fR0
effective nuclear size param
Definition: Intranuke2018.h:125
pWARN
#define pWARN
Definition: Messenger.h:60
genie::HAIntranuke2018::ElasHA
void ElasHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2018.cxx:481
genie::Intranuke2018::fDoCompoundNucleus
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
Definition: Intranuke2018.h:138
genie::Intranuke2018::ProcessEventRecord
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition: Intranuke2018.cxx:119
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:159
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::HAIntranuke2018::PnBounce
double PnBounce(void) const
Definition: HAIntranuke2018.cxx:428
genie::INukeHadroData2018::Instance
static INukeHadroData2018 * Instance(void)
Definition: INukeHadroData2018.cxx:175
genie::GHepParticle::SetLastMother
void SetLastMother(int m)
Definition: GHepParticle.h:133
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::GHepParticle::X4
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Intranuke2018::fAltOset
bool fAltOset
NuWro&#39;s table-based implementation (not recommended)
Definition: Intranuke2018.h:140
genie::Target::SetHitNucPdg
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
genie::Intranuke2018::fHadroData2018
INukeHadroData2018 * fHadroData2018
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke2018.h:116
genie::Intranuke2018::fDelRNucleon
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke2018.h:129
genie::pdg::IsNeutronOrProton
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:351
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
genie::Intranuke2018::fHadStep
double fHadStep
step size for intranuclear hadron transport
Definition: Intranuke2018.h:130
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::HAIntranuke2018::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2018.cxx:168
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
genie::HAIntranuke2018::nuclA
int nuclA
value of A for the target nucleus in hA mode
Definition: HAIntranuke2018.h:81
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::GHepParticle::A
int A(void) const
Definition: GHepParticle.cxx:184
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Intranuke2018::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke2018.h:127
genie::controls::kRjMaxIterations
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
genie::NuclearModelI::GenerateNucleon
virtual bool GenerateNucleon(const Target &) const =0
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::GHepRecord
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
genie::Intranuke2018::fUseOset
bool fUseOset
Oset model for low energy pion in hN.
Definition: Intranuke2018.h:139
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
genie::HAIntranuke2018::HadronFateHA
INukeFateHA_t HadronFateHA(const GHepParticle *p) const
Definition: HAIntranuke2018.cxx:224
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
genie::INukeHadroData2018::IntBounce
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
Definition: INukeHadroData2018.cxx:1667
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
genie::Intranuke2018::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2018.h:120
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::GHepParticle::SetPdgCode
void SetPdgCode(int c)
Definition: GHepParticle.cxx:230
genie::Intranuke2018::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2018.h:121
genie::GHepParticle::Py
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
Generated on Wed Jul 17 2024 14:46:41 for GENIEGenerator by   doxygen 1.8.5