CMS 3D CMS Logo
List of all members | Public Member Functions | Private Member Functions
FullModelReactionDynamics Class Reference

#include <FullModelReactionDynamics.h>

Public Member Functions

G4int Factorial (G4int n)
 
virtual G4double FindInelasticity ()
 
virtual G4double FindTimeDelay ()
 
 FullModelReactionDynamics ()
 
G4double GenerateNBodyEvent (const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
 
G4bool GenerateXandPt (G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool leadFlag, G4ReactionProduct &leadingStrangeParticle)
 
void NuclearReaction (G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
 
void ProduceStrangeParticlePairs (G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen, const G4ReactionProduct &modifiedOriginal, const G4DynamicParticle *originalTarget, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged)
 
void SuppressChargedPions (G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen, const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged)
 
void TwoBody (G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4DynamicParticle *originalTarget, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &targetHasChanged)
 
G4bool TwoCluster (G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool leadFlag, G4ReactionProduct &leadingStrangeParticle)
 
virtual ~FullModelReactionDynamics ()
 

Private Member Functions

void AddBlackTrackParticles (const G4double epnb, const G4int npnb, const G4double edta, const G4int ndta, const G4double sprob, const G4double kineticMinimum, const G4double kineticFactor, const G4ReactionProduct &modifiedOriginal, G4double spall, const G4Nucleus &aNucleus, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
 
void Defs1 (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
 
void MomentumCheck (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
 
G4double normal ()
 
G4int Poisson (G4double x)
 
void Rotate (const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
 

Detailed Description

Definition at line 23 of file FullModelReactionDynamics.h.

Constructor & Destructor Documentation

FullModelReactionDynamics::FullModelReactionDynamics ( )
inline

Definition at line 25 of file FullModelReactionDynamics.h.

25 {}
virtual FullModelReactionDynamics::~FullModelReactionDynamics ( )
inlinevirtual

Definition at line 27 of file FullModelReactionDynamics.h.

27 {}

Member Function Documentation

void FullModelReactionDynamics::AddBlackTrackParticles ( const G4double  epnb,
const G4int  npnb,
const G4double  edta,
const G4int  ndta,
const G4double  sprob,
const G4double  kineticMinimum,
const G4double  kineticFactor,
const G4ReactionProduct &  modifiedOriginal,
G4double  spall,
const G4Nucleus &  aNucleus,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2850 of file FullModelReactionDynamics.cc.

References funct::cos(), JetChargeProducer_cfi::exp, GeV, mps_fire::i, GetRecoTauVFromDQM_MC_cff::kk, cmsBatch::log, SiStripPI::max, MeV, p1, p2, createTree::pp, summarizeEdmComparisonLogfiles::ran, funct::sin(), and mathSSE::sqrt().

Referenced by FindTimeDelay().

2861  {
2862  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2863  //
2864  // npnb is number of proton/neutron black track particles
2865  // ndta is the number of deuterons, tritons, and alphas produced
2866  // epnb is the kinetic energy available for proton/neutron black track particles
2867  // edta is the kinetic energy available for deuteron/triton/alpha particles
2868  //
2869 
2870  G4ParticleDefinition *aProton = G4Proton::Proton();
2871  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
2872  G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
2873  G4ParticleDefinition *aTriton = G4Triton::Triton();
2874  G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
2875 
2876  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / MeV;
2877  G4double atomicWeight = targetNucleus.GetN_asInt();
2878  G4double atomicNumber = targetNucleus.GetZ_asInt();
2879 
2880  const G4double ika1 = 3.6;
2881  const G4double ika2 = 35.56;
2882  const G4double ika3 = 6.45;
2883  const G4double sp1 = 1.066;
2884 
2885  G4int i;
2886  G4double pp;
2887  // G4double totalQ = 0;
2888  G4double kinCreated = 0;
2889  G4double cfa = 0.025 * ((atomicWeight - 1.0) / 120.0) * std::exp(-(atomicWeight - 1.0) / 120.0);
2890  if (npnb > 0) // first add protons and neutrons
2891  {
2892  G4double backwardKinetic = 0.0;
2893  G4int local_npnb = npnb;
2894  for (i = 0; i < npnb; ++i)
2895  if (G4UniformRand() < sprob)
2896  local_npnb--;
2897  G4double ekin = epnb / local_npnb;
2898 
2899  for (i = 0; i < local_npnb; ++i) {
2900  G4ReactionProduct *p1 = new G4ReactionProduct();
2901  if (backwardKinetic > epnb) {
2902  delete p1;
2903  break;
2904  }
2905  G4double ran = G4UniformRand();
2906  G4double kinetic = -ekin * std::log(ran) - cfa * (1.0 + 0.5 * normal());
2907  if (kinetic < 0.0)
2908  kinetic = -0.010 * std::log(ran);
2909  backwardKinetic += kinetic;
2910  if (backwardKinetic > epnb)
2911  kinetic = std::max(kineticMinimum, epnb - (backwardKinetic - kinetic));
2912  if (G4UniformRand() > (1.0 - atomicNumber / atomicWeight))
2913  p1->SetDefinition(aProton);
2914  else
2915  p1->SetDefinition(aNeutron);
2916  vec.SetElement(vecLen, p1);
2917  ++spall;
2918  G4double cost = G4UniformRand() * 2.0 - 1.0;
2919  G4double sint = std::sqrt(std::fabs(1.0 - cost * cost));
2920  G4double phi = twopi * G4UniformRand();
2921  vec[vecLen]->SetNewlyAdded(true);
2922  vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2923  kinCreated += kinetic;
2924  pp = vec[vecLen]->GetTotalMomentum() / MeV;
2925  vec[vecLen]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2926  vecLen++;
2927  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2928  }
2929  if ((atomicWeight >= 10.0) && (ekOriginal <= 2.0 * GeV)) {
2930  G4double ekw = ekOriginal / GeV;
2931  G4int ika, kk = 0;
2932  if (ekw > 1.0)
2933  ekw *= ekw;
2934  ekw = std::max(0.1, ekw);
2935  ika = G4int(ika1 * std::exp((atomicNumber * atomicNumber / atomicWeight - ika2) / ika3) / ekw);
2936  if (ika > 0) {
2937  for (i = (vecLen - 1); i >= 0; --i) {
2938  if ((vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded()) {
2939  vec[i]->SetDefinitionAndUpdateE(aNeutron); // modified 22-Oct-97
2940  if (++kk > ika)
2941  break;
2942  }
2943  }
2944  }
2945  }
2946  }
2947  if (ndta > 0) // now, try to add deuterons, tritons and alphas
2948  {
2949  G4double backwardKinetic = 0.0;
2950  G4int local_ndta = ndta;
2951  for (i = 0; i < ndta; ++i)
2952  if (G4UniformRand() < sprob)
2953  local_ndta--;
2954  G4double ekin = edta / local_ndta;
2955 
2956  for (i = 0; i < local_ndta; ++i) {
2957  G4ReactionProduct *p2 = new G4ReactionProduct();
2958  if (backwardKinetic > edta) {
2959  delete p2;
2960  break;
2961  }
2962  G4double ran = G4UniformRand();
2963  G4double kinetic = -ekin * std::log(ran) - cfa * (1. + 0.5 * normal());
2964  if (kinetic < 0.0)
2965  kinetic = kineticFactor * std::log(ran);
2966  backwardKinetic += kinetic;
2967  if (backwardKinetic > edta)
2968  kinetic = edta - (backwardKinetic - kinetic);
2969  if (kinetic < 0.0)
2970  kinetic = kineticMinimum;
2971  G4double cost = 2.0 * G4UniformRand() - 1.0;
2972  G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
2973  G4double phi = twopi * G4UniformRand();
2974  ran = G4UniformRand();
2975  if (ran <= 0.60)
2976  p2->SetDefinition(aDeuteron);
2977  else if (ran <= 0.90)
2978  p2->SetDefinition(aTriton);
2979  else
2980  p2->SetDefinition(anAlpha);
2981  spall += p2->GetMass() / GeV * sp1;
2982  if (spall > atomicWeight) {
2983  delete p2;
2984  break;
2985  }
2986  vec.SetElement(vecLen, p2);
2987  vec[vecLen]->SetNewlyAdded(true);
2988  vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2989  kinCreated += kinetic;
2990  pp = vec[vecLen]->GetTotalMomentum() / MeV;
2991  vec[vecLen++]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2992  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2993  }
2994  }
2995  // G4double delta = epnb+edta - kinCreated;
2996 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
p2
double p2[4]
Definition: TauolaWrapper.h:90
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
p1
double p1[4]
Definition: TauolaWrapper.h:89
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
void FullModelReactionDynamics::Defs1 ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2588 of file FullModelReactionDynamics.cc.

References funct::abs(), funct::cos(), mps_fire::i, MeV, AlCaHLTBitMon_ParallelJobs::p, CombinatorialSeedGeneratorForBeamHalo_cfi::SetMomentum, funct::sin(), and mathSSE::sqrt().

Referenced by FindTimeDelay().

2592  {
2593  const G4double pjx = modifiedOriginal.GetMomentum().x() / MeV;
2594  const G4double pjy = modifiedOriginal.GetMomentum().y() / MeV;
2595  const G4double pjz = modifiedOriginal.GetMomentum().z() / MeV;
2596  const G4double p = modifiedOriginal.GetMomentum().mag() / MeV;
2597  if (pjx * pjx + pjy * pjy > 0.0) {
2598  G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2599  cost = pjz / p;
2600  sint = 0.5 * (std::sqrt(std::abs((1.0 - cost) * (1.0 + cost))) + std::sqrt(pjx * pjx + pjy * pjy) / p);
2601  if (pjy < 0.0)
2602  ph = 3 * halfpi;
2603  else
2604  ph = halfpi;
2605  if (std::abs(pjx) > 0.001 * MeV)
2606  ph = std::atan2(pjy, pjx);
2607  cosp = std::cos(ph);
2608  sinp = std::sin(ph);
2609  pix = currentParticle.GetMomentum().x() / MeV;
2610  piy = currentParticle.GetMomentum().y() / MeV;
2611  piz = currentParticle.GetMomentum().z() / MeV;
2612  currentParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2613  cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2614  -sint * pix * MeV + cost * piz * MeV);
2615  pix = targetParticle.GetMomentum().x() / MeV;
2616  piy = targetParticle.GetMomentum().y() / MeV;
2617  piz = targetParticle.GetMomentum().z() / MeV;
2618  targetParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2619  cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2620  -sint * pix * MeV + cost * piz * MeV);
2621  for (G4int i = 0; i < vecLen; ++i) {
2622  pix = vec[i]->GetMomentum().x() / MeV;
2623  piy = vec[i]->GetMomentum().y() / MeV;
2624  piz = vec[i]->GetMomentum().z() / MeV;
2625  vec[i]->SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2626  cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2627  -sint * pix * MeV + cost * piz * MeV);
2628  }
2629  } else {
2630  if (pjz < 0.0) {
2631  currentParticle.SetMomentum(-currentParticle.GetMomentum().z());
2632  targetParticle.SetMomentum(-targetParticle.GetMomentum().z());
2633  for (G4int i = 0; i < vecLen; ++i)
2634  vec[i]->SetMomentum(-vec[i]->GetMomentum().z());
2635  }
2636  }
2637 }
mps_fire.i
i
Definition: mps_fire.py:338
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
CombinatorialSeedGeneratorForBeamHalo_cfi.SetMomentum
SetMomentum
Definition: CombinatorialSeedGeneratorForBeamHalo_cfi.py:99
G4int FullModelReactionDynamics::Factorial ( G4int  n)

Definition at line 2578 of file FullModelReactionDynamics.cc.

References mps_fire::i, funct::m, min(), and mps_fire::result.

Referenced by FindTimeDelay().

2578  { // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
2579  G4int m = std::min(n, 10);
2580  G4int result = 1;
2581  if (m <= 1)
2582  return result;
2583  for (G4int i = 2; i <= m; ++i)
2584  result *= i;
2585  return result;
2586 }
mps_fire.i
i
Definition: mps_fire.py:338
mps_fire.result
result
Definition: mps_fire.py:291
min
T min(T a, T b)
Definition: MathUtil.h:58
virtual G4double FullModelReactionDynamics::FindInelasticity ( )
inlinevirtual

Definition at line 29 of file FullModelReactionDynamics.h.

29 { return 0.0; }
virtual G4double FullModelReactionDynamics::FindTimeDelay ( )
inlinevirtual

Definition at line 31 of file FullModelReactionDynamics.h.

References AddBlackTrackParticles(), Defs1(), Factorial(), GenerateNBodyEvent(), GenerateXandPt(), MomentumCheck(), gen::n, normal(), NuclearReaction(), Poisson(), ProduceStrangeParticlePairs(), Rotate(), SuppressChargedPions(), groupFilesInBlocks::temp, TwoBody(), TwoCluster(), and x.

31 { return 0.0; }
G4double FullModelReactionDynamics::GenerateNBodyEvent ( const G4double  totalEnergy,
const G4bool  constantCrossSection,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)

Definition at line 2340 of file FullModelReactionDynamics.cc.

References a, b, pfBoostedDoubleSVAK8TagInfos_cfi::beta, EnergyCorrector::c, funct::cos(), JetChargeProducer_cfi::exp, ecalTB2006H4_GenSimDigiReco_cfg::G4cerr, GeV, mps_fire::i, cmsBatch::log, ResonanceBuilder::mass, SiStripPI::max, min(), funct::pow(), summarizeEdmComparisonLogfiles::ran, alignCSCRings::s, indexGen::s2, funct::sin(), mathSSE::sqrt(), groupFilesInBlocks::temp, and mps_merge::weight.

Referenced by FindTimeDelay().

2343  {
2344  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2345  // derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
2346  // Returns the weight of the event
2347  //
2348  G4int i;
2349  const G4double expxu = 82.; // upper bound for arg. of exp
2350  const G4double expxl = -expxu; // lower bound for arg. of exp
2351  if (vecLen < 2) {
2352  G4cerr << "*** Error in FullModelReactionDynamics::GenerateNBodyEvent" << G4endl;
2353  G4cerr << " number of particles < 2" << G4endl;
2354  G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
2355  return -1.0;
2356  }
2357  G4double mass[18]; // mass of each particle
2358  G4double energy[18]; // total energy of each particle
2359  G4double pcm[3][18]; // pcm is an array with 3 rows and vecLen columns
2360  G4double totalMass = 0.0;
2361  G4double extraMass = 0;
2362  G4double sm[18];
2363 
2364  for (i = 0; i < vecLen; ++i) {
2365  mass[i] = vec[i]->GetMass() / GeV;
2366  if (vec[i]->GetSide() == -2)
2367  extraMass += vec[i]->GetMass() / GeV;
2368  vec[i]->SetMomentum(0.0, 0.0, 0.0);
2369  pcm[0][i] = 0.0; // x-momentum of i-th particle
2370  pcm[1][i] = 0.0; // y-momentum of i-th particle
2371  pcm[2][i] = 0.0; // z-momentum of i-th particle
2372  energy[i] = mass[i]; // total energy of i-th particle
2373  totalMass += mass[i];
2374  sm[i] = totalMass;
2375  }
2376  G4double totalE = totalEnergy / GeV;
2377  if (totalMass > totalE) {
2378  return -1.0;
2379  }
2380  G4double kineticEnergy = totalE - totalMass;
2381  G4double emm[18];
2382  //G4double *emm = new G4double [vecLen];
2383  emm[0] = mass[0];
2384  emm[vecLen - 1] = totalE;
2385  if (vecLen > 2) // the random numbers are sorted
2386  {
2387  G4double ran[18];
2388  for (i = 0; i < vecLen; ++i)
2389  ran[i] = G4UniformRand();
2390  for (i = 0; i < vecLen - 2; ++i) {
2391  for (G4int j = vecLen - 2; j > i; --j) {
2392  if (ran[i] > ran[j]) {
2393  G4double temp = ran[i];
2394  ran[i] = ran[j];
2395  ran[j] = temp;
2396  }
2397  }
2398  }
2399  for (i = 1; i < vecLen - 1; ++i)
2400  emm[i] = ran[i - 1] * kineticEnergy + sm[i];
2401  }
2402  // Weight is the sum of logarithms of terms instead of the product of terms
2403  G4bool lzero = true;
2404  G4double wtmax = 0.0;
2405  if (constantCrossSection) // this is KGENEV=1 in PHASP
2406  {
2407  G4double emmax = kineticEnergy + mass[0];
2408  G4double emmin = 0.0;
2409  for (i = 1; i < vecLen; ++i) {
2410  emmin += mass[i - 1];
2411  emmax += mass[i];
2412  G4double wtfc = 0.0;
2413  if (emmax * emmax > 0.0) {
2414  G4double arg = emmax * emmax +
2415  (emmin * emmin - mass[i] * mass[i]) * (emmin * emmin - mass[i] * mass[i]) / (emmax * emmax) -
2416  2.0 * (emmin * emmin + mass[i] * mass[i]);
2417  if (arg > 0.0)
2418  wtfc = 0.5 * std::sqrt(arg);
2419  }
2420  if (wtfc == 0.0) {
2421  lzero = false;
2422  break;
2423  }
2424  wtmax += std::log(wtfc);
2425  }
2426  if (lzero)
2427  wtmax = -wtmax;
2428  else
2429  wtmax = expxu;
2430  } else {
2431  // ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
2432  const G4double ffq[18] = {0.,
2433  3.141592,
2434  19.73921,
2435  62.01255,
2436  129.8788,
2437  204.0131,
2438  256.3704,
2439  268.4705,
2440  240.9780,
2441  189.2637,
2442  132.1308,
2443  83.0202,
2444  47.4210,
2445  24.8295,
2446  12.0006,
2447  5.3858,
2448  2.2560,
2449  0.8859};
2450  wtmax = std::log(std::pow(kineticEnergy, vecLen - 2) * ffq[vecLen - 1] / totalE);
2451  }
2452  lzero = true;
2453  G4double pd[50];
2454  //G4double *pd = new G4double [vecLen-1];
2455  for (i = 0; i < vecLen - 1; ++i) {
2456  pd[i] = 0.0;
2457  if (emm[i + 1] * emm[i + 1] > 0.0) {
2458  G4double arg = emm[i + 1] * emm[i + 1] +
2459  (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) * (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) /
2460  (emm[i + 1] * emm[i + 1]) -
2461  2.0 * (emm[i] * emm[i] + mass[i + 1] * mass[i + 1]);
2462  if (arg > 0.0)
2463  pd[i] = 0.5 * std::sqrt(arg);
2464  }
2465  if (pd[i] <= 0.0) // changed from == on 02 April 98
2466  lzero = false;
2467  else
2468  wtmax += std::log(pd[i]);
2469  }
2470  G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
2471  if (lzero)
2472  weight = std::exp(std::max(std::min(wtmax, expxu), expxl));
2473 
2474  G4double bang, cb, sb, s0, s1, s2, c, s, esys, a, b, gama, beta;
2475  pcm[0][0] = 0.0;
2476  pcm[1][0] = pd[0];
2477  pcm[2][0] = 0.0;
2478  for (i = 1; i < vecLen; ++i) {
2479  pcm[0][i] = 0.0;
2480  pcm[1][i] = -pd[i - 1];
2481  pcm[2][i] = 0.0;
2482  bang = twopi * G4UniformRand();
2483  cb = std::cos(bang);
2484  sb = std::sin(bang);
2485  c = 2.0 * G4UniformRand() - 1.0;
2486  s = std::sqrt(std::fabs(1.0 - c * c));
2487  if (i < vecLen - 1) {
2488  esys = std::sqrt(pd[i] * pd[i] + emm[i] * emm[i]);
2489  beta = pd[i] / esys;
2490  gama = esys / emm[i];
2491  for (G4int j = 0; j <= i; ++j) {
2492  s0 = pcm[0][j];
2493  s1 = pcm[1][j];
2494  s2 = pcm[2][j];
2495  energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2496  a = s0 * c - s1 * s; // rotation
2497  pcm[1][j] = s0 * s + s1 * c;
2498  b = pcm[2][j];
2499  pcm[0][j] = a * cb - b * sb;
2500  pcm[2][j] = a * sb + b * cb;
2501  pcm[1][j] = gama * (pcm[1][j] + beta * energy[j]);
2502  }
2503  } else {
2504  for (G4int j = 0; j <= i; ++j) {
2505  s0 = pcm[0][j];
2506  s1 = pcm[1][j];
2507  s2 = pcm[2][j];
2508  energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2509  a = s0 * c - s1 * s; // rotation
2510  pcm[1][j] = s0 * s + s1 * c;
2511  b = pcm[2][j];
2512  pcm[0][j] = a * cb - b * sb;
2513  pcm[2][j] = a * sb + b * cb;
2514  }
2515  }
2516  }
2517  for (i = 0; i < vecLen; ++i) {
2518  vec[i]->SetMomentum(pcm[0][i] * GeV, pcm[1][i] * GeV, pcm[2][i] * GeV);
2519  vec[i]->SetTotalEnergy(energy[i] * GeV);
2520  }
2521  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2522  return weight;
2523 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
pfBoostedDoubleSVAK8TagInfos_cfi.beta
beta
Definition: pfBoostedDoubleSVAK8TagInfos_cfi.py:7
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
weight
Definition: weight.py:1
funct::arg
A arg
Definition: Factorize.h:38
ecalTB2006H4_GenSimDigiReco_cfg.G4cerr
G4cerr
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:304
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
mps_merge.weight
weight
Definition: mps_merge.py:88
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
G4bool FullModelReactionDynamics::GenerateXandPt ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged,
G4bool  leadFlag,
G4ReactionProduct &  leadingStrangeParticle 
)

Definition at line 62 of file FullModelReactionDynamics.cc.

References funct::abs(), funct::cos(), diffTreeTool::diff, MillePedeFileConverter_cfg::e, stringResolutionProvider_cfi::et, cmsRelvalreport::exit, JetChargeProducer_cfi::exp, ecalTB2006H4_GenSimDigiReco_cfg::G4cerr, ecalTB2006H4_GenSimDigiReco_cfg::G4cout, GeV, mps_fire::i, checklumidiff::l, BPhysicsValidation_cfi::Lambda, cmsBatch::log, funct::m, SiStripPI::max, MeV, min(), pi, funct::pow(), createTree::pp, EnergyCorrector::pt, summarizeEdmComparisonLogfiles::ran, phase2TrackerDigitizer_cfi::SigmaZero, funct::sin(), funct::sqr(), mathSSE::sqrt(), and groupFilesInBlocks::temp.

Referenced by FullModelHadronicProcess::CalculateMomenta(), and FindTimeDelay().

73  {
74  //
75  // derived from original FORTRAN code GENXPT by H. Fesefeldt (11-Oct-1987)
76  //
77  // Generation of X- and PT- values for incident, target, and all secondary particles
78  // A simple single variable description E D3S/DP3= F(Q) with
79  // Q^2 = (M*X)^2 + PT^2 is used. Final state kinematic is produced
80  // by an FF-type iterative cascade method
81  //
82  // internal units are GeV
83  //
84  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
85 
86  // Protection in case no secondary has been created; cascades down to two-body.
87  if (vecLen == 0)
88  return false;
89 
90  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
91  // G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
92  // G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
93  G4ParticleDefinition *aProton = G4Proton::Proton();
94  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
95  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
96  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
97  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
98  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
99  G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
100  G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
101 
102  G4int i, l;
103  // G4double forVeryForward = 0.;
104  G4bool veryForward = false;
105 
106  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
107  const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
108  const G4double mOriginal = modifiedOriginal.GetMass() / GeV;
109  const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
110  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass() / GeV;
111  G4double centerofmassEnergy =
112  std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal); // GeV
113  G4double currentMass = currentParticle.GetMass() / GeV;
114  targetMass = targetParticle.GetMass() / GeV;
115  //
116  // randomize the order of the secondary particles
117  // note that the current and target particles are not affected
118  //
119  for (i = 0; i < vecLen; ++i) {
120  G4int itemp = G4int(G4UniformRand() * vecLen);
121  G4ReactionProduct pTemp = *vec[itemp];
122  *vec[itemp] = *vec[i];
123  *vec[i] = pTemp;
124  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
125  }
126 
127  if (currentMass == 0.0 && targetMass == 0.0) // annihilation
128  {
129  // no kinetic energy in target .....
130  G4double ek = currentParticle.GetKineticEnergy();
131  G4ThreeVector m = currentParticle.GetMomentum();
132  currentParticle = *vec[0];
133  targetParticle = *vec[1];
134  for (i = 0; i < (vecLen - 2); ++i)
135  *vec[i] = *vec[i + 2];
136  G4ReactionProduct *temp = vec[vecLen - 1];
137  delete temp;
138  temp = vec[vecLen - 2];
139  delete temp;
140  vecLen -= 2;
141  currentMass = currentParticle.GetMass() / GeV;
142  targetMass = targetParticle.GetMass() / GeV;
143  incidentHasChanged = true;
144  targetHasChanged = true;
145  currentParticle.SetKineticEnergy(ek);
146  currentParticle.SetMomentum(m);
147  veryForward = true;
148  }
149  const G4double atomicWeight = targetNucleus.GetN_asInt();
150  const G4double atomicNumber = targetNucleus.GetZ_asInt();
151  const G4double protonMass = aProton->GetPDGMass() / MeV;
152  if ((originalIncident->GetDefinition() == aKaonMinus || originalIncident->GetDefinition() == aKaonZeroL ||
153  originalIncident->GetDefinition() == aKaonZeroS || originalIncident->GetDefinition() == aKaonPlus) &&
154  G4UniformRand() >= 0.7) {
155  G4ReactionProduct temp = currentParticle;
156  currentParticle = targetParticle;
157  targetParticle = temp;
158  incidentHasChanged = true;
159  targetHasChanged = true;
160  currentMass = currentParticle.GetMass() / GeV;
161  targetMass = targetParticle.GetMass() / GeV;
162  }
163  const G4double afc = std::min(0.75,
164  0.312 + 0.200 * std::log(std::log(centerofmassEnergy * centerofmassEnergy)) +
165  std::pow(centerofmassEnergy * centerofmassEnergy, 1.5) / 6000.0);
166 
167  // PROBLEMET ER HER!!!
168 
169  G4double freeEnergy = centerofmassEnergy - currentMass - targetMass;
170 
171  if (freeEnergy < 0) {
172  G4cout << "Free energy < 0!" << G4endl;
173  G4cout << "E_CMS = " << centerofmassEnergy << " GeV" << G4endl;
174  G4cout << "m_curr = " << currentMass << " GeV" << G4endl;
175  G4cout << "m_orig = " << mOriginal << " GeV" << G4endl;
176  G4cout << "m_targ = " << targetMass << " GeV" << G4endl;
177  G4cout << "E_free = " << freeEnergy << " GeV" << G4endl;
178  }
179 
180  G4double forwardEnergy = freeEnergy / 2.;
181  G4int forwardCount = 1; // number of particles in forward hemisphere
182 
183  G4double backwardEnergy = freeEnergy / 2.;
184  G4int backwardCount = 1; // number of particles in backward hemisphere
185  if (veryForward) {
186  if (currentParticle.GetSide() == -1) {
187  forwardEnergy += currentMass;
188  forwardCount--;
189  backwardEnergy -= currentMass;
190  backwardCount++;
191  }
192  if (targetParticle.GetSide() != -1) {
193  backwardEnergy += targetMass;
194  backwardCount--;
195  forwardEnergy -= targetMass;
196  forwardCount++;
197  }
198  }
199  for (i = 0; i < vecLen; ++i) {
200  if (vec[i]->GetSide() == -1) {
201  ++backwardCount;
202  backwardEnergy -= vec[i]->GetMass() / GeV;
203  } else {
204  ++forwardCount;
205  forwardEnergy -= vec[i]->GetMass() / GeV;
206  }
207  }
208  //
209  // add particles from intranuclear cascade
210  // nuclearExcitationCount = number of new secondaries produced by nuclear excitation
211  // extraCount = number of nucleons within these new secondaries
212  //
213  G4double xtarg;
214  if (centerofmassEnergy < (2.0 + G4UniformRand()))
215  xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount + vecLen + 2) / 2.0;
216  else
217  xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount);
218  if (xtarg <= 0.0)
219  xtarg = 0.01;
220  G4int nuclearExcitationCount = Poisson(xtarg);
221  if (atomicWeight < 1.0001)
222  nuclearExcitationCount = 0;
223  G4int extraNucleonCount = 0;
224  G4double extraNucleonMass = 0.0;
225  if (nuclearExcitationCount > 0) {
226  const G4double nucsup[] = {1.00, 0.7, 0.5, 0.4, 0.35, 0.3};
227  const G4double psup[] = {3., 6., 20., 50., 100., 1000.};
228  G4int momentumBin = 0;
229  while ((momentumBin < 6) && (modifiedOriginal.GetTotalMomentum() / GeV > psup[momentumBin]))
230  ++momentumBin;
231  momentumBin = std::min(5, momentumBin);
232  //
233  // NOTE: in GENXPT, these new particles were given negative codes
234  // here I use NewlyAdded = true instead
235  //
236  for (i = 0; i < nuclearExcitationCount; ++i) {
237  G4ReactionProduct *pVec = new G4ReactionProduct();
238  if (G4UniformRand() < nucsup[momentumBin]) {
239  if (G4UniformRand() > 1.0 - atomicNumber / atomicWeight)
240  pVec->SetDefinition(aProton);
241  else
242  pVec->SetDefinition(aNeutron);
243  pVec->SetSide(-2); // -2 means backside nucleon
244  ++extraNucleonCount;
245  backwardEnergy += pVec->GetMass() / GeV;
246  extraNucleonMass += pVec->GetMass() / GeV;
247  } else {
248  G4double ran = G4UniformRand();
249  if (ran < 0.3181)
250  pVec->SetDefinition(aPiPlus);
251  else if (ran < 0.6819)
252  pVec->SetDefinition(aPiZero);
253  else
254  pVec->SetDefinition(aPiMinus);
255  pVec->SetSide(-1); // backside particle, but not a nucleon
256  }
257  pVec->SetNewlyAdded(true); // true is the same as IPA(i)<0
258  vec.SetElement(vecLen++, pVec);
259  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
260  backwardEnergy -= pVec->GetMass() / GeV;
261  ++backwardCount;
262  }
263  }
264  //
265  // assume conservation of kinetic energy in forward & backward hemispheres
266  //
267  G4int is, iskip;
268  while (forwardEnergy <= 0.0) // must eliminate a particle from the forward side
269  {
270  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
271  iskip = G4int(G4UniformRand() * forwardCount) + 1; // 1 <= iskip <= forwardCount
272  is = 0;
273  G4int forwardParticlesLeft = 0;
274  for (i = (vecLen - 1); i >= 0; --i) {
275  if (vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled()) {
276  forwardParticlesLeft = 1;
277  if (++is == iskip) {
278  forwardEnergy += vec[i]->GetMass() / GeV;
279  for (G4int j = i; j < (vecLen - 1); j++)
280  *vec[j] = *vec[j + 1]; // shift up
281  --forwardCount;
282  G4ReactionProduct *temp = vec[vecLen - 1];
283  delete temp;
284  if (--vecLen == 0)
285  return false; // all the secondaries have been eliminated
286  break; // --+
287  } // |
288  } // |
289  } // break goes down to here
290  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
291  if (forwardParticlesLeft == 0) {
292  forwardEnergy += currentParticle.GetMass() / GeV;
293  currentParticle.SetDefinitionAndUpdateE(targetParticle.GetDefinition());
294  targetParticle.SetDefinitionAndUpdateE(vec[0]->GetDefinition());
295  // above two lines modified 20-oct-97: were just simple equalities
296  --forwardCount;
297  for (G4int j = 0; j < (vecLen - 1); ++j)
298  *vec[j] = *vec[j + 1];
299  G4ReactionProduct *temp = vec[vecLen - 1];
300  delete temp;
301  if (--vecLen == 0)
302  return false; // all the secondaries have been eliminated
303  break;
304  }
305  }
306  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
307  while (backwardEnergy <= 0.0) // must eliminate a particle from the backward side
308  {
309  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
310  iskip = G4int(G4UniformRand() * backwardCount) + 1; // 1 <= iskip <= backwardCount
311  is = 0;
312  G4int backwardParticlesLeft = 0;
313  for (i = (vecLen - 1); i >= 0; --i) {
314  if (vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled()) {
315  backwardParticlesLeft = 1;
316  if (++is == iskip) // eliminate the i'th particle
317  {
318  if (vec[i]->GetSide() == -2) {
319  --extraNucleonCount;
320  extraNucleonMass -= vec[i]->GetMass() / GeV;
321  backwardEnergy -= vec[i]->GetTotalEnergy() / GeV;
322  }
323  backwardEnergy += vec[i]->GetTotalEnergy() / GeV;
324  for (G4int j = i; j < (vecLen - 1); ++j)
325  *vec[j] = *vec[j + 1]; // shift up
326  --backwardCount;
327  G4ReactionProduct *temp = vec[vecLen - 1];
328  delete temp;
329  if (--vecLen == 0)
330  return false; // all the secondaries have been eliminated
331  break;
332  }
333  }
334  }
335  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
336  if (backwardParticlesLeft == 0) {
337  backwardEnergy += targetParticle.GetMass() / GeV;
338  targetParticle = *vec[0];
339  --backwardCount;
340  for (G4int j = 0; j < (vecLen - 1); ++j)
341  *vec[j] = *vec[j + 1];
342  G4ReactionProduct *temp = vec[vecLen - 1];
343  delete temp;
344  if (--vecLen == 0)
345  return false; // all the secondaries have been eliminated
346  break;
347  }
348  }
349  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
350  //
351  // define initial state vectors for Lorentz transformations
352  // the pseudoParticles have non-standard masses, hence the "pseudo"
353  //
354  G4ReactionProduct pseudoParticle[10];
355  for (i = 0; i < 10; ++i)
356  pseudoParticle[i].SetZero();
357 
358  pseudoParticle[0].SetMass(mOriginal * GeV);
359  pseudoParticle[0].SetMomentum(0.0, 0.0, pOriginal * GeV);
360  pseudoParticle[0].SetTotalEnergy(std::sqrt(pOriginal * pOriginal + mOriginal * mOriginal) * GeV);
361 
362  pseudoParticle[1].SetMass(protonMass * MeV); // this could be targetMass
363  pseudoParticle[1].SetTotalEnergy(protonMass * MeV);
364 
365  pseudoParticle[3].SetMass(protonMass * (1 + extraNucleonCount) * MeV);
366  pseudoParticle[3].SetTotalEnergy(protonMass * (1 + extraNucleonCount) * MeV);
367 
368  pseudoParticle[8].SetMomentum(1.0 * GeV, 0.0, 0.0);
369 
370  pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
371  pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
372 
373  pseudoParticle[0].Lorentz(pseudoParticle[0], pseudoParticle[2]);
374  pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[2]);
375 
376  G4double dndl[20];
377  //
378  // main loop for 4-momentum generation
379  // see Pitha-report (Aachen) for a detailed description of the method
380  //
381  G4double aspar, pt, et, x, pp, pp1, rthnve, phinve, rmb, wgt;
382  G4int innerCounter, outerCounter;
383  G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
384 
385  G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
386  //
387  // process the secondary particles in reverse order
388  // the incident particle is Done after the secondaries
389  // nucleons, including the target, in the backward hemisphere are also Done later
390  //
391  G4double binl[20] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
392  1.0, 1.11, 1.25, 1.43, 1.67, 2.0, 2.5, 3.33, 5.00, 10.00};
393  G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
394  G4double totalEnergy, kineticEnergy, vecMass;
395 
396  for (i = (vecLen - 1); i >= 0; --i) {
397  G4double phi = G4UniformRand() * twopi;
398  if (vec[i]->GetNewlyAdded()) // added from intranuclear cascade
399  {
400  if (vec[i]->GetSide() == -2) // is a nucleon
401  {
402  if (backwardNucleonCount < 18) {
403  if (vec[i]->GetDefinition() == G4PionMinus::PionMinus() ||
404  vec[i]->GetDefinition() == G4PionPlus::PionPlus() || vec[i]->GetDefinition() == G4PionZero::PionZero()) {
405  for (G4int i = 0; i < vecLen; i++)
406  delete vec[i];
407  vecLen = 0;
408  throw G4HadReentrentException(
409  __FILE__,
410  __LINE__,
411  "FullModelReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
412  }
413  vec[i]->SetSide(-3);
414  ++backwardNucleonCount;
415  continue;
416  }
417  }
418  }
419  //
420  // set pt and phi values, they are changed somewhat in the iteration loop
421  // set mass parameter for lambda fragmentation model
422  //
423  vecMass = vec[i]->GetMass() / GeV;
424  G4double ran = -std::log(1.0 - G4UniformRand()) / 3.5;
425  if (vec[i]->GetSide() == -2) // backward nucleon
426  {
427  if (vec[i]->GetDefinition() == aKaonMinus || vec[i]->GetDefinition() == aKaonZeroL ||
428  vec[i]->GetDefinition() == aKaonZeroS || vec[i]->GetDefinition() == aKaonPlus ||
429  vec[i]->GetDefinition() == aPiMinus || vec[i]->GetDefinition() == aPiZero ||
430  vec[i]->GetDefinition() == aPiPlus) {
431  aspar = 0.75;
432  pt = std::sqrt(std::pow(ran, 1.7));
433  } else { // vec[i] must be a proton, neutron,
434  aspar = 0.20; // lambda, sigma, xsi, or ion
435  pt = std::sqrt(std::pow(ran, 1.2));
436  }
437  } else { // not a backward nucleon
438  if (vec[i]->GetDefinition() == aPiMinus || vec[i]->GetDefinition() == aPiZero ||
439  vec[i]->GetDefinition() == aPiPlus) {
440  aspar = 0.75;
441  pt = std::sqrt(std::pow(ran, 1.7));
442  } else if (vec[i]->GetDefinition() == aKaonMinus || vec[i]->GetDefinition() == aKaonZeroL ||
443  vec[i]->GetDefinition() == aKaonZeroS || vec[i]->GetDefinition() == aKaonPlus) {
444  aspar = 0.70;
445  pt = std::sqrt(std::pow(ran, 1.7));
446  } else { // vec[i] must be a proton, neutron,
447  aspar = 0.65; // lambda, sigma, xsi, or ion
448  pt = std::sqrt(std::pow(ran, 1.5));
449  }
450  }
451  pt = std::max(0.001, pt);
452  vec[i]->SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
453  for (G4int j = 0; j < 20; ++j)
454  binl[j] = j / (19. * pt);
455  if (vec[i]->GetSide() > 0)
456  et = pseudoParticle[0].GetTotalEnergy() / GeV;
457  else
458  et = pseudoParticle[1].GetTotalEnergy() / GeV;
459  dndl[0] = 0.0;
460  //
461  // start of outer iteration loop
462  //
463  outerCounter = 0;
464  eliminateThisParticle = true;
465  resetEnergies = true;
466  while (++outerCounter < 3) {
467  for (l = 1; l < 20; ++l) {
468  x = (binl[l] + binl[l - 1]) / 2.;
469  pt = std::max(0.001, pt);
470  if (x > 1.0 / pt)
471  dndl[l] += dndl[l - 1]; // changed from just = on 02 April 98
472  else
473  dndl[l] = et * aspar / std::sqrt(std::pow((1. + aspar * x * aspar * x), 3)) * (binl[l] - binl[l - 1]) /
474  std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
475  dndl[l - 1];
476  }
477  innerCounter = 0;
478  vec[i]->SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
479  //
480  // start of inner iteration loop
481  //
482  while (++innerCounter < 7) {
483  ran = G4UniformRand() * dndl[19];
484  l = 1;
485  while ((ran >= dndl[l]) && (l < 20))
486  l++;
487  l = std::min(19, l);
488  x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
489  if (vec[i]->GetSide() < 0)
490  x *= -1.;
491  vec[i]->SetMomentum(x * et * GeV); // set the z-momentum
492  totalEnergy = std::sqrt(x * et * x * et + pt * pt + vecMass * vecMass);
493  vec[i]->SetTotalEnergy(totalEnergy * GeV);
494  kineticEnergy = vec[i]->GetKineticEnergy() / GeV;
495  if (vec[i]->GetSide() > 0) // forward side
496  {
497  if ((forwardKinetic + kineticEnergy) < 0.95 * forwardEnergy) {
498  pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
499  forwardKinetic += kineticEnergy;
500  pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
501  pseudoParticle[6].SetMomentum(0.0); // set the z-momentum
502  phi = pseudoParticle[6].Angle(pseudoParticle[8]);
503  if (pseudoParticle[6].GetMomentum().y() / MeV < 0.0)
504  phi = twopi - phi;
505  phi += pi + normal() * pi / 12.0;
506  if (phi > twopi)
507  phi -= twopi;
508  if (phi < 0.0)
509  phi = twopi - phi;
510  outerCounter = 2; // leave outer loop
511  eliminateThisParticle = false; // don't eliminate this particle
512  resetEnergies = false;
513  break; // leave inner loop
514  }
515  if (innerCounter > 5)
516  break; // leave inner loop
517  if (backwardEnergy >= vecMass) // switch sides
518  {
519  vec[i]->SetSide(-1);
520  forwardEnergy += vecMass;
521  backwardEnergy -= vecMass;
522  ++backwardCount;
523  }
524  } else { // backward side
525  G4double xxx = 0.95 + 0.05 * extraNucleonCount / 20.0;
526  if ((backwardKinetic + kineticEnergy) < xxx * backwardEnergy) {
527  pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
528  backwardKinetic += kineticEnergy;
529  pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
530  pseudoParticle[6].SetMomentum(0.0); // set the z-momentum
531  phi = pseudoParticle[6].Angle(pseudoParticle[8]);
532  if (pseudoParticle[6].GetMomentum().y() / MeV < 0.0)
533  phi = twopi - phi;
534  phi += pi + normal() * pi / 12.0;
535  if (phi > twopi)
536  phi -= twopi;
537  if (phi < 0.0)
538  phi = twopi - phi;
539  outerCounter = 2; // leave outer loop
540  eliminateThisParticle = false; // don't eliminate this particle
541  resetEnergies = false;
542  break; // leave inner loop
543  }
544  if (innerCounter > 5)
545  break; // leave inner loop
546  if (forwardEnergy >= vecMass) // switch sides
547  {
548  vec[i]->SetSide(1);
549  forwardEnergy -= vecMass;
550  backwardEnergy += vecMass;
551  backwardCount--;
552  }
553  }
554  G4ThreeVector momentum = vec[i]->GetMomentum();
555  vec[i]->SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
556  pt *= 0.9;
557  dndl[19] *= 0.9;
558  } // closes inner loop
559  if (resetEnergies) {
560  // if we get to here, the inner loop has been Done 6 Times
561  // reset the kinetic energies of previously Done particles, if they are lighter
562  // than protons and in the forward hemisphere
563  // then continue with outer loop
564  //
565  forwardKinetic = 0.0;
566  backwardKinetic = 0.0;
567  pseudoParticle[4].SetZero();
568  pseudoParticle[5].SetZero();
569  for (l = i + 1; l < vecLen; ++l) {
570  if (vec[l]->GetSide() > 0 || vec[l]->GetDefinition() == aKaonMinus || vec[l]->GetDefinition() == aKaonZeroL ||
571  vec[l]->GetDefinition() == aKaonZeroS || vec[l]->GetDefinition() == aKaonPlus ||
572  vec[l]->GetDefinition() == aPiMinus || vec[l]->GetDefinition() == aPiZero ||
573  vec[l]->GetDefinition() == aPiPlus) {
574  G4double tempMass = vec[l]->GetMass() / MeV;
575  totalEnergy = 0.95 * vec[l]->GetTotalEnergy() / MeV + 0.05 * tempMass;
576  totalEnergy = std::max(tempMass, totalEnergy);
577  vec[l]->SetTotalEnergy(totalEnergy * MeV);
578  pp = std::sqrt(std::abs(totalEnergy * totalEnergy - tempMass * tempMass));
579  pp1 = vec[l]->GetMomentum().mag() / MeV;
580  if (pp1 < 1.0e-6 * GeV) {
581  G4double rthnve = pi * G4UniformRand();
582  G4double phinve = twopi * G4UniformRand();
583  G4double srth = std::sin(rthnve);
584  vec[l]->SetMomentum(
585  pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
586  } else {
587  vec[l]->SetMomentum(vec[l]->GetMomentum() * (pp / pp1));
588  }
589  G4double px = vec[l]->GetMomentum().x() / MeV;
590  G4double py = vec[l]->GetMomentum().y() / MeV;
591  pt = std::max(1.0, std::sqrt(px * px + py * py)) / GeV;
592  if (vec[l]->GetSide() > 0) {
593  forwardKinetic += vec[l]->GetKineticEnergy() / GeV;
594  pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
595  } else {
596  backwardKinetic += vec[l]->GetKineticEnergy() / GeV;
597  pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
598  }
599  }
600  }
601  }
602  } // closes outer loop
603 
604  if (eliminateThisParticle && vec[i]->GetMayBeKilled()) // not enough energy, eliminate this particle
605  {
606  if (vec[i]->GetSide() > 0) {
607  --forwardCount;
608  forwardEnergy += vecMass;
609  } else {
610  if (vec[i]->GetSide() == -2) {
611  --extraNucleonCount;
612  extraNucleonMass -= vecMass;
613  backwardEnergy -= vecMass;
614  }
615  --backwardCount;
616  backwardEnergy += vecMass;
617  }
618  for (G4int j = i; j < (vecLen - 1); ++j)
619  *vec[j] = *vec[j + 1]; // shift up
620  G4ReactionProduct *temp = vec[vecLen - 1];
621  delete temp;
622  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
623  if (--vecLen == 0)
624  return false; // all the secondaries have been eliminated
625  pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
626  pseudoParticle[6].SetMomentum(0.0); // set z-momentum
627  }
628  } // closes main for loop
629 
630  //
631  // for the incident particle: it was placed in the forward hemisphere
632  // set pt and phi values, they are changed somewhat in the iteration loop
633  // set mass parameter for lambda fragmentation model
634  //
635  G4double phi = G4UniformRand() * twopi;
636  G4double ran = -std::log(1.0 - G4UniformRand());
637  if (currentParticle.GetDefinition() == aPiMinus || currentParticle.GetDefinition() == aPiZero ||
638  currentParticle.GetDefinition() == aPiPlus) {
639  aspar = 0.60;
640  pt = std::sqrt(std::pow(ran / 6.0, 1.7));
641  } else if (currentParticle.GetDefinition() == aKaonMinus || currentParticle.GetDefinition() == aKaonZeroL ||
642  currentParticle.GetDefinition() == aKaonZeroS || currentParticle.GetDefinition() == aKaonPlus) {
643  aspar = 0.50;
644  pt = std::sqrt(std::pow(ran / 5.0, 1.4));
645  } else {
646  aspar = 0.40;
647  pt = std::sqrt(std::pow(ran / 4.0, 1.2));
648  }
649  for (G4int j = 0; j < 20; ++j)
650  binl[j] = j / (19. * pt);
651  currentParticle.SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
652  et = pseudoParticle[0].GetTotalEnergy() / GeV;
653  dndl[0] = 0.0;
654  vecMass = currentParticle.GetMass() / GeV;
655  for (l = 1; l < 20; ++l) {
656  x = (binl[l] + binl[l - 1]) / 2.;
657  if (x > 1.0 / pt)
658  dndl[l] += dndl[l - 1]; // changed from just = on 02 April 98
659  else
660  dndl[l] = aspar / std::sqrt(std::pow((1. + sqr(aspar * x)), 3)) * (binl[l] - binl[l - 1]) * et /
661  std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
662  dndl[l - 1];
663  }
664  ran = G4UniformRand() * dndl[19];
665  l = 1;
666  while ((ran > dndl[l]) && (l < 20))
667  l++;
668  l = std::min(19, l);
669  x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
670  currentParticle.SetMomentum(x * et * GeV); // set the z-momentum
671  if (forwardEnergy < forwardKinetic)
672  totalEnergy = vecMass + 0.04 * std::fabs(normal());
673  else
674  totalEnergy = vecMass + forwardEnergy - forwardKinetic;
675  currentParticle.SetTotalEnergy(totalEnergy * GeV);
676  pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
677  pp1 = currentParticle.GetMomentum().mag() / MeV;
678  if (pp1 < 1.0e-6 * GeV) {
679  G4double rthnve = pi * G4UniformRand();
680  G4double phinve = twopi * G4UniformRand();
681  G4double srth = std::sin(rthnve);
682  currentParticle.SetMomentum(
683  pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
684  } else {
685  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
686  }
687  pseudoParticle[4] = pseudoParticle[4] + currentParticle;
688  //
689  // this finishes the current particle
690  // now for the target particle
691  //
692  if (backwardNucleonCount < 18) {
693  targetParticle.SetSide(-3);
694  ++backwardNucleonCount;
695  } else {
696  // set pt and phi values, they are changed somewhat in the iteration loop
697  // set mass parameter for lambda fragmentation model
698  //
699  vecMass = targetParticle.GetMass() / GeV;
700  ran = -std::log(1.0 - G4UniformRand());
701  aspar = 0.40;
702  pt = std::max(0.001, std::sqrt(std::pow(ran / 4.0, 1.2)));
703  targetParticle.SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
704  for (G4int j = 0; j < 20; ++j)
705  binl[j] = (j - 1.) / (19. * pt);
706  et = pseudoParticle[1].GetTotalEnergy() / GeV;
707  dndl[0] = 0.0;
708  outerCounter = 0;
709  resetEnergies = true;
710  while (++outerCounter < 3) // start of outer iteration loop
711  {
712  for (l = 1; l < 20; ++l) {
713  x = (binl[l] + binl[l - 1]) / 2.;
714  if (x > 1.0 / pt)
715  dndl[l] += dndl[l - 1]; // changed from just = on 02 April 98
716  else
717  dndl[l] = aspar / std::sqrt(std::pow((1. + aspar * x * aspar * x), 3)) * (binl[l] - binl[l - 1]) * et /
718  std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
719  dndl[l - 1];
720  }
721  innerCounter = 0;
722  while (++innerCounter < 7) // start of inner iteration loop
723  {
724  l = 1;
725  ran = G4UniformRand() * dndl[19];
726  while ((ran >= dndl[l]) && (l < 20))
727  l++;
728  l = std::min(19, l);
729  x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
730  if (targetParticle.GetSide() < 0)
731  x *= -1.;
732  targetParticle.SetMomentum(x * et * GeV); // set the z-momentum
733  totalEnergy = std::sqrt(x * et * x * et + pt * pt + vecMass * vecMass);
734  targetParticle.SetTotalEnergy(totalEnergy * GeV);
735  if (targetParticle.GetSide() < 0) {
736  G4double xxx = 0.95 + 0.05 * extraNucleonCount / 20.0;
737  if ((backwardKinetic + totalEnergy - vecMass) < xxx * backwardEnergy) {
738  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
739  backwardKinetic += totalEnergy - vecMass;
740  pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
741  pseudoParticle[6].SetMomentum(0.0); // set z-momentum
742  outerCounter = 2; // leave outer loop
743  resetEnergies = false;
744  break; // leave inner loop
745  }
746  if (innerCounter > 5)
747  break; // leave inner loop
748  if (forwardEnergy >= vecMass) // switch sides
749  {
750  targetParticle.SetSide(1);
751  forwardEnergy -= vecMass;
752  backwardEnergy += vecMass;
753  --backwardCount;
754  }
755  G4ThreeVector momentum = targetParticle.GetMomentum();
756  targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
757  pt *= 0.9;
758  dndl[19] *= 0.9;
759  } else // target has gone to forward side
760  {
761  if (forwardEnergy < forwardKinetic)
762  totalEnergy = vecMass + 0.04 * std::fabs(normal());
763  else
764  totalEnergy = vecMass + forwardEnergy - forwardKinetic;
765  targetParticle.SetTotalEnergy(totalEnergy * GeV);
766  pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
767  pp1 = targetParticle.GetMomentum().mag() / MeV;
768  if (pp1 < 1.0e-6 * GeV) {
769  G4double rthnve = pi * G4UniformRand();
770  G4double phinve = twopi * G4UniformRand();
771  G4double srth = std::sin(rthnve);
772  targetParticle.SetMomentum(
773  pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
774  } else
775  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
776 
777  pseudoParticle[4] = pseudoParticle[4] + targetParticle;
778  outerCounter = 2; // leave outer loop
779  //eliminateThisParticle = false; // don't eliminate this particle
780  resetEnergies = false;
781  break; // leave inner loop
782  }
783  } // closes inner loop
784  if (resetEnergies) {
785  // if we get to here, the inner loop has been Done 6 Times
786  // reset the kinetic energies of previously Done particles, if they are lighter
787  // than protons and in the forward hemisphere
788  // then continue with outer loop
789 
790  forwardKinetic = backwardKinetic = 0.0;
791  pseudoParticle[4].SetZero();
792  pseudoParticle[5].SetZero();
793  for (l = 0; l < vecLen; ++l) // changed from l=1 on 02 April 98
794  {
795  if (vec[l]->GetSide() > 0 || vec[l]->GetDefinition() == aKaonMinus || vec[l]->GetDefinition() == aKaonZeroL ||
796  vec[l]->GetDefinition() == aKaonZeroS || vec[l]->GetDefinition() == aKaonPlus ||
797  vec[l]->GetDefinition() == aPiMinus || vec[l]->GetDefinition() == aPiZero ||
798  vec[l]->GetDefinition() == aPiPlus) {
799  G4double tempMass = vec[l]->GetMass() / GeV;
800  totalEnergy = std::max(tempMass, 0.95 * vec[l]->GetTotalEnergy() / GeV + 0.05 * tempMass);
801  vec[l]->SetTotalEnergy(totalEnergy * GeV);
802  pp = std::sqrt(std::abs(totalEnergy * totalEnergy - tempMass * tempMass)) * GeV;
803  pp1 = vec[l]->GetMomentum().mag() / MeV;
804  if (pp1 < 1.0e-6 * GeV) {
805  G4double rthnve = pi * G4UniformRand();
806  G4double phinve = twopi * G4UniformRand();
807  G4double srth = std::sin(rthnve);
808  vec[l]->SetMomentum(
809  pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
810  } else
811  vec[l]->SetMomentum(vec[l]->GetMomentum() * (pp / pp1));
812 
813  pt = std::max(0.001 * GeV,
814  std::sqrt(sqr(vec[l]->GetMomentum().x() / MeV) + sqr(vec[l]->GetMomentum().y() / MeV))) /
815  GeV;
816  if (vec[l]->GetSide() > 0) {
817  forwardKinetic += vec[l]->GetKineticEnergy() / GeV;
818  pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
819  } else {
820  backwardKinetic += vec[l]->GetKineticEnergy() / GeV;
821  pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
822  }
823  }
824  }
825  }
826  } // closes outer loop
827  }
828  //
829  // this finishes the target particle
830  // backward nucleons produced with a cluster model
831  //
832  pseudoParticle[6].Lorentz(pseudoParticle[3], pseudoParticle[2]);
833  pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
834  pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
835  if (backwardNucleonCount == 1) // target particle is the only backward nucleon
836  {
837  G4double ekin = std::min(backwardEnergy - backwardKinetic, centerofmassEnergy / 2.0 - protonMass / GeV);
838  if (ekin < 0.04)
839  ekin = 0.04 * std::fabs(normal());
840  vecMass = targetParticle.GetMass() / GeV;
841  totalEnergy = ekin + vecMass;
842  targetParticle.SetTotalEnergy(totalEnergy * GeV);
843  pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
844  pp1 = pseudoParticle[6].GetMomentum().mag() / MeV;
845  if (pp1 < 1.0e-6 * GeV) {
846  rthnve = pi * G4UniformRand();
847  phinve = twopi * G4UniformRand();
848  G4double srth = std::sin(rthnve);
849  targetParticle.SetMomentum(
850  pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
851  } else {
852  targetParticle.SetMomentum(pseudoParticle[6].GetMomentum() * (pp / pp1));
853  }
854  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
855  } else // more than one backward nucleon
856  {
857  const G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
858  const G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
859  // Replaced the following min function to get correct behaviour on DEC.
860  G4int tempCount = std::max(1, std::min(5, backwardNucleonCount)) - 1;
861  //cout << "backwardNucleonCount " << backwardNucleonCount << G4endl;
862  //cout << "tempCount " << tempCount << G4endl;
863  G4double rmb0 = 0.0;
864  if (targetParticle.GetSide() == -3)
865  rmb0 += targetParticle.GetMass() / GeV;
866  for (i = 0; i < vecLen; ++i) {
867  if (vec[i]->GetSide() == -3)
868  rmb0 += vec[i]->GetMass() / GeV;
869  }
870  rmb = rmb0 + std::pow(-std::log(1.0 - G4UniformRand()), cpar[tempCount]) / gpar[tempCount];
871  totalEnergy = pseudoParticle[6].GetTotalEnergy() / GeV;
872  vecMass = std::min(rmb, totalEnergy);
873  pseudoParticle[6].SetMass(vecMass * GeV);
874  pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
875  pp1 = pseudoParticle[6].GetMomentum().mag() / MeV;
876  if (pp1 < 1.0e-6 * GeV) {
877  rthnve = pi * G4UniformRand();
878  phinve = twopi * G4UniformRand();
879  G4double srth = std::sin(rthnve);
880  pseudoParticle[6].SetMomentum(
881  -pp * srth * std::cos(phinve) * MeV, -pp * srth * std::sin(phinve) * MeV, -pp * std::cos(rthnve) * MeV);
882  } else
883  pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp / pp1));
884 
885  G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV; // tempV contains the backward nucleons
886  tempV.Initialize(backwardNucleonCount);
887  G4int tempLen = 0;
888  if (targetParticle.GetSide() == -3)
889  tempV.SetElement(tempLen++, &targetParticle);
890  for (i = 0; i < vecLen; ++i) {
891  if (vec[i]->GetSide() == -3)
892  tempV.SetElement(tempLen++, vec[i]);
893  }
894  if (tempLen != backwardNucleonCount) {
895  G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
896  G4cerr << "tempLen = " << tempLen;
897  G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
898  G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
899  G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
900  for (i = 0; i < vecLen; ++i)
901  G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
902  exit(EXIT_FAILURE);
903  }
904  constantCrossSection = true;
905  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
906  if (tempLen >= 2) {
907  GenerateNBodyEvent(pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen);
908  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
909  if (targetParticle.GetSide() == -3) {
910  targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
911  // tempV contains the real stuff
912  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
913  }
914  for (i = 0; i < vecLen; ++i) {
915  if (vec[i]->GetSide() == -3) {
916  vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
917  pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
918  }
919  }
920  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
921  }
922  }
923  //
924  // Lorentz transformation in lab system
925  //
926  if (vecLen == 0)
927  return false; // all the secondaries have been eliminated
928  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
929 
930  G4int numberofFinalStateNucleons = 0;
931  if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron ||
932  currentParticle.GetDefinition() == G4SigmaMinus::SigmaMinus() ||
933  currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
934  currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
935  currentParticle.GetDefinition() == G4XiZero::XiZero() ||
936  currentParticle.GetDefinition() == G4XiMinus::XiMinus() ||
937  currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
938  currentParticle.GetDefinition() == G4Lambda::Lambda())
939  ++numberofFinalStateNucleons;
940  currentParticle.Lorentz(currentParticle, pseudoParticle[1]);
941 
942  if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron ||
943  targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() ||
944  targetParticle.GetDefinition() == G4XiMinus::XiMinus() ||
945  targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
946  targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
947  targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
948  targetParticle.GetDefinition() == G4SigmaMinus::SigmaMinus())
949  ++numberofFinalStateNucleons;
950  if (targetParticle.GetDefinition() == G4AntiProton::AntiProton())
951  --numberofFinalStateNucleons;
952  if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron())
953  --numberofFinalStateNucleons;
954  if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
955  --numberofFinalStateNucleons;
956  if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
957  --numberofFinalStateNucleons;
958  if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
959  --numberofFinalStateNucleons;
960  if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero())
961  --numberofFinalStateNucleons;
962  if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus())
963  --numberofFinalStateNucleons;
964  if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
965  --numberofFinalStateNucleons;
966  if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda())
967  --numberofFinalStateNucleons;
968  targetParticle.Lorentz(targetParticle, pseudoParticle[1]);
969 
970  for (i = 0; i < vecLen; ++i) {
971  if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron ||
972  vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() ||
973  vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
974  vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() ||
975  vec[i]->GetDefinition() == G4SigmaMinus::SigmaMinus())
976  ++numberofFinalStateNucleons;
977  if (vec[i]->GetDefinition() == G4AntiProton::AntiProton())
978  --numberofFinalStateNucleons;
979  if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron())
980  --numberofFinalStateNucleons;
981  if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
982  --numberofFinalStateNucleons;
983  if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
984  --numberofFinalStateNucleons;
985  if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
986  --numberofFinalStateNucleons;
987  if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda())
988  --numberofFinalStateNucleons;
989  if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero())
990  --numberofFinalStateNucleons;
991  if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus())
992  --numberofFinalStateNucleons;
993  if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
994  --numberofFinalStateNucleons;
995  vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
996  }
997  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
998  if (veryForward)
999  numberofFinalStateNucleons++;
1000  numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
1001  //
1002  // leadFlag will be true
1003  // iff original particle is at least as heavy as K+ and not a proton or neutron AND
1004  // if
1005  // incident particle is at least as heavy as K+ and it is not a proton or neutron
1006  // leadFlag is set to the incident particle
1007  // or
1008  // target particle is at least as heavy as K+ and it is not a proton or neutron
1009  // leadFlag is set to the target particle
1010  //
1011  G4bool leadingStrangeParticleHasChanged = true;
1012  if (leadFlag) {
1013  if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1014  leadingStrangeParticleHasChanged = false;
1015  if (leadingStrangeParticleHasChanged && (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()))
1016  leadingStrangeParticleHasChanged = false;
1017  if (leadingStrangeParticleHasChanged) {
1018  for (i = 0; i < vecLen; i++) {
1019  if (vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition()) {
1020  leadingStrangeParticleHasChanged = false;
1021  break;
1022  }
1023  }
1024  }
1025  if (leadingStrangeParticleHasChanged) {
1026  G4bool leadTest =
1027  (leadingStrangeParticle.GetDefinition() == aKaonMinus ||
1028  leadingStrangeParticle.GetDefinition() == aKaonZeroL ||
1029  leadingStrangeParticle.GetDefinition() == aKaonZeroS ||
1030  leadingStrangeParticle.GetDefinition() == aKaonPlus || leadingStrangeParticle.GetDefinition() == aPiMinus ||
1031  leadingStrangeParticle.GetDefinition() == aPiZero || leadingStrangeParticle.GetDefinition() == aPiPlus);
1032  G4bool targetTest = false;
1033 
1034  // following modified by JLC 22-Oct-97
1035 
1036  if ((leadTest && targetTest) || !(leadTest || targetTest)) // both true or both false
1037  {
1038  targetParticle.SetDefinitionAndUpdateE(leadingStrangeParticle.GetDefinition());
1039  targetHasChanged = true;
1040  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1041  } else {
1042  currentParticle.SetDefinitionAndUpdateE(leadingStrangeParticle.GetDefinition());
1043  incidentHasChanged = false;
1044  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1045  }
1046  }
1047  } // end of if( leadFlag )
1048 
1049  pseudoParticle[3].SetMomentum(0.0, 0.0, pOriginal * GeV);
1050  pseudoParticle[3].SetMass(mOriginal * GeV);
1051  pseudoParticle[3].SetTotalEnergy(std::sqrt(pOriginal * pOriginal + mOriginal * mOriginal) * GeV);
1052 
1053  const G4ParticleDefinition *aOrgDef = modifiedOriginal.GetDefinition();
1054  G4int diff = 0;
1055  if (aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron())
1056  diff = 1;
1057  if (numberofFinalStateNucleons == 1)
1058  diff = 0;
1059  pseudoParticle[4].SetMomentum(0.0, 0.0, 0.0);
1060  pseudoParticle[4].SetMass(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1061  pseudoParticle[4].SetTotalEnergy(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1062 
1063  G4double theoreticalKinetic = pseudoParticle[3].GetTotalEnergy() / MeV + pseudoParticle[4].GetTotalEnergy() / MeV -
1064  currentParticle.GetMass() / MeV - targetParticle.GetMass() / MeV;
1065 
1066  G4double simulatedKinetic = currentParticle.GetKineticEnergy() / MeV + targetParticle.GetKineticEnergy() / MeV;
1067 
1068  pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
1069  pseudoParticle[3].Lorentz(pseudoParticle[3], pseudoParticle[5]);
1070  pseudoParticle[4].Lorentz(pseudoParticle[4], pseudoParticle[5]);
1071 
1072  pseudoParticle[7].SetZero();
1073  pseudoParticle[7] = pseudoParticle[7] + currentParticle;
1074  pseudoParticle[7] = pseudoParticle[7] + targetParticle;
1075 
1076  for (i = 0; i < vecLen; ++i) {
1077  pseudoParticle[7] = pseudoParticle[7] + *vec[i];
1078  simulatedKinetic += vec[i]->GetKineticEnergy() / MeV;
1079  theoreticalKinetic -= vec[i]->GetMass() / MeV;
1080  }
1081  if (vecLen <= 16 && vecLen > 0) {
1082  // must create a new set of ReactionProducts here because GenerateNBody will
1083  // modify the momenta for the particles, and we don't want to do this
1084  //
1085  G4ReactionProduct tempR[130];
1086  tempR[0] = currentParticle;
1087  tempR[1] = targetParticle;
1088  for (i = 0; i < vecLen; ++i)
1089  tempR[i + 2] = *vec[i];
1090  G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1091  tempV.Initialize(vecLen + 2);
1092  G4int tempLen = 0;
1093  for (i = 0; i < vecLen + 2; ++i)
1094  tempV.SetElement(tempLen++, &tempR[i]);
1095  constantCrossSection = true;
1096 
1097  wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() / MeV + pseudoParticle[4].GetTotalEnergy() / MeV,
1098  constantCrossSection,
1099  tempV,
1100  tempLen);
1101  if (wgt > -.5) {
1102  theoreticalKinetic = 0.0;
1103  for (i = 0; i < tempLen; ++i) {
1104  pseudoParticle[6].Lorentz(*tempV[i], pseudoParticle[4]);
1105  theoreticalKinetic += pseudoParticle[6].GetKineticEnergy() / MeV;
1106  }
1107  }
1108  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1109  }
1110  //
1111  // Make sure, that the kinetic energies are correct
1112  //
1113  if (simulatedKinetic != 0.0) {
1114  wgt = (theoreticalKinetic) / simulatedKinetic;
1115  theoreticalKinetic = currentParticle.GetKineticEnergy() / MeV * wgt;
1116  simulatedKinetic = theoreticalKinetic;
1117  currentParticle.SetKineticEnergy(theoreticalKinetic * MeV);
1118  pp = currentParticle.GetTotalMomentum() / MeV;
1119  pp1 = currentParticle.GetMomentum().mag() / MeV;
1120  if (pp1 < 1.0e-6 * GeV) {
1121  rthnve = pi * G4UniformRand();
1122  phinve = twopi * G4UniformRand();
1123  currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1124  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1125  pp * std::cos(rthnve) * MeV);
1126  } else {
1127  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1128  }
1129  theoreticalKinetic = targetParticle.GetKineticEnergy() / MeV * wgt;
1130  targetParticle.SetKineticEnergy(theoreticalKinetic * MeV);
1131  simulatedKinetic += theoreticalKinetic;
1132  pp = targetParticle.GetTotalMomentum() / MeV;
1133  pp1 = targetParticle.GetMomentum().mag() / MeV;
1134  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1135  if (pp1 < 1.0e-6 * GeV) {
1136  rthnve = pi * G4UniformRand();
1137  phinve = twopi * G4UniformRand();
1138  targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1139  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1140  pp * std::cos(rthnve) * MeV);
1141  } else {
1142  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1143  }
1144  for (i = 0; i < vecLen; ++i) {
1145  theoreticalKinetic = vec[i]->GetKineticEnergy() / MeV * wgt;
1146  simulatedKinetic += theoreticalKinetic;
1147  vec[i]->SetKineticEnergy(theoreticalKinetic * MeV);
1148  pp = vec[i]->GetTotalMomentum() / MeV;
1149  pp1 = vec[i]->GetMomentum().mag() / MeV;
1150  if (pp1 < 1.0e-6 * GeV) {
1151  rthnve = pi * G4UniformRand();
1152  phinve = twopi * G4UniformRand();
1153  vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1154  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1155  pp * std::cos(rthnve) * MeV);
1156  } else
1157  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
1158  }
1159  }
1160  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1161  Rotate(numberofFinalStateNucleons,
1162  pseudoParticle[3].GetMomentum(),
1163  modifiedOriginal,
1164  originalIncident,
1165  targetNucleus,
1166  currentParticle,
1167  targetParticle,
1168  vec,
1169  vecLen);
1170  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1171  //
1172  // add black track particles
1173  // the total number of particles produced is restricted to 198
1174  // this may have influence on very high energies
1175  //
1176  if (atomicWeight >= 1.5) {
1177  // npnb is number of proton/neutron black track particles
1178  // ndta is the number of deuterons, tritons, and alphas produced
1179  // epnb is the kinetic energy available for proton/neutron black track particles
1180  // edta is the kinetic energy available for deuteron/triton/alpha particles
1181  //
1182  G4double epnb, edta;
1183  G4int npnb = 0;
1184  G4int ndta = 0;
1185 
1186  epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
1187  edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
1188  const G4double pnCutOff = 0.001;
1189  const G4double dtaCutOff = 0.001;
1190  const G4double kineticMinimum = 1.e-6;
1191  const G4double kineticFactor = -0.010;
1192  G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
1193  const G4double ekIncident = originalIncident->GetKineticEnergy() / GeV;
1194  if (ekIncident >= 5.0)
1195  sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0));
1196  if (epnb >= pnCutOff) {
1197  npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta));
1198  if (numberofFinalStateNucleons + npnb > atomicWeight)
1199  npnb = G4int(atomicWeight + 0.00001 - numberofFinalStateNucleons);
1200  npnb = std::min(npnb, 127 - vecLen);
1201  }
1202  if (edta >= dtaCutOff) {
1203  ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta));
1204  ndta = std::min(ndta, 127 - vecLen);
1205  }
1206  G4double spall = numberofFinalStateNucleons;
1207  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1208 
1209  AddBlackTrackParticles(epnb,
1210  npnb,
1211  edta,
1212  ndta,
1213  sprob,
1214  kineticMinimum,
1215  kineticFactor,
1216  modifiedOriginal,
1217  spall,
1218  targetNucleus,
1219  vec,
1220  vecLen);
1221 
1222  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1223  }
1224  //
1225  // calculate time delay for nuclear reactions
1226  //
1227  if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
1228  currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
1229  else
1230  currentParticle.SetTOF(1.0);
1231  return true;
1232  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1233 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
pi
const Double_t pi
Definition: trackSplitPlot.h:36
MeV
const double MeV
ecalTB2006H4_GenSimDigiReco_cfg.G4cout
G4cout
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:303
ecalTB2006H4_GenSimDigiReco_cfg.G4cerr
G4cerr
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:304
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
FullModelReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2340
FullModelReactionDynamics::AddBlackTrackParticles
void AddBlackTrackParticles(const G4double epnb, const G4int npnb, const G4double edta, const G4int ndta, const G4double sprob, const G4double kineticMinimum, const G4double kineticFactor, const G4ReactionProduct &modifiedOriginal, G4double spall, const G4Nucleus &aNucleus, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2850
FullModelReactionDynamics::Rotate
void Rotate(const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2639
stringResolutionProvider_cfi.et
et
define resolution functions of each parameter
Definition: stringResolutionProvider_cfi.py:13
funct::sqr
Square< F >::type sqr(const F &f)
Definition: Square.h:13
BPhysicsValidation_cfi.Lambda
Lambda
Definition: BPhysicsValidation_cfi.py:63
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
phase2TrackerDigitizer_cfi.SigmaZero
SigmaZero
Definition: phase2TrackerDigitizer_cfi.py:46
void FullModelReactionDynamics::MomentumCheck ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2998 of file FullModelReactionDynamics.cc.

References mps_fire::i, MeV, and mathSSE::sqrt().

Referenced by FindTimeDelay().

3002  {
3003  const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / MeV;
3004  G4double testMomentum = currentParticle.GetMomentum().mag() / MeV;
3005  G4double pMass;
3006  if (testMomentum >= pOriginal) {
3007  pMass = currentParticle.GetMass() / MeV;
3008  currentParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3009  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pOriginal / testMomentum));
3010  }
3011  testMomentum = targetParticle.GetMomentum().mag() / MeV;
3012  if (testMomentum >= pOriginal) {
3013  pMass = targetParticle.GetMass() / MeV;
3014  targetParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3015  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pOriginal / testMomentum));
3016  }
3017  for (G4int i = 0; i < vecLen; ++i) {
3018  testMomentum = vec[i]->GetMomentum().mag() / MeV;
3019  if (testMomentum >= pOriginal) {
3020  pMass = vec[i]->GetMass() / MeV;
3021  vec[i]->SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3022  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pOriginal / testMomentum));
3023  }
3024  }
3025 }
mps_fire.i
i
Definition: mps_fire.py:338
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
G4double FullModelReactionDynamics::normal ( )
private

