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::FullModelReactionDynamics ( )
inline

Definition at line 25 of file FullModelReactionDynamics.h.

25 {}

◆ ~FullModelReactionDynamics()

virtual FullModelReactionDynamics::~FullModelReactionDynamics ( )
inlinevirtual

Definition at line 27 of file FullModelReactionDynamics.h.

27 {}

Member Function Documentation

◆ AddBlackTrackParticles()

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 2848 of file FullModelReactionDynamics.cc.

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

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

◆ Defs1()

void FullModelReactionDynamics::Defs1 ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2586 of file FullModelReactionDynamics.cc.

2590  {
2591  const G4double pjx = modifiedOriginal.GetMomentum().x() / MeV;
2592  const G4double pjy = modifiedOriginal.GetMomentum().y() / MeV;
2593  const G4double pjz = modifiedOriginal.GetMomentum().z() / MeV;
2594  const G4double p = modifiedOriginal.GetMomentum().mag() / MeV;
2595  if (pjx * pjx + pjy * pjy > 0.0) {
2596  G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2597  cost = pjz / p;
2598  sint = 0.5 * (std::sqrt(std::abs((1.0 - cost) * (1.0 + cost))) + std::sqrt(pjx * pjx + pjy * pjy) / p);
2599  if (pjy < 0.0)
2600  ph = 3 * halfpi;
2601  else
2602  ph = halfpi;
2603  if (std::abs(pjx) > 0.001 * MeV)
2604  ph = std::atan2(pjy, pjx);
2605  cosp = std::cos(ph);
2606  sinp = std::sin(ph);
2607  pix = currentParticle.GetMomentum().x() / MeV;
2608  piy = currentParticle.GetMomentum().y() / MeV;
2609  piz = currentParticle.GetMomentum().z() / MeV;
2610  currentParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2611  cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2612  -sint * pix * MeV + cost * piz * MeV);
2613  pix = targetParticle.GetMomentum().x() / MeV;
2614  piy = targetParticle.GetMomentum().y() / MeV;
2615  piz = targetParticle.GetMomentum().z() / MeV;
2616  targetParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2617  cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2618  -sint * pix * MeV + cost * piz * MeV);
2619  for (G4int i = 0; i < vecLen; ++i) {
2620  pix = vec[i]->GetMomentum().x() / MeV;
2621  piy = vec[i]->GetMomentum().y() / MeV;
2622  piz = vec[i]->GetMomentum().z() / MeV;
2623  vec[i]->SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2624  cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2625  -sint * pix * MeV + cost * piz * MeV);
2626  }
2627  } else {
2628  if (pjz < 0.0) {
2629  currentParticle.SetMomentum(-currentParticle.GetMomentum().z());
2630  targetParticle.SetMomentum(-targetParticle.GetMomentum().z());
2631  for (G4int i = 0; i < vecLen; ++i)
2632  vec[i]->SetMomentum(-vec[i]->GetMomentum().z());
2633  }
2634  }
2635 }

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

◆ Factorial()

G4int FullModelReactionDynamics::Factorial ( G4int  n)

Definition at line 2576 of file FullModelReactionDynamics.cc.

2576  { // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
2577  G4int m = std::min(n, 10);
2578  G4int result = 1;
2579  if (m <= 1)
2580  return result;
2581  for (G4int i = 2; i <= m; ++i)
2582  result *= i;
2583  return result;
2584 }

References mps_fire::i, visualization-live-secondInstance_cfg::m, min(), dqmiodumpmetadata::n, and mps_fire::result.

◆ FindInelasticity()

virtual G4double FullModelReactionDynamics::FindInelasticity ( )
inlinevirtual

Definition at line 29 of file FullModelReactionDynamics.h.

29 { return 0.0; }

◆ FindTimeDelay()

virtual G4double FullModelReactionDynamics::FindTimeDelay ( )
inlinevirtual

Definition at line 31 of file FullModelReactionDynamics.h.

31 { return 0.0; }

