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

#include <FullModelReactionDynamics.h>

Public Member Functions

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

Private Member Functions

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

Detailed Description

Definition at line 23 of file FullModelReactionDynamics.h.

Constructor & Destructor Documentation

FullModelReactionDynamics::FullModelReactionDynamics ( )
inline

Definition at line 27 of file FullModelReactionDynamics.h.

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

Definition at line 29 of file FullModelReactionDynamics.h.

29 {}

Member Function Documentation

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

Definition at line 2998 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

3011  {
3012  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
3013  //
3014  // npnb is number of proton/neutron black track particles
3015  // ndta is the number of deuterons, tritons, and alphas produced
3016  // epnb is the kinetic energy available for proton/neutron black track particles
3017  // edta is the kinetic energy available for deuteron/triton/alpha particles
3018  //
3019 
3020  G4ParticleDefinition *aProton = G4Proton::Proton();
3021  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3022  G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3023  G4ParticleDefinition *aTriton = G4Triton::Triton();
3024  G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3025 
3026  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/MeV;
3027  G4double atomicWeight = targetNucleus.GetN_asInt();
3028  G4double atomicNumber = targetNucleus.GetZ_asInt();
3029 
3030  const G4double ika1 = 3.6;
3031  const G4double ika2 = 35.56;
3032  const G4double ika3 = 6.45;
3033  const G4double sp1 = 1.066;
3034 
3035  G4int i;
3036  G4double pp;
3037  // G4double totalQ = 0;
3038  G4double kinCreated = 0;
3039  G4double cfa = 0.025*((atomicWeight-1.0)/120.0) * std::exp(-(atomicWeight-1.0)/120.0);
3040  if( npnb > 0) // first add protons and neutrons
3041  {
3042  G4double backwardKinetic = 0.0;
3043  G4int local_npnb = npnb;
3044  for( i=0; i<npnb; ++i ) if( G4UniformRand() < sprob ) local_npnb--;
3045  G4double ekin = epnb/local_npnb;
3046 
3047  for( i=0; i<local_npnb; ++i )
3048  {
3049  G4ReactionProduct * p1 = new G4ReactionProduct();
3050  if( backwardKinetic > epnb )
3051  {
3052  delete p1;
3053  break;
3054  }
3055  G4double ran = G4UniformRand();
3056  G4double kinetic = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
3057  if( kinetic < 0.0 )kinetic = -0.010*std::log(ran);
3058  backwardKinetic += kinetic;
3059  if( backwardKinetic > epnb )
3060  kinetic = std::max( kineticMinimum, epnb-(backwardKinetic-kinetic) );
3061  if( G4UniformRand() > (1.0-atomicNumber/atomicWeight) )
3062  p1->SetDefinition( aProton );
3063  else
3064  p1->SetDefinition( aNeutron );
3065  vec.SetElement( vecLen, p1 );
3066  ++spall;
3067  G4double cost = G4UniformRand() * 2.0 - 1.0;
3068  G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
3069  G4double phi = twopi * G4UniformRand();
3070  vec[vecLen]->SetNewlyAdded( true );
3071  vec[vecLen]->SetKineticEnergy( kinetic*GeV );
3072  kinCreated+=kinetic;
3073  pp = vec[vecLen]->GetTotalMomentum()/MeV;
3074  vec[vecLen]->SetMomentum( pp*sint*std::sin(phi)*MeV,
3075  pp*sint*std::cos(phi)*MeV,
3076  pp*cost*MeV );
3077  vecLen++;
3078  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3079  }
3080  if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*GeV) )
3081  {
3082  G4double ekw = ekOriginal/GeV;
3083  G4int ika, kk = 0;
3084  if( ekw > 1.0 )ekw *= ekw;
3085  ekw = std::max( 0.1, ekw );
3086  ika = G4int(ika1*std::exp((atomicNumber*atomicNumber/atomicWeight-ika2)/ika3)/ekw);
3087  if( ika > 0 )
3088  {
3089  for( i=(vecLen-1); i>=0; --i )
3090  {
3091  if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
3092  {
3093  vec[i]->SetDefinitionAndUpdateE( aNeutron ); // modified 22-Oct-97
3094  if( ++kk > ika )break;
3095  }
3096  }
3097  }
3098  }
3099  }
3100  if( ndta > 0 ) // now, try to add deuterons, tritons and alphas
3101  {
3102  G4double backwardKinetic = 0.0;
3103  G4int local_ndta=ndta;
3104  for( i=0; i<ndta; ++i )if( G4UniformRand() < sprob )local_ndta--;
3105  G4double ekin = edta/local_ndta;
3106 
3107  for( i=0; i<local_ndta; ++i )
3108  {
3109  G4ReactionProduct *p2 = new G4ReactionProduct();
3110  if( backwardKinetic > edta )
3111  {
3112  delete p2;
3113  break;
3114  }
3115  G4double ran = G4UniformRand();
3116  G4double kinetic = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
3117  if( kinetic < 0.0 )kinetic = kineticFactor*std::log(ran);
3118  backwardKinetic += kinetic;
3119  if( backwardKinetic > edta )kinetic = edta-(backwardKinetic-kinetic);
3120  if( kinetic < 0.0 )kinetic = kineticMinimum;
3121  G4double cost = 2.0*G4UniformRand() - 1.0;
3122  G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
3123  G4double phi = twopi*G4UniformRand();
3124  ran = G4UniformRand();
3125  if( ran <= 0.60 )
3126  p2->SetDefinition( aDeuteron );
3127  else if( ran <= 0.90 )
3128  p2->SetDefinition( aTriton );
3129  else
3130  p2->SetDefinition( anAlpha );
3131  spall += p2->GetMass()/GeV * sp1;
3132  if( spall > atomicWeight )
3133  {
3134  delete p2;
3135  break;
3136  }
3137  vec.SetElement( vecLen, p2 );
3138  vec[vecLen]->SetNewlyAdded( true );
3139  vec[vecLen]->SetKineticEnergy( kinetic*GeV );
3140  kinCreated+=kinetic;
3141  pp = vec[vecLen]->GetTotalMomentum()/MeV;
3142  vec[vecLen++]->SetMomentum( pp*sint*std::sin(phi)*MeV,
3143  pp*sint*std::cos(phi)*MeV,
3144  pp*cost*MeV );
3145  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3146  }
3147  }
3148  // G4double delta = epnb+edta - kinCreated;
3149  }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:269
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
p2
double p2[4]
Definition: TauolaWrapper.h:90
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:83
p1
double p1[4]
Definition: TauolaWrapper.h:89
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
void FullModelReactionDynamics::Defs1 ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen 
)
private

Definition at line 2702 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

2708  {
2709  const G4double pjx = modifiedOriginal.GetMomentum().x()/MeV;
2710  const G4double pjy = modifiedOriginal.GetMomentum().y()/MeV;
2711  const G4double pjz = modifiedOriginal.GetMomentum().z()/MeV;
2712  const G4double p = modifiedOriginal.GetMomentum().mag()/MeV;
2713  if( pjx*pjx+pjy*pjy > 0.0 )
2714  {
2715  G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2716  cost = pjz/p;
2717  sint = 0.5 * ( std::sqrt(std::abs((1.0-cost)*(1.0+cost))) + std::sqrt(pjx*pjx+pjy*pjy)/p );
2718  if( pjy < 0.0 )
2719  ph = 3*halfpi;
2720  else
2721  ph = halfpi;
2722  if( std::abs( pjx ) > 0.001*MeV )ph = std::atan2(pjy,pjx);
2723  cosp = std::cos(ph);
2724  sinp = std::sin(ph);
2725  pix = currentParticle.GetMomentum().x()/MeV;
2726  piy = currentParticle.GetMomentum().y()/MeV;
2727  piz = currentParticle.GetMomentum().z()/MeV;
2728  currentParticle.SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
2729  cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
2730  -sint*pix*MeV + cost*piz*MeV );
2731  pix = targetParticle.GetMomentum().x()/MeV;
2732  piy = targetParticle.GetMomentum().y()/MeV;
2733  piz = targetParticle.GetMomentum().z()/MeV;
2734  targetParticle.SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
2735  cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
2736  -sint*pix*MeV + cost*piz*MeV );
2737  for( G4int i=0; i<vecLen; ++i )
2738  {
2739  pix = vec[i]->GetMomentum().x()/MeV;
2740  piy = vec[i]->GetMomentum().y()/MeV;
2741  piz = vec[i]->GetMomentum().z()/MeV;
2742  vec[i]->SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
2743  cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
2744  -sint*pix*MeV + cost*piz*MeV );
2745  }
2746  }
2747  else
2748  {
2749  if( pjz < 0.0 )
2750  {
2751  currentParticle.SetMomentum( -currentParticle.GetMomentum().z() );
2752  targetParticle.SetMomentum( -targetParticle.GetMomentum().z() );
2753  for( G4int i=0; i<vecLen; ++i )
2754  vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
2755  }
2756  }
2757  }
mps_fire.i
i
Definition: mps_fire.py:269
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
CombinatorialSeedGeneratorForBeamHalo_cfi.SetMomentum
SetMomentum
Definition: CombinatorialSeedGeneratorForBeamHalo_cfi.py:99
G4int FullModelReactionDynamics::Factorial ( G4int  n)

