INukeUtils2018.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Jim Dobson <j.dobson07 \at imperial.ac.uk>
8  Imperial College London
9 
10  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
11  University of Liverpool & STFC Rutherford Appleton Lab
12 
13  Aaron Meyer <asm58 \at pitt.edu>
14  Pittsburgh University
15 
16  For documentation see the corresponding header file.
17 
18  Important revisions after version 2.0.0 :
19  @ Mar 04, 2009 - JD
20  Was first added in v2.5.1. Adapted from the T2K GENIE reweighting tool.
21  @ Mar 05, 2009 - CA
22  Modified ReconstructHadronFateHA() to work with hadron+A event files in
23  addition to neutrino event files.
24  @ Sep 10, 2009 - CA
25  Added MeanFreePath(), Dist2Exit(), Dist2ExitMFP()
26  @ Sep 30, 2009 - CA
27  Added StepParticle() from Intranuke.cxx
28  @ Oct 02, 2009 - CA
29  Added test MeanFreePath_Delta().
30  @ Jul 15, 2010 - AM
31  Added common utility functions used by both hA and hN mode. Updated
32  MeanFreePath to separate proton and neutron cross sections. Added general
33  utility functions.
34  @ Jan 9, 2015 - SD, NG, TG
35  Added 2014 version of INTRANUKE codes for v2.9.0. Uses INukeHadroData2014,
36  but no changes to mean free path.
37 */
38 //____________________________________________________________________________
39 
40 #include <TLorentzVector.h>
41 #include <TMath.h>
42 #include <TSystem.h>
43 
69 #include "TComplex.h"
70 
71 using std::ostringstream;
72 using namespace genie;
73 using namespace genie::utils;
74 using namespace genie::constants;
75 using namespace genie::controls;
76 
77 //____________________________________________________________________________
79  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4,
80  double A, double Z, double nRpi, double nRnuc, const bool useOset, const bool altOset, const bool xsecNNCorr, string INukeMode)
81 {
82 // Calculate the mean free path (in fm) for a pions and nucleons in a nucleus
83 //
84 // Inputs
85 // pdgc : Hadron PDG code
86 // x4 : Hadron 4-position in the nucleus coordinate system (units: fm)
87 // p4 : Hadron 4-momentum (units: GeV)
88 // A : Nucleus atomic mass number
89 // nRpi : Controls the pion ring size in terms of de-Broglie wavelengths
90 // nRnuc: Controls the nuclepn ring size in terms of de-Broglie wavelengths
91 //
92  bool is_pion = pdgc == kPdgPiP || pdgc == kPdgPi0 || pdgc == kPdgPiM;
93  bool is_nucleon = pdgc == kPdgProton || pdgc == kPdgNeutron;
94  bool is_kaon = pdgc == kPdgKP;
95  bool is_gamma = pdgc == kPdgGamma;
96 
97  if(!is_pion && !is_nucleon && !is_kaon && !is_gamma) return 0.;
98 
99  // before getting the nuclear density at the current position
100  // check whether the nucleus has to become larger by const times the
101  // de Broglie wavelength -- that is somewhat empirical, but this
102  // is what is needed to get piA total cross sections right.
103  // The ring size is different for light nuclei (using gaus density) /
104  // heavy nuclei (using woods-saxon density).
105  // The ring size is different for pions / nucleons.
106  //
107  double momentum = p4.Vect().Mag(); // hadron momentum in GeV
108  double ring = (momentum>0) ? 1.240/momentum : 0; // de-Broglie wavelength
109 
110  if(A<=20) { ring /= 2.; }
111 
112  /*
113  if (is_pion ) { ring *= nRpi; }
114  else if (is_nucleon ) { ring *= nRnuc; }
115  else if (is_gamma || is_kaon || useOset) { ring = 0.; }
116  */
117  if(INukeMode=="hN2018")
118  {
119  if (is_pion ) { ring *= nRpi; }
120  else if (is_nucleon ) { ring *= nRnuc; }
121  else if (is_gamma || is_kaon || useOset) { ring = 0.;}
122  }
123  else
124  {
125  if (is_pion || is_kaon ) { ring *= nRpi; }
126  else if (is_nucleon ) { ring *= nRnuc; }
127  else if (is_gamma ) { ring = 0.; }
128  }
129 
130  // get the nuclear density at the current position
131  double rnow = x4.Vect().Mag();
132  double rho = A * utils::nuclear::Density(rnow,(int) A,ring);
133 
134  // the hadron+nucleon cross section will be evaluated within the range
135  // of the input spline and assumed to be const outside that range
136  //
137  double ke = (p4.Energy() - p4.M()) / units::MeV; // kinetic energy in MeV
138  ke = TMath::Max(INukeHadroData2018::fMinKinEnergy, ke);
139  ke = TMath::Min(INukeHadroData2018::fMaxKinEnergyHN, ke);
140 
141  // get total xsection for the incident hadron at its current
142  // kinetic energy
143  double sigtot = 0;
144  double ppcnt = (double) Z/ (double) A; // % of protons remaining
146 
147  if (is_pion and (INukeMode == "hN2018") and useOset and ke < 350.0)
148  sigtot = sigmaTotalOset (ke, rho, pdgc, ppcnt, altOset);
149  else if (pdgc == kPdgPiP)
150  { sigtot = fHadroData2018 -> XSecPipp_Tot() -> Evaluate(ke)*ppcnt;
151  sigtot+= fHadroData2018 -> XSecPipn_Tot() -> Evaluate(ke)*(1-ppcnt);}
152  else if (pdgc == kPdgPi0)
153  { sigtot = fHadroData2018 -> XSecPi0p_Tot() -> Evaluate(ke)*ppcnt;
154  sigtot+= fHadroData2018 -> XSecPi0n_Tot() -> Evaluate(ke)*(1-ppcnt);}
155  else if (pdgc == kPdgPiM)
156  { sigtot = fHadroData2018 -> XSecPipn_Tot() -> Evaluate(ke)*ppcnt;
157  sigtot+= fHadroData2018 -> XSecPipp_Tot() -> Evaluate(ke)*(1-ppcnt);}
158  else if (pdgc == kPdgProton)
159  {
160  sigtot = fHadroData2018 -> XSecPp_Tot() -> Evaluate(ke)*ppcnt;
161  //sigtot+= fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*(1-ppcnt);
162 
163  PDGLibrary * pLib = PDGLibrary::Instance();
164  double hc = 197.327;
165  double R0 = 1.25 * TMath::Power(A,1./3.) + 2.0 * 0.65; // should all be in units of fm
166  double Mp = pLib->Find(2212)->Mass();
167  double M = pLib->Find(pdgc)->Mass();
168  //double E = (p4.Energy() - Mp) * 1000.; // Convert GeV to MeV.
169  double E = ke;
170  if (Z*hc/137./x4.Vect().Mag() > E) // Coulomb correction (Cohen, Concepts of Nuclear Physics, pg. 259-260)
171  {
172  double z = 1.0; // charge for single proton
173  double Bc = z*Z*hc/137./R0;
174  double x = E/Bc;
175  double f = TMath::ACos(TMath::Power(x,0.5)) - TMath::Power(x*(1-x),0.5);
176  double B = 0.63*z*Z*TMath::Power((M/Mp)/E,0.5);
177  double Pc = TMath::Exp(-B*f);
178  sigtot *= Pc;
179  }
180  sigtot+= fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*(1-ppcnt);
181 
182  double E0 = TMath::Power(A,0.2)*12.;
183  if (INukeMode=="hN2018"){if(ke<E0){sigtot=0.0;}} //empirical - needed to cut off large number of low energy nucleons
184  // LOG("INukeUtils",pDEBUG) "sigtot for proton= " << sigtot << "; KE= " << ke << "; E0= " << E0;
185  }
186  else if (pdgc == kPdgNeutron)
187  {
188  sigtot = fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*ppcnt;
189  sigtot+= fHadroData2018 -> XSecNn_Tot() -> Evaluate(ke)*(1-ppcnt);
190  double E0 = TMath::Power(A,0.2)*12.;
191  if (INukeMode=="hN2018"){if(ke<E0){sigtot=0.0;}} //empirical - needed to cut off large number of low energy nucleons
192  // LOG("INukeUtils",pDEBUG) "sigtot for neutron= " << sigtot << "; KE= " << ke;
193  }
194  else if (pdgc == kPdgKP)
195  { sigtot = fHadroData2018 -> XSecKpN_Tot() -> Evaluate(ke);
196  // this factor is used to empirically get agreement with tot xs data, justified historically.
197  sigtot*=1.1;}
198  else if (pdgc == kPdgGamma)
199  { sigtot = fHadroData2018 -> XSecGamp_fs() -> Evaluate(ke)*ppcnt;
200  sigtot+= fHadroData2018 -> XSecGamn_fs() -> Evaluate(ke)*(1-ppcnt);}
201  else {
202  return 0;
203  }
204 
205  // the xsection splines in INukeHadroData return the hadron x-section in
206  // mb -> convert to fm^2
207  sigtot *= (units::mb / units::fm2);
208 
209  // avoid defective error handling
210  //if(sigtot<1E-6){sigtot=1E-6;}
211 
212  if (xsecNNCorr and is_nucleon)
213  sigtot *= INukeNucleonCorr::getInstance()->
214  getAvgCorrection (rho, A, p4.E() - PDGLibrary::Instance()->Find(pdgc)->Mass()); //uses lookup tables
215 
216  // avoid defective error handling
217  if(sigtot<1E-6){sigtot=1E-6;}
218 
219  // compute the mean free path
220  double lamda = 1. / (rho * sigtot);
221 
222  // exits if lamda is InF (if cross section is 0)
223  if( ! TMath::Finite(lamda) ) {
224  return -1;
225  }
226 
227 /*
228  LOG("INukeUtils", pDEBUG)
229  << "sig_total = " << sigtot << " fm^2, rho = " << rho
230  << " fm^-3 => mfp = " << lamda << " fm.";
231 */
232  return lamda;
233 }
234 //____________________________________________________________________________
236  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4, double A)
237 {
238 //
239 // **test**
240 //
241 
242 // Calculate the mean free path (in fm) for Delta's in a nucleus
243 //
244 // Inputs
245 // pdgc : Hadron PDG code
246 // x4 : Hadron 4-position in the nucleus coordinate system (units: fm)
247 // p4 : Hadron 4-momentum (units: GeV)
248 // A : Nucleus atomic mass number
249 //
250  bool is_deltapp = (pdgc==kPdgP33m1232_DeltaPP);
251  if(!is_deltapp) return 0.;
252 
253  // get the nuclear density at the current position
254  double rnow = x4.Vect().Mag();
255  double rho = A * utils::nuclear::Density(rnow,(int) A);
256 
257  // the Delta+N->N+N cross section will be evaluated within the range
258  // of the input spline and assumed to be const outside that range
259  double ke = (p4.Energy() - p4.M()) / units::MeV; // kinetic energy in MeV
260  ke = TMath::Max(1500., ke);
261  ke = TMath::Min( 0., ke);
262 
263  // get the Delta+N->N+N cross section
264  double sig = 0;
265  if (ke< 500) sig=20;
266  else if (ke<1000) sig=40;
267  else sig=30;
268 
269  // value is in mb -> convert to fm^2
270  sig *= (units::mb / units::fm2);
271 
272  // compute the mean free path
273  double lamda = 1. / (rho * sig);
274 
275  // exits if lamda is InF (if cross section is 0)
276  if( ! TMath::Finite(lamda) ) {
277  return -1;
278  }
279 
280  return lamda;
281 }
282 //____________________________________________________________________________
284  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4, double A, double Z,
285  double mfp_scale_factor, double nRpi, double nRnuc, double NR, double R0)
286 {
287 // Calculate the survival probability for a hadron inside a nucleus
288 //
289 // Inputs
290 // pdgc : Hadron PDG code
291 // x4 : Hadron 4-position in the nucleus coordinate system (units: fm)
292 // p4 : Hadron 4-momentum (units: GeV)
293 // A : Nucleus atomic mass number
294 // mfp_scale_factor: Tweaks the mean free path (mfp -> mfp*scale). Def: 1.0
295 // nRpi: Controls the pion ring size in terms of de-Broglie wavelengths
296 // nRnuc: Controls the nuclepn ring size in terms of de-Broglie wavelengths
297 // NR: How far away to track the hadron, in terms of the corresponding
298 // nuclear radius. Def: 3
299 // R0: R0 in R=R0*A^1/3 (units:fm). Def. 1.4
300 
301  double prob = 1.0;
302 
303  double step = 0.05; // fermi
304  double R = NR * R0 * TMath::Power(A, 1./3.);
305 
306  TVector3 dr3 = p4.Vect().Unit(); // unit vector along its direction
307  TLorentzVector dr4(dr3,0);
308 
309  LOG("INukeUtils", pDEBUG)
310  << "Calculating survival probability for hadron with PDG code = " << pdgc
311  << " and momentum = " << p4.P() << " GeV";
312  LOG("INukeUtils", pDEBUG)
313  << "mfp scale = " << mfp_scale_factor
314  << ", nRpi = " << nRpi << ", nRnuc = " << nRnuc << ", NR = " << NR
315  << ", R0 = " << R0 << " fm";
316 
317  TLorentzVector x4_curr(x4); // current position
318 
319  while(1) {
320  double rnow = x4_curr.Vect().Mag();
321  if (rnow > (R+step)) break;
322 
323  x4_curr += (step*dr4);
324  rnow = x4_curr.Vect().Mag();
325  double mfp =
326  genie::utils::intranuke2018::MeanFreePath(pdgc,x4_curr,p4,A,Z,nRpi,nRnuc);
327  double mfp_twk = mfp * mfp_scale_factor;
328 
329  double dprob = (mfp_twk>0) ? TMath::Exp(-step/mfp_twk) : 0.;
330  prob*=dprob;
331 /*
332  LOG("INukeUtils", pDEBUG)
333  << "+ step size = " << step << " fm, |r| = " << rnow << " fm, "
334  << "mfp = " << mfp_twk << "fm (nominal mfp = " << mfp << " fm): "
335  << "dPsurv = " << dprob << ", current Psurv = " << prob;
336 */
337  }
338 
339  LOG("INukeUtils", pDEBUG) << "Psurv = " << prob;
340 
341  return prob;
342 }
343 //____________________________________________________________________________
345  const TLorentzVector & x4, const TLorentzVector & p4,
346  double A, double NR, double R0)
347 {
348 // Calculate distance within a nucleus (units: fm) before we stop tracking
349 // the hadron.
350 // See previous functions for a description of inputs.
351 //
352  double R = NR * R0 * TMath::Power(A, 1./3.);
353  double step = 0.05; // fermi
354 
355  TVector3 dr3 = p4.Vect().Unit(); // unit vector along its direction
356  TLorentzVector dr4(dr3,0);
357 
358  TLorentzVector x4_curr(x4); // current position
359 
360  double d=0;
361  while(1) {
362  double rnow = x4_curr.Vect().Mag();
363  x4_curr += (step*dr4);
364  d+=step;
365  rnow = x4_curr.Vect().Mag();
366  if (rnow > R) break;
367  }
368  return d;
369 }
370 //____________________________________________________________________________
372  int pdgc, const TLorentzVector & x4, const TLorentzVector & p4,
373  double A, double Z, double NR, double R0)
374 {
375 // Calculate distance within a nucleus (expressed in terms of 'mean free
376 // paths') before we stop tracking the hadron.
377 // See previous functions for a description of inputs.
378 //
379 
380 // distance before exiting in mean free path lengths
381 //
382  double R = NR * R0 * TMath::Power(A, 1./3.);
383  double step = 0.05; // fermi
384 
385  TVector3 dr3 = p4.Vect().Unit(); // unit vector along its direction
386  TLorentzVector dr4(dr3,0);
387 
388  TLorentzVector x4_curr(x4); // current position
389 
390  double d=0;
391  double d_mfp=0;
392  while(1) {
393  double rnow = x4_curr.Vect().Mag();
394  x4_curr += (step*dr4);
395  d+=step;
396  rnow = x4_curr.Vect().Mag();
397 
398  double lambda = genie::utils::intranuke2018::MeanFreePath(pdgc,x4_curr,p4,A,Z);
399  d_mfp += (step/lambda);
400 
401  if (rnow > R) break;
402  }
403  return d_mfp;
404 }
405 //____________________________________________________________________________
407  GHepParticle * p, double step, double nuclear_radius)
408 {
409 // Steps a particle starting from its current position (in fm) and moving
410 // along the direction of its current momentum by the input step (in fm).
411 // The particle is stepped in a straight line.
412 // If a nuclear radius is set then the following check is performed:
413 // If the step is too large and takes the the particle far away from the
414 // nucleus then its position is scaled back so that the escaped particles are
415 // always within a ~1fm from the "outer nucleus surface"
416 
417 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
418  LOG("INukeUtils", pDEBUG)
419  << "Stepping particle [" << p->Name() << "] by dr = " << step << " fm";
420 #endif
421 
422  // Step particle
423  TVector3 dr = p->P4()->Vect().Unit(); // unit vector along its direction
424  dr.SetMag(step); // spatial step size
425  double dt = 0; // temporal step:
426  TLorentzVector dx4(dr,dt); // 4-vector step
427  TLorentzVector x4new = *(p->X4()) + dx4; // new position
428 
429  if(nuclear_radius > 0.) {
430  // Check position against nuclear boundary. If the particle was stepped
431  // too far away outside the nuclear boundary bring it back to within
432  // 1fm from that boundary
433  double epsilon = 1; // fm
434  double r = x4new.Vect().Mag(); // fm
435  double rmax = nuclear_radius+epsilon;
436  if(r > rmax) {
437  LOG("INukeUtils", pINFO)
438  << "Particle was stepped too far away (r = " << r << " fm)";
439  LOG("INukeUtils", pINFO)
440  << "Placing it " << epsilon
441  << " fm outside the nucleus (r' = " << rmax << " fm)";
442  double scale = rmax/r;
443  x4new *= scale;
444  }//r>rmax
445  }//nucl radius set
446 
447 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
448  LOG("INukeUtils", pDEBUG)
449  << "\n Init direction = " << print::Vec3AsString(&dr)
450  << "\n Init position (in fm,nsec) = " << print::X4AsString(p->X4())
451  << "\n Fin position (in fm,nsec) = " << print::X4AsString(&x4new);
452 #endif
453 
454  p->SetPosition(x4new);
455 }
456 
457 
458 //___________________________________________________________________________
459 // Method to handle compound nucleus considerations, preequilibrium
460 // and equilibrium
461 // Alex Bell -- 6/17/2008
463  GHepRecord * ev, GHepParticle * p,
464  int &RemnA, int &RemnZ, TLorentzVector &RemnP4,
465  bool /* DoFermi */, double /* FermiFac */,
466  const NuclearModelI* /* Nuclmodel */, double NucRmvE, EINukeMode mode)
467 {
468 
469 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
470  LOG("INukeUtils", pDEBUG)
471  << "PreEquilibrium() is invoked for a : " << p->Name()
472  << " whose kinetic energy is : " << p->KinE();
473 #endif
474 
475  // Random number generator
477  //unused PDGLibrary * pLib = PDGLibrary::Instance();
478 
479  bool allow_dup = true;
480  PDGCodeList list(allow_dup); // list of final state particles
481 
482  double ppcnt = (double) RemnZ / (double) RemnA; // % of protons left
483 
484  // figure out the final state conditions
485 
486  if(p->Pdg()==kPdgProton) list.push_back(kPdgProton);
487  else if(p->Pdg()==kPdgNeutron) list.push_back(kPdgNeutron);
488 
489  for(int i=0;i<3;i++)
490  {
491  if(rnd->RndFsi().Rndm()<ppcnt)
492  {
493  list.push_back(kPdgProton);
494  RemnZ--;
495  }
496  else list.push_back(kPdgNeutron);
497 
498  RemnA--;
499 
500  ppcnt = (double) RemnZ / (double) RemnA;
501  }
502 
503  // Add the fermi energy of the three nucleons to the phase space
504  /*
505  if(DoFermi)
506  {
507  Target target(ev->TargetNucleus()->Pdg());
508  TVector3 pBuf = p->P4()->Vect();
509  double mBuf = p->Mass();
510  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
511  TLorentzVector tSum(pBuf,eBuf);
512  double mSum = 0.0;
513  vector<int>::const_iterator pdg_iter;
514  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
515  {
516  target.SetHitNucPdg(*pdg_iter);
517  Nuclmodel->GenerateNucleon(target);
518  mBuf = pLib->Find(*pdg_iter)->Mass();
519  mSum += mBuf;
520  pBuf = FermiFac * Nuclmodel->Momentum3();
521  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
522  tSum += TLorentzVector(pBuf,eBuf);
523  RemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
524  }
525  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
526  p->SetMomentum(dP4);
527  }
528  */
529  // do the phase space decay & save all f/s particles to the event record
530  bool success = genie::utils::intranuke2018::PhaseSpaceDecay(ev,p,list,RemnP4,NucRmvE,mode);
531  if(success) LOG("INukeUtils2018",pINFO) << "Successful phase space decay for pre-equilibrium nucleus FSI event";
532  else
533  {
535  exception.SetReason("Phase space generation of pre-equilibrium nucleus final state failed, details above");
536  throw exception;
537  }
538 
539  int p_loc = 0;
540  while(p_loc<ev->GetEntries())
541  {
542  GHepParticle * p_ref = ev->Particle(p_loc);
543  if(!p->ComparePdgCodes(p_ref)) p_loc++;
544  else
545  {
546  if(!p->CompareStatusCodes(p_ref)) p_loc++;
547  else
548  {
549  if(!p->CompareMomentum(p_ref)) p_loc++;
550  else break;
551  }
552  }
553  }
554 
555 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
556  LOG("INukeUtils", pDEBUG)
557  << "Particle at: " << p_loc;
558 #endif
559 
560  // find the appropriate daughter
561  vector<int> * descendants = ev->GetStableDescendants(p_loc);
562 
563  int loc = p_loc + 1;
564  int f_loc = p_loc + 1;
565  double energy = ev->Particle(loc)->E();
566 
567 /* // (1) least energetic
568  double min_en = energy;
569 
570  for(unsigned int j=0;j<descendants->size();j++)
571  {
572  loc = (*descendants)[j];
573  energy = ev->Particle(loc)->E();
574  if(energy<min_en)
575  {
576  f_loc = loc;
577  min_en = energy;
578  }
579  }
580 */
581  // (2) most energetic
582  double max_en = energy;
583 
584  for(unsigned int j=0;j<descendants->size();j++)
585  {
586  loc = (*descendants)[j];
587  energy = ev->Particle(loc)->E();
588  if(energy>max_en)
589  {
590  f_loc = loc;
591  max_en = energy;
592  }
593  }
594 
595  // (3) 1st particle
596  // ...just use the defaulted f_loc
597 
598  delete descendants;
599 
600  // change particle status for decaying particle - take out as test
601  //ev->Particle(f_loc)->SetStatus(kIStIntermediateState);
602  // decay a clone particle
603  GHepParticle * t = new GHepParticle(*(ev->Particle(f_loc)));
604  t->SetFirstMother(f_loc);
605  //next statement was in Alex Bell's original code - PreEq, then Equilibrium using particle with highest energy. Note it gets IST=kIStIntermediateState.
606  //genie::utils::intranuke2018::Equilibrium(ev,t,RemnA,RemnZ,RemnP4,DoFermi,FermiFac,Nuclmodel,NucRmvE,mode);
607 
608  delete t;
609 }
610 //___________________________________________________________________________
611 // Method to handle Equilibrium reaction
612 // Alex Bell -- 6/17/2008
614  GHepRecord * ev, GHepParticle * p,
615  int &RemnA, int &RemnZ, TLorentzVector &RemnP4,
616  bool /* DoFermi */, double /* FermiFac */,
617  const NuclearModelI* /* Nuclmodel */, double NucRmvE, EINukeMode mode)
618 {
619 
620 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
621  LOG("INukeUtils", pDEBUG)
622  << "Equilibrium() is invoked for a : " << p->Name()
623  << " whose kinetic energy is : " << p->KinE();
624 #endif
625 
626  // Random number generator
628  //usused PDGLibrary * pLib = PDGLibrary::Instance();
629 
630  bool allow_dup = true;
631  PDGCodeList list(allow_dup); // list of final state particles
632 
633  // % of protons left
634  double ppcnt = (double) RemnZ / (double) RemnA;
635 
636  // figure out the final state conditions
637 
638  if(p->Pdg()==kPdgProton) list.push_back(kPdgProton);
639  else if(p->Pdg()==kPdgNeutron) list.push_back(kPdgNeutron);
640 
641  //add additonal particles to stack
642  for(int i=0;i<4;i++)
643  {
644  if(rnd->RndFsi().Rndm()<ppcnt)
645  {
646  list.push_back(kPdgProton);
647  RemnZ--;
648  }
649  else list.push_back(kPdgNeutron);
650 
651  RemnA--;
652 
653  ppcnt = (double) RemnZ / (double) RemnA;
654  }
655 
656 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
657  LOG("INukeUtils2018", pDEBUG)
658  << "Remnant nucleus (A,Z) = (" << RemnA << ", " << RemnZ << ")";
659 #endif
660 
661  // Add the fermi energy of the three nucleons to the phase space
662  /* if(DoFermi)
663  {
664  Target target(ev->TargetNucleus()->Pdg());
665  TVector3 pBuf = p->P4()->Vect();
666  double mBuf = p->Mass();
667  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
668  TLorentzVector tSum(pBuf,eBuf);
669  double mSum = 0.0;
670  vector<int>::const_iterator pdg_iter;
671  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
672  {
673  target.SetHitNucPdg(*pdg_iter);
674  Nuclmodel->GenerateNucleon(target);
675  mBuf = pLib->Find(*pdg_iter)->Mass();
676  mSum += mBuf;
677  pBuf = FermiFac * Nuclmodel->Momentum3();
678  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
679  tSum += TLorentzVector(pBuf,eBuf);
680  RemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
681  }
682  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
683  p->SetMomentum(dP4);
684  }
685  */
686  // do the phase space decay & save all f/s particles to the record
687  bool success = genie::utils::intranuke2018::PhaseSpaceDecay(ev,p,list,RemnP4,NucRmvE,mode);
688  if (success) LOG("INukeUtils",pINFO) << "successful equilibrium interaction";
689  else
690  {
692  exception.SetReason("Phase space generation of compound nucleus final state failed, details above");
693  throw exception;
694  }
695 
696 }
697 
698 
699 //___________________________________________________________________________
701  GHepRecord* ev, int pcode, int tcode, int scode, int s2code, double C3CM,
702  GHepParticle* p, GHepParticle* t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode)
703 {
704  // Aaron Meyer (10/29/09)
705  // Adapted from kinematics in other function calls
706  //
707  // C3CM is the cosine of the scattering angle, calculated before calling
708  // p and t are the output particles, must be allocated before calling
709  // pcode,tcode,scode,s2code are initial and final particle PDG codes in scattering
710  // return value used for error checking
711 
712  // Kinematic variables
713 
714  double M1, /* M2, */ M3, M4; // rest energies, in GeV
715  double E3L, P3L, E4L, P4L;
716  TVector3 tP1L, tPtot, tbeta, tbetadir, tTrans, tVect;
717  TVector3 tP1zCM, tP2zCM, tP3L, tP4L;
718 
719  // Library instance for reference
720  PDGLibrary * pLib = PDGLibrary::Instance();
721 
722  // random number generator
723  //RandomGen * rnd = RandomGen::Instance();
724 
725  // handle fermi target
726  Target target(ev->TargetNucleus()->Pdg());
727 
728  // get mass for particles
729  M1 = pLib->Find(pcode)->Mass();
730  // usused // M2 = pLib->Find(tcode)->Mass();
731  M3 = pLib->Find(scode)->Mass();
732  M4 = pLib->Find(s2code)->Mass();
733 
734  // get lab energy and momenta and assign to 4 vectors
735  TLorentzVector t4P1L = *p->P4();
736  TLorentzVector t4P2L = *t->P4();
737 
738  // binding energy
739  double bindE = 0.025; // empirical choice, might need to be improved
740  //double bindE = 0.0;
741 
742  LOG("TwoBodyCollision",pNOTICE) << "M1 = " << t4P1L.M() << " , M2 = " << t4P2L.M();
743  LOG("TwoBodyCollision",pNOTICE) << "E1 = " << t4P1L.E() << " , E2 = " << t4P2L.E();
744 
745  if ( (p->Energy()-p->Mass()) < bindE ){bindE = 0.;} // if the probe's energy is less than the binding energy, set the binding energy to 0.
746 
747  // don't use BE unless kinetic energy >> BE.
748  if((pcode==2112||pcode==2212)&&(t4P1L.E()-M1)<.1) bindE = 0.0;
749  if((pcode==211||pcode==-211||pcode==111)&&(t4P1L.E()-M1)<.08) bindE = 0.0;
750  if((pcode==321)&&(t4P1L.E()-M1)<.1) bindE = 0.0;
751 
752  // carry out scattering
753  TLorentzVector t4P3L, t4P4L;
754  if (!TwoBodyKinematics(M3,M4,t4P1L,t4P2L,t4P3L,t4P4L,C3CM,RemnP4,bindE))
755  {
756  P3L = t4P3L.Vect().Mag();
757  P4L = t4P4L.Vect().Mag();
758  E3L = t4P3L.E();
759  E4L = t4P4L.E();
760 
761  LOG("TwoBodyCollision",pNOTICE)
762  << "TwoBodyKinematics fails: C3CM, P3 = " << C3CM << " "
763  << P3L << " " << E3L << "\n" << " P4 = "
764  << P4L << " " << E4L ;
765  return false; //covers all possiblities for now
766  }
767 
768  // error checking
769  P3L = t4P3L.Vect().Mag();
770  P4L = t4P4L.Vect().Mag();
771  E3L = t4P3L.E();
772  E4L = t4P4L.E();
773  LOG("INukeUtils",pINFO)
774  << "C3CM, P3 = " << C3CM << " "
775  << P3L << " " << E3L << "\n" << " P4 = "
776  << P4L << " " << E4L ;
777 
778  // handle very low momentum particles
779  if(!(TMath::Finite(P3L)) || P3L<.001)
780  {
781  LOG("INukeUtils",pINFO)
782  << "Particle 3 momentum small or non-finite: " << P3L
783  << "\n" << "--> Assigning .001 as new momentum";
784  P3L = .001;
785  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
786  }
787  if(!(TMath::Finite(P4L)) || P4L<.001)
788  {
789  LOG("INukeUtils",pINFO)
790  << "Particle 4 momentum small or non-finite: " << P4L
791  << "\n" << "--> Assigning .001 as new momentum";
792  P4L = .001;
793  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
794  }
795 
796  // if this is going to be on in the future, remember to not apply PB for Oset
797  /*// pauli blocking turn off for now to better match data
798  // if(P3L<fFermiMomentum && pdg::IsNeutronOrProton(scode) ||
799  // P4L<fFermiMomentum && pdg::IsNeutronOrProton(s2code) )
800  if(P3L<.25 && pdg::IsNeutronOrProton(scode) ||
801  P4L<.25 && pdg::IsNeutronOrProton(s2code) )
802  {
803  LOG("INukeUtils",pNOTICE)<< "TwoBodyCollision failed: Pauli blocking";
804  p->SetStatus(kIStHadronInTheNucleus);
805  RemnP4 -= TLorentzVector(0,0,0,bindE);
806  return false;
807  }*/
808 
809  // update remnant nucleus
810  RemnP4 -= t4P2L;
811  LOG("INukeUtils",pINFO)
812  << "t4P2L= " << t4P2L.E() << " " << t4P2L.Z()
813  << " RemnP4= " << RemnP4.E() << " " << RemnP4.Z() ;
814  if (tcode==kPdgProton) {RemnZ--;RemnA--;}
815  else if(tcode==kPdgNeutron) RemnA--;
816 
817  // create t particle w/ appropriate momenta, code, and status
818  // Set target's mom to be the mom of the hadron that was cloned
819  t->SetFirstMother(p->FirstMother());
820  t->SetLastMother(p->LastMother());
821 
822  // adjust p to reflect scattering
823  p->SetPdgCode(scode);
824  p->SetMomentum(t4P3L);
825 
826  t->SetPdgCode(s2code);
827  t->SetMomentum(t4P4L);
828 
829  if (mode==kIMdHN)
830  {
833  }
834  else
835  {
838  }
839  LOG("INukeUtils",pINFO) << "Successful 2 body collision";
840  return true;
841 
842 }
843 //___________________________________________________________________________
845  double M3, double M4, TLorentzVector t4P1L, TLorentzVector t4P2L,
846  TLorentzVector &t4P3L, TLorentzVector &t4P4L, double C3CM, TLorentzVector &RemnP4, double bindE)
847 {
848  // Aaron Meyer (05/17/10)
849  // Adapted from kinematics in other function calls
850  //
851  // Outgoing particle masses M3,M4
852  // Scatters particles according to normal two body collisions
853  //
854  // bindE is the binding energy (GeV) of a particle removed from the nucleus (default 0)
855  // For nonzero binding energy, remove the binding energy from the total energy,
856  // then put both of the particles back on mass shell by shifting momentum/energy
857  // of target
858  // Momentum only shifted in the direction parallel to the probe's motion
859  //
860  // Rotates final transverse component of momentum by a random angle from 0->2pi
861  // Return value for error checking
862  // Gives outgoing 4-momenta of particles 3 and 4 (t4P3L, t4P4L respectively)
863  //
864  // All 4-momenta should be on mass shell
865 
866  double E1L, E2L, P1L, P2L, E3L, P3L;
867  double beta, gm; // speed and gamma for CM frame in Lab
868  double S3CM; // sin of scattering angle
869  double PHI3;
870  double E1CM, E2CM, E3CM, P3CM;//, E4CM, P4CM;
871  double P3zL, P3tL;//, P4zL, P4tL;
872  double Et;
873  double theta1, theta2, theta5, P1zL, P2zL, P1tL, P2tL;
874  TVector3 tbeta, tbetadir, tTrans, tVect;
875  TVector3 tP1zCM, tP2zCM, vP3L;
876  TLorentzVector t4P1buf, t4P2buf, t4Ptot;
877 
878  // Library instance for reference
879  //PDGLibrary * pLib = PDGLibrary::Instance();
880 
881  // random number generator
883 
884  // error checking
885  if (C3CM < -1. || C3CM > 1.) return false;
886 
887  // calculate sine from scattering angle
888  S3CM = TMath::Sqrt(1.0 - C3CM*C3CM);
889 
890  // fill buffers
891  t4P1buf = t4P1L;
892  t4P2buf = t4P2L;
893 
894  // get lab energy and momenta
895  E1L = t4P1buf.E();
896  P1L = t4P1buf.P();
897  E2L = t4P2buf.E();
898  P2L = t4P2buf.P();
899  t4Ptot = t4P1buf + t4P2buf;
900 
901  LOG("INukeUtils",pINFO) <<"M1 "<<t4P1L.M()<< ", M2 "<<t4P2L.M();
902  LOG("INukeUtils",pINFO) <<"E1L "<<E1L<< ", E1CM "<<E1CM;
903  LOG("INukeUtils",pINFO) <<"bindE = " << bindE;
904 
905  // binding energy
906  if (bindE!=0)
907  {
908 
909  E1L -= bindE;
910 
911  if (E1L+E2L < M3+M4)
912  {
913  LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Forbidden by binding energy";
914  LOG("INukeUtils",pNOTICE) <<"E1L, E2L, M3, M4 : "<< E1L <<", "<< E2L <<", "<< M3 <<", "<< M4;
915  t4P3L.SetPxPyPzE(0,0,0,0);
916  t4P4L.SetPxPyPzE(0,0,0,0);
917  return false;
918  }
919  }
920 
921  // calculate beta and gamma
922  tbeta = t4Ptot.Vect() * (1.0 / (E1L + E2L));
923  tbetadir = tbeta.Unit();
924  beta = tbeta.Mag();
925  gm = 1.0 / TMath::Sqrt(1.0 - beta*beta);
926 
927  // get angle and component info
928  theta1 = t4P1buf.Angle(tbeta);
929  theta2 = t4P2buf.Angle(tbeta);
930  P1zL = P1L*TMath::Cos(theta1);
931  P2zL = P2L*TMath::Cos(theta2);
932  P1tL = TMath::Sqrt(P1L*P1L - P1zL*P1zL);
933  P2tL = -TMath::Sqrt(P2L*P2L - P2zL*P2zL);
934  tVect.SetXYZ(1,0,0);
935  if(TMath::Abs((tVect - tbetadir).Mag())<.01) tVect.SetXYZ(0,1,0);
936  theta5 = tVect.Angle(tbetadir);
937  tTrans = (tVect - TMath::Cos(theta5)*tbetadir).Unit();
938 
939  // boost to CM frame to get scattered particle momenta
940  E1CM = gm*E1L - gm*beta*P1zL;
941  tP1zCM = gm*P1zL*tbetadir - gm*tbeta*E1L;
942  E2CM = gm*E2L - gm*beta*P2zL;
943  tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
944  Et = E1CM + E2CM;
945 //-------
946 
947  LOG("INukeUtils",pINFO) <<"M1 "<<t4P1L.M()<< ", M2 "<<t4P2L.M();
948  LOG("INukeUtils",pINFO) <<"E1L "<<E1L<< ", E1CM "<<E1CM;
949  LOG("INukeUtils",pINFO) <<"P1zL "<<P1zL<<", P1zCM "<<tP1zCM.Mag()<<", P1tL "<<P1tL;
950  LOG("INukeUtils",pINFO) <<"E2L "<<E2L<< ", E2CM "<<E2CM;
951  LOG("INukeUtils",pINFO) <<"P2zL "<<P2zL<<", P2zCM "<<tP2zCM.Mag()<<", P2tL "<<P2tL;
952  LOG("INukeUtils",pINFO) <<"C3CM "<<C3CM;
953 
954 //-------
955  E3CM = (Et*Et + M3*M3 - M4*M4) / (2.0 * Et);
956 
957  // check to see if decay is viable
958  if(E3CM*E3CM - M3*M3<0 || E3CM<0 || Et<0)
959  {
960  if (Et<0) LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Total energy is negative";
961  if (E3CM<M3) LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Scattered Particle 3 CM energy is too small";
962  if (E3CM*E3CM - M3*M3<0) LOG("INukeUtils",pNOTICE) <<"TwoBodyKinematics Failed: Scattered Particle 3 CM momentum is nonreal";
963  t4P3L.SetPxPyPzE(0,0,0,0);
964  t4P4L.SetPxPyPzE(0,0,0,0);
965  return false;
966  }
967  P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
968 
969  // boost back to lab
970  P3zL = gm*beta*E3CM + gm*P3CM*C3CM;
971  P3tL = P3CM*S3CM;
972 
973  P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
974  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
975 
976  //-------
977 
978  double E4CM = Et-E3CM;
979  double P4zL = gm*beta*E4CM - gm*P3CM*C3CM;
980  double P4tL = -1.*P3tL;
981  double P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
982  double E4L = TMath::Sqrt(P4L*P4L + M4*M4);
983 
984  LOG("INukeUtils",pINFO) <<"M3 "<< M3 << ", M4 "<< M4;
985  LOG("INukeUtils",pINFO) <<"E3L "<<E3L<< ", E3CM "<<E3CM;
986  LOG("INukeUtils",pINFO) <<"P3zL "<<P3zL<<", P3tL "<<P3tL;
987  LOG("INukeUtils",pINFO) <<"C3L "<<P3zL/P3L;
988  LOG("INukeUtils",pINFO) <<"Check:";
989  LOG("INukeUtils",pINFO) <<"E4L "<<E4L<< ", E4CM "<<E4CM;
990  LOG("INukeUtils",pINFO) <<"P4zL "<<P4zL<<", P4tL "<<P4tL;
991  LOG("INukeUtils",pINFO) <<"P4L "<<P4L;
992  LOG("INukeUtils",pINFO) <<"C4L "<<P4zL/P4L;
993 
994  double echeck = E1L + E2L - (E3L + E4L);
995  double pzcheck = P1zL+ P2zL - (P3zL + P4zL);
996  double ptcheck = P1tL+ P2tL - (P3tL + P4tL);
997 
998  LOG("INukeUtils",pINFO) <<"Check 4-momentum conservation - Energy "<<echeck<<", z momentum "<<pzcheck << ", transverse momentum " << ptcheck ;
999 
1000  // -------
1001 
1002  // handle very low momentum particles
1003  if(!(TMath::Finite(P3L)) || P3L<.001)
1004  {
1005  LOG("INukeUtils",pINFO)
1006  << "Particle 3 momentum small or non-finite: " << P3L
1007  << "\n" << "--> Assigning .001 as new momentum";
1008  P3tL = 0;
1009  P3zL = .001;
1010  P3L = .001;
1011  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1012  }
1013 
1014  // get random phi angle, distributed uniformally in 360 deg
1015  PHI3 = 2 * kPi * rnd->RndFsi().Rndm();
1016 
1017  vP3L = P3zL*tbetadir + P3tL*tTrans;
1018  vP3L.Rotate(PHI3,tbetadir);
1019 
1020  t4P3L.SetVect(vP3L);
1021  t4P3L.SetE(E3L);
1022 
1023  t4P4L = t4P1buf + t4P2buf - t4P3L;
1024  t4P4L-= TLorentzVector(0,0,0,bindE);
1025  /*LOG("INukeUtils",pINFO) <<"GENIE:";
1026  LOG("INukeUtils",pINFO) <<"E4L "<<t4P4L.E();
1027  LOG("INukeUtils",pINFO) <<"P4zL "<<t4P4L.Vect()*tbetadir<<", P4tL "<<-1.*TMath::Sqrt(t4P4L.Vect().Mag2()-TMath::Power(t4P4L.Vect()*tbetadir,2.));
1028  LOG("INukeUtils",pINFO) <<"P4L "<<t4P4L.Vect().Mag();
1029  LOG("INukeUtils",pINFO) <<"C4L "<<t4P4L.Vect()*tbetadir/t4P4L.Vect().Mag(); */
1030 
1031  if(t4P4L.Mag2()<0 || t4P4L.E()<0)
1032  {
1033  LOG("INukeUtils",pNOTICE)<<"TwoBodyKinematics Failed: Target mass or energy is negative";
1034  t4P3L.SetPxPyPzE(0,0,0,0);
1035  t4P4L.SetPxPyPzE(0,0,0,0);
1036  return false;
1037  }
1038 
1039  if (bindE!=0) RemnP4 += TLorentzVector(0,0,0,bindE);
1040  return true;
1041 }
1042 //___________________________________________________________________________
1044  GHepRecord* ev, GHepParticle* p, int tcode, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3,
1045  bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
1046 {
1047 
1048  // Aaron Meyer (7/15/10)
1049  //
1050  // Kinematics used in utils::intranuke2018::PionProduction
1051  // Handles the kinematics of three body scattering
1052  //
1053  // s1,s2,and s3 are pointers to particles with the PDG code that needs to be scattered
1054  // the last four variables are for Fermi momentum and pauli blocking, default will use none of them
1055 
1056  // kinematic variables
1057  double M1, M2, M3, M4, M5; // rest energies, in GeV
1058  double P1L, P2L, P3L, P4L, P5L;
1059  double E1L, E2L, E3L, E4L, E5L;
1060  double E1CM, E2CM, P3tL;
1061  double PizL, PitL, PiL, EiL;
1062  double EiCM, P4CM2, E4CM2, E5CM2, P3CM, E3CM;
1063  double beta, gm, beta2, gm2;
1064  double P3zL, P4zL, P4tL, P5zL, P5tL;
1065  double Et, M, theta1, theta2;
1066  double P1zL, P2zL;
1067  double theta3, theta4, phi3, phi4, theta5;
1068  TVector3 tP2L, tP1L, tPtot, tbeta, tbetadir, tTrans, tP4L, tP5L;
1069  TVector3 tP1zCM, tP2zCM, tP3L, tPiL, tbeta2, tbetadir2, tVect, tTrans2;
1070 
1071  // Library instance for reference
1072  PDGLibrary * pLib = PDGLibrary::Instance();
1073 
1074  // random number generator
1076 
1077  M1 = pLib->Find(p->Pdg())->Mass();
1078  M2 = pLib->Find(tcode)->Mass();
1079  M3 = pLib->Find(s1->Pdg())->Mass();
1080  M4 = pLib->Find(s2->Pdg())->Mass();
1081  M5 = pLib->Find(s3->Pdg())->Mass();
1082 
1083  // set up fermi target
1084  Target target(ev->TargetNucleus()->Pdg());
1085 
1086  // handle fermi momentum
1087  if(DoFermi)
1088  {
1089  target.SetHitNucPdg(tcode);
1090  Nuclmodel->GenerateNucleon(target);
1091  tP2L = FermiFac * Nuclmodel->Momentum3();
1092  P2L = tP2L.Mag();
1093  E2L = TMath::Sqrt(tP2L.Mag2() + M2*M2);
1094  }
1095  else
1096  {
1097  tP2L.SetXYZ(0.0, 0.0, 0.0);
1098  P2L = 0;
1099  E2L = M2;
1100  }
1101 
1102  // first sequence, handle 4th and 5th products as composite
1103  E1L = p->E();
1104 
1105  P1L = TMath::Sqrt(E1L*E1L - M1*M1);
1106  tP1L = p->P4()->Vect();
1107  tPtot = tP1L + tP2L;
1108 
1109  tbeta = tPtot * (1.0 / (E1L + E2L));
1110  tbetadir = tbeta.Unit();
1111  beta = tbeta.Mag();
1112  gm = 1.0 / TMath::Sqrt(1.0 - beta*beta);
1113 
1114  theta1 = tP1L.Angle(tbeta);
1115  theta2 = tP2L.Angle(tbeta);
1116  P1zL = P1L*TMath::Cos(theta1);
1117  P2zL = P2L*TMath::Cos(theta2);
1118  tVect.SetXYZ(1,0,0);
1119  if(TMath::Abs((tVect - tbetadir).Mag())<.01) tVect.SetXYZ(0,1,0);
1120  theta5 = tVect.Angle(tbetadir);
1121  tTrans = (tVect - TMath::Cos(theta5)*tbetadir).Unit();
1122 
1123  E1CM = gm*E1L - gm*beta*P1zL;
1124  tP1zCM = gm*P1zL*tbetadir - gm*tbeta*E1L;
1125  E2CM = gm*E2L - gm*beta*P2zL;
1126  tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
1127  Et = E1CM + E2CM;
1128  M = (rnd->RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
1129  E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
1130  EiCM = Et - E3CM;
1131  if(E3CM*E3CM - M3*M3<0)
1132  {
1133  LOG("INukeUtils",pNOTICE)
1134  << "PionProduction P3 has non-real momentum - retry kinematics";
1135  LOG("INukeUtils",pNOTICE) << "Energy, masses of 3 fs particales:"
1136  << E3CM << " " << M3 << " " << " " << M4 << " " << M5;
1138  exception.SetReason("PionProduction particle 3 has non-real momentum");
1139  throw exception;
1140  return false;
1141  }
1142  P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
1143 
1144  theta3 = kPi * rnd->RndFsi().Rndm();
1145  theta4 = kPi * rnd->RndFsi().Rndm();
1146  phi3 = 2*kPi * rnd->RndFsi().Rndm();
1147  phi4 = 2*kPi * rnd->RndFsi().Rndm();
1148 
1149  P3zL = gm*beta*E3CM + gm*P3CM*TMath::Cos(theta3);
1150  P3tL = P3CM*TMath::Sin(theta3);
1151  PizL = gm*beta*EiCM - gm*P3CM*TMath::Cos(theta3);
1152  PitL = -P3CM*TMath::Sin(theta3);
1153 
1154  P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
1155  PiL = TMath::Sqrt(PizL*PizL + PitL*PitL);
1156  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1157  EiL = TMath::Sqrt(PiL*PiL + M*M);
1158 
1159  // handle very low momentum particles
1160  if(!(TMath::Finite(P3L)) || P3L < .001)
1161  {
1162  LOG("INukeUtils",pINFO)
1163  << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
1164  << "\n" << "--> Assigning .001 as new momentum";
1165  P3tL = 0;
1166  P3zL = .001;
1167  P3L = .001;
1168  E3L = TMath::Sqrt(P3L*P3L + M3*M3);
1169  }
1170 
1171  tP3L = P3zL*tbetadir + P3tL*tTrans;
1172  tPiL = PizL*tbetadir + PitL*tTrans;
1173  tP3L.Rotate(phi3,tbetadir);
1174  tPiL.Rotate(phi3,tbetadir);
1175 
1176  // second sequence, handle formally composite particles 4 and 5
1177  tbeta2 = tPiL * (1.0 / EiL);
1178  tbetadir2 = tbeta2.Unit();
1179  beta2 = tbeta2.Mag();
1180  gm2 = 1.0 / TMath::Sqrt(1.0 - beta2*beta2);
1181 
1182  E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
1183  E5CM2 = M - E4CM2;
1184  P4CM2 = TMath::Sqrt(E4CM2*E4CM2 - M4*M4);
1185 
1186  tVect.SetXYZ(1,0,0);
1187  if(TMath::Abs((tVect - tbetadir2).Mag())<.01) tVect.SetXYZ(0,1,0);
1188  theta5 = tVect.Angle(tbetadir2);
1189  tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).Unit();
1190 
1191  P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*TMath::Cos(theta4);
1192  P4tL = P4CM2*TMath::Sin(theta4);
1193  P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*TMath::Cos(theta4);
1194  P5tL = - P4tL;
1195 
1196  P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
1197  P5L = TMath::Sqrt(P5zL*P5zL + P5tL*P5tL);
1198  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
1199  E5L = TMath::Sqrt(P5L*P5L + M5*M5);
1200 
1201  // handle very low momentum particles
1202  if(!(TMath::Finite(P4L)) || P4L < .001)
1203  {
1204  LOG("INukeUtils",pINFO)
1205  << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
1206  << "\n" << "--> Assigning .001 as new momentum";
1207  P4tL = 0;
1208  P4zL = .001;
1209  P4L = .001;
1210  E4L = TMath::Sqrt(P4L*P4L + M4*M4);
1211  }
1212  if(!(TMath::Finite(P5L)) || P5L < .001)
1213  {
1214  LOG("INukeUtils",pINFO)
1215  << "Particle 5 " << M5 << " momentum small or non-finite: " << P5L
1216  << "\n" << "--> Assigning .001 as new momentum";
1217  P5tL = 0;
1218  P5zL = .001;
1219  P5L = .001;
1220  E5L = TMath::Sqrt(P5L*P5L + M5*M5);
1221  }
1222 
1223  tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
1224  tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
1225  tP4L.Rotate(phi4,tbetadir2);
1226  tP5L.Rotate(phi4,tbetadir2);
1227 
1228  // pauli blocking
1229  if(P3L < FermiMomentum || ( pdg::IsNeutronOrProton(s2->Pdg()) && P4L < FermiMomentum ) )
1230  {
1231  LOG("INukeUtils",pNOTICE)
1232  << "PionProduction fails because of Pauli blocking - retry kinematics";
1234  exception.SetReason("PionProduction final state not determined");
1235  throw exception;
1236  return false;
1237  }
1238 
1239  // create scattered particles w/ appropriate momenta, code, and status
1240  // Set moms to be the moms of the hadron that was cloned
1241  s1->SetFirstMother(p->FirstMother());
1242  s2->SetFirstMother(p->FirstMother());
1243  s3->SetFirstMother(p->FirstMother());
1244  s1->SetLastMother(p->LastMother());
1245  s2->SetLastMother(p->LastMother());
1246  s3->SetLastMother(p->LastMother());
1247 
1248  TLorentzVector(tP3L,E3L);
1249  TLorentzVector(tP4L,E4L);
1250  TLorentzVector(tP5L,E5L);
1251 
1252  s1->SetMomentum(TLorentzVector(tP3L,E3L));
1253  s2->SetMomentum(TLorentzVector(tP4L,E4L));
1254  s3->SetMomentum(TLorentzVector(tP5L,E5L));
1255  int mode = kIMdHA;
1256  LOG ("INukeUtils",pDEBUG) << "in Pi Prod, mode = " << mode;
1257  if (mode==kIMdHN)
1258  {
1262  }
1263  else
1264  {
1268  }
1269  return true;
1270 }
1271 //___________________________________________________________________________
1273  GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ,
1274  TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
1275 {
1276  // Aaron Meyer (7/15/2010)
1277  //
1278  // Handles pion production reactions in both hA and hN mode
1279  // Calculates fundamental cross sections from fit functions
1280  // Uses isospin relations to determine the rest of cross sections
1281  //
1282  // p is the probe particle
1283  // s1, s2, and s3 are the particles produced in the reaction
1284  // must set the status and add particles to the event record after returning from this method
1285  // return value for error checking
1286 
1287 
1288  // random number generator
1290 
1291  // library reference
1292  PDGLibrary * pLib = PDGLibrary::Instance();
1293 
1294  bool ptarg = false;
1295  int pcode = p->Pdg();
1296 
1297  int p1code = p->Pdg();
1298  // Randomly determine target and 1st product baryons
1299  int p3code = 0, p4code = 0, p5code = 0;
1300 
1301  //
1302  // Uses a fit curve log(sigma) = a - b/(T_pi - c) for pions
1303  // Fit parameters determined by Roman Tacik (4/3/09)
1304  // pi- & p cross sections are assumed to be the same as pi+ & n
1305  //
1306  // Fit curve for nucleons:
1307  // sigma = a*(1+b*e^(-c*(eta-d)^2))*(1-e^(-(f*eta)^g))*(1-e^(-h/eta^2))
1308  // 7 parameters (a,b,c,d,f,g,h)
1309  // eta is maximum kinematically allowed momentum of the pion, normalized by the mass
1310  // Uses isotopic spin decomposition of total cross sections
1311  //
1312 
1313  if ((p1code==kPdgPi0)||(p1code==kPdgPiP)||(p1code==kPdgPiM)) {
1314 
1315  double kine = 1000*p->KinE();
1316 
1317  // Determine cross sections
1318 
1319  // pion
1320  // pi- & p
1321  // -> pi0 & pi0 & n
1322  // a = 8.82; b = 573.2; c = 107.3;
1323  double xsec2pi0n = TMath::Max(0.,TMath::Exp(8.82 - (573.2/(kine-107.3))));
1324  // -> pi- & pi+ & n
1325  // a = 11.06; b = 985.9; c = 88.2;
1326  double xsecpippimn = TMath::Max(0.,TMath::Exp(11.06 - (985.9/(kine-88.2))));
1327  // -> pi- & pi0 & p
1328  // a = 9.58; b = 1229.4; c = 60.5;
1329  double xsecpimpi0p = TMath::Max(0.,TMath::Exp(9.58 - (1229.4/(kine-60.5))));
1330  double totpimp = xsec2pi0n + xsecpippimn + xsecpimpi0p;
1331 
1332 
1333  // pi+ & p
1334  // -> pi+ & pi+ & n
1335  // a = 5.64; b = 222.6; c = 150.0;
1336  double xsec2pipn = TMath::Max(0.,TMath::Exp(5.64 - (222.6/(kine-150.))));
1337  // -> pi+ & pi0 & p
1338  // a = 7.95; b = 852.6; c = 77.8;
1339  double xsecpippi0p = TMath::Max(0.,TMath::Exp(7.95 - (852.6/(kine-77.8))));
1340  double totpipp = xsec2pipn + xsecpippi0p;
1341 
1342  if (totpimp<=0 && totpipp<=0) {
1343  LOG("INukeUtils",pNOTICE) << "InelasticHN called below threshold energy";
1345  ev->AddParticle(*p);
1346  return false;
1347  }
1348 
1349  double xsecp, xsecn;
1350  switch (p1code) {
1351  case kPdgPi0: xsecp = 0.5 * (totpimp + totpipp); xsecn = xsecp; break;
1352  case kPdgPiP: xsecp = totpipp; xsecn = totpimp; break;
1353  case kPdgPiM: xsecp = totpimp; xsecn = totpipp; break;
1354  default:
1355  LOG("INukeUtils",pWARN) << "InelasticHN cannot handle probe: "
1356  << PDGLibrary::Instance()->Find(p1code)->GetName();
1358  exception.SetReason("PionProduction final state not determined");
1359  throw exception;
1360  return false;
1361  break;
1362  }
1363 
1364  // Normalize cross sections by Z or A-Z
1365 
1366  xsecp *= RemnZ;
1367  xsecn *= RemnA-RemnZ;
1368 
1369  // determine target
1370 
1371  double rand = rnd->RndFsi().Rndm() * (xsecp + xsecn);
1372  if (rand < xsecp) // proton target
1373  { rand /= RemnZ; ptarg = true;}
1374  else // neutron target
1375  { rand -= xsecp; rand /= RemnA-RemnZ; ptarg = false;}
1376 
1377 
1378  // determine final state
1379 
1380  if (((ptarg==true)&&(p1code==kPdgPiP))
1381  || ((ptarg==false)&&(p1code==kPdgPiM)))
1382  {
1383  if (rand < xsec2pipn) // pi+ & pi+ & n final state
1384  {
1385  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1386  p4code = p1code;
1387  p5code = p4code;
1388  }
1389  else { // pi+ & pi0 & p final state
1390  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1391  p4code = p1code;
1392  p5code = kPdgPi0;
1393  }
1394  }
1395  else if (((ptarg==false)&&(p1code==kPdgPiP))
1396  || ((ptarg==true)&&(p1code==kPdgPiM)))
1397  {
1398  if (rand < xsec2pi0n) // pi0 & pi0 & n final state
1399  {
1400  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1401  p4code = kPdgPi0;
1402  p5code = p4code;
1403  }
1404  else if (rand < (xsec2pi0n + xsecpippimn)) // pi+ & pi- & n final state
1405  {
1406  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1407  p4code = p1code;
1408  p5code = ((p1code==kPdgPiP) ? kPdgPiM : kPdgPiP);
1409  }
1410  else // pi0 & pi- & p final state
1411  {
1412  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1413  p4code = p1code;
1414  p5code = kPdgPi0;
1415  }
1416  }
1417  else if (p1code==kPdgPi0)
1418  {
1419  rand = rnd->RndFsi().Rndm();
1420  if (rand < 191./270.)
1421  { // pi+ & pi- & p final state
1422  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1423  p4code = kPdgPiP;
1424  p5code = kPdgPiM;
1425  }
1426  else if (rand < 7./135.)
1427  { // pi0 & pi0 & p final state
1428  p3code = (ptarg ? kPdgProton : kPdgNeutron);
1429  p4code = kPdgPi0;
1430  p5code = p4code;
1431  }
1432  else
1433  { // pi+ & pi0 & n final state
1434  p3code = (ptarg ? kPdgNeutron : kPdgProton);
1435  p4code = (ptarg ? kPdgPiP : kPdgPiM);
1436  p5code = kPdgPi0;
1437  }
1438  }
1439  else // unhandled
1440  {
1441  LOG("INukeUtils",pNOTICE) << "Pi production final state unable to be determined, picode, ptarg = " <<PDGLibrary::Instance()->Find(p1code)->GetName() << " " << PDGLibrary::Instance()->Find(ptarg)->GetName();
1443  exception.SetReason("PionProduction final state not determined");
1444  throw exception;
1445  return false;
1446  }
1447 
1448  } else if(p1code==kPdgProton||p1code==kPdgNeutron) //nucleon probes
1449  {
1450 
1451  double tote = p->Energy();
1452  double pMass = pLib->Find(2212)->Mass();
1453  double nMass = pLib->Find(2112)->Mass();
1454  double etapp2ppPi0 =
1455  utils::intranuke2018::CalculateEta(pMass,tote,pMass,pMass+pMass,pLib->Find(111)->Mass());
1456  double etapp2pnPip =
1457  utils::intranuke2018::CalculateEta(pLib->Find(p1code)->Mass(),tote,((p1code==kPdgProton)?pMass:nMass),
1458  pMass+nMass,pLib->Find(211)->Mass());
1459  double etapn2nnPip =
1460  utils::intranuke2018::CalculateEta(pMass,tote,nMass,nMass+nMass,pLib->Find(211)->Mass());
1461  double etapn2ppPim =
1462  utils::intranuke2018::CalculateEta(pMass,tote,nMass,pMass+pMass,pLib->Find(211)->Mass());
1463 
1464  if ((etapp2ppPi0<=0.)&&(etapp2pnPip<=0.)&&(etapn2nnPip<=0.)&&(etapn2ppPim<=0.)) { // below threshold
1465  LOG("INukeUtils",pNOTICE) << "PionProduction() called below threshold energy";
1467  exception.SetReason("PionProduction final state not possible - below threshold");
1468  throw exception;
1469  return false;
1470  }
1471 
1472  // calculate cross sections
1473  double xsecppPi0=0,xsecpnPiP=0,xsecnnPiP=0,xsecppPiM=0;
1474  if (etapp2ppPi0>0){
1475  xsecppPi0 = 4511*(1-.91*TMath::Exp(-TMath::Power((etapp2ppPi0-.705),2)));
1476  xsecppPi0 *= (1-TMath::Exp(-TMath::Power((.556*etapp2ppPi0),3.5)));
1477  xsecppPi0 *= (1-TMath::Exp(-56.897/(etapp2ppPi0*etapp2ppPi0)));
1478  xsecppPi0 = TMath::Max(0.,xsecppPi0);}
1479 
1480  if (etapp2pnPip>0){
1481  xsecpnPiP = 18840*(1-.808*TMath::Exp(-TMath::Power((etapp2pnPip-.371),2)));
1482  xsecpnPiP *= (1-TMath::Exp(-TMath::Power((.568*etapp2pnPip),3.2)));
1483  xsecpnPiP *= (1-TMath::Exp(-39.818/(etapp2pnPip*etapp2pnPip)));
1484  xsecpnPiP = TMath::Max(0.,xsecpnPiP);}
1485 
1486  if (etapn2nnPip>0){
1487  xsecnnPiP = 7670*(1-.479*TMath::Exp(-TMath::Power((etapn2nnPip-.947),2)));
1488  xsecnnPiP *= (1-TMath::Exp(-TMath::Power((.35*etapn2nnPip),3.2)));
1489  xsecnnPiP *= (1-TMath::Exp(-71.53/(etapn2nnPip*etapn2nnPip)));
1490  xsecnnPiP = TMath::Max(0.,xsecnnPiP);}
1491 
1492  if (etapn2ppPim>0){
1493  xsecppPiM = 7670*(1-.479*TMath::Exp(-TMath::Power((etapn2ppPim-.947),2)));
1494  xsecppPiM *= (1-TMath::Exp(-TMath::Power((.35*etapn2ppPim),3.2)));
1495  xsecppPiM *= (1-TMath::Exp(-71.53/(etapn2ppPim*etapn2ppPim)));
1496  xsecppPiM = TMath::Max(0.,xsecppPiM);}
1497 
1498  // double sigma11 = xsecppPi0;
1499  double sigma10 = TMath::Max(0.,xsecpnPiP - xsecppPi0); // Fundamental cross sections
1500  double sigma01 = TMath::Max(0.,xsecppPiM + xsecnnPiP - xsecppPi0);
1501 
1502  double xsecpnPi0 = .5*(sigma10 + sigma01);
1503  xsecpnPi0 = TMath::Max(xsecpnPi0,0.);
1504 
1505  LOG("INukeUtils",pDEBUG) << '\n' << "Cross section values: "<<'\n'
1506  << xsecppPi0 << " PP pi0" <<'\n'
1507  << xsecpnPiP << " PN pi+" <<'\n'
1508  << xsecnnPiP << " NN pi+" <<'\n'
1509  << xsecpnPi0 << " PN pi0";
1510 
1511  double xsecp=0,xsecn=0;
1512  switch (p1code) {
1513  case kPdgProton: xsecp=xsecppPi0+xsecpnPiP; xsecn=xsecppPiM+xsecnnPiP+xsecpnPi0; break;
1514  case kPdgNeutron: xsecp=xsecppPiM+xsecnnPiP+xsecpnPi0; xsecn=xsecppPi0+xsecpnPiP; break;
1515  default:
1516  LOG("INukeUtils",pWARN) << "InelasticHN cannot handle probe: "
1517  << PDGLibrary::Instance()->Find(p1code)->GetName();
1518  return false;
1519  break;
1520  }
1521 
1522  // Normalize cross sections by Z or (A-Z)
1523 
1524  xsecp *= RemnZ;
1525  xsecn *= RemnA-RemnZ;
1526 
1527  // determine target
1528 
1529  double rand = rnd->RndFsi().Rndm() * (xsecp + xsecn);
1530  if (rand < xsecp) // proton target
1531  { rand /= RemnZ; ptarg = true;}
1532  else // neutron target
1533  { rand -= xsecp; rand /= RemnA-RemnZ; ptarg = false;}
1534 
1535  if(p1code==kPdgProton) // Cross sections not explicitly given are calculated from isospin relations
1536  {
1537  if(ptarg)
1538  {
1539  if (rand<xsecppPi0) {p3code=kPdgProton; p4code=kPdgProton; p5code=kPdgPi0;}
1540  else {p3code=kPdgProton; p4code=kPdgNeutron; p5code=kPdgPiP;}
1541  }
1542  else
1543  {
1544  if (rand<xsecnnPiP) {p3code=kPdgNeutron; p4code=kPdgNeutron; p5code=kPdgPiP;}
1545  else if (rand<xsecppPiM+xsecnnPiP) {p3code=kPdgProton; p4code=kPdgProton; p5code=kPdgPiM;}
1546  else {p3code=kPdgProton; p4code=kPdgNeutron; p5code=kPdgPi0;}
1547  }
1548  }
1549  else if(p1code==kPdgNeutron)
1550  {
1551  if(ptarg)
1552  {
1553  if (rand<xsecnnPiP) {p3code=kPdgNeutron; p4code=kPdgNeutron; p5code=kPdgPiP;}
1554  else if (rand<xsecppPiM+xsecnnPiP) {p3code=kPdgProton; p4code=kPdgProton; p5code=kPdgPiM;}
1555  else {p3code=kPdgProton; p4code=kPdgNeutron; p5code=kPdgPi0;}
1556  }
1557  else
1558  {
1559  if (rand<xsecpnPiP) {p3code=kPdgNeutron; p4code=kPdgProton; p5code=kPdgPiM;}
1560  else {p3code=kPdgNeutron; p4code=kPdgNeutron; p5code=kPdgPi0;}
1561  }
1562  }
1563  }
1564  else
1565  {
1566  LOG("INukeUtils",pWARN)
1567  << "Unable to handle probe (=" << p1code << ") in InelasticHN()";
1568  return false;
1569  }
1570 
1571  // determine if reaction is allowed
1572  if ( RemnA < 1 )
1573  {
1574  LOG("INukeUtils",pNOTICE) << "PionProduction() failed : no nucleons to produce pions";
1575  return false;
1576  }
1577  else if ( RemnZ + ((pcode==kPdgProton || pcode==kPdgPiP)?1:0) - ((pcode==kPdgPiM)?1:0)
1578  < ((p3code==kPdgProton || p3code==kPdgPiP)?1:0) - ((p3code==kPdgPiM)?1:0)
1579  + ((p4code==kPdgProton || p4code==kPdgPiP)?1:0) - ((p4code==kPdgPiM)?1:0)
1580  + ((p5code==kPdgProton || p5code==kPdgPiP)?1:0) - ((p5code==kPdgPiM)?1:0) )
1581  {
1582  LOG("INukeUtils",pNOTICE) << "PionProduction() failed : too few protons in nucleus";
1584  exception.SetReason("PionProduction fails - too few protons available");
1585  throw exception;
1586  return false;
1587  }
1588 
1589  s1->SetPdgCode(p3code);
1590  s2->SetPdgCode(p4code);
1591  s3->SetPdgCode(p5code);
1592 
1594  ev,p,(ptarg?kPdgProton:kPdgNeutron),s1,s2,s3,DoFermi,FermiFac,FermiMomentum,Nuclmodel))
1595  {
1596  // okay, handle remnants and return true
1597  // assumes first particle is always the nucleon,
1598  // second can be either nucleon or pion
1599  // last always pion
1600  if (pcode==kPdgProton || pcode==kPdgPiP) RemnZ++;
1601  if (pcode==kPdgPiM) RemnZ--;
1602  if (pdg::IsPion(pcode)) RemnA--;
1603  if (pdg::IsProton(p3code)) RemnZ--;
1604  if (pdg::IsNeutronOrProton(p4code)) RemnA--;
1605  if (p4code==kPdgPiP || p4code==kPdgProton) RemnZ--;
1606  if (p4code==kPdgPiM) RemnZ++;
1607  if (p5code==kPdgPiP) RemnZ--;
1608  if (p5code==kPdgPiM) RemnZ++;
1609 
1610  LOG("INukeUtils",pDEBUG) << "Remnant (A,Z) = (" <<RemnA<<','<<RemnZ<<')';
1611 
1612  RemnP4 -= *s1->P4() + *s2->P4() + *s3->P4() - *p->P4();
1613  return true;
1614  }
1615  else {
1617  exception.SetReason("PionProduction final state not determined");
1618  throw exception;
1619  return false;
1620  }
1621 }
1622 //___________________________________________________________________________
1623 double genie::utils::intranuke2018::CalculateEta(double Minc, double nrg, double Mtarg,
1624  double Mtwopart, double Mpi)
1625 {
1626  //Aaron Meyer (1/20/2010)
1627 
1628  //Used to calculate the maximum kinematically allowed CM frame pion momentum
1629  // ke in MeV, eta normalized by pion mass
1630  // approximated by taking two ejected nucleons to be one particle of the same mass
1631  //For pion cross sections, in utils::intranuke2018::PionProduction
1632 
1633  //LOG("INukeUtils",pDEBUG) << "Input values: "<<Minc<<' '<<nrg<<' '<<Mtarg<<' '<<Mtwopart<<' '<<Mpi;
1634  double Scm = Minc*Minc + Mtarg*Mtarg + 2*Mtarg*nrg;
1635  double eta = 0;
1636  eta= TMath::Power(Scm,2) + TMath::Power(Mtwopart,4) + TMath::Power(Mpi,4);
1637  eta-= 2*TMath::Power(Mtwopart*Mpi,2);
1638  eta-= 2*Scm*TMath::Power(Mtwopart,2);
1639  eta-= 2*Scm*TMath::Power(Mpi,2);
1640  eta = TMath::Power(eta/Scm,1./2.);
1641  eta/= (2*Mpi);
1642 
1643  eta=TMath::Max(eta,0.);
1644  return eta;
1645 }
1646 //___________________________________________________________________________
1647 // Generic Phase Space Decay methods
1648 //___________________________________________________________________________
1650  GHepRecord* ev, GHepParticle* p, const PDGCodeList & pdgv,
1651  TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode)
1652 {
1653 // General method decaying the input particle system 'pdgv' with available 4-p
1654 // given by 'pd'. The decayed system is used to populate the input TMCParticle
1655 // array starting from the slot 'offset'.
1656 //
1657  LOG("INukeUtils", pINFO) << "*** Performing a Phase Space Decay";
1658  assert(pdgv.size() > 1);
1659 
1660  LOG("INukeUtils",pINFO) << "probe mass: M = " << p->Mass();
1661 
1662  // Get the decay product masses & names
1663 
1664  ostringstream state_sstream;
1665  state_sstream << "( ";
1666  vector<int>::const_iterator pdg_iter;
1667  int i = 0;
1668  double * mass = new double[pdgv.size()];
1669  double mass_sum = 0;
1670  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1671  int pdgc = *pdg_iter;
1672  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
1673  string nm = PDGLibrary::Instance()->Find(pdgc)->GetName();
1674  mass[i++] = m;
1675  mass_sum += m;
1676  state_sstream << nm << " ";
1677  }
1678  state_sstream << ")";
1679 
1680  TLorentzVector * pd = p->GetP4(); // incident particle 4p
1681 
1682  bool is_nuc = pdg::IsNeutronOrProton(p->Pdg());
1683  bool is_kaon = p->Pdg()==kPdgKP || p->Pdg()==kPdgKM;
1684  // not used // bool is_pion = p->Pdg()==kPdgPiP || p->Pdg()==kPdgPi0 || p->Pdg()==kPdgPiM;
1685  // update available energy -> init (mass + kinetic) + sum of f/s masses
1686  // for pion only. Probe mass not available for nucleon, kaon
1687  double availE = pd->Energy() + mass_sum;
1688  if(is_nuc||is_kaon) availE -= p->Mass();
1689  pd->SetE(availE);
1690 
1691  LOG("INukeUtils",pNOTICE)
1692  << "size, mass_sum, availE, pd mass, energy = " << pdgv.size() << " "
1693  << mass_sum << " " << availE << " " << p->Mass() << " " << p->Energy() ;
1694 
1695  // compute the 4p transfer to the hadronic blob
1696  double dE = mass_sum;
1697  if(is_nuc||is_kaon) dE -= p->Mass();
1698  TLorentzVector premnsub(0,0,0,dE);
1699  RemnP4 -= premnsub;
1700 
1701  LOG("INukeUtils", pINFO)
1702  << "Final state = " << state_sstream.str() << " has N = " << pdgv.size()
1703  << " particles / total mass = " << mass_sum;
1704  LOG("INukeUtils", pINFO)
1705  << "Composite system p4 = " << utils::print::P4AsString(pd);
1706 
1707  // Set the decay
1708  TGenPhaseSpace GenPhaseSpace;
1709  bool permitted = GenPhaseSpace.SetDecay(*pd, pdgv.size(), mass);
1710  if(!permitted) {
1711  LOG("INukeUtils", pERROR)
1712  << " *** Phase space decay is not permitted \n"
1713  << " Total particle mass = " << mass_sum << "\n"
1714  << " Decaying system p4 = " << utils::print::P4AsString(pd);
1715 
1716  // clean-up and return
1717  RemnP4 += premnsub;
1718  delete [] mass;
1719  delete pd;
1720  return false;
1721  }
1722 
1723  // The decay is permitted - add the incident particle at the event record
1724  // and mark is as 'Nucleon Cluster Target' (used to be confusing 'Decayed State')
1725  p->SetStatus(kIStNucleonClusterTarget); //kIStDecayedState);
1727  ev->AddParticle(*p);
1728  // Get the maximum weight
1729  double wmax = -1;
1730  for(int k=0; k<200; k++) {
1731  double w = GenPhaseSpace.Generate();
1732  wmax = TMath::Max(wmax,w);
1733  }
1734  assert(wmax>0);
1735 
1736  LOG("INukeUtils", pINFO)
1737  << "Max phase space gen. weight @ current hadronic interaction: " << wmax;
1738 
1739  // Generate an unweighted decay
1740 
1742  wmax *= 1.2;
1743 
1744  bool accept_decay=false;
1745  unsigned int itry=0;
1746 
1747  while(!accept_decay)
1748  {
1749  itry++;
1750 
1751  if(itry>kMaxUnweightDecayIterations) {
1752  // report, clean-up and return
1753  LOG("INukeUtils", pNOTICE)
1754  << "Couldn't generate an unweighted phase space decay after "
1755  << itry << " attempts";
1756  delete [] mass;
1757  delete pd;
1758  return false;
1759  }
1760 
1761  double w = GenPhaseSpace.Generate();
1762  double gw = wmax * rnd->RndFsi().Rndm();
1763 
1764  if(w > wmax) {
1765  LOG("INukeUtils", pNOTICE)
1766  << "Decay weight = " << w << " > max decay weight = " << wmax;
1767  }
1768 
1769  LOG("INukeUtils", pNOTICE) << "Decay weight = " << w << " / R = " << gw;
1770  accept_decay = (gw<=w);
1771  }
1772 
1773  // Insert final state products into the event record
1774  // - the particles are added as daughters of the decayed state
1775  // - the particles are marked as final stable state (in hA mode)
1776  i=0;
1777  int mom = ev->ParticlePosition(p);
1778  LOG("INukeUtils", pNOTICE) << "mother index = " << mom;
1781 
1782  TLorentzVector * v4 = p->GetX4();
1783 
1784  double checkpx = p->Px();
1785  double checkpy = p->Py();
1786  double checkpz = p->Pz();
1787  double checkE = p->E();
1788 
1789  //LOG("INukeUtils", PNOTICE)
1790 
1791  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1792 
1793  //-- current PDG code
1794  int pdgc = *pdg_iter;
1795  bool isnuc = pdg::IsNeutronOrProton(pdgc);
1796 
1797  //-- get the 4-momentum of the i-th final state particle
1798  TLorentzVector * p4fin = GenPhaseSpace.GetDecay(i++);
1799 
1800  //-- intranuke no longer throws "bindinos" but adds all the energy
1801  // not going at a simulated f/s particle at a "hadronic blob"
1802  // representing the remnant system: do the binding energy subtraction
1803  // here & update the remnant hadronic system 4p
1804  double M = PDGLibrary::Instance()->Find(pdgc)->Mass();
1805  double En = p4fin->Energy();
1806 
1807  double KE = En-M;
1808 
1809  //double KE = En;
1810  //if(is_pion) KE -= M;
1811 
1812  double dE_leftover = TMath::Min(NucRmvE, KE);
1813  KE -= dE_leftover;
1814  En = KE+M;
1815  double pmag_old = p4fin->P();
1816  double pmag_new = TMath::Sqrt(TMath::Max(0.,En*En-M*M));
1817  double scale = pmag_new / pmag_old;
1818  double pxn = scale * p4fin->Px();
1819  double pyn = scale * p4fin->Py();
1820  double pzn = scale * p4fin->Pz();
1821 
1822  TLorentzVector p4n(pxn,pyn,pzn,En);
1823  // LOG("INukeUtils", pNOTICE) << "Px = " << pxn << " Py = " << pyn
1824  // << " Pz = " << pzn << " E = " << KE;
1825  checkpx -= pxn;
1826  checkpy -= pyn;
1827  checkpz -= pzn;
1828  checkE -= KE;
1829 
1830  if (mode==kIMdHA &&
1831  (pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) )
1832  {
1833  if (p4n.Vect().Mag()>=0.001)
1834  {
1835  GHepParticle new_particle(pdgc, ist_pi, mom,-1,-1,-1, p4n, *v4);
1836  ev->AddParticle(new_particle);
1837  }
1838  else
1839  {
1840  // Momentum too small, assign a non-zero momentum to the particle
1841  // Conserve momentum with the remnant nucleus
1842 
1843  LOG("INukeUtils", pINFO)<<"Momentum too small; assigning 0.001 as new momentum";
1844 
1845  double phi = 2*kPi*rnd->RndFsi().Rndm();
1846  double omega = 2*rnd->RndFsi().Rndm();
1847  // throw number against solid angle for uniform distribution
1848 
1849  double E4n = TMath::Sqrt(0.001*0.001+M*M);
1850  p4n.SetPxPyPzE(0.001,0,0,E4n);
1851  p4n.Rotate(TMath::ACos(1-omega),TVector3(0,0,1));
1852  p4n.Rotate(phi,TVector3(1,0,0));
1853 
1854  RemnP4 -= (p4n - TLorentzVector(0,0,0,M));
1855 
1856  GHepParticle new_particle(pdgc, ist, mom,-1,-1,-1, p4n, *v4);
1857  ev->AddParticle(new_particle);
1858  }
1859  }
1860  else
1861  {
1862  GHepParticle new_particle(pdgc, ist, mom,-1,-1,-1, p4n, *v4);
1863 
1864  if(isnuc) new_particle.SetRemovalEnergy(0.);
1865  ev->AddParticle(new_particle);
1866  }
1867 
1868  double dpx = (1-scale)*p4fin->Px();
1869  double dpy = (1-scale)*p4fin->Py();
1870  double dpz = (1-scale)*p4fin->Pz();
1871  TLorentzVector premnadd(dpx,dpy,dpz,dE_leftover);
1872  RemnP4 += premnadd;
1873  }
1874  //LOG("INukeUtils", pNOTICE) << "TEST: " << p->Mass();
1875  LOG("INukeUtils", pNOTICE) << "check conservation: Px = " << checkpx << " Py = " << checkpy
1876  << " Pz = " << checkpz << " E = " << checkE;
1877 
1878  // Clean-up
1879  delete [] mass;
1880  delete pd;
1881  delete v4;
1882 
1883  return true;
1884 }
1885 
1887  const double &pionKineticEnergy,
1888  const double &density,
1889  const int &pionPDG,
1890  const double &protonFraction,
1891  const bool &isTableChosen
1892  )
1893 {
1894  // ------ OsetCrossSection init (only first time function is called) ------ //
1895  static INukeOset *iNukeOset = NULL;
1896 
1897  if (iNukeOset == NULL)
1898  {
1899  if (isTableChosen)
1900  {
1901  // set directory with data on first call
1902  static const std::string dataDir = (gSystem->Getenv("GINUKEHADRONDATA")) ?
1903  string(gSystem->Getenv("GINUKEHADRONDATA")) :
1904  string(gSystem->Getenv("GENIE")) +
1905  string("/data/evgen/intranuke/");
1906  // set file with Oset table on first call
1907  static const std::string dataFile = dataDir + "tot_xsec/"
1908  "intranuke-xsection-pi+N-Oset.dat";
1909  // initialize OsetCrossSection on first call
1910  iNukeOset = new INukeOsetTable (dataFile.c_str());
1911  }
1912  else iNukeOset = new INukeOsetFormula();
1913  }
1914  // ------ OsetCrossSection init (only first time function is called) ------ //
1915 
1916  // set up Oset class (assign pion Tk, nuclear density etc)
1917  iNukeOset->setupOset (density, pionKineticEnergy, pionPDG, protonFraction);
1918 
1919  return iNukeOset->getTotalCrossSection();
1920 
1921 }
double MeanFreePath_Delta(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A)
Mean free path (Delta++ test)
const double kPi
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:91
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
void SetFirstMother(int m)
Definition: GHepParticle.h:133
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:289
Basic constants.
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
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 int ParticlePosition(int pdg, GHepStatus_t i, int start=0) const
Definition: GHepRecord.cxx:181
TLorentzVector * GetX4(void) const
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
Float_t mfp
Definition: plot.C:20
const XML_Char * target
Definition: expat.h:268
static INukeNucleonCorr * getInstance()
get single instance of INukeNucleonCorr; create if necessary
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
double Density(double r, int A, double ring=0.)
Table-based implementation of Oset model.
enum genie::EGHepStatus GHepStatus_t
void StepParticle(GHepParticle *p, double step, double nuclear_radius=-1.)
Step particle.
const char * p
Definition: xmltok.h:285
static const double MeV
Definition: Units.h:130
bool ThreeBodyKinematics(GHepRecord *ev, GHepParticle *p, int tcode, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI *Nuclmodel=(const NuclearModelI *) 0)
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)
static constexpr Double_t nm
Definition: Munits.h:133
double ProbSurvival(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double mfp_scale_factor=1.0, double nRpi=0.5, double nRnuc=1.0, double NR=3, double R0=1.4)
Hadron survival probability.
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double Mass(Resonance_t res)
resonance mass (GeV)
EINukeMode
Definition: INukeMode.h:29
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
Double_t beta
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
cout<< t1-> GetEntries()
Definition: plottest35.C:29
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double dE
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
double Energy(void) const
Get energy.
Definition: GHepParticle.h:93
A list of PDG codes.
Definition: PDGCodeList.h:33
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
bool CompareMomentum(const GHepParticle *p) const
int LastMother(void) const
Definition: GHepParticle.h:68
double Dist2Exit(const TLorentzVector &x4, const TLorentzVector &p4, double A, double NR=3, double R0=1.4)
Distance to exit.
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
int FirstMother(void) const
Definition: GHepParticle.h:67
string Name(void) const
Name that corresponds to the PDG code.
const int kPdgCompNuclCluster
Definition: PDGCodes.h:194
Definition: Cand.cxx:23
Double_t scale
Definition: plot.C:25
void SetPosition(const TLorentzVector &v4)
double getTotalCrossSection() const
return total = (qel+cex+abs) cross section
Definition: INukeOset.h:35
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
const double R0
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const double NR
Float_t E
Definition: plot.C:20
const int kPdgKM
Definition: PDGCodes.h:150
const int kPdgGamma
Definition: PDGCodes.h:166
TLorentzVector * GetP4(void) const
bool ComparePdgCodes(const GHepParticle *p) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
double energy
Definition: plottest35.C:25
const int kPdgKP
Definition: PDGCodes.h:149
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils.cxx:754
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Float_t d
Definition: plot.C:236
bool TwoBodyCollision(GHepRecord *ev, int pcode, int tcode, int scode, int s2code, double C3CM, GHepParticle *p, GHepParticle *t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode=kIMdHA)
Intranuke utility functions.
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
#define R(x)
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
virtual vector< int > * GetStableDescendants(int position) const
Definition: GHepRecord.cxx:218
const double j
Definition: BetheBloch.cxx:29
def success(message)
Definition: log.py:5
double CalculateEta(double Minc, double ke, double Mtarg, double Mtwopart, double Mpi)
#define pINFO
Definition: Messenger.h:63
TVector3 Unit() const
double MeanFreePath(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double nRpi=0.5, double nRnuc=1.0, const bool useOset=false, const bool altOset=false, const bool xsecNNCorr=false, string INukeMode="XX2018")
Mean free path (pions, nucleons)
void Equilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Misc GENIE control constants.
z
Definition: test.py:28
#define pWARN
Definition: Messenger.h:61
virtual void setupOset(const double &density, const double &pionTk, const int &pionPDG, const double &protonFraction)=0
use to set up Oset class (assign pion Tk, nuclear density etc)
bool CompareStatusCodes(const GHepParticle *p) const
void SetRemovalEnergy(double Erm)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double A
Definition: Units.h:82
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static INukeHadroData2018 * Instance(void)
void SetLastMother(int m)
Definition: GHepParticle.h:134
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
double epsilon
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:314
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
float Mag() const
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
const int kPdgPiM
Definition: PDGCodes.h:136
virtual bool GenerateNucleon(const Target &) const =0
static const double fm2
Definition: Units.h:84
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const int kPdgProton
Definition: PDGCodes.h:65
tuple density
Definition: barite.py:33
#define pNOTICE
Definition: Messenger.h:62
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
static const double mb
Definition: Units.h:87
double sigmaTotalOset(const double &pionKineticEnergy, const double &density, const int &pionPDG, const double &protonFraction, const bool &isTableChosen=true)
Float_t w
Definition: plot.C:20
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
double Dist2ExitMFP(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double NR=3, double R0=1.4)
Distance to exit.
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:87
TGeoManager * gm
Definition: make_fe_box.C:4
Root of GENIE utility namespaces.
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
#define pDEBUG
Definition: Messenger.h:64
void SetPdgCode(int c)
double Py(void) const
Get Py.
Definition: GHepParticle.h:90
Formula-based implementation of Oset model.