◆ GenerateNBodyEvent()

G4double FullModelReactionDynamics::GenerateNBodyEvent ( const G4double  totalEnergy,
const G4bool  constantCrossSection,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)

Definition at line 2338 of file FullModelReactionDynamics.cc.

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

References a, b, HLT_FULL_cff::beta, c, funct::cos(), HCALHighEnergyHPDFilter_cfi::energy, JetChargeProducer_cfi::exp, ecalTB2006H4_GenSimDigiReco_cfg::G4cerr, GeV, mps_fire::i, dqmiolumiharvest::j, dqm-mbProfile::log, EgHLTOffHistBins_cfi::mass, SiStripPI::max, min(), funct::pow(), summarizeEdmComparisonLogfiles::ran, alignCSCRings::s, funct::sin(), mathSSE::sqrt(), groupFilesInBlocks::temp, and mps_merge::weight.

◆ GenerateXandPt()

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 60 of file FullModelReactionDynamics.cc.

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

References funct::abs(), funct::cos(), change_name::diff, MillePedeFileConverter_cfg::e, EgHLTOffHistBins_cfi::et, beamvalidation::exit(), JetChargeProducer_cfi::exp, ecalTB2006H4_GenSimDigiReco_cfg::G4cerr, ecalTB2006H4_GenSimDigiReco_cfg::G4cout, GeV, mps_fire::i, dqmiolumiharvest::j, cmsLHEtoEOSManager::l, BPhysicsValidation_cfi::Lambda, dqm-mbProfile::log, visualization-live-secondInstance_cfg::m, SiStripPI::max, MeV, min(), pi, funct::pow(), createTree::pp, DiDispStaMuonMonitor_cfi::pt, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, summarizeEdmComparisonLogfiles::ran, phase2TrackerDigitizer_cfi::SigmaZero, funct::sin(), sqr(), mathSSE::sqrt(), groupFilesInBlocks::temp, and testFFTrackingParticleHits::xxx.

Referenced by FullModelHadronicProcess::CalculateMomenta().

◆ MomentumCheck()

void FullModelReactionDynamics::MomentumCheck ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2996 of file FullModelReactionDynamics.cc.

3000  {
3001  const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / MeV;
3002  G4double testMomentum = currentParticle.GetMomentum().mag() / MeV;
3003  G4double pMass;
3004  if (testMomentum >= pOriginal) {
3005  pMass = currentParticle.GetMass() / MeV;
3006  currentParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3007  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pOriginal / testMomentum));
3008  }
3009  testMomentum = targetParticle.GetMomentum().mag() / MeV;
3010  if (testMomentum >= pOriginal) {
3011  pMass = targetParticle.GetMass() / MeV;
3012  targetParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3013  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pOriginal / testMomentum));
3014  }
3015  for (G4int i = 0; i < vecLen; ++i) {
3016  testMomentum = vec[i]->GetMomentum().mag() / MeV;
3017  if (testMomentum >= pOriginal) {
3018  pMass = vec[i]->GetMass() / MeV;
3019  vec[i]->SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3020  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pOriginal / testMomentum));
3021  }
3022  }
3023 }

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

◆ normal()

G4double FullModelReactionDynamics::normal ( )
private

Definition at line 2523 of file FullModelReactionDynamics.cc.

2523  {
2524  G4double ran = -6.0;
2525  for (G4int i = 0; i < 12; ++i)
2526  ran += G4UniformRand();
2527  return ran;
2528 }

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

◆ NuclearReaction()

void FullModelReactionDynamics::NuclearReaction ( G4FastVector< G4ReactionProduct, 4 > &  vec,
G4int &  vecLen,
const G4HadProjectile *  originalIncident,
const G4Nucleus &  aNucleus,
const G4double  theAtomicMass,
const G4double *  massVec 
)

Definition at line 3392 of file FullModelReactionDynamics.cc.

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