Definition at line 2525 of file FullModelReactionDynamics.cc.

References mps_fire::i, and summarizeEdmComparisonLogfiles::ran.

Referenced by FindTimeDelay().

2525  {
2526  G4double ran = -6.0;
2527  for (G4int i = 0; i < 12; ++i)
2528  ran += G4UniformRand();
2529  return ran;
2530 }
mps_fire.i
i
Definition: mps_fire.py:338
void FullModelReactionDynamics::NuclearReaction ( G4FastVector< G4ReactionProduct, 4 > &  vec,
G4int &  vecLen,
const G4HadProjectile *  originalIncident,
const G4Nucleus &  aNucleus,
const G4double  theAtomicMass,
const G4double *  massVec 
)

Definition at line 3394 of file FullModelReactionDynamics.cc.

References patCaloMETCorrections_cff::A, funct::cos(), Gamma, GeV, mps_fire::i, ke, SiStripPI::max, MeV, min(), nt, AlCaHLTBitMon_ParallelJobs::p, createTree::pp, summarizeEdmComparisonLogfiles::ran, funct::sin(), and findQualityFiles::v.

Referenced by FindTimeDelay().

3399  {
3400  // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
3401  //
3402  // Nuclear reaction kinematics at low energies
3403  //
3404  G4ParticleDefinition *aGamma = G4Gamma::Gamma();
3405  G4ParticleDefinition *aProton = G4Proton::Proton();
3406  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3407  G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3408  G4ParticleDefinition *aTriton = G4Triton::Triton();
3409  G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3410 
3411  const G4double aProtonMass = aProton->GetPDGMass() / MeV;
3412  const G4double aNeutronMass = aNeutron->GetPDGMass() / MeV;
3413  const G4double aDeuteronMass = aDeuteron->GetPDGMass() / MeV;
3414  const G4double aTritonMass = aTriton->GetPDGMass() / MeV;
3415  const G4double anAlphaMass = anAlpha->GetPDGMass() / MeV;
3416 
3417  G4ReactionProduct currentParticle;
3418  currentParticle = *originalIncident;
3419  //
3420  // Set beam particle, take kinetic energy of current particle as the
3421  // fundamental quantity. Due to the difficult kinematic, all masses have to
3422  // be assigned the best measured values
3423  //
3424  G4double p = currentParticle.GetTotalMomentum();
3425  G4double pp = currentParticle.GetMomentum().mag();
3426  if (pp <= 0.001 * MeV) {
3427  G4double phinve = twopi * G4UniformRand();
3428  G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3429  currentParticle.SetMomentum(
3430  p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3431  } else
3432  currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3433  //
3434  // calculate Q-value of reactions
3435  //
3436  G4double currentKinetic = currentParticle.GetKineticEnergy() / MeV;
3437  G4double currentMass = currentParticle.GetDefinition()->GetPDGMass() / MeV;
3438  G4double qv = currentKinetic + theAtomicMass + currentMass;
3439 
3440  G4double qval[9];
3441  qval[0] = qv - mass[0];
3442  qval[1] = qv - mass[1] - aNeutronMass;
3443  qval[2] = qv - mass[2] - aProtonMass;
3444  qval[3] = qv - mass[3] - aDeuteronMass;
3445  qval[4] = qv - mass[4] - aTritonMass;
3446  qval[5] = qv - mass[5] - anAlphaMass;
3447  qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3448  qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3449  qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3450 
3451  if (currentParticle.GetDefinition() == aNeutron) {
3452  const G4double A = targetNucleus.GetN_asInt(); // atomic weight
3453  if (G4UniformRand() > ((A - 1.0) / 230.0) * ((A - 1.0) / 230.0))
3454  qval[0] = 0.0;
3455  if (G4UniformRand() >= currentKinetic / 7.9254 * A)
3456  qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3457  } else
3458  qval[0] = 0.0;
3459 
3460  G4int i;
3461  qv = 0.0;
3462  for (i = 0; i < 9; ++i) {
3463  if (mass[i] < 500.0 * MeV)
3464  qval[i] = 0.0;
3465  if (qval[i] < 0.0)
3466  qval[i] = 0.0;
3467  qv += qval[i];
3468  }
3469  G4double qv1 = 0.0;
3470  G4double ran = G4UniformRand();
3471  G4int index;
3472  for (index = 0; index < 9; ++index) {
3473  if (qval[index] > 0.0) {
3474  qv1 += qval[index] / qv;
3475  if (ran <= qv1)
3476  break;
3477  }
3478  }
3479  if (index == 9) // loop continued to the end
3480  {
3481  throw G4HadronicException(
3482  __FILE__,
3483  __LINE__,
3484  "FullModelReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3485  }
3486  G4double ke = currentParticle.GetKineticEnergy() / GeV;
3487  G4int nt = 2;
3488  if ((index >= 6) || (G4UniformRand() < std::min(0.5, ke * 10.0)))
3489  nt = 3;
3490 
3491  G4ReactionProduct **v = new G4ReactionProduct *[3];
3492  v[0] = new G4ReactionProduct;
3493  v[1] = new G4ReactionProduct;
3494  v[2] = new G4ReactionProduct;
3495 
3496  v[0]->SetMass(mass[index] * MeV);
3497  switch (index) {
3498  case 0:
3499  v[1]->SetDefinition(aGamma);
3500  v[2]->SetDefinition(aGamma);
3501  break;
3502  case 1:
3503  v[1]->SetDefinition(aNeutron);
3504  v[2]->SetDefinition(aGamma);
3505  break;
3506  case 2:
3507  v[1]->SetDefinition(aProton);
3508  v[2]->SetDefinition(aGamma);
3509  break;
3510  case 3:
3511  v[1]->SetDefinition(aDeuteron);
3512  v[2]->SetDefinition(aGamma);
3513  break;
3514  case 4:
3515  v[1]->SetDefinition(aTriton);
3516  v[2]->SetDefinition(aGamma);
3517  break;
3518  case 5:
3519  v[1]->SetDefinition(anAlpha);
3520  v[2]->SetDefinition(aGamma);
3521  break;
3522  case 6:
3523  v[1]->SetDefinition(aNeutron);
3524  v[2]->SetDefinition(aNeutron);
3525  break;
3526  case 7:
3527  v[1]->SetDefinition(aNeutron);
3528  v[2]->SetDefinition(aProton);
3529  break;
3530  case 8:
3531  v[1]->SetDefinition(aProton);
3532  v[2]->SetDefinition(aProton);
3533  break;
3534  }
3535  //
3536  // calculate centre of mass energy
3537  //
3538  G4ReactionProduct pseudo1;
3539  pseudo1.SetMass(theAtomicMass * MeV);
3540  pseudo1.SetTotalEnergy(theAtomicMass * MeV);
3541  G4ReactionProduct pseudo2 = currentParticle + pseudo1;
3542  pseudo2.SetMomentum(pseudo2.GetMomentum() * (-1.0));
3543  //
3544  // use phase space routine in centre of mass system
3545  //
3546  G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
3547  tempV.Initialize(nt);
3548  G4int tempLen = 0;
3549  tempV.SetElement(tempLen++, v[0]);
3550  tempV.SetElement(tempLen++, v[1]);
3551  if (nt == 3)
3552  tempV.SetElement(tempLen++, v[2]);
3553  G4bool constantCrossSection = true;
3554  GenerateNBodyEvent(pseudo2.GetMass() / MeV, constantCrossSection, tempV, tempLen);
3555  v[0]->Lorentz(*v[0], pseudo2);
3556  v[1]->Lorentz(*v[1], pseudo2);
3557  if (nt == 3)
3558  v[2]->Lorentz(*v[2], pseudo2);
3559 
3560  G4bool particleIsDefined = false;
3561  if (v[0]->GetMass() / MeV - aProtonMass < 0.1) {
3562  v[0]->SetDefinition(aProton);
3563  particleIsDefined = true;
3564  } else if (v[0]->GetMass() / MeV - aNeutronMass < 0.1) {
3565  v[0]->SetDefinition(aNeutron);
3566  particleIsDefined = true;
3567  } else if (v[0]->GetMass() / MeV - aDeuteronMass < 0.1) {
3568  v[0]->SetDefinition(aDeuteron);
3569  particleIsDefined = true;
3570  } else if (v[0]->GetMass() / MeV - aTritonMass < 0.1) {
3571  v[0]->SetDefinition(aTriton);
3572  particleIsDefined = true;
3573  } else if (v[0]->GetMass() / MeV - anAlphaMass < 0.1) {
3574  v[0]->SetDefinition(anAlpha);
3575  particleIsDefined = true;
3576  }
3577  currentParticle.SetKineticEnergy(std::max(0.001, currentParticle.GetKineticEnergy() / MeV));
3578  p = currentParticle.GetTotalMomentum();
3579  pp = currentParticle.GetMomentum().mag();
3580  if (pp <= 0.001 * MeV) {
3581  G4double phinve = twopi * G4UniformRand();
3582  G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3583  currentParticle.SetMomentum(
3584  p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3585  } else
3586  currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3587 
3588  if (particleIsDefined) {
3589  v[0]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[0]->GetKineticEnergy() / MeV));
3590  p = v[0]->GetTotalMomentum();
3591  pp = v[0]->GetMomentum().mag();
3592  if (pp <= 0.001 * MeV) {
3593  G4double phinve = twopi * G4UniformRand();
3594  G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3595  v[0]->SetMomentum(
3596  p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3597  } else
3598  v[0]->SetMomentum(v[0]->GetMomentum() * (p / pp));
3599  }
3600  if ((v[1]->GetDefinition() == aDeuteron) || (v[1]->GetDefinition() == aTriton) || (v[1]->GetDefinition() == anAlpha))
3601  v[1]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[1]->GetKineticEnergy() / MeV));
3602  else
3603  v[1]->SetKineticEnergy(std::max(0.001, v[1]->GetKineticEnergy() / MeV));
3604 
3605  p = v[1]->GetTotalMomentum();
3606  pp = v[1]->GetMomentum().mag();
3607  if (pp <= 0.001 * MeV) {
3608  G4double phinve = twopi * G4UniformRand();
3609  G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3610  v[1]->SetMomentum(
3611  p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3612  } else
3613  v[1]->SetMomentum(v[1]->GetMomentum() * (p / pp));
3614 
3615  if (nt == 3) {
3616  if ((v[2]->GetDefinition() == aDeuteron) || (v[2]->GetDefinition() == aTriton) ||
3617  (v[2]->GetDefinition() == anAlpha))
3618  v[2]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[2]->GetKineticEnergy() / MeV));
3619  else
3620  v[2]->SetKineticEnergy(std::max(0.001, v[2]->GetKineticEnergy() / MeV));
3621 
3622  p = v[2]->GetTotalMomentum();
3623  pp = v[2]->GetMomentum().mag();
3624  if (pp <= 0.001 * MeV) {
3625  G4double phinve = twopi * G4UniformRand();
3626  G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3627  v[2]->SetMomentum(
3628  p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3629  } else
3630  v[2]->SetMomentum(v[2]->GetMomentum() * (p / pp));
3631  }
3632  G4int del;
3633  for (del = 0; del < vecLen; del++)
3634  delete vec[del];
3635  vecLen = 0;
3636  if (particleIsDefined) {
3637  vec.SetElement(vecLen++, v[0]);
3638  } else {
3639  delete v[0];
3640  }
3641  vec.SetElement(vecLen++, v[1]);
3642  if (nt == 3) {
3643  vec.SetElement(vecLen++, v[2]);
3644  } else {
3645  delete v[2];
3646  }
3647  delete[] v;
3648  return;
3649 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MeV
const double MeV
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
FullModelReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2340
nt
int nt
Definition: AMPTWrapper.h:32
ke
int ke
Definition: CascadeWrapper.h:14
Gamma
dbl * Gamma
Definition: mlp_gen.cc:38
patCaloMETCorrections_cff.A
A
Definition: patCaloMETCorrections_cff.py:43
G4int FullModelReactionDynamics::Poisson ( G4double  x)
private

Definition at line 2532 of file FullModelReactionDynamics.cc.

References JetChargeProducer_cfi::exp, mps_fire::i, cmsBatch::log, SiStripPI::max, p1, p2, p3, funct::pow(), alignCSCRings::r, summarizeEdmComparisonLogfiles::ran, findQualityFiles::rr, and mathSSE::sqrt().

Referenced by FindTimeDelay().

2533 {
2534  G4int iran;
2535  G4double ran;
2536 
2537  if (x > 9.9) // use normal distribution with sigma^2 = <x>
2538  iran = static_cast<G4int>(std::max(0.0, x + normal() * std::sqrt(x)));
2539  else {
2540  G4int mm = G4int(5.0 * x);
2541  if (mm <= 0) // for very small x try iran=1,2,3
2542  {
2543  G4double p1 = x * std::exp(-x);
2544  G4double p2 = x * p1 / 2.0;
2545  G4double p3 = x * p2 / 3.0;
2546  ran = G4UniformRand();
2547  if (ran < p3)
2548  iran = 3;
2549  else if (ran < p2) // this is original Geisha, it should be ran < p2+p3
2550  iran = 2;
2551  else if (ran < p1) // should be ran < p1+p2+p3
2552  iran = 1;
2553  else
2554  iran = 0;
2555  } else {
2556  iran = 0;
2557  G4double r = std::exp(-x);
2558  ran = G4UniformRand();
2559  if (ran > r) {
2560  G4double rrr;
2561  G4double rr = r;
2562  for (G4int i = 1; i <= mm; ++i) {
2563  iran++;
2564  if (i > 5) // Stirling's formula for large numbers
2565  rrr = std::exp(i * std::log(x) - (i + 0.5) * std::log((G4double)i) + i - 0.9189385);
2566  else
2567  rrr = std::pow(x, i) / Factorial(i);
2568  rr += r * rrr;
2569  if (ran <= rr)
2570  break;
2571  }
2572  }
2573  }
2574  }
2575  return iran;
2576 }
mps_fire.i
i
Definition: mps_fire.py:338
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
p2
double p2[4]
Definition: TauolaWrapper.h:90
p1
double p1[4]
Definition: TauolaWrapper.h:89
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
p3
double p3[4]
Definition: TauolaWrapper.h:91
void FullModelReactionDynamics::ProduceStrangeParticlePairs ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
const G4ReactionProduct &  modifiedOriginal,
const G4DynamicParticle *  originalTarget,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged 
)

Definition at line 3027 of file FullModelReactionDynamics.cc.

References JetChargeProducer_cfi::exp, GeV, mps_fire::i, BPhysicsValidation_cfi::Lambda, cmsBatch::log, SiStripPI::max, p1, summarizeEdmComparisonLogfiles::ran, phase2TrackerDigitizer_cfi::SigmaZero, and mathSSE::sqrt().

Referenced by FindTimeDelay().

3034  {
3035  // derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
3036  //
3037  // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
3038  // K+ Y0, K0 Y+, K0 Y-
3039  // For antibaryon induced reactions half of the cross sections KB YB
3040  // pairs are produced. Charge is not conserved, no experimental data available
3041  // for exclusive reactions, therefore some average behaviour assumed.
3042  // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
3043  //
3044  if (vecLen == 0)
3045  return;
3046  //
3047  // the following protects against annihilation processes
3048  //
3049  if (currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0)
3050  return;
3051 
3052  const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
3053  const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass() / GeV;
3054  G4double targetMass = originalTarget->GetDefinition()->GetPDGMass() / GeV;
3055  G4double centerofmassEnergy =
3056  std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal); // GeV
3057  G4double currentMass = currentParticle.GetMass() / GeV;
3058  G4double availableEnergy = centerofmassEnergy - (targetMass + currentMass);
3059  if (availableEnergy <= 1.0)
3060  return;
3061 
3062  G4ParticleDefinition *aProton = G4Proton::Proton();
3063  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
3064  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3065  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
3066  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
3067  G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
3068  G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
3069  G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
3070  G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
3071  G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
3072  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
3073  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
3074  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
3075  G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
3076  G4ParticleDefinition *aLambda = G4Lambda::Lambda();
3077  G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
3078 
3079  const G4double protonMass = aProton->GetPDGMass() / GeV;
3080  const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass() / GeV;
3081  //
3082  // determine the center of mass energy bin
3083  //
3084  const G4double avrs[] = {3., 4., 5., 6., 7., 8., 9., 10., 20., 30., 40., 50.};
3085 
3086  G4int ibin, i3, i4;
3087  G4double avk, avy, avn, ran;
3088  G4int i = 1;
3089  while ((i < 12) && (centerofmassEnergy > avrs[i]))
3090  ++i;
3091  if (i == 12)
3092  ibin = 11;
3093  else
3094  ibin = i;
3095  //
3096  // the fortran code chooses a random replacement of produced kaons
3097  // but does not take into account charge conservation
3098  //
3099  if (vecLen == 1) // we know that vecLen > 0
3100  {
3101  i3 = 0;
3102  i4 = 1; // note that we will be adding a new secondary particle in this case only
3103  } else // otherwise 0 <= i3,i4 < vecLen
3104  {
3105  G4double ran = G4UniformRand();
3106  while (ran == 1.0)
3107  ran = G4UniformRand();
3108  i4 = i3 = G4int(vecLen * ran);
3109  while (i3 == i4) {
3110  ran = G4UniformRand();
3111  while (ran == 1.0)
3112  ran = G4UniformRand();
3113  i4 = G4int(vecLen * ran);
3114  }
3115  }
3116  //
3117  // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
3118  //
3119  const G4double avkkb[] = {0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075, 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573};
3120  const G4double avky[] = {0.005, 0.03, 0.064, 0.095, 0.115, 0.13, 0.145, 0.155, 0.20, 0.205, 0.210, 0.212};
3121  const G4double avnnb[] = {0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02, 0.04, 0.05, 0.12, 0.15, 0.18, 0.20};
3122 
3123  avk = (std::log(avkkb[ibin]) - std::log(avkkb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3124  (avrs[ibin] - avrs[ibin - 1]) +
3125  std::log(avkkb[ibin - 1]);
3126  avk = std::exp(avk);
3127 
3128  avy = (std::log(avky[ibin]) - std::log(avky[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3129  (avrs[ibin] - avrs[ibin - 1]) +
3130  std::log(avky[ibin - 1]);
3131  avy = std::exp(avy);
3132 
3133  avn = (std::log(avnnb[ibin]) - std::log(avnnb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3134  (avrs[ibin] - avrs[ibin - 1]) +
3135  std::log(avnnb[ibin - 1]);
3136  avn = std::exp(avn);
3137 
3138  if (avk + avy + avn <= 0.0)
3139  return;
3140 
3141  if (currentMass < protonMass)
3142  avy /= 2.0;
3143  if (targetMass < protonMass)
3144  avy = 0.0;
3145  avy += avk + avn;
3146  avk += avn;
3147  ran = G4UniformRand();
3148  if (ran < avn) {
3149  if (availableEnergy < 2.0)
3150  return;
3151  if (vecLen == 1) // add a new secondary
3152  {
3153  G4ReactionProduct *p1 = new G4ReactionProduct;
3154  if (G4UniformRand() < 0.5) {
3155  vec[0]->SetDefinition(aNeutron);
3156  p1->SetDefinition(anAntiNeutron);
3157  (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3158  vec[0]->SetMayBeKilled(false);
3159  p1->SetMayBeKilled(false);
3160  } else {
3161  vec[0]->SetDefinition(aProton);
3162  p1->SetDefinition(anAntiProton);
3163  (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3164  vec[0]->SetMayBeKilled(false);
3165  p1->SetMayBeKilled(false);
3166  }
3167  vec.SetElement(vecLen++, p1);
3168  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3169  } else { // replace two secondaries
3170  if (G4UniformRand() < 0.5) {
3171  vec[i3]->SetDefinition(aNeutron);
3172  vec[i4]->SetDefinition(anAntiNeutron);
3173  vec[i3]->SetMayBeKilled(false);
3174  vec[i4]->SetMayBeKilled(false);
3175  } else {
3176  vec[i3]->SetDefinition(aProton);
3177  vec[i4]->SetDefinition(anAntiProton);
3178  vec[i3]->SetMayBeKilled(false);
3179  vec[i4]->SetMayBeKilled(false);
3180  }
3181  }
3182  } else if (ran < avk) {
3183  if (availableEnergy < 1.0)
3184  return;
3185 
3186  const G4double kkb[] = {0.2500, 0.3750, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8750, 1.000};
3187  const G4int ipakkb1[] = {10, 10, 10, 11, 11, 12, 12, 11, 12};
3188  const G4int ipakkb2[] = {13, 11, 12, 11, 12, 11, 12, 13, 13};
3189  ran = G4UniformRand();
3190  i = 0;
3191  while ((i < 9) && (ran >= kkb[i]))
3192  ++i;
3193  if (i == 9)
3194  return;
3195  //
3196  // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
3197  // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
3198  //
3199  switch (ipakkb1[i]) {
3200  case 10:
3201  vec[i3]->SetDefinition(aKaonPlus);
3202  vec[i3]->SetMayBeKilled(false);
3203  break;
3204  case 11:
3205  vec[i3]->SetDefinition(aKaonZS);
3206  vec[i3]->SetMayBeKilled(false);
3207  break;
3208  case 12:
3209  vec[i3]->SetDefinition(aKaonZL);
3210  vec[i3]->SetMayBeKilled(false);
3211  break;
3212  }
3213  if (vecLen == 1) // add a secondary
3214  {
3215  G4ReactionProduct *p1 = new G4ReactionProduct;
3216  switch (ipakkb2[i]) {
3217  case 11:
3218  p1->SetDefinition(aKaonZS);
3219  p1->SetMayBeKilled(false);
3220  break;
3221  case 12:
3222  p1->SetDefinition(aKaonZL);
3223  p1->SetMayBeKilled(false);
3224  break;
3225  case 13:
3226  p1->SetDefinition(aKaonMinus);
3227  p1->SetMayBeKilled(false);
3228  break;
3229  }
3230  (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3231  vec.SetElement(vecLen++, p1);
3232  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3233  } else // replace
3234  {
3235  switch (ipakkb2[i]) {
3236  case 11:
3237  vec[i4]->SetDefinition(aKaonZS);
3238  vec[i4]->SetMayBeKilled(false);
3239  break;
3240  case 12:
3241  vec[i4]->SetDefinition(aKaonZL);
3242  vec[i4]->SetMayBeKilled(false);
3243  break;
3244  case 13:
3245  vec[i4]->SetDefinition(aKaonMinus);
3246  vec[i4]->SetMayBeKilled(false);
3247  break;
3248  }
3249  }
3250  } else if (ran < avy) {
3251  if (availableEnergy < 1.6)
3252  return;
3253 
3254  const G4double ky[] = {0.200, 0.300, 0.400, 0.550, 0.625, 0.700, 0.800, 0.850, 0.900, 0.950, 0.975, 1.000};
3255  const G4int ipaky1[] = {18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22};
3256  const G4int ipaky2[] = {10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12};
3257  const G4int ipakyb1[] = {19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25};
3258  const G4int ipakyb2[] = {13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11};
3259  ran = G4UniformRand();
3260  i = 0;
3261  while ((i < 12) && (ran > ky[i]))
3262  ++i;
3263  if (i == 12)
3264  return;
3265  if ((currentMass < protonMass) || (G4UniformRand() < 0.5)) {
3266  // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
3267  // 0 + 0 0 0 0 + + + 0 + 0
3268  //
3269  // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
3270  // 0 + 0 0 0 0 - + - 0 - 0
3271  switch (ipaky1[i]) {
3272  case 18:
3273  targetParticle.SetDefinition(aLambda);
3274  break;
3275  case 20:
3276  targetParticle.SetDefinition(aSigmaPlus);
3277  break;
3278  case 21:
3279  targetParticle.SetDefinition(aSigmaZero);
3280  break;
3281  case 22:
3282  targetParticle.SetDefinition(aSigmaMinus);
3283  break;
3284  }
3285  targetHasChanged = true;
3286  switch (ipaky2[i]) {
3287  case 10:
3288  vec[i3]->SetDefinition(aKaonPlus);
3289  vec[i3]->SetMayBeKilled(false);
3290  break;
3291  case 11:
3292  vec[i3]->SetDefinition(aKaonZS);
3293  vec[i3]->SetMayBeKilled(false);
3294  break;
3295  case 12:
3296  vec[i3]->SetDefinition(aKaonZL);
3297  vec[i3]->SetMayBeKilled(false);
3298  break;
3299  }
3300  } else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
3301  {
3302  if ((currentParticle.GetDefinition() == anAntiProton) || (currentParticle.GetDefinition() == anAntiNeutron) ||
3303  (currentParticle.GetDefinition() == anAntiLambda) || (currentMass > sigmaMinusMass)) {
3304  switch (ipakyb1[i]) {
3305  case 19:
3306  currentParticle.SetDefinitionAndUpdateE(anAntiLambda);
3307  break;
3308  case 23:
3309  currentParticle.SetDefinitionAndUpdateE(anAntiSigmaPlus);
3310  break;
3311  case 24:
3312  currentParticle.SetDefinitionAndUpdateE(anAntiSigmaZero);
3313  break;
3314  case 25:
3315  currentParticle.SetDefinitionAndUpdateE(anAntiSigmaMinus);
3316  break;
3317  }
3318  incidentHasChanged = true;
3319  switch (ipakyb2[i]) {
3320  case 11:
3321  vec[i3]->SetDefinition(aKaonZS);
3322  vec[i3]->SetMayBeKilled(false);
3323  break;
3324  case 12:
3325  vec[i3]->SetDefinition(aKaonZL);
3326  vec[i3]->SetMayBeKilled(false);
3327  break;
3328  case 13:
3329  vec[i3]->SetDefinition(aKaonMinus);
3330  vec[i3]->SetMayBeKilled(false);
3331  break;
3332  }
3333  } else {
3334  switch (ipaky1[i]) {
3335  case 18:
3336  currentParticle.SetDefinitionAndUpdateE(aLambda);
3337  break;
3338  case 20:
3339  currentParticle.SetDefinitionAndUpdateE(aSigmaPlus);
3340  break;
3341  case 21:
3342  currentParticle.SetDefinitionAndUpdateE(aSigmaZero);
3343  break;
3344  case 22:
3345  currentParticle.SetDefinitionAndUpdateE(aSigmaMinus);
3346  break;
3347  }
3348  incidentHasChanged = true;
3349  switch (ipaky2[i]) {
3350  case 10:
3351  vec[i3]->SetDefinition(aKaonPlus);
3352  vec[i3]->SetMayBeKilled(false);
3353  break;
3354  case 11:
3355  vec[i3]->SetDefinition(aKaonZS);
3356  vec[i3]->SetMayBeKilled(false);
3357  break;
3358  case 12:
3359  vec[i3]->SetDefinition(aKaonZL);
3360  vec[i3]->SetMayBeKilled(false);
3361  break;
3362  }
3363  }
3364  }
3365  } else
3366  return;
3367  //
3368  // check the available energy
3369  // if there is not enough energy for kkb/ky pair production
3370  // then reduce the number of secondary particles
3371  // NOTE:
3372  // the number of secondaries may have been changed
3373  // the incident and/or target particles may have changed
3374  // charge conservation is ignored (as well as strangness conservation)
3375  //
3376  currentMass = currentParticle.GetMass() / GeV;
3377  targetMass = targetParticle.GetMass() / GeV;
3378 
3379  G4double energyCheck = centerofmassEnergy - (currentMass + targetMass);
3380  for (i = 0; i < vecLen; ++i) {
3381  energyCheck -= vec[i]->GetMass() / GeV;
3382  if (energyCheck < 0.0) // chop off the secondary List
3383  {
3384  vecLen = std::max(0, --i); // looks like a memory leak @@@@@@@@@@@@
3385  G4int j;
3386  for (j = i; j < vecLen; j++)
3387  delete vec[j];
3388  break;
3389  }
3390  }
3391  return;
3392 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
p1
double p1[4]
Definition: TauolaWrapper.h:89
BPhysicsValidation_cfi.Lambda
Lambda
Definition: BPhysicsValidation_cfi.py:63
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
phase2TrackerDigitizer_cfi.SigmaZero
SigmaZero
Definition: phase2TrackerDigitizer_cfi.py:46
void FullModelReactionDynamics::Rotate ( const G4double  numberofFinalStateNucleons,
const G4ThreeVector &  temp,
const G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
const G4Nucleus &  targetNucleus,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2639 of file FullModelReactionDynamics.cc.

References funct::cos(), MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, fermi, GeV, mps_fire::i, cuy::ii, cmsBatch::log, SiStripPI::max, MeV, AlCaHLTBitMon_ParallelJobs::p, pi, createTree::pp, diffTwoXMLs::r1, diffTwoXMLs::r2, idealTransformation::rotation, funct::sin(), mathSSE::sqrt(), and groupFilesInBlocks::temp.

Referenced by FindTimeDelay().

2648  {
2649  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2650  //
2651  // Rotate in direction of z-axis, this does disturb in some way our
2652  // inclusive distributions, but it is necessary for momentum conservation
2653  //
2654  const G4double atomicWeight = targetNucleus.GetN_asInt();
2655  const G4double logWeight = std::log(atomicWeight);
2656 
2657  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2658  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2659  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2660 
2661  G4int i;
2662  G4ThreeVector pseudoParticle[4];
2663  for (i = 0; i < 4; ++i)
2664  pseudoParticle[i].set(0, 0, 0);
2665  pseudoParticle[0] = currentParticle.GetMomentum() + targetParticle.GetMomentum();
2666  for (i = 0; i < vecLen; ++i)
2667  pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
2668  //
2669  // Some smearing in transverse direction from Fermi motion
2670  //
2671  G4float pp, pp1;
2672  G4double alekw, p, rthnve, phinve;
2673  G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
2674 
2675  r1 = twopi * G4UniformRand();
2676  r2 = G4UniformRand();
2677  a1 = std::sqrt(-2.0 * std::log(r2));
2678  ran1 = a1 * std::sin(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2679  ran2 = a1 * std::cos(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2680  G4ThreeVector fermi(ran1, ran2, 0);
2681 
2682  pseudoParticle[0] = pseudoParticle[0] + fermi; // all particles + fermi
2683  pseudoParticle[2] = temp; // original in cms system
2684  pseudoParticle[3] = pseudoParticle[0];
2685 
2686  pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
2687  G4double rotation = 2. * pi * G4UniformRand();
2688  pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
2689  pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
2690  for (G4int ii = 1; ii <= 3; ii++) {
2691  p = pseudoParticle[ii].mag();
2692  if (p == 0.0)
2693  pseudoParticle[ii] = G4ThreeVector(0.0, 0.0, 0.0);
2694  else
2695  pseudoParticle[ii] = pseudoParticle[ii] * (1. / p);
2696  }
2697 
2698  pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
2699  pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
2700  pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
2701  currentParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2702 
2703  pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
2704  pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
2705  pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
2706  targetParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2707 
2708  for (i = 0; i < vecLen; ++i) {
2709  pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
2710  pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
2711  pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
2712  vec[i]->SetMomentum(pxTemp, pyTemp, pzTemp);
2713  }
2714  //
2715  // Rotate in direction of primary particle, subtract binding energies
2716  // and make some further corrections if required
2717  //
2718  Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2719  G4double ekin;
2720  G4double dekin = 0.0;
2721  G4double ek1 = 0.0;
2722  G4int npions = 0;
2723  if (atomicWeight >= 1.5) // self-absorption in heavy molecules
2724  {
2725  // corrections for single particle spectra (shower particles)
2726  //
2727  const G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00};
2728  const G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65};
2729  alekw = std::log(originalIncident->GetKineticEnergy() / GeV);
2730  exh = 1.0;
2731  if (alekw > alem[0]) // get energy bin
2732  {
2733  exh = val0[6];
2734  for (G4int j = 1; j < 7; ++j) {
2735  if (alekw < alem[j]) // use linear interpolation/extrapolation
2736  {
2737  G4double rcnve = (val0[j] - val0[j - 1]) / (alem[j] - alem[j - 1]);
2738  exh = rcnve * alekw + val0[j - 1] - rcnve * alem[j - 1];
2739  break;
2740  }
2741  }
2742  exh = 1.0 - exh;
2743  }
2744  const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2745  ekin = currentParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2746  ekin = std::max(1.0e-6, ekin);
2747  xxh = 1.0;
2748  dekin += ekin * (1.0 - xxh);
2749  ekin *= xxh;
2750  currentParticle.SetKineticEnergy(ekin * GeV);
2751  pp = currentParticle.GetTotalMomentum() / MeV;
2752  pp1 = currentParticle.GetMomentum().mag() / MeV;
2753  if (pp1 < 0.001 * MeV) {
2754  rthnve = pi * G4UniformRand();
2755  phinve = twopi * G4UniformRand();
2756  currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2757  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2758  pp * std::cos(rthnve) * MeV);
2759  } else
2760  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2761  ekin = targetParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2762  ekin = std::max(1.0e-6, ekin);
2763  xxh = 1.0;
2764  if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2765  (targetParticle.GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2766  xxh = exh;
2767  dekin += ekin * (1.0 - xxh);
2768  ekin *= xxh;
2769  if ((targetParticle.GetDefinition() == aPiPlus) || (targetParticle.GetDefinition() == aPiZero) ||
2770  (targetParticle.GetDefinition() == aPiMinus)) {
2771  ++npions;
2772  ek1 += ekin;
2773  }
2774  targetParticle.SetKineticEnergy(ekin * GeV);
2775  pp = targetParticle.GetTotalMomentum() / MeV;
2776  pp1 = targetParticle.GetMomentum().mag() / MeV;
2777  if (pp1 < 0.001 * MeV) {
2778  rthnve = pi * G4UniformRand();
2779  phinve = twopi * G4UniformRand();
2780  targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2781  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2782  pp * std::cos(rthnve) * MeV);
2783  } else
2784  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2785  for (i = 0; i < vecLen; ++i) {
2786  ekin = vec[i]->GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2787  ekin = std::max(1.0e-6, ekin);
2788  xxh = 1.0;
2789  if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2790  (vec[i]->GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2791  xxh = exh;
2792  dekin += ekin * (1.0 - xxh);
2793  ekin *= xxh;
2794  if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2795  (vec[i]->GetDefinition() == aPiMinus)) {
2796  ++npions;
2797  ek1 += ekin;
2798  }
2799  vec[i]->SetKineticEnergy(ekin * GeV);
2800  pp = vec[i]->GetTotalMomentum() / MeV;
2801  pp1 = vec[i]->GetMomentum().mag() / MeV;
2802  if (pp1 < 0.001 * MeV) {
2803  rthnve = pi * G4UniformRand();
2804  phinve = twopi * G4UniformRand();
2805  vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2806  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2807  pp * std::cos(rthnve) * MeV);
2808  } else
2809  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2810  }
2811  }
2812  if ((ek1 != 0.0) && (npions > 0)) {
2813  dekin = 1.0 + dekin / ek1;
2814  //
2815  // first do the incident particle
2816  //
2817  if ((currentParticle.GetDefinition() == aPiPlus) || (currentParticle.GetDefinition() == aPiZero) ||
2818  (currentParticle.GetDefinition() == aPiMinus)) {
2819  currentParticle.SetKineticEnergy(std::max(0.001 * MeV, dekin * currentParticle.GetKineticEnergy()));
2820  pp = currentParticle.GetTotalMomentum() / MeV;
2821  pp1 = currentParticle.GetMomentum().mag() / MeV;
2822  if (pp1 < 0.001) {
2823  rthnve = pi * G4UniformRand();
2824  phinve = twopi * G4UniformRand();
2825  currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2826  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2827  pp * std::cos(rthnve) * MeV);
2828  } else
2829  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2830  }
2831  for (i = 0; i < vecLen; ++i) {
2832  if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2833  (vec[i]->GetDefinition() == aPiMinus)) {
2834  vec[i]->SetKineticEnergy(std::max(0.001 * MeV, dekin * vec[i]->GetKineticEnergy()));
2835  pp = vec[i]->GetTotalMomentum() / MeV;
2836  pp1 = vec[i]->GetMomentum().mag() / MeV;
2837  if (pp1 < 0.001) {
2838  rthnve = pi * G4UniformRand();
2839  phinve = twopi * G4UniformRand();
2840  vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2841  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2842  pp * std::cos(rthnve) * MeV);
2843  } else
2844  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2845  }
2846  }
2847  }
2848 }
GeV
const double GeV
Definition: MathUtil.h:16
FullModelReactionDynamics::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2588
mps_fire.i
i
Definition: mps_fire.py:338
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
pi
const Double_t pi
Definition: trackSplitPlot.h:36
fermi
const double fermi
Definition: MathUtil.h:17
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
cuy.ii
ii
Definition: cuy.py:590
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
void FullModelReactionDynamics::SuppressChargedPions ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged 
)

Definition at line 1235 of file FullModelReactionDynamics.cc.

References GeV, and mps_fire::i.

Referenced by FullModelHadronicProcess::CalculateMomenta(), and FindTimeDelay().

1242  {
1243  // this code was originally in the fortran code TWOCLU
1244  //
1245  // suppress charged pions, for various reasons
1246  //
1247  const G4double atomicWeight = targetNucleus.GetN_asInt();
1248  const G4double atomicNumber = targetNucleus.GetZ_asInt();
1249  const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / GeV;
1250 
1251  // G4ParticleDefinition *aGamma = G4Gamma::Gamma();
1252  G4ParticleDefinition *aProton = G4Proton::Proton();
1253  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
1254  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1255  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
1256  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1257  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1258 
1259  const G4bool antiTest =
1260  modifiedOriginal.GetDefinition() != anAntiProton && modifiedOriginal.GetDefinition() != anAntiNeutron;
1261  if (antiTest &&
1262  (
1263  // currentParticle.GetDefinition() == aGamma ||
1264  currentParticle.GetDefinition() == aPiPlus || currentParticle.GetDefinition() == aPiMinus) &&
1265  (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) {
1266  if (G4UniformRand() > atomicNumber / atomicWeight)
1267  currentParticle.SetDefinitionAndUpdateE(aNeutron);
1268  else
1269  currentParticle.SetDefinitionAndUpdateE(aProton);
1270  incidentHasChanged = true;
1271  }
1272  for (G4int i = 0; i < vecLen; ++i) {
1273  if (antiTest &&
1274  (
1275  // vec[i]->GetDefinition() == aGamma ||
1276  vec[i]->GetDefinition() == aPiPlus || vec[i]->GetDefinition() == aPiMinus) &&
1277  (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) {
1278  if (G4UniformRand() > atomicNumber / atomicWeight)
1279  vec[i]->SetDefinitionAndUpdateE(aNeutron);
1280  else
1281  vec[i]->SetDefinitionAndUpdateE(aProton);
1282  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1283  }
1284  }
1285  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1286 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
void FullModelReactionDynamics::TwoBody ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4DynamicParticle *  originalTarget,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  targetHasChanged 
)

Definition at line 2069 of file FullModelReactionDynamics.cc.

References b, funct::cos(), JetChargeProducer_cfi::exp, GeV, mps_fire::i, cmsBatch::log, SiStripPI::max, MeV, min(), packedPFCandidateRefMixer_cfi::pf, createTree::pp, funct::sin(), and mathSSE::sqrt().

Referenced by FullModelHadronicProcess::CalculateMomenta(), and FindTimeDelay().

2076  {
2077  // G4cout<<"TwoBody called"<<G4endl;
2078  //
2079  // derived from original FORTRAN code TWOB by H. Fesefeldt (15-Sep-1987)
2080  //
2081  // Generation of momenta for elastic and quasi-elastic 2 body reactions
2082  //
2083  // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
2084  // The b values are parametrizations from experimental data.
2085  // Not available values are taken from those of similar reactions.
2086  //
2087  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2088  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2089  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2090  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
2091  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
2092  G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
2093  G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
2094 
2095  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2096  static const G4double expxu = 82.; // upper bound for arg. of exp
2097  static const G4double expxl = -expxu; // lower bound for arg. of exp
2098 
2099  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
2100  const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
2101  G4double currentMass = currentParticle.GetMass() / GeV;
2102  G4double targetMass = targetParticle.GetMass() / GeV;
2103  const G4double atomicWeight = targetNucleus.GetN_asInt();
2104  // G4cout<<"Atomic weight is found to be: "<<atomicWeight<<G4endl;
2105  G4double etCurrent = currentParticle.GetTotalEnergy() / GeV;
2106  G4double pCurrent = currentParticle.GetTotalMomentum() / GeV;
2107 
2108  G4double cmEnergy =
2109  std::sqrt(currentMass * currentMass + targetMass * targetMass + 2.0 * targetMass * etCurrent); // in GeV
2110 
2111  //if( (pOriginal < 0.1) ||
2112  // (centerofmassEnergy < 0.01) ) // 2-body scattering not possible
2113  // Continue with original particle, but spend the nuclear evaporation energy
2114  // targetParticle.SetMass( 0.0 ); // flag that the target doesn't exist
2115  //else // Two-body scattering is possible
2116 
2117  if ((pCurrent < 0.1) || (cmEnergy < 0.01)) // 2-body scattering not possible
2118  {
2119  targetParticle.SetMass(0.0); // flag that the target particle doesn't exist
2120  } else {
2121  // moved this if-block to a later stage, i.e. to the assignment of the scattering angle
2122  // @@@@@ double-check.
2123  // if( targetParticle.GetDefinition() == aKaonMinus ||
2124  // targetParticle.GetDefinition() == aKaonZeroL ||
2125  // targetParticle.GetDefinition() == aKaonZeroS ||
2126  // targetParticle.GetDefinition() == aKaonPlus ||
2127  // targetParticle.GetDefinition() == aPiMinus ||
2128  // targetParticle.GetDefinition() == aPiZero ||
2129  // targetParticle.GetDefinition() == aPiPlus )
2130  // {
2131  // if( G4UniformRand() < 0.5 )
2132  // targetParticle.SetDefinitionAndUpdateE( aNeutron );
2133  // else
2134  // targetParticle.SetDefinitionAndUpdateE( aProton );
2135  // targetHasChanged = true;
2136  // targetMass = targetParticle.GetMass()/GeV;
2137  // }
2138  //
2139  // Set masses and momenta for final state particles
2140  //
2141  /*
2142  G4cout<<"Check 0"<<G4endl;
2143  G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()/GeV<<G4endl;
2144  G4cout<<"target mass: "<<targetParticle.GetMass()/GeV<<G4endl;
2145  G4cout<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
2146  G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()/GeV<<G4endl;
2147  G4cout<<"current mass: "<<currentParticle.GetMass()/GeV<<G4endl;
2148  G4cout<<currentParticle.GetDefinition()->GetParticleName()<<G4endl;
2149  */
2150  G4double pf = cmEnergy * cmEnergy + targetMass * targetMass - currentMass * currentMass;
2151  pf = pf * pf - 4 * cmEnergy * cmEnergy * targetMass * targetMass;
2152  // G4cout << "pf: " << pf<< G4endl;
2153 
2154  if (pf <= 0.) // 0.001 )
2155  {
2156  for (G4int i = 0; i < vecLen; i++)
2157  delete vec[i];
2158  vecLen = 0;
2159  throw G4HadronicException(__FILE__, __LINE__, "FullModelReactionDynamics::TwoBody: pf is too small ");
2160  }
2161 
2162  pf = std::sqrt(pf) / (2.0 * cmEnergy);
2163  //
2164  // Set beam and target in centre of mass system
2165  //
2166  G4ReactionProduct pseudoParticle[3];
2167  pseudoParticle[0].SetMass(currentMass * GeV);
2168  pseudoParticle[0].SetTotalEnergy(etCurrent * GeV);
2169  pseudoParticle[0].SetMomentum(0.0, 0.0, pCurrent * GeV);
2170 
2171  pseudoParticle[1].SetMomentum(0.0, 0.0, 0.0);
2172  pseudoParticle[1].SetMass(targetMass * GeV);
2173  pseudoParticle[1].SetKineticEnergy(0.0);
2174  //
2175  // Transform into centre of mass system
2176  //
2177  pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
2178  pseudoParticle[0].Lorentz(pseudoParticle[0], pseudoParticle[2]);
2179  pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[2]);
2180  //
2181  // Set final state masses and energies in centre of mass system
2182  //
2183  currentParticle.SetTotalEnergy(std::sqrt(pf * pf + currentMass * currentMass) * GeV);
2184  targetParticle.SetTotalEnergy(std::sqrt(pf * pf + targetMass * targetMass) * GeV);
2185  //
2186  // Set |t| and |tmin|
2187  //
2188  const G4double cb = 0.01;
2189  const G4double b1 = 4.225;
2190  const G4double b2 = 1.795;
2191  //
2192  // Calculate slope b for elastic scattering on proton/neutron
2193  //
2194  G4double b = std::max(cb, b1 + b2 * std::log(pOriginal));
2195  // G4double b = std::max( cb, b1+b2*std::log(ptemp) );
2196  G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag() / GeV;
2197 
2198  G4double exindt = -1.0;
2199  exindt += std::exp(std::max(-btrang, expxl));
2200  //
2201  // Calculate sqr(std::sin(teta/2.) and std::cos(teta), set azimuth angle phi
2202  //
2203  G4double ctet = 1.0 + 2 * std::log(1.0 + G4UniformRand() * exindt) / btrang;
2204  if (std::fabs(ctet) > 1.0)
2205  ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
2206  G4double stet = std::sqrt((1.0 - ctet) * (1.0 + ctet));
2207  G4double phi = twopi * G4UniformRand();
2208  //
2209  // Calculate final state momenta in centre of mass system
2210  //
2211  if (targetParticle.GetDefinition() == aKaonMinus || targetParticle.GetDefinition() == aKaonZeroL ||
2212  targetParticle.GetDefinition() == aKaonZeroS || targetParticle.GetDefinition() == aKaonPlus ||
2213  targetParticle.GetDefinition() == aPiMinus || targetParticle.GetDefinition() == aPiZero ||
2214  targetParticle.GetDefinition() == aPiPlus) {
2215  currentParticle.SetMomentum(-pf * stet * std::sin(phi) * GeV, -pf * stet * std::cos(phi) * GeV, -pf * ctet * GeV);
2216  } else {
2217  currentParticle.SetMomentum(pf * stet * std::sin(phi) * GeV, pf * stet * std::cos(phi) * GeV, pf * ctet * GeV);
2218  }
2219  targetParticle.SetMomentum(currentParticle.GetMomentum() * (-1.0));
2220  //
2221  // Transform into lab system
2222  //
2223  currentParticle.Lorentz(currentParticle, pseudoParticle[1]);
2224  targetParticle.Lorentz(targetParticle, pseudoParticle[1]);
2225 
2226  /*
2227  G4cout<<"Check 1"<<G4endl;
2228  G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
2229  G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
2230  G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
2231  G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
2232  */
2233 
2234  Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2235 
2236  G4double pp, pp1, ekin;
2237  if (atomicWeight >= 1.5) {
2238  const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2239  pp1 = currentParticle.GetMomentum().mag() / MeV;
2240  if (pp1 >= 1.0) {
2241  ekin = currentParticle.GetKineticEnergy() / MeV - cfa * (1.0 + 0.5 * normal()) * GeV;
2242  ekin = std::max(0.0001 * GeV, ekin);
2243  currentParticle.SetKineticEnergy(ekin * MeV);
2244  pp = currentParticle.GetTotalMomentum() / MeV;
2245  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2246  }
2247  pp1 = targetParticle.GetMomentum().mag() / MeV;
2248  if (pp1 >= 1.0) {
2249  ekin = targetParticle.GetKineticEnergy() / MeV - cfa * (1.0 + normal() / 2.) * GeV;
2250  ekin = std::max(0.0001 * GeV, ekin);
2251  targetParticle.SetKineticEnergy(ekin * MeV);
2252  pp = targetParticle.GetTotalMomentum() / MeV;
2253  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2254  }
2255  }
2256  }
2257  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2258  if (atomicWeight >= 1.5) {
2259  // Add black track particles
2260  // the procedure is somewhat different than in TwoCluster and GenerateXandPt.
2261  // The reason is that we have to also simulate the nuclear reactions
2262  // at low energies like a(h,p)b, a(h,p p)b, a(h,n)b etc.
2263  //
2264  // npnb is number of proton/neutron black track particles
2265  // ndta is the number of deuterons, tritons, and alphas produced
2266  // epnb is the kinetic energy available for proton/neutron black track particles
2267  // edta is the kinetic energy available for deuteron/triton/alpha particles
2268  //
2269  G4double epnb, edta;
2270  G4int npnb = 0, ndta = 0;
2271 
2272  epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
2273  edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
2274  const G4double pnCutOff = 0.0001; // GeV
2275  const G4double dtaCutOff = 0.0001; // GeV
2276  const G4double kineticMinimum = 0.0001;
2277  const G4double kineticFactor = -0.010;
2278  G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
2279  if (epnb >= pnCutOff) {
2280  npnb = Poisson(epnb / 0.02);
2281  /*
2282  G4cout<<"A couple of Poisson numbers:"<<G4endl;
2283  for (int n=0;n!=10;n++) G4cout<<Poisson(epnb/0.02)<<", ";
2284  G4cout<<G4endl;
2285  */
2286  if (npnb > atomicWeight)
2287  npnb = G4int(atomicWeight);
2288  if ((epnb > pnCutOff) && (npnb <= 0))
2289  npnb = 1;
2290  npnb = std::min(npnb, 127 - vecLen);
2291  }
2292  if (edta >= dtaCutOff) {
2293  ndta = G4int(2.0 * std::log(atomicWeight));
2294  ndta = std::min(ndta, 127 - vecLen);
2295  }
2296  G4double spall = 0.0;
2297  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2298 
2299  /*
2300  G4cout<<"Check 2"<<G4endl;
2301  G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
2302  G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
2303  G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
2304  G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
2305 
2306  G4cout<<"------------------------------------------------------------------------"<<G4endl;
2307  G4cout<<"Atomic weight: "<<atomicWeight<<G4endl;
2308  G4cout<<"number of proton/neutron black track particles: "<<npnb<<G4endl;
2309  G4cout<<"number of deuterons, tritons, and alphas produced: "<<ndta <<G4endl;
2310  G4cout<<"kinetic energy available for proton/neutron black track particles: "<<epnb/GeV<<" GeV" <<G4endl;
2311  G4cout<<"kinetic energy available for deuteron/triton/alpha particles: "<<edta/GeV <<" GeV"<<G4endl;
2312  G4cout<<"------------------------------------------------------------------------"<<G4endl;
2313  */
2314 
2315  AddBlackTrackParticles(epnb,
2316  npnb,
2317  edta,
2318  ndta,
2319  sprob,
2320  kineticMinimum,
2321  kineticFactor,
2322  modifiedOriginal,
2323  spall,
2324  targetNucleus,
2325  vec,
2326  vecLen);
2327 
2328  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2329  }
2330  //
2331  // calculate time delay for nuclear reactions
2332  //
2333  if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
2334  currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
2335  else
2336  currentParticle.SetTOF(1.0);
2337  return;
2338 }
GeV
const double GeV
Definition: MathUtil.h:16
FullModelReactionDynamics::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2588
mps_fire.i
i
Definition: mps_fire.py:338
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
FullModelReactionDynamics::AddBlackTrackParticles
void AddBlackTrackParticles(const G4double epnb, const G4int npnb, const G4double edta, const G4int ndta, const G4double sprob, const G4double kineticMinimum, const G4double kineticFactor, const G4ReactionProduct &modifiedOriginal, G4double spall, const G4Nucleus &aNucleus, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2850
b
double b
Definition: hdecay.h:120
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
G4bool FullModelReactionDynamics::TwoCluster ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged,
G4bool  leadFlag,
G4ReactionProduct &  leadingStrangeParticle 
)

Definition at line 1288 of file FullModelReactionDynamics.cc.

References a, funct::abs(), funct::cos(), diffTreeTool::diff, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, GeV, mps_fire::i, BPhysicsValidation_cfi::Lambda, cmsBatch::log, funct::m, SiStripPI::max, MeV, min(), packedPFCandidateRefMixer_cfi::pf, pi, funct::pow(), createTree::pp, summarizeEdmComparisonLogfiles::ran, phase2TrackerDigitizer_cfi::SigmaZero, funct::sin(), mathSSE::sqrt(), protons_cff::t, and groupFilesInBlocks::temp.

Referenced by FullModelHadronicProcess::CalculateMomenta(), and FindTimeDelay().

1299  {
1300  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1301  // derived from original FORTRAN code TWOCLU by H. Fesefeldt (11-Oct-1987)
1302  //
1303  // Generation of X- and PT- values for incident, target, and all secondary particles
1304  //
1305  // A simple two cluster model is used.
1306  // This should be sufficient for low energy interactions.
1307  //
1308 
1309  // + debugging
1310  // raise(SIGSEGV);
1311  // - debugging
1312 
1313  G4int i;
1314  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1315  G4ParticleDefinition *aProton = G4Proton::Proton();
1316  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1317  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1318  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
1319  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
1320  const G4double protonMass = aProton->GetPDGMass() / MeV;
1321  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
1322  const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
1323  const G4double mOriginal = modifiedOriginal.GetMass() / GeV;
1324  const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
1325  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass() / GeV;
1326  G4double centerofmassEnergy =
1327  std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal); // GeV
1328  G4double currentMass = currentParticle.GetMass() / GeV;
1329  targetMass = targetParticle.GetMass() / GeV;
1330 
1331  if (currentMass == 0.0 && targetMass == 0.0) {
1332  G4double ek = currentParticle.GetKineticEnergy();
1333  G4ThreeVector m = currentParticle.GetMomentum();
1334  currentParticle = *vec[0];
1335  targetParticle = *vec[1];
1336  for (i = 0; i < (vecLen - 2); ++i)
1337  *vec[i] = *vec[i + 2];
1338  if (vecLen < 2) {
1339  for (G4int i = 0; i < vecLen; i++)
1340  delete vec[i];
1341  vecLen = 0;
1342  throw G4HadReentrentException(
1343  __FILE__, __LINE__, "FullModelReactionDynamics::TwoCluster: Negative number of particles");
1344  }
1345  delete vec[vecLen - 1];
1346  delete vec[vecLen - 2];
1347  vecLen -= 2;
1348  incidentHasChanged = true;
1349  targetHasChanged = true;
1350  currentParticle.SetKineticEnergy(ek);
1351  currentParticle.SetMomentum(m);
1352  }
1353  const G4double atomicWeight = targetNucleus.GetN_asInt();
1354  const G4double atomicNumber = targetNucleus.GetZ_asInt();
1355  //
1356  // particles have been distributed in forward and backward hemispheres
1357  // in center of mass system of the hadron nucleon interaction
1358  //
1359  // incident is always in forward hemisphere
1360  G4int forwardCount = 1; // number of particles in forward hemisphere
1361  currentParticle.SetSide(1);
1362  G4double forwardMass = currentParticle.GetMass() / GeV;
1363  G4double cMass = forwardMass;
1364 
1365  // target is always in backward hemisphere
1366  G4int backwardCount = 1; // number of particles in backward hemisphere
1367  G4int backwardNucleonCount = 1; // number of nucleons in backward hemisphere
1368  targetParticle.SetSide(-1);
1369  G4double backwardMass = targetParticle.GetMass() / GeV;
1370  G4double bMass = backwardMass;
1371 
1372  for (i = 0; i < vecLen; ++i) {
1373  if (vec[i]->GetSide() < 0)
1374  vec[i]->SetSide(-1); // added by JLC, 2Jul97
1375  // to take care of the case where vec has been preprocessed by GenerateXandPt
1376  // and some of them have been set to -2 or -3
1377  if (vec[i]->GetSide() == -1) {
1378  ++backwardCount;
1379  backwardMass += vec[i]->GetMass() / GeV;
1380  } else {
1381  ++forwardCount;
1382  forwardMass += vec[i]->GetMass() / GeV;
1383  }
1384  }
1385  //
1386  // nucleons and some pions from intranuclear cascade
1387  //
1388  G4double term1 = std::log(centerofmassEnergy * centerofmassEnergy);
1389  if (term1 < 0)
1390  term1 = 0.0001; // making sure xtarg<0;
1391  const G4double afc = 0.312 + 0.2 * std::log(term1);
1392  G4double xtarg;
1393  if (centerofmassEnergy < 2.0 + G4UniformRand()) // added +2 below, JLC 4Jul97
1394  xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount + vecLen + 2) / 2.0;
1395  else
1396  xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount);
1397  if (xtarg <= 0.0)
1398  xtarg = 0.01;
1399  G4int nuclearExcitationCount = Poisson(xtarg);
1400  if (atomicWeight < 1.0001)
1401  nuclearExcitationCount = 0;
1402  G4int extraNucleonCount = 0;
1403  G4double extraMass = 0.0;
1404  G4double extraNucleonMass = 0.0;
1405  if (nuclearExcitationCount > 0) {
1406  G4int momentumBin = std::min(4, G4int(pOriginal / 3.0));
1407  const G4double nucsup[] = {1.0, 0.8, 0.6, 0.5, 0.4};
1408  //
1409  // NOTE: in TWOCLU, these new particles were given negative codes
1410  // here we use NewlyAdded = true instead
1411  //
1412  for (i = 0; i < nuclearExcitationCount; ++i) {
1413  G4ReactionProduct *pVec = new G4ReactionProduct();
1414  if (G4UniformRand() < nucsup[momentumBin]) // add proton or neutron
1415  {
1416  if (G4UniformRand() > 1.0 - atomicNumber / atomicWeight)
1417  // HPW: looks like a gheisha bug
1418  pVec->SetDefinition(aProton);
1419  else
1420  pVec->SetDefinition(aNeutron);
1421  ++backwardNucleonCount;
1422  ++extraNucleonCount;
1423  extraNucleonMass += pVec->GetMass() / GeV;
1424  } else { // add a pion
1425  G4double ran = G4UniformRand();
1426  if (ran < 0.3181)
1427  pVec->SetDefinition(aPiPlus);
1428  else if (ran < 0.6819)
1429  pVec->SetDefinition(aPiZero);
1430  else
1431  pVec->SetDefinition(aPiMinus);
1432  }
1433  pVec->SetSide(-2); // backside particle
1434  extraMass += pVec->GetMass() / GeV;
1435  pVec->SetNewlyAdded(true);
1436  vec.SetElement(vecLen++, pVec);
1437  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1438  }
1439  }
1440  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1441  G4double forwardEnergy = (centerofmassEnergy - cMass - bMass) / 2.0 + cMass - forwardMass;
1442  G4double backwardEnergy = (centerofmassEnergy - cMass - bMass) / 2.0 + bMass - backwardMass;
1443  G4double eAvailable = centerofmassEnergy - (forwardMass + backwardMass);
1444  G4bool secondaryDeleted;
1445  G4double pMass;
1446  while (eAvailable <= 0.0) // must eliminate a particle
1447  {
1448  secondaryDeleted = false;
1449  for (i = (vecLen - 1); i >= 0; --i) {
1450  if (vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled()) {
1451  pMass = vec[i]->GetMass() / GeV;
1452  for (G4int j = i; j < (vecLen - 1); ++j)
1453  *vec[j] = *vec[j + 1]; // shift up
1454  --forwardCount;
1455  forwardEnergy += pMass;
1456  forwardMass -= pMass;
1457  secondaryDeleted = true;
1458  break;
1459  } else if (vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled()) {
1460  pMass = vec[i]->GetMass() / GeV;
1461  for (G4int j = i; j < (vecLen - 1); ++j)
1462  *vec[j] = *vec[j + 1]; // shift up
1463  --backwardCount;
1464  backwardEnergy += pMass;
1465  backwardMass -= pMass;
1466  secondaryDeleted = true;
1467  break;
1468  }
1469  } // breaks go down to here
1470  if (secondaryDeleted) {
1471  G4ReactionProduct *temp = vec[vecLen - 1];
1472  delete temp;
1473  --vecLen;
1474  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1475  } else {
1476  if (vecLen == 0) {
1477  return false; // all secondaries have been eliminated
1478  }
1479  if (targetParticle.GetSide() == -1) {
1480  pMass = targetParticle.GetMass() / GeV;
1481  targetParticle = *vec[0];
1482  for (G4int j = 0; j < (vecLen - 1); ++j)
1483  *vec[j] = *vec[j + 1]; // shift up
1484  --backwardCount;
1485  backwardEnergy += pMass;
1486  backwardMass -= pMass;
1487  secondaryDeleted = true;
1488  } else if (targetParticle.GetSide() == 1) {
1489  pMass = targetParticle.GetMass() / GeV;
1490  targetParticle = *vec[0];
1491  for (G4int j = 0; j < (vecLen - 1); ++j)
1492  *vec[j] = *vec[j + 1]; // shift up
1493  --forwardCount;
1494  forwardEnergy += pMass;
1495  forwardMass -= pMass;
1496  secondaryDeleted = true;
1497  }
1498  if (secondaryDeleted) {
1499  G4ReactionProduct *temp = vec[vecLen - 1];
1500  delete temp;
1501  --vecLen;
1502  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1503  } else {
1504  if (currentParticle.GetSide() == -1) {
1505  pMass = currentParticle.GetMass() / GeV;
1506  currentParticle = *vec[0];
1507  for (G4int j = 0; j < (vecLen - 1); ++j)
1508  *vec[j] = *vec[j + 1]; // shift up
1509  --backwardCount;
1510  backwardEnergy += pMass;
1511  backwardMass -= pMass;
1512  secondaryDeleted = true;
1513  } else if (currentParticle.GetSide() == 1) {
1514  pMass = currentParticle.GetMass() / GeV;
1515  currentParticle = *vec[0];
1516  for (G4int j = 0; j < (vecLen - 1); ++j)
1517  *vec[j] = *vec[j + 1]; // shift up
1518  --forwardCount; //This line can cause infinite loop
1519  forwardEnergy += pMass;
1520  forwardMass -= pMass;
1521  secondaryDeleted = true;
1522  }
1523  if (secondaryDeleted) {
1524  G4ReactionProduct *temp = vec[vecLen - 1];
1525  delete temp;
1526  --vecLen;
1527  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1528  } else
1529  break;
1530  }
1531  }
1532  eAvailable = centerofmassEnergy - (forwardMass + backwardMass);
1533  }
1534  //
1535  // This is the start of the TwoCluster function
1536  // Choose masses for the 3 clusters:
1537  // forward cluster
1538  // backward meson cluster
1539  // backward nucleon cluster
1540  //
1541  G4double rmc = 0.0, rmd = 0.0; // rme = 0.0;
1542  const G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
1543  const G4double gpar[] = {2.6, 2.6, 1.8, 1.30, 1.20};
1544 
1545  if (forwardCount == 0)
1546  return false;
1547 
1548  if (forwardCount == 1)
1549  rmc = forwardMass;
1550  else {
1551  G4int ntc = std::max(1, std::min(5, forwardCount)) - 1; // check if offset by 1 @@
1552  rmc = forwardMass + std::pow(-std::log(1.0 - G4UniformRand()), cpar[ntc]) / gpar[ntc];
1553  }
1554  if (backwardCount == 1)
1555  rmd = backwardMass;
1556  else {
1557  G4int ntc = std::max(1, std::min(5, backwardCount)) - 1; // check, if offfset by 1 @@
1558  rmd = backwardMass + std::pow(-std::log(1.0 - G4UniformRand()), cpar[ntc]) / gpar[ntc];
1559  }
1560  while (rmc + rmd > centerofmassEnergy) {
1561  if ((rmc <= forwardMass) && (rmd <= backwardMass)) {
1562  G4double temp = 0.999 * centerofmassEnergy / (rmc + rmd);
1563  rmc *= temp;
1564  rmd *= temp;
1565  } else {
1566  rmc = 0.1 * forwardMass + 0.9 * rmc;
1567  rmd = 0.1 * backwardMass + 0.9 * rmd;
1568  }
1569  }
1570  //
1571  // Set beam, target of first interaction in centre of mass system
1572  //
1573  G4ReactionProduct pseudoParticle[8];
1574  for (i = 0; i < 8; ++i)
1575  pseudoParticle[i].SetZero();
1576 
1577  pseudoParticle[1].SetMass(mOriginal * GeV);
1578  pseudoParticle[1].SetTotalEnergy(etOriginal * GeV);
1579  pseudoParticle[1].SetMomentum(0.0, 0.0, pOriginal * GeV);
1580 
1581  pseudoParticle[2].SetMass(protonMass * MeV);
1582  pseudoParticle[2].SetTotalEnergy(protonMass * MeV);
1583  pseudoParticle[2].SetMomentum(0.0, 0.0, 0.0);
1584  //
1585  // transform into centre of mass system
1586  //
1587  pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
1588  pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[0]);
1589  pseudoParticle[2].Lorentz(pseudoParticle[2], pseudoParticle[0]);
1590 
1591  const G4double pfMin = 0.0001;
1592  G4double pf = (centerofmassEnergy * centerofmassEnergy + rmd * rmd - rmc * rmc);
1593  pf *= pf;
1594  pf -= 4 * centerofmassEnergy * centerofmassEnergy * rmd * rmd;
1595  pf = std::sqrt(std::max(pf, pfMin)) / (2.0 * centerofmassEnergy);
1596  //
1597  // set final state masses and energies in centre of mass system
1598  //
1599  pseudoParticle[3].SetMass(rmc * GeV);
1600  pseudoParticle[3].SetTotalEnergy(std::sqrt(pf * pf + rmc * rmc) * GeV);
1601 
1602  pseudoParticle[4].SetMass(rmd * GeV);
1603  pseudoParticle[4].SetTotalEnergy(std::sqrt(pf * pf + rmd * rmd) * GeV);
1604  //
1605  // set |T| and |TMIN|
1606  //
1607  const G4double bMin = 0.01;
1608  const G4double b1 = 4.0;
1609  const G4double b2 = 1.6;
1610  G4double t = std::log(1.0 - G4UniformRand()) / std::max(bMin, b1 + b2 * std::log(pOriginal));
1611  G4double t1 = pseudoParticle[1].GetTotalEnergy() / GeV - pseudoParticle[3].GetTotalEnergy() / GeV;
1612  G4double pin = pseudoParticle[1].GetMomentum().mag() / GeV;
1613  G4double tacmin = t1 * t1 - (pin - pf) * (pin - pf);
1614  //
1615  // calculate (std::sin(teta/2.)^2 and std::cos(teta), set azimuth angle phi
1616  //
1617  const G4double smallValue = 1.0e-10;
1618  G4double dumnve = 4.0 * pin * pf;
1619  if (dumnve == 0.0)
1620  dumnve = smallValue;
1621  G4double ctet = std::max(-1.0, std::min(1.0, 1.0 + 2.0 * (t - tacmin) / dumnve));
1622  dumnve = std::max(0.0, 1.0 - ctet * ctet);
1623  G4double stet = std::sqrt(dumnve);
1624  G4double phi = G4UniformRand() * twopi;
1625  //
1626  // calculate final state momenta in centre of mass system
1627  //
1628  pseudoParticle[3].SetMomentum(pf * stet * std::sin(phi) * GeV, pf * stet * std::cos(phi) * GeV, pf * ctet * GeV);
1629  pseudoParticle[4].SetMomentum(pseudoParticle[3].GetMomentum() * (-1.0));
1630  //
1631  // simulate backward nucleon cluster in lab. system and transform in cms
1632  //
1633  G4double pp, pp1, rthnve, phinve;
1634  if (nuclearExcitationCount > 0) {
1635  const G4double ga = 1.2;
1636  G4double ekit1 = 0.04;
1637  G4double ekit2 = 0.6;
1638  if (ekOriginal <= 5.0) {
1639  ekit1 *= ekOriginal * ekOriginal / 25.0;
1640  ekit2 *= ekOriginal * ekOriginal / 25.0;
1641  }
1642  const G4double a = (1.0 - ga) / (std::pow(ekit2, (1.0 - ga)) - std::pow(ekit1, (1.0 - ga)));
1643  for (i = 0; i < vecLen; ++i) {
1644  if (vec[i]->GetSide() == -2) {
1645  G4double kineticE =
1646  std::pow((G4UniformRand() * (1.0 - ga) / a + std::pow(ekit1, (1.0 - ga))), (1.0 / (1.0 - ga)));
1647  vec[i]->SetKineticEnergy(kineticE * GeV);
1648  G4double vMass = vec[i]->GetMass() / MeV;
1649  G4double totalE = kineticE + vMass;
1650  pp = std::sqrt(std::abs(totalE * totalE - vMass * vMass));
1651  G4double cost = std::min(1.0, std::max(-1.0, std::log(2.23 * G4UniformRand() + 0.383) / 0.96));
1652  G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
1653  phi = twopi * G4UniformRand();
1654  vec[i]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
1655  vec[i]->Lorentz(*vec[i], pseudoParticle[0]);
1656  }
1657  }
1658  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1659  }
1660  //
1661  // fragmentation of forward cluster and backward meson cluster
1662  //
1663  currentParticle.SetMomentum(pseudoParticle[3].GetMomentum());
1664  currentParticle.SetTotalEnergy(pseudoParticle[3].GetTotalEnergy());
1665 
1666  targetParticle.SetMomentum(pseudoParticle[4].GetMomentum());
1667  targetParticle.SetTotalEnergy(pseudoParticle[4].GetTotalEnergy());
1668 
1669  pseudoParticle[5].SetMomentum(pseudoParticle[3].GetMomentum() * (-1.0));
1670  pseudoParticle[5].SetMass(pseudoParticle[3].GetMass());
1671  pseudoParticle[5].SetTotalEnergy(pseudoParticle[3].GetTotalEnergy());
1672 
1673  pseudoParticle[6].SetMomentum(pseudoParticle[4].GetMomentum() * (-1.0));
1674  pseudoParticle[6].SetMass(pseudoParticle[4].GetMass());
1675  pseudoParticle[6].SetTotalEnergy(pseudoParticle[4].GetTotalEnergy());
1676 
1677  G4double wgt;
1678  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1679  if (forwardCount > 1) // tempV will contain the forward particles
1680  {
1681  G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1682  tempV.Initialize(forwardCount);
1683  G4bool constantCrossSection = true;
1684  G4int tempLen = 0;
1685  if (currentParticle.GetSide() == 1)
1686  tempV.SetElement(tempLen++, &currentParticle);
1687  if (targetParticle.GetSide() == 1)
1688  tempV.SetElement(tempLen++, &targetParticle);
1689  for (i = 0; i < vecLen; ++i) {
1690  if (vec[i]->GetSide() == 1) {
1691  if (tempLen < 18)
1692  tempV.SetElement(tempLen++, vec[i]);
1693  else {
1694  vec[i]->SetSide(-1);
1695  continue;
1696  }
1697  }
1698  }
1699  if (tempLen >= 2) {
1700  GenerateNBodyEvent(pseudoParticle[3].GetMass() / MeV, constantCrossSection, tempV, tempLen);
1701  if (currentParticle.GetSide() == 1)
1702  currentParticle.Lorentz(currentParticle, pseudoParticle[5]);
1703  if (targetParticle.GetSide() == 1)
1704  targetParticle.Lorentz(targetParticle, pseudoParticle[5]);
1705  for (i = 0; i < vecLen; ++i) {
1706  if (vec[i]->GetSide() == 1)
1707  vec[i]->Lorentz(*vec[i], pseudoParticle[5]);
1708  }
1709  }
1710  }
1711  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1712  if (backwardCount > 1) // tempV will contain the backward particles,
1713  { // but not those created from the intranuclear cascade
1714  G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1715  tempV.Initialize(backwardCount);
1716  G4bool constantCrossSection = true;
1717  G4int tempLen = 0;
1718  if (currentParticle.GetSide() == -1)
1719  tempV.SetElement(tempLen++, &currentParticle);
1720  if (targetParticle.GetSide() == -1)
1721  tempV.SetElement(tempLen++, &targetParticle);
1722  for (i = 0; i < vecLen; ++i) {
1723  if (vec[i]->GetSide() == -1) {
1724  if (tempLen < 18)
1725  tempV.SetElement(tempLen++, vec[i]);
1726  else {
1727  vec[i]->SetSide(-2);
1728  vec[i]->SetKineticEnergy(0.0);
1729  vec[i]->SetMomentum(0.0, 0.0, 0.0);
1730  continue;
1731  }
1732  }
1733  }
1734  if (tempLen >= 2) {
1735  GenerateNBodyEvent(pseudoParticle[4].GetMass() / MeV, constantCrossSection, tempV, tempLen);
1736  if (currentParticle.GetSide() == -1)
1737  currentParticle.Lorentz(currentParticle, pseudoParticle[6]);
1738  if (targetParticle.GetSide() == -1)
1739  targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
1740  for (i = 0; i < vecLen; ++i) {
1741  if (vec[i]->GetSide() == -1)
1742  vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
1743  }
1744  }
1745  }
1746  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1747  //
1748  // Lorentz transformation in lab system
1749  //
1750  G4int numberofFinalStateNucleons = 0;
1751  if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron ||
1752  currentParticle.GetDefinition() == aSigmaMinus || currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
1753  currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
1754  currentParticle.GetDefinition() == G4XiZero::XiZero() ||
1755  currentParticle.GetDefinition() == G4XiMinus::XiMinus() ||
1756  currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1757  currentParticle.GetDefinition() == G4Lambda::Lambda())
1758  ++numberofFinalStateNucleons;
1759  currentParticle.Lorentz(currentParticle, pseudoParticle[2]);
1760  if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron ||
1761  targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() ||
1762  targetParticle.GetDefinition() == G4XiMinus::XiMinus() ||
1763  targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1764  targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
1765  targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || targetParticle.GetDefinition() == aSigmaMinus)
1766  ++numberofFinalStateNucleons;
1767  if (targetParticle.GetDefinition() == G4AntiProton::AntiProton())
1768  --numberofFinalStateNucleons;
1769  if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron())
1770  --numberofFinalStateNucleons;
1771  if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
1772  --numberofFinalStateNucleons;
1773  if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
1774  --numberofFinalStateNucleons;
1775  if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
1776  --numberofFinalStateNucleons;
1777  if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero())
1778  --numberofFinalStateNucleons;
1779  if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus())
1780  --numberofFinalStateNucleons;
1781  if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
1782  --numberofFinalStateNucleons;
1783  if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda())
1784  --numberofFinalStateNucleons;
1785  targetParticle.Lorentz(targetParticle, pseudoParticle[2]);
1786  for (i = 0; i < vecLen; ++i) {
1787  if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron ||
1788  vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() ||
1789  vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1790  vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() ||
1791  vec[i]->GetDefinition() == aSigmaMinus)
1792  ++numberofFinalStateNucleons;
1793  if (vec[i]->GetDefinition() == G4AntiProton::AntiProton())
1794  --numberofFinalStateNucleons;
1795  if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron())
1796  --numberofFinalStateNucleons;
1797  if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
1798  --numberofFinalStateNucleons;
1799  if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
1800  --numberofFinalStateNucleons;
1801  if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
1802  --numberofFinalStateNucleons;
1803  if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda())
1804  --numberofFinalStateNucleons;
1805  if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero())
1806  --numberofFinalStateNucleons;
1807  if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus())
1808  --numberofFinalStateNucleons;
1809  if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
1810  --numberofFinalStateNucleons;
1811  vec[i]->Lorentz(*vec[i], pseudoParticle[2]);
1812  }
1813  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1814  numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
1815  //
1816  // sometimes the leading strange particle is lost, set it back
1817  //
1818  G4bool dum = true;
1819  if (leadFlag) {
1820  // leadFlag will be true
1821  // iff original particle is at least as heavy as K+ and not a proton or neutron AND
1822  // if
1823  // incident particle is at least as heavy as K+ and it is not a proton or neutron
1824  // leadFlag is set to the incident particle
1825  // or
1826  // target particle is at least as heavy as K+ and it is not a proton or neutron
1827  // leadFlag is set to the target particle
1828  //
1829  if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1830  dum = false;
1831  else if (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1832  dum = false;
1833  else {
1834  for (i = 0; i < vecLen; ++i) {
1835  if (vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition()) {
1836  dum = false;
1837  break;
1838  }
1839  }
1840  }
1841  if (dum) {
1842  G4double leadMass = leadingStrangeParticle.GetMass() / MeV;
1843  G4double ekin;
1844  if (((leadMass < protonMass) && (targetParticle.GetMass() / MeV < protonMass)) ||
1845  ((leadMass >= protonMass) && (targetParticle.GetMass() / MeV >= protonMass))) {
1846  ekin = targetParticle.GetKineticEnergy() / GeV;
1847  pp1 = targetParticle.GetMomentum().mag() / MeV; // old momentum
1848  targetParticle.SetDefinition(leadingStrangeParticle.GetDefinition());
1849  targetParticle.SetKineticEnergy(ekin * GeV);
1850  pp = targetParticle.GetTotalMomentum() / MeV; // new momentum
1851  if (pp1 < 1.0e-3) {
1852  rthnve = pi * G4UniformRand();
1853  phinve = twopi * G4UniformRand();
1854  targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1855  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1856  pp * std::cos(rthnve) * MeV);
1857  } else
1858  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1859 
1860  targetHasChanged = true;
1861  } else {
1862  ekin = currentParticle.GetKineticEnergy() / GeV;
1863  pp1 = currentParticle.GetMomentum().mag() / MeV;
1864  currentParticle.SetDefinition(leadingStrangeParticle.GetDefinition());
1865  currentParticle.SetKineticEnergy(ekin * GeV);
1866  pp = currentParticle.GetTotalMomentum() / MeV;
1867  if (pp1 < 1.0e-3) {
1868  rthnve = pi * G4UniformRand();
1869  phinve = twopi * G4UniformRand();
1870  currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1871  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1872  pp * std::cos(rthnve) * MeV);
1873  } else
1874  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1875 
1876  incidentHasChanged = true;
1877  }
1878  }
1879  } // end of if( leadFlag )
1880  //
1881  // for various reasons, the energy balance is not sufficient,
1882  // check that, energy balance, angle of final system, etc.
1883  //
1884  pseudoParticle[4].SetMass(mOriginal * GeV);
1885  pseudoParticle[4].SetTotalEnergy(etOriginal * GeV);
1886  pseudoParticle[4].SetMomentum(0.0, 0.0, pOriginal * GeV);
1887 
1888  const G4ParticleDefinition *aOrgDef = modifiedOriginal.GetDefinition();
1889  G4int diff = 0;
1890  if (aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron())
1891  diff = 1;
1892  if (numberofFinalStateNucleons == 1)
1893  diff = 0;
1894  pseudoParticle[5].SetMomentum(0.0, 0.0, 0.0);
1895  pseudoParticle[5].SetMass(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1896  pseudoParticle[5].SetTotalEnergy(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1897 
1898  // G4double ekin0 = pseudoParticle[4].GetKineticEnergy()/GeV;
1899  G4double theoreticalKinetic = pseudoParticle[4].GetTotalEnergy() / GeV + pseudoParticle[5].GetTotalEnergy() / GeV;
1900 
1901  pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
1902  pseudoParticle[4].Lorentz(pseudoParticle[4], pseudoParticle[6]);
1903  pseudoParticle[5].Lorentz(pseudoParticle[5], pseudoParticle[6]);
1904 
1905  if (vecLen < 16) {
1906  G4ReactionProduct tempR[130];
1907  //G4ReactionProduct *tempR = new G4ReactionProduct [vecLen+2];
1908  tempR[0] = currentParticle;
1909  tempR[1] = targetParticle;
1910  for (i = 0; i < vecLen; ++i)
1911  tempR[i + 2] = *vec[i];
1912 
1913  G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1914  tempV.Initialize(vecLen + 2);
1915  G4bool constantCrossSection = true;
1916  G4int tempLen = 0;
1917  for (i = 0; i < vecLen + 2; ++i)
1918  tempV.SetElement(tempLen++, &tempR[i]);
1919 
1920  if (tempLen >= 2) {
1921  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1922  GenerateNBodyEvent(pseudoParticle[4].GetTotalEnergy() / MeV + pseudoParticle[5].GetTotalEnergy() / MeV,
1923  constantCrossSection,
1924  tempV,
1925  tempLen);
1926  theoreticalKinetic = 0.0;
1927  for (i = 0; i < vecLen + 2; ++i) {
1928  pseudoParticle[7].SetMomentum(tempV[i]->GetMomentum());
1929  pseudoParticle[7].SetMass(tempV[i]->GetMass());
1930  pseudoParticle[7].SetTotalEnergy(tempV[i]->GetTotalEnergy());
1931  pseudoParticle[7].Lorentz(pseudoParticle[7], pseudoParticle[5]);
1932  theoreticalKinetic += pseudoParticle[7].GetKineticEnergy() / GeV;
1933  }
1934  }
1935  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1936  //delete [] tempR;
1937  } else {
1938  theoreticalKinetic -= (currentParticle.GetMass() / GeV + targetParticle.GetMass() / GeV);
1939  for (i = 0; i < vecLen; ++i)
1940  theoreticalKinetic -= vec[i]->GetMass() / GeV;
1941  }
1942  G4double simulatedKinetic = currentParticle.GetKineticEnergy() / GeV + targetParticle.GetKineticEnergy() / GeV;
1943  for (i = 0; i < vecLen; ++i)
1944  simulatedKinetic += vec[i]->GetKineticEnergy() / GeV;
1945  //
1946  // make sure that kinetic energies are correct
1947  // the backward nucleon cluster is not produced within proper kinematics!!!
1948  //
1949 
1950  if (simulatedKinetic != 0.0) {
1951  wgt = (theoreticalKinetic) / simulatedKinetic;
1952  currentParticle.SetKineticEnergy(wgt * currentParticle.GetKineticEnergy());
1953  pp = currentParticle.GetTotalMomentum() / MeV;
1954  pp1 = currentParticle.GetMomentum().mag() / MeV;
1955  if (pp1 < 0.001 * MeV) {
1956  rthnve = pi * G4UniformRand();
1957  phinve = twopi * G4UniformRand();
1958  currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1959  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1960  pp * std::cos(rthnve) * MeV);
1961  } else
1962  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1963 
1964  targetParticle.SetKineticEnergy(wgt * targetParticle.GetKineticEnergy());
1965  pp = targetParticle.GetTotalMomentum() / MeV;
1966  pp1 = targetParticle.GetMomentum().mag() / MeV;
1967  if (pp1 < 0.001 * MeV) {
1968  rthnve = pi * G4UniformRand();
1969  phinve = twopi * G4UniformRand();
1970  targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1971  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1972  pp * std::cos(rthnve) * MeV);
1973  } else
1974  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1975 
1976  for (i = 0; i < vecLen; ++i) {
1977  vec[i]->SetKineticEnergy(wgt * vec[i]->GetKineticEnergy());
1978  pp = vec[i]->GetTotalMomentum() / MeV;
1979  pp1 = vec[i]->GetMomentum().mag() / MeV;
1980  if (pp1 < 0.001) {
1981  rthnve = pi * G4UniformRand();
1982  phinve = twopi * G4UniformRand();
1983  vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1984  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1985  pp * std::cos(rthnve) * MeV);
1986  } else
1987  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
1988  }
1989  }
1990  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1991  Rotate(numberofFinalStateNucleons,
1992  pseudoParticle[4].GetMomentum(),
1993  modifiedOriginal,
1994  originalIncident,
1995  targetNucleus,
1996  currentParticle,
1997  targetParticle,
1998  vec,
1999  vecLen);
2000  //
2001  // add black track particles
2002  // the total number of particles produced is restricted to 198
2003  // this may have influence on very high energies
2004  //
2005  if (atomicWeight >= 1.5) {
2006  // npnb is number of proton/neutron black track particles
2007  // ndta is the number of deuterons, tritons, and alphas produced
2008  // epnb is the kinetic energy available for proton/neutron black track particles
2009  // edta is the kinetic energy available for deuteron/triton/alpha particles
2010  //
2011  G4double epnb, edta;
2012  G4int npnb = 0;
2013  G4int ndta = 0;
2014 
2015  epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
2016  edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
2017  const G4double pnCutOff = 0.001; // GeV
2018  const G4double dtaCutOff = 0.001; // GeV
2019  const G4double kineticMinimum = 1.e-6;
2020  const G4double kineticFactor = -0.005;
2021 
2022  G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
2023  const G4double ekIncident = originalIncident->GetKineticEnergy() / GeV;
2024  if (ekIncident >= 5.0)
2025  sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0));
2026 
2027  if (epnb >= pnCutOff) {
2028  npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta));
2029  if (numberofFinalStateNucleons + npnb > atomicWeight)
2030  npnb = G4int(atomicWeight - numberofFinalStateNucleons);
2031  npnb = std::min(npnb, 127 - vecLen);
2032  }
2033  if (edta >= dtaCutOff) {
2034  ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta));
2035  ndta = std::min(ndta, 127 - vecLen);
2036  }
2037  G4double spall = numberofFinalStateNucleons;
2038  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2039 
2040  AddBlackTrackParticles(epnb,
2041  npnb,
2042  edta,
2043  ndta,
2044  sprob,
2045  kineticMinimum,
2046  kineticFactor,
2047  modifiedOriginal,
2048  spall,
2049  targetNucleus,
2050  vec,
2051  vecLen);
2052 
2053  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2054  }
2055  //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
2056  // MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2057  //
2058  // calculate time delay for nuclear reactions
2059  //
2060  if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
2061  currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
2062  else
2063  currentParticle.SetTOF(1.0);
2064 
2065  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2066  return true;
2067 }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:338
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
pi
const Double_t pi
Definition: trackSplitPlot.h:36
MeV
const double MeV
protons_cff.t
t
Definition: protons_cff.py:34
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
FullModelReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2340
FullModelReactionDynamics::AddBlackTrackParticles
void AddBlackTrackParticles(const G4double epnb, const G4int npnb, const G4double edta, const G4int ndta, const G4double sprob, const G4double kineticMinimum, const G4double kineticFactor, const G4ReactionProduct &modifiedOriginal, G4double spall, const G4Nucleus &aNucleus, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2850
FullModelReactionDynamics::Rotate
void Rotate(const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2639
a
double a
Definition: hdecay.h:121
BPhysicsValidation_cfi.Lambda
Lambda
Definition: BPhysicsValidation_cfi.py:63
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
phase2TrackerDigitizer_cfi.SigmaZero
SigmaZero
Definition: phase2TrackerDigitizer_cfi.py:46
Generated for CMSSW Reference Manual by doxygen 1.8.11