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