References VtxSmearedParameters_cfi::Alpha, funct::cos(), GeV, mps_fire::i, ke, EgHLTOffHistBins_cfi::mass, SiStripPI::max, MeV, min(), nt, AlCaHLTBitMon_ParallelJobs::p, createTree::pp, summarizeEdmComparisonLogfiles::ran, funct::sin(), and findQualityFiles::v.

◆ Poisson()

G4int FullModelReactionDynamics::Poisson ( G4double  x)
private

Definition at line 2530 of file FullModelReactionDynamics.cc.

2531 {
2532  G4int iran;
2533  G4double ran;
2534 
2535  if (x > 9.9) // use normal distribution with sigma^2 = <x>
2536  iran = static_cast<G4int>(std::max(0.0, x + normal() * std::sqrt(x)));
2537  else {
2538  G4int mm = G4int(5.0 * x);
2539  if (mm <= 0) // for very small x try iran=1,2,3
2540  {
2541  G4double p1 = x * std::exp(-x);
2542  G4double p2 = x * p1 / 2.0;
2543  G4double p3 = x * p2 / 3.0;
2544  ran = G4UniformRand();
2545  if (ran < p3)
2546  iran = 3;
2547  else if (ran < p2) // this is original Geisha, it should be ran < p2+p3
2548  iran = 2;
2549  else if (ran < p1) // should be ran < p1+p2+p3
2550  iran = 1;
2551  else
2552  iran = 0;
2553  } else {
2554  iran = 0;
2555  G4double r = std::exp(-x);
2556  ran = G4UniformRand();
2557  if (ran > r) {
2558  G4double rrr;
2559  G4double rr = r;
2560  for (G4int i = 1; i <= mm; ++i) {
2561  iran++;
2562  if (i > 5) // Stirling's formula for large numbers
2563  rrr = std::exp(i * std::log(x) - (i + 0.5) * std::log((G4double)i) + i - 0.9189385);
2564  else
2565  rrr = std::pow(x, i) / Factorial(i);
2566  rr += r * rrr;
2567  if (ran <= rr)
2568  break;
2569  }
2570  }
2571  }
2572  }
2573  return iran;
2574 }

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

◆ ProduceStrangeParticlePairs()

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 3025 of file FullModelReactionDynamics.cc.

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

References JetChargeProducer_cfi::exp, GeV, mps_fire::i, testProducerWithPsetDescEmpty_cfi::i3, dqmiolumiharvest::j, BPhysicsValidation_cfi::Lambda, dqm-mbProfile::log, SiStripPI::max, p1, summarizeEdmComparisonLogfiles::ran, phase2TrackerDigitizer_cfi::SigmaZero, and mathSSE::sqrt().

◆ Rotate()

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 2637 of file FullModelReactionDynamics.cc.

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

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

◆ SuppressChargedPions()

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 1232 of file FullModelReactionDynamics.cc.

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

References GeV, and mps_fire::i.

Referenced by FullModelHadronicProcess::CalculateMomenta().

◆ TwoBody()

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 2067 of file FullModelReactionDynamics.cc.

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

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

Referenced by FullModelHadronicProcess::CalculateMomenta().

◆ TwoCluster()

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 1285 of file FullModelReactionDynamics.cc.

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

References a, funct::abs(), b1, b2, amptDefault_cfi::bMin, funct::cos(), change_name::diff, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, GeV, mps_fire::i, dqmiolumiharvest::j, BPhysicsValidation_cfi::Lambda, dqm-mbProfile::log, visualization-live-secondInstance_cfg::m, SiStripPI::max, MeV, min(), packedPFCandidateRefMixer_cfi::pf, pi, funct::pow(), createTree::pp, summarizeEdmComparisonLogfiles::ran, phase2TrackerDigitizer_cfi::SigmaZero, funct::sin(), mathSSE::sqrt(), submitPVValidationJobs::t, RandomServiceHelper::t1, and groupFilesInBlocks::temp.

Referenced by FullModelHadronicProcess::CalculateMomenta().