Definition at line 2693 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

2694  { // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
2695  G4int m = std::min(n,10);
2696  G4int result = 1;
2697  if( m <= 1 )return result;
2698  for( G4int i=2; i<=m; ++i )result *= i;
2699  return result;
2700  }
mps_fire.i
i
Definition: mps_fire.py:269
mps_fire.result
result
Definition: mps_fire.py:226
min
T min(T a, T b)
Definition: MathUtil.h:58
virtual G4double FullModelReactionDynamics::FindInelasticity ( )
inlinevirtual

Definition at line 31 of file FullModelReactionDynamics.h.

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

Definition at line 34 of file FullModelReactionDynamics.h.

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

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

Definition at line 2447 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

2452  {
2453  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2454  // derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
2455  // Returns the weight of the event
2456  //
2457  G4int i;
2458  const G4double expxu = 82.; // upper bound for arg. of exp
2459  const G4double expxl = -expxu; // lower bound for arg. of exp
2460  if( vecLen < 2 )
2461  {
2462  G4cerr << "*** Error in FullModelReactionDynamics::GenerateNBodyEvent" << G4endl;
2463  G4cerr << " number of particles < 2" << G4endl;
2464  G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
2465  return -1.0;
2466  }
2467  G4double mass[18]; // mass of each particle
2468  G4double energy[18]; // total energy of each particle
2469  G4double pcm[3][18]; // pcm is an array with 3 rows and vecLen columns
2470  G4double totalMass = 0.0;
2471  G4double extraMass = 0;
2472  G4double sm[18];
2473 
2474  for( i=0; i<vecLen; ++i )
2475  {
2476  mass[i] = vec[i]->GetMass()/GeV;
2477  if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
2478  vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
2479  pcm[0][i] = 0.0; // x-momentum of i-th particle
2480  pcm[1][i] = 0.0; // y-momentum of i-th particle
2481  pcm[2][i] = 0.0; // z-momentum of i-th particle
2482  energy[i] = mass[i]; // total energy of i-th particle
2483  totalMass += mass[i];
2484  sm[i] = totalMass;
2485  }
2486  G4double totalE = totalEnergy/GeV;
2487  if( totalMass > totalE )
2488  {
2489  return -1.0;
2490  }
2491  G4double kineticEnergy = totalE - totalMass;
2492  G4double emm[18];
2493  //G4double *emm = new G4double [vecLen];
2494  emm[0] = mass[0];
2495  emm[vecLen-1] = totalE;
2496  if( vecLen > 2 ) // the random numbers are sorted
2497  {
2498  G4double ran[18];
2499  for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
2500  for( i=0; i<vecLen-2; ++i )
2501  {
2502  for( G4int j=vecLen-2; j>i; --j )
2503  {
2504  if( ran[i] > ran[j] )
2505  {
2506  G4double temp = ran[i];
2507  ran[i] = ran[j];
2508  ran[j] = temp;
2509  }
2510  }
2511  }
2512  for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
2513  }
2514  // Weight is the sum of logarithms of terms instead of the product of terms
2515  G4bool lzero = true;
2516  G4double wtmax = 0.0;
2517  if( constantCrossSection ) // this is KGENEV=1 in PHASP
2518  {
2519  G4double emmax = kineticEnergy + mass[0];
2520  G4double emmin = 0.0;
2521  for( i=1; i<vecLen; ++i )
2522  {
2523  emmin += mass[i-1];
2524  emmax += mass[i];
2525  G4double wtfc = 0.0;
2526  if( emmax*emmax > 0.0 )
2527  {
2528  G4double arg = emmax*emmax
2529  + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
2530  - 2.0*(emmin*emmin+mass[i]*mass[i]);
2531  if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
2532  }
2533  if( wtfc == 0.0 )
2534  {
2535  lzero = false;
2536  break;
2537  }
2538  wtmax += std::log( wtfc );
2539  }
2540  if( lzero )
2541  wtmax = -wtmax;
2542  else
2543  wtmax = expxu;
2544  }
2545  else
2546  {
2547  // ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
2548  const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
2549  256.3704, 268.4705, 240.9780, 189.2637,
2550  132.1308, 83.0202, 47.4210, 24.8295,
2551  12.0006, 5.3858, 2.2560, 0.8859 };
2552  wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
2553  }
2554  lzero = true;
2555  G4double pd[50];
2556  //G4double *pd = new G4double [vecLen-1];
2557  for( i=0; i<vecLen-1; ++i )
2558  {
2559  pd[i] = 0.0;
2560  if( emm[i+1]*emm[i+1] > 0.0 )
2561  {
2562  G4double arg = emm[i+1]*emm[i+1]
2563  + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
2564  /(emm[i+1]*emm[i+1])
2565  - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
2566  if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
2567  }
2568  if( pd[i] <= 0.0 ) // changed from == on 02 April 98
2569  lzero = false;
2570  else
2571  wtmax += std::log( pd[i] );
2572  }
2573  G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
2574  if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
2575 
2576  G4double bang, cb, sb, s0, s1, s2, c, s, esys, a, b, gama, beta;
2577  pcm[0][0] = 0.0;
2578  pcm[1][0] = pd[0];
2579  pcm[2][0] = 0.0;
2580  for( i=1; i<vecLen; ++i )
2581  {
2582  pcm[0][i] = 0.0;
2583  pcm[1][i] = -pd[i-1];
2584  pcm[2][i] = 0.0;
2585  bang = twopi*G4UniformRand();
2586  cb = std::cos(bang);
2587  sb = std::sin(bang);
2588  c = 2.0*G4UniformRand() - 1.0;
2589  s = std::sqrt( std::fabs( 1.0-c*c ) );
2590  if( i < vecLen-1 )
2591  {
2592  esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
2593  beta = pd[i]/esys;
2594  gama = esys/emm[i];
2595  for( G4int j=0; j<=i; ++j )
2596  {
2597  s0 = pcm[0][j];
2598  s1 = pcm[1][j];
2599  s2 = pcm[2][j];
2600  energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
2601  a = s0*c - s1*s; // rotation
2602  pcm[1][j] = s0*s + s1*c;
2603  b = pcm[2][j];
2604  pcm[0][j] = a*cb - b*sb;
2605  pcm[2][j] = a*sb + b*cb;
2606  pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
2607  }
2608  }
2609  else
2610  {
2611  for( G4int j=0; j<=i; ++j )
2612  {
2613  s0 = pcm[0][j];
2614  s1 = pcm[1][j];
2615  s2 = pcm[2][j];
2616  energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
2617  a = s0*c - s1*s; // rotation
2618  pcm[1][j] = s0*s + s1*c;
2619  b = pcm[2][j];
2620  pcm[0][j] = a*cb - b*sb;
2621  pcm[2][j] = a*sb + b*cb;
2622  }
2623  }
2624  }
2625  for( i=0; i<vecLen; ++i )
2626  {
2627  vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
2628  vec[i]->SetTotalEnergy( energy[i]*GeV );
2629  }
2630  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2631  return weight;
2632  }
beta
const double beta
Definition: CosmicMuonParameters.h:69
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:269
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
weight
Definition: weight.py:1
funct::arg
A arg
Definition: Factorize.h:37
ecalTB2006H4_GenSimDigiReco_cfg.G4cerr
G4cerr
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:303
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
mps_merge.weight
weight
Definition: mps_merge.py:86
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
G4bool FullModelReactionDynamics::GenerateXandPt ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged,
G4bool  leadFlag,
G4ReactionProduct &  leadingStrangeParticle 
)