change_name.diff
diff
Definition: change_name.py:13
FullModelReactionDynamics::Poisson
G4int Poisson(G4double x)
Definition: FullModelReactionDynamics.cc:2530
DDAxes::y
VtxSmearedParameters_cfi.Alpha
Alpha
Definition: VtxSmearedParameters_cfi.py:115
testProducerWithPsetDescEmpty_cfi.i3
i3
Definition: testProducerWithPsetDescEmpty_cfi.py:47
mps_fire.i
i
Definition: mps_fire.py:428
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
nt
int nt
Definition: AMPTWrapper.h:42
summarizeEdmComparisonLogfiles.ran
ran
Definition: summarizeEdmComparisonLogfiles.py:113
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
min
T min(T a, T b)
Definition: MathUtil.h:58
HLT_FULL_cff.beta
beta
Definition: HLT_FULL_cff.py:8651
findQualityFiles.rr
string rr
Definition: findQualityFiles.py:185
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
mps_merge.weight
weight
Definition: mps_merge.py:88
testFFTrackingParticleHits.xxx
xxx
Definition: testFFTrackingParticleHits.py:31
sqr
int sqr(const T &t)
Definition: pfalgo_common_ref.h:9
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:2848
MeV
const double MeV
CombinatorialSeedGeneratorForBeamHalo_cfi.SetMomentum
SetMomentum
Definition: CombinatorialSeedGeneratorForBeamHalo_cfi.py:99
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:2637
hgcal_conditions::parameters
Definition: HGCConditions.h:86
DDAxes::x
findQualityFiles.v
v
Definition: findQualityFiles.py:179
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
FullModelReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2338
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
FullModelReactionDynamics::normal
G4double normal()
Definition: FullModelReactionDynamics.cc:2523
fermi
const double fermi
Definition: MathUtil.h:17
BPhysicsValidation_cfi.Lambda
Lambda
Definition: BPhysicsValidation_cfi.py:63
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
alignCSCRings.s
s
Definition: alignCSCRings.py:92
ecalTB2006H4_GenSimDigiReco_cfg.G4cerr
G4cerr
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:286
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
phase2TrackerDigitizer_cfi.SigmaZero
SigmaZero
Definition: phase2TrackerDigitizer_cfi.py:28
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:78
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
p2
double p2[4]
Definition: TauolaWrapper.h:90
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
b
double b
Definition: hdecay.h:118
idealTransformation.rotation
dictionary rotation
Definition: idealTransformation.py:1
a
double a
Definition: hdecay.h:119
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
FullModelReactionDynamics::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2586
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
ke
int ke
Definition: CascadeWrapper.h:13
GeV
const double GeV
Definition: MathUtil.h:16
diffTwoXMLs.r2
r2
Definition: diffTwoXMLs.py:73
p1
double p1[4]
Definition: TauolaWrapper.h:89
EgHLTOffHistBins_cfi.et
et
Definition: EgHLTOffHistBins_cfi.py:8
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:204
FullModelReactionDynamics::Factorial
G4int Factorial(G4int n)
Definition: FullModelReactionDynamics.cc:2576
alignCSCRings.r
r
Definition: alignCSCRings.py:93
DDAxes::phi
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
diffTwoXMLs.r1
r1
Definition: diffTwoXMLs.py:53
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
p3
double p3[4]
Definition: TauolaWrapper.h:91
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
funct::arg
A arg
Definition: Factorize.h:31
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
mps_fire.result
result
Definition: mps_fire.py:311
pi
const Double_t pi
Definition: trackSplitPlot.h:36
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
beamvalidation.exit
def exit(msg="")
Definition: beamvalidation.py:52
createTree.pp
pp
Definition: createTree.py:17
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:56
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
ecalTB2006H4_GenSimDigiReco_cfg.G4cout
G4cout
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:285
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
A
cuy.ii
ii
Definition: cuy.py:589
amptDefault_cfi.bMin
bMin
Definition: amptDefault_cfi.py:19
weight
Definition: weight.py:1
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16