Definition at line 62 of file FullModelReactionDynamics.cc.

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

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

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

Definition at line 3151 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

3157  {
3158  const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/MeV;
3159  G4double testMomentum = currentParticle.GetMomentum().mag()/MeV;
3160  G4double pMass;
3161  if( testMomentum >= pOriginal )
3162  {
3163  pMass = currentParticle.GetMass()/MeV;
3164  currentParticle.SetTotalEnergy(
3165  std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
3166  currentParticle.SetMomentum(
3167  currentParticle.GetMomentum() * (pOriginal/testMomentum) );
3168  }
3169  testMomentum = targetParticle.GetMomentum().mag()/MeV;
3170  if( testMomentum >= pOriginal )
3171  {
3172  pMass = targetParticle.GetMass()/MeV;
3173  targetParticle.SetTotalEnergy(
3174  std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
3175  targetParticle.SetMomentum(
3176  targetParticle.GetMomentum() * (pOriginal/testMomentum) );
3177  }
3178  for( G4int i=0; i<vecLen; ++i )
3179  {
3180  testMomentum = vec[i]->GetMomentum().mag()/MeV;
3181  if( testMomentum >= pOriginal )
3182  {
3183  pMass = vec[i]->GetMass()/MeV;
3184  vec[i]->SetTotalEnergy(
3185  std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
3186  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
3187  }
3188  }
3189  }
mps_fire.i
i
Definition: mps_fire.py:269
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
G4double FullModelReactionDynamics::normal ( )
private

Definition at line 2635 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

2636  {
2637  G4double ran = -6.0;
2638  for( G4int i=0; i<12; ++i )ran += G4UniformRand();
2639  return ran;
2640  }
mps_fire.i
i
Definition: mps_fire.py:269
void FullModelReactionDynamics::NuclearReaction ( G4FastVector< G4ReactionProduct, 4 > &  vec,
G4int &  vecLen,
const G4HadProjectile *  originalIncident,
const G4Nucleus &  aNucleus,
const G4double  theAtomicMass,
const G4double *  massVec 
)

Definition at line 3580 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

3587  {
3588  // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
3589  //
3590  // Nuclear reaction kinematics at low energies
3591  //
3592  G4ParticleDefinition *aGamma = G4Gamma::Gamma();
3593  G4ParticleDefinition *aProton = G4Proton::Proton();
3594  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3595  G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3596  G4ParticleDefinition *aTriton = G4Triton::Triton();
3597  G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3598 
3599  const G4double aProtonMass = aProton->GetPDGMass()/MeV;
3600  const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
3601  const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
3602  const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
3603  const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;
3604 
3605  G4ReactionProduct currentParticle;
3606  currentParticle = *originalIncident;
3607  //
3608  // Set beam particle, take kinetic energy of current particle as the
3609  // fundamental quantity. Due to the difficult kinematic, all masses have to
3610  // be assigned the best measured values
3611  //
3612  G4double p = currentParticle.GetTotalMomentum();
3613  G4double pp = currentParticle.GetMomentum().mag();
3614  if( pp <= 0.001*MeV )
3615  {
3616  G4double phinve = twopi*G4UniformRand();
3617  G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
3618  currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3619  p*std::sin(rthnve)*std::sin(phinve),
3620  p*std::cos(rthnve) );
3621  }
3622  else
3623  currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
3624  //
3625  // calculate Q-value of reactions
3626  //
3627  G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
3628  G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
3629  G4double qv = currentKinetic + theAtomicMass + currentMass;
3630 
3631  G4double qval[9];
3632  qval[0] = qv - mass[0];
3633  qval[1] = qv - mass[1] - aNeutronMass;
3634  qval[2] = qv - mass[2] - aProtonMass;
3635  qval[3] = qv - mass[3] - aDeuteronMass;
3636  qval[4] = qv - mass[4] - aTritonMass;
3637  qval[5] = qv - mass[5] - anAlphaMass;
3638  qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3639  qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3640  qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3641 
3642  if( currentParticle.GetDefinition() == aNeutron )
3643  {
3644  const G4double A = targetNucleus.GetN_asInt(); // atomic weight
3645  if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
3646  qval[0] = 0.0;
3647  if( G4UniformRand() >= currentKinetic/7.9254*A )
3648  qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3649  }
3650  else
3651  qval[0] = 0.0;
3652 
3653  G4int i;
3654  qv = 0.0;
3655  for( i=0; i<9; ++i )
3656  {
3657  if( mass[i] < 500.0*MeV )qval[i] = 0.0;
3658  if( qval[i] < 0.0 )qval[i] = 0.0;
3659  qv += qval[i];
3660  }
3661  G4double qv1 = 0.0;
3662  G4double ran = G4UniformRand();
3663  G4int index;
3664  for( index=0; index<9; ++index )
3665  {
3666  if( qval[index] > 0.0 )
3667  {
3668  qv1 += qval[index]/qv;
3669  if( ran <= qv1 )break;
3670  }
3671  }
3672  if( index == 9 ) // loop continued to the end
3673  {
3674  throw G4HadronicException(__FILE__, __LINE__,
3675  "FullModelReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3676  }
3677  G4double ke = currentParticle.GetKineticEnergy()/GeV;
3678  G4int nt = 2;
3679  if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
3680 
3681  G4ReactionProduct **v = new G4ReactionProduct * [3];
3682  v[0] = new G4ReactionProduct;
3683  v[1] = new G4ReactionProduct;
3684  v[2] = new G4ReactionProduct;
3685 
3686  v[0]->SetMass( mass[index]*MeV );
3687  switch( index )
3688  {
3689  case 0:
3690  v[1]->SetDefinition( aGamma );
3691  v[2]->SetDefinition( aGamma );
3692  break;
3693  case 1:
3694  v[1]->SetDefinition( aNeutron );
3695  v[2]->SetDefinition( aGamma );
3696  break;
3697  case 2:
3698  v[1]->SetDefinition( aProton );
3699  v[2]->SetDefinition( aGamma );
3700  break;
3701  case 3:
3702  v[1]->SetDefinition( aDeuteron );
3703  v[2]->SetDefinition( aGamma );
3704  break;
3705  case 4:
3706  v[1]->SetDefinition( aTriton );
3707  v[2]->SetDefinition( aGamma );
3708  break;
3709  case 5:
3710  v[1]->SetDefinition( anAlpha );
3711  v[2]->SetDefinition( aGamma );
3712  break;
3713  case 6:
3714  v[1]->SetDefinition( aNeutron );
3715  v[2]->SetDefinition( aNeutron );
3716  break;
3717  case 7:
3718  v[1]->SetDefinition( aNeutron );
3719  v[2]->SetDefinition( aProton );
3720  break;
3721  case 8:
3722  v[1]->SetDefinition( aProton );
3723  v[2]->SetDefinition( aProton );
3724  break;
3725  }
3726  //
3727  // calculate centre of mass energy
3728  //
3729  G4ReactionProduct pseudo1;
3730  pseudo1.SetMass( theAtomicMass*MeV );
3731  pseudo1.SetTotalEnergy( theAtomicMass*MeV );
3732  G4ReactionProduct pseudo2 = currentParticle + pseudo1;
3733  pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
3734  //
3735  // use phase space routine in centre of mass system
3736  //
3737  G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> tempV;
3738  tempV.Initialize( nt );
3739  G4int tempLen = 0;
3740  tempV.SetElement( tempLen++, v[0] );
3741  tempV.SetElement( tempLen++, v[1] );
3742  if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
3743  G4bool constantCrossSection = true;
3744  GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
3745  v[0]->Lorentz( *v[0], pseudo2 );
3746  v[1]->Lorentz( *v[1], pseudo2 );
3747  if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
3748 
3749  G4bool particleIsDefined = false;
3750  if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
3751  {
3752  v[0]->SetDefinition( aProton );
3753  particleIsDefined = true;
3754  }
3755  else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
3756  {
3757  v[0]->SetDefinition( aNeutron );
3758  particleIsDefined = true;
3759  }
3760  else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
3761  {
3762  v[0]->SetDefinition( aDeuteron );
3763  particleIsDefined = true;
3764  }
3765  else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
3766  {
3767  v[0]->SetDefinition( aTriton );
3768  particleIsDefined = true;
3769  }
3770  else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
3771  {
3772  v[0]->SetDefinition( anAlpha );
3773  particleIsDefined = true;
3774  }
3775  currentParticle.SetKineticEnergy(
3776  std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
3777  p = currentParticle.GetTotalMomentum();
3778  pp = currentParticle.GetMomentum().mag();
3779  if( pp <= 0.001*MeV )
3780  {
3781  G4double phinve = twopi*G4UniformRand();
3782  G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
3783  currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3784  p*std::sin(rthnve)*std::sin(phinve),
3785  p*std::cos(rthnve) );
3786  }
3787  else
3788  currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
3789 
3790  if( particleIsDefined )
3791  {
3792  v[0]->SetKineticEnergy(
3793  std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
3794  p = v[0]->GetTotalMomentum();
3795  pp = v[0]->GetMomentum().mag();
3796  if( pp <= 0.001*MeV )
3797  {
3798  G4double phinve = twopi*G4UniformRand();
3799  G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
3800  v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3801  p*std::sin(rthnve)*std::sin(phinve),
3802  p*std::cos(rthnve) );
3803  }
3804  else
3805  v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
3806  }
3807  if( (v[1]->GetDefinition() == aDeuteron) ||
3808  (v[1]->GetDefinition() == aTriton) ||
3809  (v[1]->GetDefinition() == anAlpha) )
3810  v[1]->SetKineticEnergy(
3811  std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
3812  else
3813  v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
3814 
3815  p = v[1]->GetTotalMomentum();
3816  pp = v[1]->GetMomentum().mag();
3817  if( pp <= 0.001*MeV )
3818  {
3819  G4double phinve = twopi*G4UniformRand();
3820  G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
3821  v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3822  p*std::sin(rthnve)*std::sin(phinve),
3823  p*std::cos(rthnve) );
3824  }
3825  else
3826  v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
3827 
3828  if( nt == 3 )
3829  {
3830  if( (v[2]->GetDefinition() == aDeuteron) ||
3831  (v[2]->GetDefinition() == aTriton) ||
3832  (v[2]->GetDefinition() == anAlpha) )
3833  v[2]->SetKineticEnergy(
3834  std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
3835  else
3836  v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
3837 
3838  p = v[2]->GetTotalMomentum();
3839  pp = v[2]->GetMomentum().mag();
3840  if( pp <= 0.001*MeV )
3841  {
3842  G4double phinve = twopi*G4UniformRand();
3843  G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
3844  v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3845  p*std::sin(rthnve)*std::sin(phinve),
3846  p*std::cos(rthnve) );
3847  }
3848  else
3849  v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
3850  }
3851  G4int del;
3852  for(del=0; del<vecLen; del++) delete vec[del];
3853  vecLen = 0;
3854  if( particleIsDefined )
3855  {
3856  vec.SetElement( vecLen++, v[0] );
3857  }
3858  else
3859  {
3860  delete v[0];
3861  }
3862  vec.SetElement( vecLen++, v[1] );
3863  if( nt == 3 )
3864  {
3865  vec.SetElement( vecLen++, v[2] );
3866  }
3867  else
3868  {
3869  delete v[2];
3870  }
3871  delete [] v;
3872  return;
3873  }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:269
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MeV
const double MeV
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
FullModelReactionDynamics::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2447
nt
int nt
Definition: AMPTWrapper.h:32
Gamma
dbl * Gamma
Definition: mlp_gen.cc:38
patCaloMETCorrections_cff.A
A
Definition: patCaloMETCorrections_cff.py:43
G4int FullModelReactionDynamics::Poisson ( G4double  x)
private

Definition at line 2643 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

2644  {
2645  G4int iran;
2646  G4double ran;
2647 
2648  if( x > 9.9 ) // use normal distribution with sigma^2 = <x>
2649  iran = static_cast<G4int>(std::max( 0.0, x+normal()*std::sqrt(x) ) );
2650  else {
2651  G4int mm = G4int(5.0*x);
2652  if( mm <= 0 ) // for very small x try iran=1,2,3
2653  {
2654  G4double p1 = x*std::exp(-x);
2655  G4double p2 = x*p1/2.0;
2656  G4double p3 = x*p2/3.0;
2657  ran = G4UniformRand();
2658  if( ran < p3 )
2659  iran = 3;
2660  else if( ran < p2 ) // this is original Geisha, it should be ran < p2+p3
2661  iran = 2;
2662  else if( ran < p1 ) // should be ran < p1+p2+p3
2663  iran = 1;
2664  else
2665  iran = 0;
2666  }
2667  else
2668  {
2669  iran = 0;
2670  G4double r = std::exp(-x);
2671  ran = G4UniformRand();
2672  if( ran > r )
2673  {
2674  G4double rrr;
2675  G4double rr = r;
2676  for( G4int i=1; i<=mm; ++i )
2677  {
2678  iran++;
2679  if( i > 5 ) // Stirling's formula for large numbers
2680  rrr = std::exp(i*std::log(x)-(i+0.5)*std::log((G4double)i)+i-0.9189385);
2681  else
2682  rrr = std::pow(x,i)/Factorial(i);
2683  rr += r*rrr;
2684  if( ran <= rr )break;
2685  }
2686  }
2687  }
2688  }
2689  return iran;
2690  }
mps_fire.i
i
Definition: mps_fire.py:269
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
p2
double p2[4]
Definition: TauolaWrapper.h:90
p1
double p1[4]
Definition: TauolaWrapper.h:89
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
p3
double p3[4]
Definition: TauolaWrapper.h:91
void FullModelReactionDynamics::ProduceStrangeParticlePairs ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
const G4ReactionProduct &  modifiedOriginal,
const G4DynamicParticle *  originalTarget,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged 
)

Definition at line 3191 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

3200  {
3201  // derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
3202  //
3203  // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
3204  // K+ Y0, K0 Y+, K0 Y-
3205  // For antibaryon induced reactions half of the cross sections KB YB
3206  // pairs are produced. Charge is not conserved, no experimental data available
3207  // for exclusive reactions, therefore some average behaviour assumed.
3208  // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
3209  //
3210  if( vecLen == 0 )return;
3211  //
3212  // the following protects against annihilation processes
3213  //
3214  if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return;
3215 
3216  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
3217  const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
3218  G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
3219  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
3220  targetMass*targetMass +
3221  2.0*targetMass*etOriginal ); // GeV
3222  G4double currentMass = currentParticle.GetMass()/GeV;
3223  G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
3224  if( availableEnergy <= 1.0 )return;
3225 
3226  G4ParticleDefinition *aProton = G4Proton::Proton();
3227  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
3228  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3229  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
3230  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
3231  G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
3232  G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
3233  G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
3234  G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
3235  G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
3236  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
3237  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
3238  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
3239  G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
3240  G4ParticleDefinition *aLambda = G4Lambda::Lambda();
3241  G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
3242 
3243  const G4double protonMass = aProton->GetPDGMass()/GeV;
3244  const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
3245  //
3246  // determine the center of mass energy bin
3247  //
3248  const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
3249 
3250  G4int ibin, i3, i4;
3251  G4double avk, avy, avn, ran;
3252  G4int i = 1;
3253  while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
3254  if( i == 12 )
3255  ibin = 11;
3256  else
3257  ibin = i;
3258  //
3259  // the fortran code chooses a random replacement of produced kaons
3260  // but does not take into account charge conservation
3261  //
3262  if( vecLen == 1 ) // we know that vecLen > 0
3263  {
3264  i3 = 0;
3265  i4 = 1; // note that we will be adding a new secondary particle in this case only
3266  }
3267  else // otherwise 0 <= i3,i4 < vecLen
3268  {
3269  G4double ran = G4UniformRand();
3270  while( ran == 1.0 )ran = G4UniformRand();
3271  i4 = i3 = G4int( vecLen*ran );
3272  while( i3 == i4 )
3273  {
3274  ran = G4UniformRand();
3275  while( ran == 1.0 )ran = G4UniformRand();
3276  i4 = G4int( vecLen*ran );
3277  }
3278  }
3279  //
3280  // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
3281  //
3282  const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
3283  0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
3284  const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
3285  0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
3286  const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
3287  0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
3288 
3289  avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3290  /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
3291  avk = std::exp(avk);
3292 
3293  avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3294  /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
3295  avy = std::exp(avy);
3296 
3297  avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3298  /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
3299  avn = std::exp(avn);
3300 
3301  if( avk+avy+avn <= 0.0 )return;
3302 
3303  if( currentMass < protonMass )avy /= 2.0;
3304  if( targetMass < protonMass )avy = 0.0;
3305  avy += avk+avn;
3306  avk += avn;
3307  ran = G4UniformRand();
3308  if( ran < avn )
3309  {
3310  if( availableEnergy < 2.0 )return;
3311  if( vecLen == 1 ) // add a new secondary
3312  {
3313  G4ReactionProduct *p1 = new G4ReactionProduct;
3314  if( G4UniformRand() < 0.5 )
3315  {
3316  vec[0]->SetDefinition( aNeutron );
3317  p1->SetDefinition( anAntiNeutron );
3318  (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
3319  vec[0]->SetMayBeKilled(false);
3320  p1->SetMayBeKilled(false);
3321  }
3322  else
3323  {
3324  vec[0]->SetDefinition( aProton );
3325  p1->SetDefinition( anAntiProton );
3326  (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
3327  vec[0]->SetMayBeKilled(false);
3328  p1->SetMayBeKilled(false);
3329  }
3330  vec.SetElement( vecLen++, p1 );
3331  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3332  }
3333  else
3334  { // replace two secondaries
3335  if( G4UniformRand() < 0.5 )
3336  {
3337  vec[i3]->SetDefinition( aNeutron );
3338  vec[i4]->SetDefinition( anAntiNeutron );
3339  vec[i3]->SetMayBeKilled(false);
3340  vec[i4]->SetMayBeKilled(false);
3341  }
3342  else
3343  {
3344  vec[i3]->SetDefinition( aProton );
3345  vec[i4]->SetDefinition( anAntiProton );
3346  vec[i3]->SetMayBeKilled(false);
3347  vec[i4]->SetMayBeKilled(false);
3348  }
3349  }
3350  }
3351  else if( ran < avk )
3352  {
3353  if( availableEnergy < 1.0 )return;
3354 
3355  const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
3356  0.6875, 0.7500, 0.8750, 1.000 };
3357  const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
3358  const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
3359  ran = G4UniformRand();
3360  i = 0;
3361  while( (i<9) && (ran>=kkb[i]) )++i;
3362  if( i == 9 )return;
3363  //
3364  // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
3365  // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
3366  //
3367  switch( ipakkb1[i] )
3368  {
3369  case 10:
3370  vec[i3]->SetDefinition( aKaonPlus );
3371  vec[i3]->SetMayBeKilled(false);
3372  break;
3373  case 11:
3374  vec[i3]->SetDefinition( aKaonZS );
3375  vec[i3]->SetMayBeKilled(false);
3376  break;
3377  case 12:
3378  vec[i3]->SetDefinition( aKaonZL );
3379  vec[i3]->SetMayBeKilled(false);
3380  break;
3381  }
3382  if( vecLen == 1 ) // add a secondary
3383  {
3384  G4ReactionProduct *p1 = new G4ReactionProduct;
3385  switch( ipakkb2[i] )
3386  {
3387  case 11:
3388  p1->SetDefinition( aKaonZS );
3389  p1->SetMayBeKilled(false);
3390  break;
3391  case 12:
3392  p1->SetDefinition( aKaonZL );
3393  p1->SetMayBeKilled(false);
3394  break;
3395  case 13:
3396  p1->SetDefinition( aKaonMinus );
3397  p1->SetMayBeKilled(false);
3398  break;
3399  }
3400  (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
3401  vec.SetElement( vecLen++, p1 );
3402  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3403  }
3404  else // replace
3405  {
3406  switch( ipakkb2[i] )
3407  {
3408  case 11:
3409  vec[i4]->SetDefinition( aKaonZS );
3410  vec[i4]->SetMayBeKilled(false);
3411  break;
3412  case 12:
3413  vec[i4]->SetDefinition( aKaonZL );
3414  vec[i4]->SetMayBeKilled(false);
3415  break;
3416  case 13:
3417  vec[i4]->SetDefinition( aKaonMinus );
3418  vec[i4]->SetMayBeKilled(false);
3419  break;
3420  }
3421  }
3422  }
3423  else if( ran < avy )
3424  {
3425  if( availableEnergy < 1.6 )return;
3426 
3427  const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
3428  0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
3429  const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
3430  const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
3431  const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
3432  const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
3433  ran = G4UniformRand();
3434  i = 0;
3435  while( (i<12) && (ran>ky[i]) )++i;
3436  if( i == 12 )return;
3437  if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
3438  {
3439  // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
3440  // 0 + 0 0 0 0 + + + 0 + 0
3441  //
3442  // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
3443  // 0 + 0 0 0 0 - + - 0 - 0
3444  switch( ipaky1[i] )
3445  {
3446  case 18:
3447  targetParticle.SetDefinition( aLambda );
3448  break;
3449  case 20:
3450  targetParticle.SetDefinition( aSigmaPlus );
3451  break;
3452  case 21:
3453  targetParticle.SetDefinition( aSigmaZero );
3454  break;
3455  case 22:
3456  targetParticle.SetDefinition( aSigmaMinus );
3457  break;
3458  }
3459  targetHasChanged = true;
3460  switch( ipaky2[i] )
3461  {
3462  case 10:
3463  vec[i3]->SetDefinition( aKaonPlus );
3464  vec[i3]->SetMayBeKilled(false);
3465  break;
3466  case 11:
3467  vec[i3]->SetDefinition( aKaonZS );
3468  vec[i3]->SetMayBeKilled(false);
3469  break;
3470  case 12:
3471  vec[i3]->SetDefinition( aKaonZL );
3472  vec[i3]->SetMayBeKilled(false);
3473  break;
3474  }
3475  }
3476  else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
3477  {
3478  if( (currentParticle.GetDefinition() == anAntiProton) ||
3479  (currentParticle.GetDefinition() == anAntiNeutron) ||
3480  (currentParticle.GetDefinition() == anAntiLambda) ||
3481  (currentMass > sigmaMinusMass) )
3482  {
3483  switch( ipakyb1[i] )
3484  {
3485  case 19:
3486  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
3487  break;
3488  case 23:
3489  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
3490  break;
3491  case 24:
3492  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
3493  break;
3494  case 25:
3495  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
3496  break;
3497  }
3498  incidentHasChanged = true;
3499  switch( ipakyb2[i] )
3500  {
3501  case 11:
3502  vec[i3]->SetDefinition( aKaonZS );
3503  vec[i3]->SetMayBeKilled(false);
3504  break;
3505  case 12:
3506  vec[i3]->SetDefinition( aKaonZL );
3507  vec[i3]->SetMayBeKilled(false);
3508  break;
3509  case 13:
3510  vec[i3]->SetDefinition( aKaonMinus );
3511  vec[i3]->SetMayBeKilled(false);
3512  break;
3513  }
3514  }
3515  else
3516  {
3517  switch( ipaky1[i] )
3518  {
3519  case 18:
3520  currentParticle.SetDefinitionAndUpdateE( aLambda );
3521  break;
3522  case 20:
3523  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
3524  break;
3525  case 21:
3526  currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
3527  break;
3528  case 22:
3529  currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
3530  break;
3531  }
3532  incidentHasChanged = true;
3533  switch( ipaky2[i] )
3534  {
3535  case 10:
3536  vec[i3]->SetDefinition( aKaonPlus );
3537  vec[i3]->SetMayBeKilled(false);
3538  break;
3539  case 11:
3540  vec[i3]->SetDefinition( aKaonZS );
3541  vec[i3]->SetMayBeKilled(false);
3542  break;
3543  case 12:
3544  vec[i3]->SetDefinition( aKaonZL );
3545  vec[i3]->SetMayBeKilled(false);
3546  break;
3547  }
3548  }
3549  }
3550  }
3551  else return;
3552  //
3553  // check the available energy
3554  // if there is not enough energy for kkb/ky pair production
3555  // then reduce the number of secondary particles
3556  // NOTE:
3557  // the number of secondaries may have been changed
3558  // the incident and/or target particles may have changed
3559  // charge conservation is ignored (as well as strangness conservation)
3560  //
3561  currentMass = currentParticle.GetMass()/GeV;
3562  targetMass = targetParticle.GetMass()/GeV;
3563 
3564  G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
3565  for( i=0; i<vecLen; ++i )
3566  {
3567  energyCheck -= vec[i]->GetMass()/GeV;
3568  if( energyCheck < 0.0 ) // chop off the secondary List
3569  {
3570  vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
3571  G4int j;
3572  for(j=i; j<vecLen; j++) delete vec[j];
3573  break;
3574  }
3575  }
3576  return;
3577  }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:269
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
p1
double p1[4]
Definition: TauolaWrapper.h:89
BPhysicsValidation_cfi.Lambda
Lambda
Definition: BPhysicsValidation_cfi.py:63
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
phase2TrackerDigitizer_cfi.SigmaZero
SigmaZero
Definition: phase2TrackerDigitizer_cfi.py:41
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 2759 of file FullModelReactionDynamics.cc.

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

Referenced by FindTimeDelay().

2769  {
2770  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2771  //
2772  // Rotate in direction of z-axis, this does disturb in some way our
2773  // inclusive distributions, but it is necessary for momentum conservation
2774  //
2775  const G4double atomicWeight = targetNucleus.GetN_asInt();
2776  const G4double logWeight = std::log(atomicWeight);
2777 
2778  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2779  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2780  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2781 
2782  G4int i;
2783  G4ThreeVector pseudoParticle[4];
2784  for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
2785  pseudoParticle[0] = currentParticle.GetMomentum()
2786  + targetParticle.GetMomentum();
2787  for( i=0; i<vecLen; ++i )
2788  pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
2789  //
2790  // Some smearing in transverse direction from Fermi motion
2791  //
2792  G4float pp, pp1;
2793  G4double alekw, p, rthnve, phinve;
2794  G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
2795 
2796  r1 = twopi*G4UniformRand();
2797  r2 = G4UniformRand();
2798  a1 = std::sqrt(-2.0*std::log(r2));
2799  ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*GeV;
2800  ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*GeV;
2801  G4ThreeVector fermi(ran1, ran2, 0);
2802 
2803  pseudoParticle[0] = pseudoParticle[0]+fermi; // all particles + fermi
2804  pseudoParticle[2] = temp; // original in cms system
2805  pseudoParticle[3] = pseudoParticle[0];
2806 
2807  pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
2808  G4double rotation = 2.*pi*G4UniformRand();
2809  pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
2810  pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
2811  for(G4int ii=1; ii<=3; ii++)
2812  {
2813  p = pseudoParticle[ii].mag();
2814  if( p == 0.0 )
2815  pseudoParticle[ii]= G4ThreeVector( 0.0, 0.0, 0.0 );
2816  else
2817  pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
2818  }
2819 
2820  pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
2821  pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
2822  pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
2823  currentParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
2824 
2825  pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
2826  pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
2827  pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
2828  targetParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
2829 
2830  for( i=0; i<vecLen; ++i )
2831  {
2832  pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
2833  pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
2834  pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
2835  vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
2836  }
2837  //
2838  // Rotate in direction of primary particle, subtract binding energies
2839  // and make some further corrections if required
2840  //
2841  Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2842  G4double ekin;
2843  G4double dekin = 0.0;
2844  G4double ek1 = 0.0;
2845  G4int npions = 0;
2846  if( atomicWeight >= 1.5 ) // self-absorption in heavy molecules
2847  {
2848  // corrections for single particle spectra (shower particles)
2849  //
2850  const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
2851  const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
2852  alekw = std::log( originalIncident->GetKineticEnergy()/GeV );
2853  exh = 1.0;
2854  if( alekw > alem[0] ) // get energy bin
2855  {
2856  exh = val0[6];
2857  for( G4int j=1; j<7; ++j )
2858  {
2859  if( alekw < alem[j] ) // use linear interpolation/extrapolation
2860  {
2861  G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
2862  exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
2863  break;
2864  }
2865  }
2866  exh = 1.0 - exh;
2867  }
2868  const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
2869  ekin = currentParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
2870  ekin = std::max( 1.0e-6, ekin );
2871  xxh = 1.0;
2872  dekin += ekin*(1.0-xxh);
2873  ekin *= xxh;
2874  currentParticle.SetKineticEnergy( ekin*GeV );
2875  pp = currentParticle.GetTotalMomentum()/MeV;
2876  pp1 = currentParticle.GetMomentum().mag()/MeV;
2877  if( pp1 < 0.001*MeV )
2878  {
2879  rthnve = pi*G4UniformRand();
2880  phinve = twopi*G4UniformRand();
2881  currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV,
2882  pp*std::sin(rthnve)*std::sin(phinve)*MeV,
2883  pp*std::cos(rthnve)*MeV );
2884  }
2885  else
2886  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
2887  ekin = targetParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
2888  ekin = std::max( 1.0e-6, ekin );
2889  xxh = 1.0;
2890  if( ( (modifiedOriginal.GetDefinition() == aPiPlus) ||
2891  (modifiedOriginal.GetDefinition() == aPiMinus) ) &&
2892  (targetParticle.GetDefinition() == aPiZero) &&
2893  (G4UniformRand() < logWeight) )xxh = exh;
2894  dekin += ekin*(1.0-xxh);
2895  ekin *= xxh;
2896  if( (targetParticle.GetDefinition() == aPiPlus) ||
2897  (targetParticle.GetDefinition() == aPiZero) ||
2898  (targetParticle.GetDefinition() == aPiMinus) )
2899  {
2900  ++npions;
2901  ek1 += ekin;
2902  }
2903  targetParticle.SetKineticEnergy( ekin*GeV );
2904  pp = targetParticle.GetTotalMomentum()/MeV;
2905  pp1 = targetParticle.GetMomentum().mag()/MeV;
2906  if( pp1 < 0.001*MeV )
2907  {
2908  rthnve = pi*G4UniformRand();
2909  phinve = twopi*G4UniformRand();
2910  targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV,
2911  pp*std::sin(rthnve)*std::sin(phinve)*MeV,
2912  pp*std::cos(rthnve)*MeV );
2913  }
2914  else
2915  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
2916  for( i=0; i<vecLen; ++i )
2917  {
2918  ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
2919  ekin = std::max( 1.0e-6, ekin );
2920  xxh = 1.0;
2921  if( ( (modifiedOriginal.GetDefinition() == aPiPlus) ||
2922  (modifiedOriginal.GetDefinition() == aPiMinus) ) &&
2923  (vec[i]->GetDefinition() == aPiZero) &&
2924  (G4UniformRand() < logWeight) )xxh = exh;
2925  dekin += ekin*(1.0-xxh);
2926  ekin *= xxh;
2927  if( (vec[i]->GetDefinition() == aPiPlus) ||
2928  (vec[i]->GetDefinition() == aPiZero) ||
2929  (vec[i]->GetDefinition() == aPiMinus) )
2930  {
2931  ++npions;
2932  ek1 += ekin;
2933  }
2934  vec[i]->SetKineticEnergy( ekin*GeV );
2935  pp = vec[i]->GetTotalMomentum()/MeV;
2936  pp1 = vec[i]->GetMomentum().mag()/MeV;
2937  if( pp1 < 0.001*MeV )
2938  {
2939  rthnve = pi*G4UniformRand();
2940  phinve = twopi*G4UniformRand();
2941  vec[i]->SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV,
2942  pp*std::sin(rthnve)*std::sin(phinve)*MeV,
2943  pp*std::cos(rthnve)*MeV );
2944  }
2945  else
2946  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
2947  }
2948  }
2949  if( (ek1 != 0.0) && (npions > 0) )
2950  {
2951  dekin = 1.0 + dekin/ek1;
2952  //
2953  // first do the incident particle
2954  //
2955  if( (currentParticle.GetDefinition() == aPiPlus) ||
2956  (currentParticle.GetDefinition() == aPiZero) ||
2957  (currentParticle.GetDefinition() == aPiMinus) )
2958  {
2959  currentParticle.SetKineticEnergy(
2960  std::max( 0.001*MeV, dekin*currentParticle.GetKineticEnergy() ) );
2961  pp = currentParticle.GetTotalMomentum()/MeV;
2962  pp1 = currentParticle.GetMomentum().mag()/MeV;
2963  if( pp1 < 0.001 )
2964  {
2965  rthnve = pi*G4UniformRand();
2966  phinve = twopi*G4UniformRand();
2967  currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV,
2968  pp*std::sin(rthnve)*std::sin(phinve)*MeV,
2969  pp*std::cos(rthnve)*MeV );
2970  }
2971  else
2972  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
2973  }
2974  for( i=0; i<vecLen; ++i )
2975  {
2976  if( (vec[i]->GetDefinition() == aPiPlus) ||
2977  (vec[i]->GetDefinition() == aPiZero) ||
2978  (vec[i]->GetDefinition() == aPiMinus) )
2979  {
2980  vec[i]->SetKineticEnergy( std::max( 0.001*MeV, dekin*vec[i]->GetKineticEnergy() ) );
2981  pp = vec[i]->GetTotalMomentum()/MeV;
2982  pp1 = vec[i]->GetMomentum().mag()/MeV;
2983  if( pp1 < 0.001 )
2984  {
2985  rthnve = pi*G4UniformRand();
2986  phinve = twopi*G4UniformRand();
2987  vec[i]->SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV,
2988  pp*std::sin(rthnve)*std::sin(phinve)*MeV,
2989  pp*std::cos(rthnve)*MeV );
2990  }
2991  else
2992  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
2993  }
2994  }
2995  }
2996  }
GeV
const double GeV
Definition: MathUtil.h:16
FullModelReactionDynamics::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2702
mps_fire.i
i
Definition: mps_fire.py:269
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
pi
const Double_t pi
Definition: trackSplitPlot.h:36
fermi
const double fermi
Definition: MathUtil.h:17
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
cuy.ii
ii
Definition: cuy.py:588
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
void FullModelReactionDynamics::SuppressChargedPions ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged 
)

Definition at line 1290 of file FullModelReactionDynamics.cc.

References GeV, and mps_fire::i.

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

1299  {
1300  // this code was originally in the fortran code TWOCLU
1301  //
1302  // suppress charged pions, for various reasons
1303  //
1304  const G4double atomicWeight = targetNucleus.GetN_asInt();
1305  const G4double atomicNumber = targetNucleus.GetZ_asInt();
1306  const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
1307 
1308  // G4ParticleDefinition *aGamma = G4Gamma::Gamma();
1309  G4ParticleDefinition *aProton = G4Proton::Proton();
1310  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
1311  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1312  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
1313  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1314  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1315 
1316  const G4bool antiTest =
1317  modifiedOriginal.GetDefinition() != anAntiProton &&
1318  modifiedOriginal.GetDefinition() != anAntiNeutron;
1319  if( antiTest && (
1320  // currentParticle.GetDefinition() == aGamma ||
1321  currentParticle.GetDefinition() == aPiPlus ||
1322  currentParticle.GetDefinition() == aPiMinus ) &&
1323  ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
1324  ( G4UniformRand() <= atomicWeight/300.0 ) )
1325  {
1326  if( G4UniformRand() > atomicNumber/atomicWeight )
1327  currentParticle.SetDefinitionAndUpdateE( aNeutron );
1328  else
1329  currentParticle.SetDefinitionAndUpdateE( aProton );
1330  incidentHasChanged = true;
1331  }
1332  for( G4int i=0; i<vecLen; ++i )
1333  {
1334  if( antiTest && (
1335  // vec[i]->GetDefinition() == aGamma ||
1336  vec[i]->GetDefinition() == aPiPlus ||
1337  vec[i]->GetDefinition() == aPiMinus ) &&
1338  ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
1339  ( G4UniformRand() <= atomicWeight/300.0 ) )
1340  {
1341  if( G4UniformRand() > atomicNumber/atomicWeight )
1342  vec[i]->SetDefinitionAndUpdateE( aNeutron );
1343  else
1344  vec[i]->SetDefinitionAndUpdateE( aProton );
1345  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1346  }
1347  }
1348  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1349  }
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:269
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 2167 of file FullModelReactionDynamics.cc.

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

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

2176  {
2177  // G4cout<<"TwoBody called"<<G4endl;
2178  //
2179  // derived from original FORTRAN code TWOB by H. Fesefeldt (15-Sep-1987)
2180  //
2181  // Generation of momenta for elastic and quasi-elastic 2 body reactions
2182  //
2183  // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
2184  // The b values are parametrizations from experimental data.
2185  // Not available values are taken from those of similar reactions.
2186  //
2187  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2188  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2189  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2190  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
2191  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
2192  G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
2193  G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
2194 
2195  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2196  static const G4double expxu = 82.; // upper bound for arg. of exp
2197  static const G4double expxl = -expxu; // lower bound for arg. of exp
2198 
2199  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
2200  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
2201  G4double currentMass = currentParticle.GetMass()/GeV;
2202  G4double targetMass = targetParticle.GetMass()/GeV;
2203  const G4double atomicWeight = targetNucleus.GetN_asInt();
2204  // G4cout<<"Atomic weight is found to be: "<<atomicWeight<<G4endl;
2205  G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
2206  G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;
2207 
2208  G4double cmEnergy = std::sqrt( currentMass*currentMass +
2209  targetMass*targetMass +
2210  2.0*targetMass*etCurrent ); // in GeV
2211 
2212  //if( (pOriginal < 0.1) ||
2213  // (centerofmassEnergy < 0.01) ) // 2-body scattering not possible
2214  // Continue with original particle, but spend the nuclear evaporation energy
2215  // targetParticle.SetMass( 0.0 ); // flag that the target doesn't exist
2216  //else // Two-body scattering is possible
2217 
2218  if( (pCurrent < 0.1) || (cmEnergy < 0.01) ) // 2-body scattering not possible
2219  {
2220  targetParticle.SetMass( 0.0 ); // flag that the target particle doesn't exist
2221  }
2222  else
2223  {
2224 // moved this if-block to a later stage, i.e. to the assignment of the scattering angle
2225 // @@@@@ double-check.
2226 // if( targetParticle.GetDefinition() == aKaonMinus ||
2227 // targetParticle.GetDefinition() == aKaonZeroL ||
2228 // targetParticle.GetDefinition() == aKaonZeroS ||
2229 // targetParticle.GetDefinition() == aKaonPlus ||
2230 // targetParticle.GetDefinition() == aPiMinus ||
2231 // targetParticle.GetDefinition() == aPiZero ||
2232 // targetParticle.GetDefinition() == aPiPlus )
2233 // {
2234 // if( G4UniformRand() < 0.5 )
2235 // targetParticle.SetDefinitionAndUpdateE( aNeutron );
2236 // else
2237 // targetParticle.SetDefinitionAndUpdateE( aProton );
2238 // targetHasChanged = true;
2239 // targetMass = targetParticle.GetMass()/GeV;
2240 // }
2241  //
2242  // Set masses and momenta for final state particles
2243  //
2244  /*
2245  G4cout<<"Check 0"<<G4endl;
2246  G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()/GeV<<G4endl;
2247  G4cout<<"target mass: "<<targetParticle.GetMass()/GeV<<G4endl;
2248  G4cout<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
2249  G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()/GeV<<G4endl;
2250  G4cout<<"current mass: "<<currentParticle.GetMass()/GeV<<G4endl;
2251  G4cout<<currentParticle.GetDefinition()->GetParticleName()<<G4endl;
2252  */
2253  G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
2254  pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
2255  // G4cout << "pf: " << pf<< G4endl;
2256 
2257  if( pf <= 0.)// 0.001 )
2258  {
2259  for(G4int i=0; i<vecLen; i++) delete vec[i];
2260  vecLen = 0;
2261  throw G4HadronicException(__FILE__, __LINE__, "FullModelReactionDynamics::TwoBody: pf is too small ");
2262  }
2263 
2264  pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
2265  //
2266  // Set beam and target in centre of mass system
2267  //
2268  G4ReactionProduct pseudoParticle[3];
2269  pseudoParticle[0].SetMass( currentMass*GeV );
2270  pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
2271  pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
2272 
2273  pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
2274  pseudoParticle[1].SetMass( targetMass*GeV );
2275  pseudoParticle[1].SetKineticEnergy( 0.0 );
2276  //
2277  // Transform into centre of mass system
2278  //
2279  pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
2280  pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
2281  pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
2282  //
2283  // Set final state masses and energies in centre of mass system
2284  //
2285  currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
2286  targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
2287  //
2288  // Set |t| and |tmin|
2289  //
2290  const G4double cb = 0.01;
2291  const G4double b1 = 4.225;
2292  const G4double b2 = 1.795;
2293  //
2294  // Calculate slope b for elastic scattering on proton/neutron
2295  //
2296  G4double b = std::max( cb, b1+b2*std::log(pOriginal) );
2297  // G4double b = std::max( cb, b1+b2*std::log(ptemp) );
2298  G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
2299 
2300  G4double exindt = -1.0;
2301  exindt += std::exp(std::max(-btrang,expxl));
2302  //
2303  // Calculate sqr(std::sin(teta/2.) and std::cos(teta), set azimuth angle phi
2304  //
2305  G4double ctet = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) / btrang;
2306  if( std::fabs(ctet) > 1.0 )ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
2307  G4double stet = std::sqrt( (1.0-ctet)*(1.0+ctet) );
2308  G4double phi = twopi * G4UniformRand();
2309  //
2310  // Calculate final state momenta in centre of mass system
2311  //
2312  if(
2313  targetParticle.GetDefinition() == aKaonMinus ||
2314  targetParticle.GetDefinition() == aKaonZeroL ||
2315  targetParticle.GetDefinition() == aKaonZeroS ||
2316  targetParticle.GetDefinition() == aKaonPlus ||
2317  targetParticle.GetDefinition() == aPiMinus ||
2318  targetParticle.GetDefinition() == aPiZero ||
2319  targetParticle.GetDefinition() == aPiPlus )
2320  {
2321  currentParticle.SetMomentum( -pf*stet*std::sin(phi)*GeV,
2322  -pf*stet*std::cos(phi)*GeV,
2323  -pf*ctet*GeV );
2324  }
2325  else
2326  {
2327  currentParticle.SetMomentum( pf*stet*std::sin(phi)*GeV,
2328  pf*stet*std::cos(phi)*GeV,
2329  pf*ctet*GeV );
2330  }
2331  targetParticle.SetMomentum( currentParticle.GetMomentum() * (-1.0) );
2332  //
2333  // Transform into lab system
2334  //
2335  currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
2336  targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
2337 
2338  /*
2339  G4cout<<"Check 1"<<G4endl;
2340  G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
2341  G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
2342  G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
2343  G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
2344  */
2345 
2346  Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2347 
2348  G4double pp, pp1, ekin;
2349  if( atomicWeight >= 1.5 )
2350  {
2351  const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
2352  pp1 = currentParticle.GetMomentum().mag()/MeV;
2353  if( pp1 >= 1.0 )
2354  {
2355  ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
2356  ekin = std::max( 0.0001*GeV, ekin );
2357  currentParticle.SetKineticEnergy( ekin*MeV );
2358  pp = currentParticle.GetTotalMomentum()/MeV;
2359  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
2360  }
2361  pp1 = targetParticle.GetMomentum().mag()/MeV;
2362  if( pp1 >= 1.0 )
2363  {
2364  ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
2365  ekin = std::max( 0.0001*GeV, ekin );
2366  targetParticle.SetKineticEnergy( ekin*MeV );
2367  pp = targetParticle.GetTotalMomentum()/MeV;
2368  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
2369  }
2370  }
2371  }
2372  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2373  if( atomicWeight >= 1.5 )
2374  {
2375  // Add black track particles
2376  // the procedure is somewhat different than in TwoCluster and GenerateXandPt.
2377  // The reason is that we have to also simulate the nuclear reactions
2378  // at low energies like a(h,p)b, a(h,p p)b, a(h,n)b etc.
2379  //
2380  // npnb is number of proton/neutron black track particles
2381  // ndta is the number of deuterons, tritons, and alphas produced
2382  // epnb is the kinetic energy available for proton/neutron black track particles
2383  // edta is the kinetic energy available for deuteron/triton/alpha particles
2384  //
2385  G4double epnb, edta;
2386  G4int npnb=0, ndta=0;
2387 
2388  epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
2389  edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
2390  const G4double pnCutOff = 0.0001; // GeV
2391  const G4double dtaCutOff = 0.0001; // GeV
2392  const G4double kineticMinimum = 0.0001;
2393  const G4double kineticFactor = -0.010;
2394  G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
2395  if( epnb >= pnCutOff )
2396  {
2397  npnb = Poisson( epnb/0.02 );
2398  /*
2399  G4cout<<"A couple of Poisson numbers:"<<G4endl;
2400  for (int n=0;n!=10;n++) G4cout<<Poisson(epnb/0.02)<<", ";
2401  G4cout<<G4endl;
2402  */
2403  if( npnb > atomicWeight )npnb = G4int(atomicWeight);
2404  if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
2405  npnb = std::min( npnb, 127-vecLen );
2406  }
2407  if( edta >= dtaCutOff )
2408  {
2409  ndta = G4int(2.0 * std::log(atomicWeight));
2410  ndta = std::min( ndta, 127-vecLen );
2411  }
2412  G4double spall = 0.0;
2413  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2414 
2415  /*
2416  G4cout<<"Check 2"<<G4endl;
2417  G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
2418  G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
2419  G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
2420  G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
2421 
2422  G4cout<<"------------------------------------------------------------------------"<<G4endl;
2423  G4cout<<"Atomic weight: "<<atomicWeight<<G4endl;
2424  G4cout<<"number of proton/neutron black track particles: "<<npnb<<G4endl;
2425  G4cout<<"number of deuterons, tritons, and alphas produced: "<<ndta <<G4endl;
2426  G4cout<<"kinetic energy available for proton/neutron black track particles: "<<epnb/GeV<<" GeV" <<G4endl;
2427  G4cout<<"kinetic energy available for deuteron/triton/alpha particles: "<<edta/GeV <<" GeV"<<G4endl;
2428  G4cout<<"------------------------------------------------------------------------"<<G4endl;
2429  */
2430 
2431  AddBlackTrackParticles( epnb, npnb, edta, ndta, sprob, kineticMinimum, kineticFactor,
2432  modifiedOriginal, spall, targetNucleus,
2433  vec, vecLen );
2434 
2435  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2436  }
2437  //
2438  // calculate time delay for nuclear reactions
2439  //
2440  if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
2441  currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
2442  else
2443  currentParticle.SetTOF( 1.0 );
2444  return;
2445  }
GeV
const double GeV
Definition: MathUtil.h:16
FullModelReactionDynamics::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2702
mps_fire.i
i
Definition: mps_fire.py:269
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
MeV
const double MeV
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
FullModelReactionDynamics::AddBlackTrackParticles
void AddBlackTrackParticles(const G4double epnb, const G4int npnb, const G4double edta, const G4int ndta, const G4double sprob, const G4double kineticMinimum, const G4double kineticFactor, const G4ReactionProduct &modifiedOriginal, G4double spall, const G4Nucleus &aNucleus, G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &vec, G4int &vecLen)
Definition: FullModelReactionDynamics.cc:2998
b
double b
Definition: hdecay.h:120
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
G4bool FullModelReactionDynamics::TwoCluster ( G4FastVector< G4ReactionProduct, MYGHADLISTSIZE > &  vec,
G4int &  vecLen,
G4ReactionProduct &  modifiedOriginal,
const G4HadProjectile *  originalIncident,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
const G4Nucleus &  targetNucleus,
G4bool &  incidentHasChanged,
G4bool &  targetHasChanged,
G4bool  leadFlag,
G4ReactionProduct &  leadingStrangeParticle 
)

Definition at line 1351 of file FullModelReactionDynamics.cc.

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

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

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