G4ReactionDynamics Class Reference

#include <G4ReactionDynamics.hh>

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

---

Detailed Description

Definition at line 48 of file G4ReactionDynamics.hh.

---

Constructor & Destructor Documentation

G4ReactionDynamics::G4ReactionDynamics

(

)

[inline]

Definition at line 52 of file G4ReactionDynamics.hh.

{}

virtual G4ReactionDynamics::~G4ReactionDynamics

(

)

[inline, virtual]

Definition at line 54 of file G4ReactionDynamics.hh.

{}

---

Member Function Documentation

G4int G4ReactionDynamics::Factorial

(

G4int

n

)

Definition at line 2743 of file G4ReactionDynamics.cc.

   {   // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
     G4int mn = std::min(n,10);
     G4int result = 1;
     if( mn <= 1 )return result;
     for( G4int i=2; i<=mn; ++i )result *= i;
     return result;
   }

virtual G4double G4ReactionDynamics::FindInelasticity

(

)

[inline, virtual]

Definition at line 56 of file G4ReactionDynamics.hh.

    { return 0.0; }

virtual G4double G4ReactionDynamics::FindTimeDelay

(

)

[inline, virtual]

Definition at line 59 of file G4ReactionDynamics.hh.

    { return 0.0; }

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

Definition at line 2491 of file G4ReactionDynamics.cc.

References G4cerr, G4endl, G4UniformRand, G4InuclParticleNames::s0, and G4InuclParticleNames::sm.

Referenced by GenerateXandPt(), NuclearReaction(), and TwoCluster().

  {
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    // derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
    // Returns the weight of the event
    //
    G4int i;
    const G4double expxu =  82.;           // upper bound for arg. of exp
    const G4double expxl = -expxu;         // lower bound for arg. of exp
    if( vecLen < 2 )
    {
      G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
      G4cerr << "    number of particles < 2" << G4endl;
      G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
      return -1.0;
    }
    G4double mass[18];    // mass of each particle
    G4double energy[18];  // total energy of each particle
    G4double pcm[3][18];           // pcm is an array with 3 rows and vecLen columns
    
    G4double totalMass = 0.0;
    G4double extraMass = 0;
    G4double sm[18];
    
    for( i=0; i<vecLen; ++i )
    {
      mass[i] = vec[i]->GetMass()/GeV;
      if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
      vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
      pcm[0][i] = 0.0;      // x-momentum of i-th particle
      pcm[1][i] = 0.0;      // y-momentum of i-th particle
      pcm[2][i] = 0.0;      // z-momentum of i-th particle
      energy[i] = mass[i];  // total energy of i-th particle
      totalMass += mass[i];
      sm[i] = totalMass;
    }
    G4double totalE = totalEnergy/GeV;
    if( totalMass > totalE )
    {
      //G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
      //G4cerr << "    total mass (" << totalMass*GeV << "MeV) > total energy ("
      //     << totalEnergy << "MeV)" << G4endl;
      totalE = totalMass;
      return -1.0;
    }
    G4double kineticEnergy = totalE - totalMass;
    G4double emm[18];
    //G4double *emm = new G4double [vecLen];
    emm[0] = mass[0];
    emm[vecLen-1] = totalE;
    if( vecLen > 2 )          // the random numbers are sorted
    {
      G4double ran[18];
      for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
      for( i=0; i<vecLen-2; ++i )
      {
        for( G4int j=vecLen-2; j>i; --j )
        {
          if( ran[i] > ran[j] )
          {
            G4double temp = ran[i];
            ran[i] = ran[j];
            ran[j] = temp;
          }
        }
      }
      for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
    }
    //   Weight is the sum of logarithms of terms instead of the product of terms
    G4bool lzero = true;    
    G4double wtmax = 0.0;
    if( constantCrossSection )     // this is KGENEV=1 in PHASP
    {
      G4double emmax = kineticEnergy + mass[0];
      G4double emmin = 0.0;
      for( i=1; i<vecLen; ++i )
      {
        emmin += mass[i-1];
        emmax += mass[i];
        G4double wtfc = 0.0;
        if( emmax*emmax > 0.0 )
        {
          G4double arg = emmax*emmax
            + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
            - 2.0*(emmin*emmin+mass[i]*mass[i]);
          if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
        }
        if( wtfc == 0.0 )
        {
          lzero = false;
          break;
        }
        wtmax += std::log( wtfc );
      }
      if( lzero )
        wtmax = -wtmax;
      else
        wtmax = expxu;
    }
    else
    {
      //   ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
      const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
                                 256.3704, 268.4705, 240.9780, 189.2637,
                                 132.1308,  83.0202,  47.4210,  24.8295,
                                 12.0006,   5.3858,   2.2560,   0.8859 };
      wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
    }
    lzero = true;
    G4double pd[50];
    //G4double *pd = new G4double [vecLen-1];
    for( i=0; i<vecLen-1; ++i )
    {
      pd[i] = 0.0;
      if( emm[i+1]*emm[i+1] > 0.0 )
      {
        G4double arg = emm[i+1]*emm[i+1]
          + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
            /(emm[i+1]*emm[i+1])
          - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
        if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
      }
      if( pd[i] <= 0.0 )    //  changed from  ==  on 02 April 98
        lzero = false;
      else
        wtmax += std::log( pd[i] );
    }
    G4double weight = 0.0;           // weight is returned by GenerateNBodyEvent
    if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
    
    G4double bang, cb, sb, s0, s1, s2, c, ss, esys, a, b, gama, beta;
    pcm[0][0] = 0.0;
    pcm[1][0] = pd[0];
    pcm[2][0] = 0.0;
    for( i=1; i<vecLen; ++i )
    {
      pcm[0][i] = 0.0;
      pcm[1][i] = -pd[i-1];
      pcm[2][i] = 0.0;
      bang = twopi*G4UniformRand();
      cb = std::cos(bang);
      sb = std::sin(bang);
      c = 2.0*G4UniformRand() - 1.0;
      ss = std::sqrt( std::fabs( 1.0-c*c ) );
      if( i < vecLen-1 )
      {
        esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
        beta = pd[i]/esys;
        gama = esys/emm[i];
        for( G4int j=0; j<=i; ++j )
        {
          s0 = pcm[0][j];
          s1 = pcm[1][j];
          s2 = pcm[2][j];
          energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
          a = s0*c - s1*ss;                           //  rotation
          pcm[1][j] = s0*ss + s1*c;
          b = pcm[2][j];
          pcm[0][j] = a*cb - b*sb;
          pcm[2][j] = a*sb + b*cb;
          pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
        }
      }
      else
      {
        for( G4int j=0; j<=i; ++j )
        {
          s0 = pcm[0][j];
          s1 = pcm[1][j];
          s2 = pcm[2][j];
          energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
          a = s0*c - s1*ss;                           //  rotation
          pcm[1][j] = s0*ss + s1*c;
          b = pcm[2][j];
          pcm[0][j] = a*cb - b*sb;
          pcm[2][j] = a*sb + b*cb;
        }
      }
    }
    for( i=0; i<vecLen; ++i )
    {
      vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
      vec[i]->SetTotalEnergy( energy[i]*GeV );
    }

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    return weight;
  }

G4bool G4ReactionDynamics::GenerateXandPt	(	G4FastVector< G4ReactionProduct, GHADLISTSIZE > &	vec,
		G4int &	vecLen,
		G4ReactionProduct &	modifiedOriginal,
		const G4HadProjectile *	originalIncident,
		G4ReactionProduct &	currentParticle,
		G4ReactionProduct &	targetParticle,
		const G4DynamicParticle *	originalTarget,
		const G4Nucleus &	targetNucleus,
		G4bool &	incidentHasChanged,
		G4bool &	targetHasChanged,
		G4bool	leadFlag,
		G4ReactionProduct &	leadingStrangeParticle
	)

Definition at line 105 of file G4ReactionDynamics.cc.

References FatalException, G4cerr, G4endl, G4Exception(), G4UniformRand, GenerateNBodyEvent(), G4Nucleus::GetA_asInt(), G4Nucleus::GetAnnihilationDTABlackTrackEnergy(), G4Nucleus::GetAnnihilationPNBlackTrackEnergy(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4Nucleus::GetN_asInt(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetSide(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4FastVector< Type, N >::Initialize(), G4Lambda::Lambda(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), G4INCL::Math::pi, G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4InuclParticleNames::pp, G4Proton::Proton(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetSide(), G4ReactionProduct::SetTotalEnergy(), G4ReactionProduct::SetZero(), and sqr().

Referenced by G4InelasticInteraction::CalculateMomenta().

  {
    // 
    // derived from original FORTRAN code GENXPT by H. Fesefeldt (11-Oct-1987)
    //
    //  Generation of X- and PT- values for incident, target, and all secondary particles 
    //  A simple single variable description E D3S/DP3= F(Q) with
    //  Q^2 = (M*X)^2 + PT^2 is used. Final state kinematic is produced
    //  by an FF-type iterative cascade method
    //
    //  internal units are GeV
    //
    // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    
    // Protection in case no secondary has been created; cascades down to two-body.
    if(vecLen == 0) return false;

    G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
    G4ParticleDefinition *aProton = G4Proton::Proton();
    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
    G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
    G4ParticleDefinition *aPiZero = G4PionZero::PionZero();

    G4int i, l;
    G4bool veryForward = false;
    
    const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
    const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
    const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
    const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
    G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
    G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
                                        targetMass*targetMass +
                                        2.0*targetMass*etOriginal );  // GeV
    G4double currentMass = currentParticle.GetMass()/GeV;
    targetMass = targetParticle.GetMass()/GeV;
    //
    //  randomize the order of the secondary particles
    //  note that the current and target particles are not affected
    //
    for( i=0; i<vecLen; ++i )
    {
      G4int itemp = G4int( G4UniformRand()*vecLen );
      G4ReactionProduct pTemp = *vec[itemp];
      *vec[itemp] = *vec[i];
      *vec[i] = pTemp;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }

    if( currentMass == 0.0 && targetMass == 0.0 )
    {
      // Target and projectile have annihilated.  Replace them with the first 
      // two secondaries in the list.  Current particle KE is maintained.
 
      G4double ek = currentParticle.GetKineticEnergy();
      G4ThreeVector mom = currentParticle.GetMomentum();
      currentParticle = *vec[0];
      targetParticle = *vec[1];
      for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
      G4ReactionProduct *temp = vec[vecLen-1];
      delete temp;
      temp = vec[vecLen-2];
      delete temp;
      vecLen -= 2;
      currentMass = currentParticle.GetMass()/GeV;
      targetMass = targetParticle.GetMass()/GeV;
      incidentHasChanged = true;
      targetHasChanged = true;
      currentParticle.SetKineticEnergy( ek );
      currentParticle.SetMomentum( mom );
      veryForward = true;
    }
    const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
    const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
    const G4double protonMass = aProton->GetPDGMass()/MeV;

    if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
        && G4UniformRand() >= 0.7) {
      G4ReactionProduct temp = currentParticle;
      currentParticle = targetParticle;
      targetParticle = temp;
      incidentHasChanged = true;
      targetHasChanged = true;
      currentMass = currentParticle.GetMass()/GeV;
      targetMass = targetParticle.GetMass()/GeV;
    }
    const G4double afc = std::min( 0.75,
          0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
          std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
    
    G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
    G4double forwardEnergy = freeEnergy/2.;
    G4int forwardCount = 1;         // number of particles in forward hemisphere
    
    G4double backwardEnergy = freeEnergy/2.;
    G4int backwardCount = 1;        // number of particles in backward hemisphere

    if(veryForward)
    {
      if(currentParticle.GetSide()==-1)
      {
        forwardEnergy += currentMass;
        forwardCount --;
        backwardEnergy -= currentMass;
        backwardCount ++;
      }
      if(targetParticle.GetSide()!=-1)
      {
        backwardEnergy += targetMass;
        backwardCount --;
        forwardEnergy -= targetMass;
        forwardCount ++;
      }
    }

    for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetSide() == -1 )
      {
        ++backwardCount;
        backwardEnergy -= vec[i]->GetMass()/GeV;
      } else {
        ++forwardCount;
        forwardEnergy -= vec[i]->GetMass()/GeV;
      }
    }
    //
    //  Add particles from intranuclear cascade.
    //  nuclearExcitationCount = number of new secondaries produced by nuclear excitation
    //  extraCount = number of nucleons within these new secondaries
    //
    G4double xtarg;
    if( centerofmassEnergy < (2.0+G4UniformRand()) )
      xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
    else
      xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
    if( xtarg <= 0.0 )xtarg = 0.01;
    G4int nuclearExcitationCount = Poisson( xtarg );
    if(atomicWeight<1.0001) nuclearExcitationCount = 0;
    G4int extraNucleonCount = 0;
    G4double extraNucleonMass = 0.0;
    if( nuclearExcitationCount > 0 )
    {
      const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
      const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
      G4int momentumBin = 0;
      while( (momentumBin < 6) &&
             (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
        ++momentumBin;
      momentumBin = std::min( 5, momentumBin );
      //
      //  NOTE: in GENXPT, these new particles were given negative codes
      //        here I use  NewlyAdded = true  instead
      //

      for( i=0; i<nuclearExcitationCount; ++i )
      {
        G4ReactionProduct * pVec = new G4ReactionProduct();
        if( G4UniformRand() < nucsup[momentumBin] )
        {
          if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
            pVec->SetDefinition( aProton );
          else
            pVec->SetDefinition( aNeutron );
          pVec->SetSide( -2 );                // -2 means backside nucleon
          ++extraNucleonCount;
          backwardEnergy += pVec->GetMass()/GeV;
          extraNucleonMass += pVec->GetMass()/GeV;
        }
        else
        {
          G4double ran = G4UniformRand();
          if( ran < 0.3181 )
            pVec->SetDefinition( aPiPlus );
          else if( ran < 0.6819 )
            pVec->SetDefinition( aPiZero );
          else
            pVec->SetDefinition( aPiMinus );
          pVec->SetSide( -1 );                // backside particle, but not a nucleon
        }
        pVec->SetNewlyAdded( true );                // true is the same as IPA(i)<0
        vec.SetElement( vecLen++, pVec );    
        // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        backwardEnergy -= pVec->GetMass()/GeV;
        ++backwardCount;
      }
    }
    //
    //  assume conservation of kinetic energy in forward & backward hemispheres
    //
    G4int is, iskip;
    while( forwardEnergy <= 0.0 )  // must eliminate a particle from the forward side
    {
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      iskip = G4int(G4UniformRand()*forwardCount) + 1; // 1 <= iskip <= forwardCount
      is = 0;
      G4int forwardParticlesLeft = 0;
      for( i=(vecLen-1); i>=0; --i )
      {
        if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
        {
          forwardParticlesLeft = 1;  
          if( ++is == iskip )
          { 
            forwardEnergy += vec[i]->GetMass()/GeV;
            for( G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1];    // shift up
            --forwardCount;
            delete vec[vecLen-1];
            if( --vecLen == 0 )return false;  // all the secondaries have been eliminated
            break;  // --+
          }         //   |
        }           //   |
      }             // break goes down to here
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      if( forwardParticlesLeft == 0 )
      {
        G4int iremove = -1;
        for (i = 0; i < vecLen; i++) {
          if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
            iremove = i;
            break;
          }
        }
        if (iremove == -1) {
          for (i = 0; i < vecLen; i++) {
            if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
              iremove = i;
              break;
            }
          }
        }
        if (iremove == -1) iremove = 0;

        forwardEnergy += vec[iremove]->GetMass()/GeV;
        if (vec[iremove]->GetSide() > 0) --forwardCount;
 
        for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
        delete vec[vecLen-1];
        vecLen--;
        if (vecLen == 0) return false;  // all secondaries have been eliminated
        break;
      }
    } // while

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    while( backwardEnergy <= 0.0 )  // must eliminate a particle from the backward side
    {
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      iskip = G4int(G4UniformRand()*backwardCount) + 1; // 1 <= iskip <= backwardCount
      is = 0;
      G4int backwardParticlesLeft = 0;
      for( i=(vecLen-1); i>=0; --i )
      {
        if( vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled())
        {
          backwardParticlesLeft = 1;
          if( ++is == iskip )        // eliminate the i'th particle
          {
            if( vec[i]->GetSide() == -2 )
            {
              --extraNucleonCount;
              extraNucleonMass -= vec[i]->GetMass()/GeV;
              backwardEnergy -= vec[i]->GetTotalEnergy()/GeV;
            }
            backwardEnergy += vec[i]->GetTotalEnergy()/GeV;
            for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];   // shift up
            --backwardCount;
            delete vec[vecLen-1];
            if( --vecLen == 0 )return false;  // all the secondaries have been eliminated
            break;
          }
        }
      }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      if( backwardParticlesLeft == 0 )
      {
        G4int iremove = -1;
        for (i = 0; i < vecLen; i++) {
          if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
            iremove = i;
            break;
          }
        }
        if (iremove == -1) {
          for (i = 0; i < vecLen; i++) {
            if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
              iremove = i;
              break;
            }
          }
        }
        if (iremove == -1) iremove = 0;
 
        backwardEnergy += vec[iremove]->GetMass()/GeV;
        if (vec[iremove]->GetSide() > 0) --backwardCount;
 
        for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
        delete vec[vecLen-1];
        vecLen--;
        if (vecLen == 0) return false;  // all secondaries have been eliminated
        break;
      }
    } // while

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    //
    //  define initial state vectors for Lorentz transformations
    //  the pseudoParticles have non-standard masses, hence the "pseudo"
    //
    G4ReactionProduct pseudoParticle[10];
    for( i=0; i<10; ++i )pseudoParticle[i].SetZero();
    
    pseudoParticle[0].SetMass( mOriginal*GeV );
    pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
    pseudoParticle[0].SetTotalEnergy(
     std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
    
    pseudoParticle[1].SetMass( protonMass*MeV );        // this could be targetMass
    pseudoParticle[1].SetTotalEnergy( protonMass*MeV );
    
    pseudoParticle[3].SetMass( protonMass*(1+extraNucleonCount)*MeV );
    pseudoParticle[3].SetTotalEnergy( protonMass*(1+extraNucleonCount)*MeV );
    
    pseudoParticle[8].SetMomentum( 1.0*GeV, 0.0, 0.0 );
    
    pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
    pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
    
    pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
    pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
    
    G4double dndl[20];
    //
    //  main loop for 4-momentum generation
    //  see Pitha-report (Aachen) for a detailed description of the method
    //
    G4double aspar, pt, et, x, pp, pp1, rmb, wgt;
    G4int    innerCounter, outerCounter;
    G4bool   eliminateThisParticle, resetEnergies, constantCrossSection;
    
    G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
    //
    // process the secondary particles in reverse order
    // the incident particle is Done after the secondaries
    // nucleons, including the target, in the backward hemisphere are also Done later
    //
    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,
                         1.43,1.67,2.0,2.5,3.33,5.00,10.00};
    G4int backwardNucleonCount = 0;       // number of nucleons in backward hemisphere
    G4double totalEnergy, kineticEnergy, vecMass;

    for( i=(vecLen-1); i>=0; --i )
    {
      G4double phi = G4UniformRand()*twopi;
      if( vec[i]->GetNewlyAdded() )           // added from intranuclear cascade
      {
        if( vec[i]->GetSide() == -2 )         //  is a nucleon
        {
          if( backwardNucleonCount < 18 )
          {
            if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
              for(G4int j=0; j<vecLen; j++) delete vec[j];
              vecLen = 0;
              throw G4HadReentrentException(__FILE__, __LINE__,
              "G4ReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
            }
            vec[i]->SetSide( -3 );
            ++backwardNucleonCount;
            continue;
          }
        }
      }
      //
      //  set pt and phi values, they are changed somewhat in the iteration loop
      //  set mass parameter for lambda fragmentation model
      //
      vecMass = vec[i]->GetMass()/GeV;
      G4double ran = -std::log(1.0-G4UniformRand())/3.5;
      if( vec[i]->GetSide() == -2 )   // backward nucleon
      {
        if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon" ||
            vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
          aspar = 0.75;
          pt = std::sqrt( std::pow( ran, 1.7 ) );
        } else {                          // vec[i] must be a proton, neutron,
          aspar = 0.20;                   //  lambda, sigma, xsi, or ion
          pt = std::sqrt( std::pow( ran, 1.2 ) );
        }

      } else {                          // not a backward nucleon

        if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
          aspar = 0.75;
          pt = std::sqrt( std::pow( ran, 1.7 ) );
        } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
          aspar = 0.70;
          pt = std::sqrt( std::pow( ran, 1.7 ) );
        } else {                        // vec[i] must be a proton, neutron, 
          aspar = 0.65;                 //  lambda, sigma, xsi, or ion
          pt = std::sqrt( std::pow( ran, 1.5 ) );
        }
      }
      pt = std::max( 0.001, pt );
      vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
      for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
      if( vec[i]->GetSide() > 0 )
        et = pseudoParticle[0].GetTotalEnergy()/GeV;
      else
        et = pseudoParticle[1].GetTotalEnergy()/GeV;
      dndl[0] = 0.0;
      //
      //   start of outer iteration loop
      //
      outerCounter = 0;
      eliminateThisParticle = true;
      resetEnergies = true;
      while( ++outerCounter < 3 )
      {
        for( l=1; l<20; ++l )
        {
          x = (binl[l]+binl[l-1])/2.;
          pt = std::max( 0.001, pt );
          if( x > 1.0/pt )
            dndl[l] += dndl[l-1];   //  changed from just  =  on 02 April 98
          else
            dndl[l] = et * aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) )
              * (binl[l]-binl[l-1]) / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass )
              + dndl[l-1];
        }
        innerCounter = 0;
        vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
        //
        //   start of inner iteration loop
        //
        while( ++innerCounter < 7 )
        {
          ran = G4UniformRand()*dndl[19];
          l = 1;
          while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
          l = std::min( 19, l );
          x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
          if( vec[i]->GetSide() < 0 )x *= -1.;
          vec[i]->SetMomentum( x*et*GeV );              // set the z-momentum
          totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
          vec[i]->SetTotalEnergy( totalEnergy*GeV );
          kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
          if( vec[i]->GetSide() > 0 )                            // forward side
          {
            if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy )
            {
              pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
              forwardKinetic += kineticEnergy;
              pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
              pseudoParticle[6].SetMomentum( 0.0 );           // set the z-momentum
              phi = pseudoParticle[6].Angle( pseudoParticle[8] );
              if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
              phi += pi + normal()*pi/12.0;
              if( phi > twopi )phi -= twopi;
              if( phi < 0.0 )phi = twopi - phi;
              outerCounter = 2;                     // leave outer loop
              eliminateThisParticle = false;        // don't eliminate this particle
              resetEnergies = false;
              break;                                // leave inner loop
            }
            if( innerCounter > 5 )break;           // leave inner loop
            if( backwardEnergy >= vecMass )  // switch sides
            {
              vec[i]->SetSide( -1 );
              forwardEnergy += vecMass;
              backwardEnergy -= vecMass;
              ++backwardCount;
            }
          } else {                                                 // backward side
           if( extraNucleonCount > 19 )   // commented out to duplicate ?bug? in GENXPT
             x = 0.999;
           G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
           if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy )
            {
              pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
              backwardKinetic += kineticEnergy;
              pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
              pseudoParticle[6].SetMomentum( 0.0 );           // set the z-momentum
              phi = pseudoParticle[6].Angle( pseudoParticle[8] );
              if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
              phi += pi + normal() * pi / 12.0;
              if( phi > twopi )phi -= twopi;
              if( phi < 0.0 )phi = twopi - phi;
              outerCounter = 2;                    // leave outer loop
              eliminateThisParticle = false;       // don't eliminate this particle
              resetEnergies = false;
              break;                               // leave inner loop
            }
            if( innerCounter > 5 )break;           // leave inner loop
            if( forwardEnergy >= vecMass )  // switch sides
            {
              vec[i]->SetSide( 1 );
              forwardEnergy -= vecMass;
              backwardEnergy += vecMass;
              backwardCount--;
            }
          }
          G4ThreeVector momentum = vec[i]->GetMomentum();
          vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
          pt *= 0.9;
          dndl[19] *= 0.9;
        }                       // closes inner loop
        if( resetEnergies )
        {
          //   if we get to here, the inner loop has been Done 6 Times
          //   reset the kinetic energies of previously Done particles, if they are lighter
          //    than protons and in the forward hemisphere
          //   then continue with outer loop
          //
          forwardKinetic = 0.0;
          backwardKinetic = 0.0;
          pseudoParticle[4].SetZero();
          pseudoParticle[5].SetZero();
          for( l=i+1; l<vecLen; ++l )
          {
            if (vec[l]->GetSide() > 0 ||
                vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
                vec[l]->GetDefinition()->GetParticleSubType() == "pi") {

              G4double tempMass = vec[l]->GetMass()/MeV;
              totalEnergy = 0.95*vec[l]->GetTotalEnergy()/MeV + 0.05*tempMass;
              totalEnergy = std::max( tempMass, totalEnergy );
              vec[l]->SetTotalEnergy( totalEnergy*MeV );
              pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
              pp1 = vec[l]->GetMomentum().mag()/MeV;
              if( pp1 < 1.0e-6*GeV )
              {
                G4double costheta = 2.*G4UniformRand() - 1.;
                G4double sintheta = std::sqrt(1. - costheta*costheta);
                G4double phi2 = twopi*G4UniformRand();
                vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                     pp*sintheta*std::sin(phi2)*MeV,
                                     pp*costheta*MeV ) ;
              } else {
                vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
              }
              G4double px = vec[l]->GetMomentum().x()/MeV;
              G4double py = vec[l]->GetMomentum().y()/MeV;
              pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/GeV;
              if( vec[l]->GetSide() > 0 )
              {
                forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
                pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
              } else {
                backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
                pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
              }
            } // if pi, K or forward
          } // for l
        } // if resetEnergies
      } // closes outer loop
              
      if( eliminateThisParticle && vec[i]->GetMayBeKilled())  // not enough energy, eliminate this particle
      {
        if( vec[i]->GetSide() > 0 )
        {
          --forwardCount;
          forwardEnergy += vecMass;
        } else {
          if( vec[i]->GetSide() == -2 )
          {
            --extraNucleonCount;
            extraNucleonMass -= vecMass;
            backwardEnergy -= vecMass;
          }
          --backwardCount;
          backwardEnergy += vecMass;
        }
        for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
        G4ReactionProduct *temp = vec[vecLen-1];
        delete temp;
        // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        if( --vecLen == 0 )return false;  // all the secondaries have been eliminated
        pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
        pseudoParticle[6].SetMomentum( 0.0 );                 // set z-momentum
        phi = pseudoParticle[6].Angle( pseudoParticle[8] );
        if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
        phi += pi + normal() * pi / 12.0;
        if( phi > twopi )phi -= twopi;
        if( phi < 0.0 )phi = twopi - phi;
      }
    }   // closes main for loop

    //
    //  for the incident particle:  it was placed in the forward hemisphere
    //   set pt and phi values, they are changed somewhat in the iteration loop
    //   set mass parameter for lambda fragmentation model
    //
    G4double phi = G4UniformRand()*twopi;
    G4double ran = -std::log(1.0-G4UniformRand());
    if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
      aspar = 0.60;
      pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
    } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
      aspar = 0.50;
      pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
    } else {
      aspar = 0.40;
      pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
    }

    for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
    currentParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
    et = pseudoParticle[0].GetTotalEnergy()/GeV;
    dndl[0] = 0.0;
    vecMass = currentParticle.GetMass()/GeV;
    for( l=1; l<20; ++l )
    {
      x = (binl[l]+binl[l-1])/2.;
      if( x > 1.0/pt )
        dndl[l] += dndl[l-1];   //  changed from just  =   on 02 April 98
      else
        dndl[l] = aspar/std::sqrt( std::pow((1.+sqr(aspar*x)),3) ) *
          (binl[l]-binl[l-1]) * et / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass ) +
          dndl[l-1];
    }
    ran = G4UniformRand()*dndl[19];
    l = 1;
    while( (ran>dndl[l]) && (l<20) )l++;
    l = std::min( 19, l );
    x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
    currentParticle.SetMomentum( x*et*GeV );                 // set the z-momentum
    if( forwardEnergy < forwardKinetic )
      totalEnergy = vecMass + 0.04*std::fabs(normal());
    else
      totalEnergy = vecMass + forwardEnergy - forwardKinetic;
    currentParticle.SetTotalEnergy( totalEnergy*GeV );
    pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
    pp1 = currentParticle.GetMomentum().mag()/MeV;
    if( pp1 < 1.0e-6*GeV )
    {
      G4double costheta = 2.*G4UniformRand() - 1.;
      G4double sintheta = std::sqrt(1. - costheta*costheta);
      G4double phi2 = twopi*G4UniformRand();
      currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                   pp*sintheta*std::sin(phi2)*MeV,
                                   pp*costheta*MeV ) ;
    } else {
      currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
    }
    pseudoParticle[4] = pseudoParticle[4] + currentParticle;

    //
    // Current particle now finished
    //
    // Begin target particle
    //

    if( backwardNucleonCount < 18 )
    {
      targetParticle.SetSide( -3 );
      ++backwardNucleonCount;
    }
    else
    {
      //  set pt and phi values, they are changed somewhat in the iteration loop
      //  set mass parameter for lambda fragmentation model
      //
      vecMass = targetParticle.GetMass()/GeV;
      ran = -std::log(1.0-G4UniformRand());
      aspar = 0.40;
      pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
      targetParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
      for( G4int j=0; j<20; ++j )binl[j] = (j-1.)/(19.*pt);
      et = pseudoParticle[1].GetTotalEnergy()/GeV;
      dndl[0] = 0.0;
      outerCounter = 0;
      eliminateThisParticle = true;     // should never eliminate the target particle
      resetEnergies = true;
      while( ++outerCounter < 3 )     // start of outer iteration loop
      {
        for( l=1; l<20; ++l )
        {
          x = (binl[l]+binl[l-1])/2.;
          if( x > 1.0/pt )
            dndl[l] += dndl[l-1];   // changed from just  =  on 02 April 98
          else
            dndl[l] = aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) ) *
              (binl[l]-binl[l-1])*et / std::sqrt( pt*x*et*pt*x*et+pt*pt+vecMass*vecMass ) +
              dndl[l-1];
        }
        innerCounter = 0;
        while( ++innerCounter < 7 )    // start of inner iteration loop
        {
          l = 1;
          ran = G4UniformRand()*dndl[19];
          while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
          l = std::min( 19, l );
          x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
          if( targetParticle.GetSide() < 0 )x *= -1.;
          targetParticle.SetMomentum( x*et*GeV );                // set the z-momentum
          totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
          targetParticle.SetTotalEnergy( totalEnergy*GeV );
          if( targetParticle.GetSide() < 0 )
          {
            if( extraNucleonCount > 19 )x=0.999;
            G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
            if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy )
            {
              pseudoParticle[5] = pseudoParticle[5] + targetParticle;
              backwardKinetic += totalEnergy - vecMass;
              pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
              pseudoParticle[6].SetMomentum( 0.0 );                      // set z-momentum
              phi = pseudoParticle[6].Angle( pseudoParticle[8] );
              if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
              phi += pi + normal() * pi / 12.0;
              if( phi > twopi )phi -= twopi;
              if( phi < 0.0 )phi = twopi - phi;
              outerCounter = 2;                    // leave outer loop
              eliminateThisParticle = false;       // don't eliminate this particle
              resetEnergies = false;
              break;                               // leave inner loop
            }
            if( innerCounter > 5 )break;           // leave inner loop
            if( forwardEnergy >= vecMass )  // switch sides
            {
              targetParticle.SetSide( 1 );
              forwardEnergy -= vecMass;
              backwardEnergy += vecMass;
              --backwardCount;
            }
            G4ThreeVector momentum = targetParticle.GetMomentum();
            targetParticle.SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
            pt *= 0.9;
            dndl[19] *= 0.9;
          }
          else                    // target has gone to forward side
          {
            if( forwardEnergy < forwardKinetic )
              totalEnergy = vecMass + 0.04*std::fabs(normal());
            else
              totalEnergy = vecMass + forwardEnergy - forwardKinetic;
            targetParticle.SetTotalEnergy( totalEnergy*GeV );
            pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
            pp1 = targetParticle.GetMomentum().mag()/MeV;
            if( pp1 < 1.0e-6*GeV )
            {
              G4double costheta = 2.*G4UniformRand() - 1.;
              G4double sintheta = std::sqrt(1. - costheta*costheta);
              G4double phi2 = twopi*G4UniformRand();
              targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                          pp*sintheta*std::sin(phi2)*MeV,
                                          pp*costheta*MeV ) ;
            }
            else
              targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );

            pseudoParticle[4] = pseudoParticle[4] + targetParticle;
            outerCounter = 2;                    // leave outer loop
            eliminateThisParticle = false;       // don't eliminate this particle
            resetEnergies = false;
            break;                               // leave inner loop
          }
        }    // closes inner loop
        if( resetEnergies )
        {
          //   if we get to here, the inner loop has been Done 6 Times
          //   reset the kinetic energies of previously Done particles, if they are lighter
          //    than protons and in the forward hemisphere
          //   then continue with outer loop
          
          forwardKinetic = backwardKinetic = 0.0;
          pseudoParticle[4].SetZero();
          pseudoParticle[5].SetZero();
          for( l=0; l<vecLen; ++l ) // changed from l=1  on 02 April 98
          {
            if (vec[l]->GetSide() > 0 ||
                vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
                vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
              G4double tempMass = vec[l]->GetMass()/GeV;
              totalEnergy =
                std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/GeV + 0.05*tempMass );
              vec[l]->SetTotalEnergy( totalEnergy*GeV );
              pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*GeV;
              pp1 = vec[l]->GetMomentum().mag()/MeV;
              if( pp1 < 1.0e-6*GeV )
              {
                G4double costheta = 2.*G4UniformRand() - 1.;
                G4double sintheta = std::sqrt(1. - costheta*costheta);
                G4double phi2 = twopi*G4UniformRand();
                vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                     pp*sintheta*std::sin(phi2)*MeV,
                                     pp*costheta*MeV ) ;
              }
              else
                vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );

              pt = std::max( 0.001*GeV, std::sqrt( sqr(vec[l]->GetMomentum().x()/MeV) +
                                         sqr(vec[l]->GetMomentum().y()/MeV) ) )/GeV;
              if( vec[l]->GetSide() > 0)
              {
                forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
                pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
              } else {
                backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
                pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
              }
            } // if pi, K or forward
          } // for l
        } // if (resetEnergies)
      } // closes outer loop
    }

    //
    // Target particle finished.
    //
    // Now produce backward nucleons with a cluster model
    //
    pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
    pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
    pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
    if( backwardNucleonCount == 1 )  // target particle is the only backward nucleon
    {
      G4double ekin =
        std::min( backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV );

      if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
      vecMass = targetParticle.GetMass()/GeV;
      totalEnergy = ekin+vecMass;
      targetParticle.SetTotalEnergy( totalEnergy*GeV );
      pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
      pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
      if( pp1 < 1.0e-6*GeV )
      { 
        G4double costheta = 2.*G4UniformRand() - 1.;
        G4double sintheta = std::sqrt(1. - costheta*costheta);
        G4double phi2 = twopi*G4UniformRand();
        targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                    pp*sintheta*std::sin(phi2)*MeV,
                                    pp*costheta*MeV ) ;
      } else {
        targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1) );
      }
      pseudoParticle[5] = pseudoParticle[5] + targetParticle;
    }
    else  // more than one backward nucleon
    {
      const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
      const G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
      // Replaced the following min function to get correct behaviour on DEC.
      // G4int tempCount = std::min( 5, backwardNucleonCount ) - 1;
      G4int tempCount;
      if (backwardNucleonCount < 5)
        {
          tempCount = backwardNucleonCount;
        }
      else
        {
          tempCount = 5;
        }
      tempCount--;
      //cout << "backwardNucleonCount " << backwardNucleonCount << G4endl;
      //cout << "tempCount " << tempCount << G4endl;
      G4double rmb0 = 0.0;
      if( targetParticle.GetSide() == -3 )
        rmb0 += targetParticle.GetMass()/GeV;
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == -3 )rmb0 += vec[i]->GetMass()/GeV;
      }
      rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()),cpar[tempCount]) / gpar[tempCount];
      totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
      vecMass = std::min( rmb, totalEnergy );
      pseudoParticle[6].SetMass( vecMass*GeV );
      pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
      pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
      if( pp1 < 1.0e-6*GeV )
      {
        G4double costheta = 2.*G4UniformRand() - 1.;
        G4double sintheta = std::sqrt(1. - costheta*costheta);
        G4double phi2 = twopi*G4UniformRand();
        pseudoParticle[6].SetMomentum( -pp*sintheta*std::cos(phi2)*MeV,
                                       -pp*sintheta*std::sin(phi2)*MeV,
                                       -pp*costheta*MeV ) ;
      }
      else
        pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
      
      G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;  // tempV contains the backward nucleons
      tempV.Initialize( backwardNucleonCount );
      G4int tempLen = 0;
      if( targetParticle.GetSide() == -3 )tempV.SetElement( tempLen++, &targetParticle );
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == -3 )tempV.SetElement( tempLen++, vec[i] );
      }
      if( tempLen != backwardNucleonCount )
      {
        G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
        G4cerr << "tempLen = " << tempLen;
        G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
        G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
        G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
        for( i=0; i<vecLen; ++i )
          G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
        G4Exception("G4ReactionDynamics::GenerateXandPt", "601", 
                    FatalException, "Mismatch in nucleon count");
      }
      constantCrossSection = true;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      if( tempLen >= 2 )
      {
        wgt = GenerateNBodyEvent(
         pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen );
         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        if( targetParticle.GetSide() == -3 )
        {
          targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
          // tempV contains the real stuff
          pseudoParticle[5] = pseudoParticle[5] + targetParticle;
        }
        for( i=0; i<vecLen; ++i )
        {
          if( vec[i]->GetSide() == -3 )
          {
            vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
            pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
          }
        }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      }
    }
    //
    //  Lorentz transformation in lab system
    //
    if( vecLen == 0 )return false;  // all the secondaries have been eliminated
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    
    currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
    targetParticle.Lorentz( targetParticle, pseudoParticle[1] );    
    for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[1] );

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    //
    // leadFlag will be true
    //  iff original particle is at least as heavy as K+ and not a proton or 
    //  neutron AND if incident particle is at least as heavy as K+ and it is 
    //  not a proton or neutron leadFlag is set to the incident particle
    //  or
    //  target particle is at least as heavy as K+ and it is not a proton or 
    //  neutron leadFlag is set to the target particle
    //
    G4bool leadingStrangeParticleHasChanged = true;
    if( leadFlag )
    {
      if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
        leadingStrangeParticleHasChanged = false;
      if( leadingStrangeParticleHasChanged &&
          ( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() ) )
        leadingStrangeParticleHasChanged = false;
      if( leadingStrangeParticleHasChanged )
      {
        for( i=0; i<vecLen; i++ )
        {
          if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
          {
            leadingStrangeParticleHasChanged = false;
            break;
          }
        }
      }
      if( leadingStrangeParticleHasChanged )
      {
        G4bool leadTest = 
          (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
           leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
        G4bool targetTest =
          (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
           targetParticle.GetDefinition()->GetParticleSubType() == "pi");
        
        // following modified by JLC 22-Oct-97
          
        if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
        {
          targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
          targetHasChanged = true;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
        else
        {
          currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
          incidentHasChanged = false;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
      }
    }   // end of if( leadFlag )

    // Get number of final state nucleons and nucleons remaining in
    // target nucleus
    
    std::pair<G4int, G4int> finalStateNucleons = 
      GetFinalStateNucleons(originalTarget, vec, vecLen);

    G4int protonsInFinalState = finalStateNucleons.first;
    G4int neutronsInFinalState = finalStateNucleons.second;

    G4int numberofFinalStateNucleons = 
      protonsInFinalState + neutronsInFinalState;

    if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
        targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
        originalIncident->GetDefinition()->GetPDGMass() < 
                                   G4Lambda::Lambda()->GetPDGMass())
      numberofFinalStateNucleons++;

    numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);

    G4int PinNucleus = std::max(0, 
      targetNucleus.GetZ_asInt() - protonsInFinalState);
    G4int NinNucleus = std::max(0,
      targetNucleus.GetN_asInt() - neutronsInFinalState);

    pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
    pseudoParticle[3].SetMass( mOriginal*GeV );
    pseudoParticle[3].SetTotalEnergy(
     std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
    
    G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
    G4int diff = 0;
    if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() )  diff = 1;
    if(numberofFinalStateNucleons == 1) diff = 0;
    pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
    pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV );
    pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV );
    
    G4double theoreticalKinetic =
      pseudoParticle[3].GetTotalEnergy()/MeV +
      pseudoParticle[4].GetTotalEnergy()/MeV -
      currentParticle.GetMass()/MeV -
      targetParticle.GetMass()/MeV;
    
    G4double simulatedKinetic =
      currentParticle.GetKineticEnergy()/MeV + targetParticle.GetKineticEnergy()/MeV;
    
    pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
    pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
    pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
    
    pseudoParticle[7].SetZero();
    pseudoParticle[7] = pseudoParticle[7] + currentParticle;
    pseudoParticle[7] = pseudoParticle[7] + targetParticle;

    for( i=0; i<vecLen; ++i )
    {
      pseudoParticle[7] = pseudoParticle[7] + *vec[i];
      simulatedKinetic += vec[i]->GetKineticEnergy()/MeV;
      theoreticalKinetic -= vec[i]->GetMass()/MeV;
    }

    if( vecLen <= 16 && vecLen > 0 )
    {
      // must create a new set of ReactionProducts here because GenerateNBody will
      // modify the momenta for the particles, and we don't want to do this
      //
      G4ReactionProduct tempR[130];
      tempR[0] = currentParticle;
      tempR[1] = targetParticle;
      for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
      G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
      tempV.Initialize( vecLen+2 );
      G4int tempLen = 0;
      for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
      constantCrossSection = true;

      wgt = GenerateNBodyEvent( pseudoParticle[3].GetTotalEnergy()/MeV+
                                pseudoParticle[4].GetTotalEnergy()/MeV,
                                constantCrossSection, tempV, tempLen );
      if (wgt == -1) {
        G4double Qvalue = 0;
        for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
        wgt = GenerateNBodyEvent( Qvalue/MeV,
                                  constantCrossSection, tempV, tempLen );
      }
      if(wgt>-.5)
      {
        theoreticalKinetic = 0.0;
        for( i=0; i<tempLen; ++i )
        {
          pseudoParticle[6].Lorentz( *tempV[i], pseudoParticle[4] );
          theoreticalKinetic += pseudoParticle[6].GetKineticEnergy()/MeV;
        }
      }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    //
    //  Make sure, that the kinetic energies are correct
    //
    if( simulatedKinetic != 0.0 )
    {
      wgt = (theoreticalKinetic)/simulatedKinetic;
      theoreticalKinetic = currentParticle.GetKineticEnergy()/MeV * wgt;
      simulatedKinetic = theoreticalKinetic;
      currentParticle.SetKineticEnergy( theoreticalKinetic*MeV );
      pp = currentParticle.GetTotalMomentum()/MeV;
      pp1 = currentParticle.GetMomentum().mag()/MeV;
      if( pp1 < 1.0e-6*GeV )
      {
        G4double costheta = 2.*G4UniformRand() - 1.;
        G4double sintheta = std::sqrt(1. - costheta*costheta);
        G4double phi2 = twopi*G4UniformRand();
        currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                     pp*sintheta*std::sin(phi2)*MeV,
                                     pp*costheta*MeV ) ;
      }
      else
      {
        currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
      }
      theoreticalKinetic = targetParticle.GetKineticEnergy()/MeV * wgt;
      targetParticle.SetKineticEnergy( theoreticalKinetic*MeV );
      simulatedKinetic += theoreticalKinetic;
      pp = targetParticle.GetTotalMomentum()/MeV;
      pp1 = targetParticle.GetMomentum().mag()/MeV;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      if( pp1 < 1.0e-6*GeV )
      {
        G4double costheta = 2.*G4UniformRand() - 1.;
        G4double sintheta = std::sqrt(1. - costheta*costheta);
        G4double phi2 = twopi*G4UniformRand();
        targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                                    pp*sintheta*std::sin(phi2)*MeV,
                                    pp*costheta*MeV ) ;
      } else {
        targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
      }
      for( i=0; i<vecLen; ++i )
      {
        theoreticalKinetic = vec[i]->GetKineticEnergy()/MeV * wgt;
        simulatedKinetic += theoreticalKinetic;
        vec[i]->SetKineticEnergy( theoreticalKinetic*MeV );
        pp = vec[i]->GetTotalMomentum()/MeV;
        pp1 = vec[i]->GetMomentum().mag()/MeV;
        if( pp1 < 1.0e-6*GeV )
        {
          G4double costheta = 2.*G4UniformRand() - 1.;
          G4double sintheta = std::sqrt(1. - costheta*costheta);
          G4double phi2 = twopi*G4UniformRand();
          vec[i]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
                               pp*sintheta*std::sin(phi2)*MeV,
                               pp*costheta*MeV ) ;
        }
        else
          vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);

    Rotate( numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
            modifiedOriginal, originalIncident, targetNucleus,
            currentParticle, targetParticle, vec, vecLen );
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    //
    // add black track particles
    // the total number of particles produced is restricted to 198
    // this may have influence on very high energies
    //
    if( atomicWeight >= 1.5 )
    {
      // npnb is number of proton/neutron black track particles
      // ndta is the number of deuterons, tritons, and alphas produced
      // epnb is the kinetic energy available for proton/neutron black track particles
      // edta is the kinetic energy available for deuteron/triton/alpha particles
      //
      G4int npnb = 0;
      G4int ndta = 0;
      
      G4double epnb, edta;
      if (veryForward) {
        epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
        edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
      } else {
        epnb = targetNucleus.GetPNBlackTrackEnergy();
        edta = targetNucleus.GetDTABlackTrackEnergy();
      }

      const G4double pnCutOff = 0.001;
      const G4double dtaCutOff = 0.001;
      const G4double kineticMinimum = 1.e-6;
      const G4double kineticFactor = -0.010;
      G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
      const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
      if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
      if( epnb >= pnCutOff )
      {
        npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
        if( numberofFinalStateNucleons + npnb > atomicWeight )
          npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
        npnb = std::min( npnb, 127-vecLen );
      }
      if( edta >= dtaCutOff )
      {
        ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
        ndta = std::min( ndta, 127-vecLen );
      }
      if (npnb == 0 && ndta == 0) npnb = 1;

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);

      AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum, 
                             kineticFactor, modifiedOriginal,
                             PinNucleus, NinNucleus, targetNucleus,
                             vec, vecLen);

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
    //  MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
    //
    //  calculate time delay for nuclear reactions
    //
    if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
      currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
    else
      currentParticle.SetTOF( 1.0 );
    return true;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  }

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

Definition at line 3777 of file G4ReactionDynamics.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4UniformRand, G4Gamma::Gamma(), GenerateNBodyEvent(), G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalMomentum(), G4FastVector< Type, N >::Initialize(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTotalEnergy(), and G4Triton::Triton().

Referenced by G4LETritonInelastic::ApplyYourself(), G4LEDeuteronInelastic::ApplyYourself(), and G4LEAlphaInelastic::ApplyYourself().

  {
    // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
    //
    // Nuclear reaction kinematics at low energies
    //
    G4ParticleDefinition *aGamma = G4Gamma::Gamma();
    G4ParticleDefinition *aProton = G4Proton::Proton();
    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
    G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
    G4ParticleDefinition *aTriton = G4Triton::Triton();
    G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
    
    const G4double aProtonMass = aProton->GetPDGMass()/MeV;
    const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
    const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
    const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
    const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;

    G4ReactionProduct currentParticle;
    currentParticle = *originalIncident;
    //
    // Set beam particle, take kinetic energy of current particle as the
    // fundamental quantity.  Due to the difficult kinematic, all masses have to
    // be assigned the best measured values
    //
    G4double p = currentParticle.GetTotalMomentum();
    G4double pp = currentParticle.GetMomentum().mag();
    if( pp <= 0.001*MeV )
    {
      G4double phinve = twopi*G4UniformRand();
      G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
      currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
                                   p*std::sin(rthnve)*std::sin(phinve),
                                   p*std::cos(rthnve) );
    }
    else
      currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
    //
    // calculate Q-value of reactions
    //
    G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
    G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
    G4double qv = currentKinetic + theAtomicMass + currentMass;
    
    G4double qval[9];
    qval[0] = qv - mass[0];
    qval[1] = qv - mass[1] - aNeutronMass;
    qval[2] = qv - mass[2] - aProtonMass;
    qval[3] = qv - mass[3] - aDeuteronMass;
    qval[4] = qv - mass[4] - aTritonMass;
    qval[5] = qv - mass[5] - anAlphaMass;
    qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
    qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
    qval[8] = qv - mass[8] - aProtonMass  - aProtonMass;
    
    if( currentParticle.GetDefinition() == aNeutron )
    {
      const G4double A = G4double(targetNucleus.GetA_asInt());    // atomic weight
      if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
        qval[0] = 0.0;
      if( G4UniformRand() >= currentKinetic/7.9254*A )
        qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
    }
    else
      qval[0] = 0.0;
    
    G4int i;
    qv = 0.0;
    for( i=0; i<9; ++i )
    {
      if( mass[i] < 500.0*MeV )qval[i] = 0.0;
      if( qval[i] < 0.0 )qval[i] = 0.0;
      qv += qval[i];
    }
    G4double qv1 = 0.0;
    G4double ran = G4UniformRand();
    G4int index;
    for( index=0; index<9; ++index )
    {
      if( qval[index] > 0.0 )
      {
        qv1 += qval[index]/qv;
        if( ran <= qv1 )break;
      }
    }
    if( index == 9 )  // loop continued to the end
    {
      throw G4HadronicException(__FILE__, __LINE__,
           "G4ReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
    }
    G4double ke = currentParticle.GetKineticEnergy()/GeV;
    G4int nt = 2;
    if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
    
    G4ReactionProduct **v = new G4ReactionProduct * [3];
    v[0] =  new G4ReactionProduct;
    v[1] =  new G4ReactionProduct;
    v[2] =  new G4ReactionProduct;
    
    v[0]->SetMass( mass[index]*MeV );
    switch( index )
    {
     case 0:
       v[1]->SetDefinition( aGamma );
       v[2]->SetDefinition( aGamma );
       break;
     case 1:
       v[1]->SetDefinition( aNeutron );
       v[2]->SetDefinition( aGamma );
       break;
     case 2:
       v[1]->SetDefinition( aProton );
       v[2]->SetDefinition( aGamma );
       break;
     case 3:
       v[1]->SetDefinition( aDeuteron );
       v[2]->SetDefinition( aGamma );
       break;
     case 4:
       v[1]->SetDefinition( aTriton );
       v[2]->SetDefinition( aGamma );
       break;
     case 5:
       v[1]->SetDefinition( anAlpha );
       v[2]->SetDefinition( aGamma );
       break;
     case 6:
       v[1]->SetDefinition( aNeutron );
       v[2]->SetDefinition( aNeutron );
       break;
     case 7:
       v[1]->SetDefinition( aNeutron );
       v[2]->SetDefinition( aProton );
       break;
     case 8:
       v[1]->SetDefinition( aProton );
       v[2]->SetDefinition( aProton );
       break;
    }
    //
    // calculate centre of mass energy
    //
    G4ReactionProduct pseudo1;
    pseudo1.SetMass( theAtomicMass*MeV );
    pseudo1.SetTotalEnergy( theAtomicMass*MeV );
    G4ReactionProduct pseudo2 = currentParticle + pseudo1;
    pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
    //
    // use phase space routine in centre of mass system
    //
    G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
    tempV.Initialize( nt );
    G4int tempLen = 0;
    tempV.SetElement( tempLen++, v[0] );
    tempV.SetElement( tempLen++, v[1] );
    if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
    G4bool constantCrossSection = true;
    GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
    v[0]->Lorentz( *v[0], pseudo2 );
    v[1]->Lorentz( *v[1], pseudo2 );
    if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
    
    G4bool particleIsDefined = false;
    if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
    {
      v[0]->SetDefinition( aProton );
      particleIsDefined = true;
    }
    else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
    {
      v[0]->SetDefinition( aNeutron );
      particleIsDefined = true;
    }
    else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
    {
      v[0]->SetDefinition( aDeuteron );
      particleIsDefined = true;
    }
    else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
    {
      v[0]->SetDefinition( aTriton );
      particleIsDefined = true;
    }
    else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
    {
      v[0]->SetDefinition( anAlpha );
      particleIsDefined = true;
    }
    currentParticle.SetKineticEnergy(
     std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
    p = currentParticle.GetTotalMomentum();
    pp = currentParticle.GetMomentum().mag();
    if( pp <= 0.001*MeV )
    {
      G4double phinve = twopi*G4UniformRand();
      G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
      currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
                                   p*std::sin(rthnve)*std::sin(phinve),
                                   p*std::cos(rthnve) );
    }
    else
      currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
    
    if( particleIsDefined )
    {
      v[0]->SetKineticEnergy(
       std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
      p = v[0]->GetTotalMomentum();
      pp = v[0]->GetMomentum().mag();
      if( pp <= 0.001*MeV )
      {
        G4double phinve = twopi*G4UniformRand();
        G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
        v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
                          p*std::sin(rthnve)*std::sin(phinve),
                          p*std::cos(rthnve) );
      }
      else
        v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
    }
    if( (v[1]->GetDefinition() == aDeuteron) ||
        (v[1]->GetDefinition() == aTriton)   ||
        (v[1]->GetDefinition() == anAlpha) ) 
      v[1]->SetKineticEnergy(
       std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
    else
      v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
    
    p = v[1]->GetTotalMomentum();
    pp = v[1]->GetMomentum().mag();
    if( pp <= 0.001*MeV )
    {
      G4double phinve = twopi*G4UniformRand();
      G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
      v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
                        p*std::sin(rthnve)*std::sin(phinve),
                        p*std::cos(rthnve) );
    }
    else
      v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
    
    if( nt == 3 )
    {
      if( (v[2]->GetDefinition() == aDeuteron) ||
          (v[2]->GetDefinition() == aTriton)   ||
          (v[2]->GetDefinition() == anAlpha) ) 
        v[2]->SetKineticEnergy(
         std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
      else
        v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
      
      p = v[2]->GetTotalMomentum();
      pp = v[2]->GetMomentum().mag();
      if( pp <= 0.001*MeV )
      {
        G4double phinve = twopi*G4UniformRand();
        G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
        v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
                          p*std::sin(rthnve)*std::sin(phinve),
                          p*std::cos(rthnve) );
      }
      else
        v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
    }
    G4int del;
    for(del=0; del<vecLen; del++) delete vec[del];
    vecLen = 0;
    if( particleIsDefined )
    {
      vec.SetElement( vecLen++, v[0] );
    }
    else
    {
      delete v[0];
    }
    vec.SetElement( vecLen++, v[1] );
    if( nt == 3 )
    {
      vec.SetElement( vecLen++, v[2] );
    }
    else
    {
      delete v[2];
    }
    delete [] v;
    return;
  }

void G4ReactionDynamics::ProduceStrangeParticlePairs	(	G4FastVector< G4ReactionProduct, GHADLISTSIZE > &	vec,
		G4int &	vecLen,
		const G4ReactionProduct &	modifiedOriginal,
		const G4DynamicParticle *	originalTarget,
		G4ReactionProduct &	currentParticle,
		G4ReactionProduct &	targetParticle,
		G4bool &	incidentHasChanged,
		G4bool &	targetHasChanged
	)

Definition at line 3386 of file G4ReactionDynamics.cc.

References G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4UniformRand, G4DynamicParticle::GetDefinition(), G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4Neutron::Neutron(), G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetMayBeKilled(), G4ReactionProduct::SetSide(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), and G4SigmaZero::SigmaZero().

Referenced by G4InelasticInteraction::CalculateMomenta().

  {
    // derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
    //
    // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
    //                            K+ Y0, K0 Y+,  K0 Y-
    // For antibaryon induced reactions half of the cross sections KB YB
    // pairs are produced.  Charge is not conserved, no experimental data available
    // for exclusive reactions, therefore some average behaviour assumed.
    // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
    //
    if( vecLen == 0 )return;
    //
    // the following protects against annihilation processes
    //
    if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return;
    
    const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
    const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
    G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
    G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
                                        targetMass*targetMass +
                                        2.0*targetMass*etOriginal );  // GeV
    G4double currentMass = currentParticle.GetMass()/GeV;
    G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
    if( availableEnergy <= 1.0 )return;
    
    G4ParticleDefinition *aProton = G4Proton::Proton();
    G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
    G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
    G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
    G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
    G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
    G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
    G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
    G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
    G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
    G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
    G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
    G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
    G4ParticleDefinition *aLambda = G4Lambda::Lambda();
    G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
    
    const G4double protonMass = aProton->GetPDGMass()/GeV;
    const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
    //
    // determine the center of mass energy bin
    //
    const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};

    G4int ibin, i3, i4;
    G4double avk, avy, avn, ran;
    G4int i = 1;
    while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
    if( i == 12 )
      ibin = 11;
    else
      ibin = i;
    //
    // the fortran code chooses a random replacement of produced kaons
    //  but does not take into account charge conservation 
    //
    if( vecLen == 1 )  // we know that vecLen > 0
    {
      i3 = 0;
      i4 = 1;   // note that we will be adding a new secondary particle in this case only
    }
    else               // otherwise  0 <= i3,i4 < vecLen
    {
      ran = G4UniformRand();
      while( ran == 1.0 )ran = G4UniformRand();
      i4 = i3 = G4int( vecLen*ran );
      while( i3 == i4 )
      {
        ran = G4UniformRand();
        while( ran == 1.0 )ran = G4UniformRand();
        i4 = G4int( vecLen*ran );
      }
    }
    //
    // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
    //
    const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
                               0.0975, 0.123, 0.28,  0.398,  0.495,  0.573 };
    const G4double avky[] = { 0.005, 0.03,  0.064, 0.095, 0.115, 0.13,
                              0.145, 0.155, 0.20,  0.205, 0.210, 0.212 };
    const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
                               0.04,    0.05,   0.12,   0.15,   0.18, 0.20 };
    
    avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
      /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
    avk = std::exp(avk);
    
    avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
      /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
    avy = std::exp(avy);
    
    avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
      /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
    avn = std::exp(avn);
    
    if( avk+avy+avn <= 0.0 )return;
    
    if( currentMass < protonMass )avy /= 2.0;
    if( targetMass < protonMass )avy = 0.0;
    avy += avk+avn;
    avk += avn;
    ran = G4UniformRand();
    if(  ran < avn )
    {
      if( availableEnergy < 2.0 )return;
      if( vecLen == 1 )                              // add a new secondary
      {
        G4ReactionProduct *p1 = new G4ReactionProduct;
        if( G4UniformRand() < 0.5 )
        {
          vec[0]->SetDefinition( aNeutron );
          p1->SetDefinition( anAntiNeutron );
          (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
          vec[0]->SetMayBeKilled(false);
          p1->SetMayBeKilled(false);
        }
        else
        {
          vec[0]->SetDefinition( aProton );
          p1->SetDefinition( anAntiProton );
          (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
          vec[0]->SetMayBeKilled(false);
          p1->SetMayBeKilled(false);
        }
        vec.SetElement( vecLen++, p1 );
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      }
      else
      {                                             // replace two secondaries
        if( G4UniformRand() < 0.5 )
        {
          vec[i3]->SetDefinition( aNeutron );
          vec[i4]->SetDefinition( anAntiNeutron );
          vec[i3]->SetMayBeKilled(false);
          vec[i4]->SetMayBeKilled(false);
        }
        else
        {
          vec[i3]->SetDefinition( aProton );
          vec[i4]->SetDefinition( anAntiProton );
          vec[i3]->SetMayBeKilled(false);
          vec[i4]->SetMayBeKilled(false);
        }
      }
    }
    else if( ran < avk )
    {
      if( availableEnergy < 1.0 )return;
      
      const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
                               0.6875, 0.7500, 0.8750, 1.000 };
      const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
      const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
      ran = G4UniformRand();
      i = 0;
      while( (i<9) && (ran>=kkb[i]) )++i;
      if( i == 9 )return;
      //
      // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
      // charge       +  -   +  0   +  0   0  0   0  0   0  0   0  0   0  -   0  -
      //
      switch( ipakkb1[i] )
      {
       case 10:
         vec[i3]->SetDefinition( aKaonPlus );
         vec[i3]->SetMayBeKilled(false);
         break;
       case 11:
         vec[i3]->SetDefinition( aKaonZS );
         vec[i3]->SetMayBeKilled(false);
         break;
       case 12:
         vec[i3]->SetDefinition( aKaonZL );
         vec[i3]->SetMayBeKilled(false);
         break;
      }
      if( vecLen == 1 )                          // add a secondary
      {
        G4ReactionProduct *p1 = new G4ReactionProduct;
        switch( ipakkb2[i] )
        {
         case 11:
           p1->SetDefinition( aKaonZS );
           p1->SetMayBeKilled(false);
           break;
         case 12:
           p1->SetDefinition( aKaonZL );
           p1->SetMayBeKilled(false);
           break;
         case 13:
           p1->SetDefinition( aKaonMinus );
           p1->SetMayBeKilled(false);
           break;
        }
        (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
        vec.SetElement( vecLen++, p1 );
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      }
      else                                        // replace
      {
        switch( ipakkb2[i] )
        {
         case 11:
           vec[i4]->SetDefinition( aKaonZS );
           vec[i4]->SetMayBeKilled(false);
           break;
         case 12:
           vec[i4]->SetDefinition( aKaonZL );
           vec[i4]->SetMayBeKilled(false);
           break;
         case 13:
           vec[i4]->SetDefinition( aKaonMinus );
           vec[i4]->SetMayBeKilled(false);
           break;
        }
      }
    }
    else if( ran < avy )
    {
      if( availableEnergy < 1.6 )return;
      
      const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
                              0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
      const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
      const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
      const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
      const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
      ran = G4UniformRand();
      i = 0;
      while( (i<12) && (ran>ky[i]) )++i;
      if( i == 12 )return;
      if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
      {
        // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
        //             0  +   0  0   0  0   +  +   +  0   +  0
        //
        //             21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
        //             0  +   0  0   0  0   -  +   -  0   -  0
        switch( ipaky1[i] )
        {
         case 18:
           targetParticle.SetDefinition( aLambda );
           break;
         case 20:
           targetParticle.SetDefinition( aSigmaPlus );
           break;
         case 21:
           targetParticle.SetDefinition( aSigmaZero );
           break;
         case 22:
           targetParticle.SetDefinition( aSigmaMinus );
           break;
        }
        targetHasChanged = true;
        switch( ipaky2[i] )
        {
         case 10:
           vec[i3]->SetDefinition( aKaonPlus ); 
           vec[i3]->SetMayBeKilled(false);
           break;
         case 11:
           vec[i3]->SetDefinition( aKaonZS );
           vec[i3]->SetMayBeKilled(false);
           break;
         case 12:
           vec[i3]->SetDefinition( aKaonZL );
           vec[i3]->SetMayBeKilled(false);
           break;
        }
      }
      else  // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
      {
        // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
        //              24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
        if( (currentParticle.GetDefinition() == anAntiProton) ||
            (currentParticle.GetDefinition() == anAntiNeutron) ||
            (currentParticle.GetDefinition() == anAntiLambda) ||
            (currentMass > sigmaMinusMass) )
        {
          switch( ipakyb1[i] )
          {
           case 19:
             currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
             break;
           case 23:
             currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
             break;
           case 24:
             currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
             break;
           case 25:
             currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
             break;
          }
          incidentHasChanged = true;
          switch( ipakyb2[i] )
          {
           case 11:
             vec[i3]->SetDefinition( aKaonZS ); 
             vec[i3]->SetMayBeKilled(false);
             break;
           case 12:
             vec[i3]->SetDefinition( aKaonZL );
             vec[i3]->SetMayBeKilled(false);
             break;
           case 13:
             vec[i3]->SetDefinition( aKaonMinus );
             vec[i3]->SetMayBeKilled(false);
             break;
          }
        }
        else
        {
          switch( ipaky1[i] )
          {
           case 18:
             currentParticle.SetDefinitionAndUpdateE( aLambda );
             break;
           case 20:
             currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
             break;
           case 21:
             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
             break;
           case 22:
             currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
             break;
          }
          incidentHasChanged = true;
          switch( ipaky2[i] )
          {
           case 10:
             vec[i3]->SetDefinition( aKaonPlus ); 
             vec[i3]->SetMayBeKilled(false);
             break;
           case 11:
             vec[i3]->SetDefinition( aKaonZS );
             vec[i3]->SetMayBeKilled(false);
             break;
           case 12:
             vec[i3]->SetDefinition( aKaonZL );
             vec[i3]->SetMayBeKilled(false);
             break;
          }
        }
      }
    }
    else return;
    //
    //  check the available energy
    //   if there is not enough energy for kkb/ky pair production
    //   then reduce the number of secondary particles 
    //  NOTE:
    //        the number of secondaries may have been changed
    //        the incident and/or target particles may have changed
    //        charge conservation is ignored (as well as strangness conservation)
    //
    currentMass = currentParticle.GetMass()/GeV;
    targetMass = targetParticle.GetMass()/GeV;
    
    G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
    for( i=0; i<vecLen; ++i )
    {
      energyCheck -= vec[i]->GetMass()/GeV;
      if( energyCheck < 0.0 )      // chop off the secondary List
      {
        vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
        G4int j;
        for(j=i; j<vecLen; j++) delete vec[j];
        break;
      }
    }
    return;
  }

void G4ReactionDynamics::SuppressChargedPions	(	G4FastVector< G4ReactionProduct, GHADLISTSIZE > &	vec,
		G4int &	vecLen,
		const G4ReactionProduct &	modifiedOriginal,
		G4ReactionProduct &	currentParticle,
		G4ReactionProduct &	targetParticle,
		const G4Nucleus &	targetNucleus,
		G4bool &	incidentHasChanged,
		G4bool &	targetHasChanged
	)

Definition at line 1339 of file G4ReactionDynamics.cc.

References G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiOmegaMinus::AntiOmegaMinus(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiXiMinus::AntiXiMinus(), G4AntiXiZero::AntiXiZero(), G4UniformRand, G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), and G4ReactionProduct::SetDefinitionAndUpdateE().

Referenced by G4InelasticInteraction::CalculateMomenta().

  {
    // this code was originally in the fortran code TWOCLU
    //
    // suppress charged pions, for various reasons
    //
    G4double mOriginal = modifiedOriginal.GetMass()/GeV;
    G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
    G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
    G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
                                   2.0*targetMass*etOriginal ); 
    G4double eAvailable = cmEnergy - mOriginal - targetMass;
    for (G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;

    const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
    const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
    const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
    
    G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
    G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
    G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
    G4ParticleDefinition *aProton = G4Proton::Proton();
    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
    G4double piMass = aPiPlus->GetPDGMass()/GeV;
    G4double nucleonMass = aNeutron->GetPDGMass()/GeV;
    
    const G4bool antiTest =
      modifiedOriginal.GetDefinition() != G4AntiProton::AntiProton() &&
      modifiedOriginal.GetDefinition() != G4AntiNeutron::AntiNeutron() &&
      modifiedOriginal.GetDefinition() != G4AntiLambda::AntiLambda() &&
      modifiedOriginal.GetDefinition() != G4AntiSigmaPlus::AntiSigmaPlus() &&
      modifiedOriginal.GetDefinition() != G4AntiSigmaMinus::AntiSigmaMinus() &&
      modifiedOriginal.GetDefinition() != G4AntiXiZero::AntiXiZero() &&
      modifiedOriginal.GetDefinition() != G4AntiXiMinus::AntiXiMinus() &&
      modifiedOriginal.GetDefinition() != G4AntiOmegaMinus::AntiOmegaMinus();

    if( antiTest && (
          currentParticle.GetDefinition() == aPiPlus ||
          currentParticle.GetDefinition() == aPiZero ||
          currentParticle.GetDefinition() == aPiMinus ) &&
        ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
        ( G4UniformRand() <= atomicWeight/300.0 ) )
    {
      if (eAvailable > nucleonMass - piMass) {
        if( G4UniformRand() > atomicNumber/atomicWeight )
          currentParticle.SetDefinitionAndUpdateE( aNeutron );
        else
          currentParticle.SetDefinitionAndUpdateE( aProton );
        incidentHasChanged = true;
      }
    }
    if( antiTest && (
          targetParticle.GetDefinition() == aPiPlus ||
          targetParticle.GetDefinition() == aPiZero ||
          targetParticle.GetDefinition() == aPiMinus ) &&
        ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
        ( G4UniformRand() <= atomicWeight/300.0 ) )
    {
      if (eAvailable > nucleonMass - piMass) {
        if( G4UniformRand() > atomicNumber/atomicWeight )
          targetParticle.SetDefinitionAndUpdateE( aNeutron );
        else
          targetParticle.SetDefinitionAndUpdateE( aProton );
        targetHasChanged = true;
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    for( G4int i=0; i<vecLen; ++i )
    {
      if( antiTest && (
            vec[i]->GetDefinition() == aPiPlus ||
            vec[i]->GetDefinition() == aPiZero ||
            vec[i]->GetDefinition() == aPiMinus ) &&
          ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
          ( G4UniformRand() <= atomicWeight/300.0 ) )
      {
        if (eAvailable > nucleonMass - piMass) {
          if( G4UniformRand() > atomicNumber/atomicWeight )
            vec[i]->SetDefinitionAndUpdateE( aNeutron );
          else
            vec[i]->SetDefinitionAndUpdateE( aProton );
        }
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  }

void G4ReactionDynamics::TwoBody	(	G4FastVector< G4ReactionProduct, GHADLISTSIZE > &	vec,
		G4int &	vecLen,
		G4ReactionProduct &	modifiedOriginal,
		const G4DynamicParticle *	originalTarget,
		G4ReactionProduct &	currentParticle,
		G4ReactionProduct &	targetParticle,
		const G4Nucleus &	targetNucleus,
		G4bool &	targetHasChanged
	)

Definition at line 2243 of file G4ReactionDynamics.cc.

References G4UniformRand, G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4Nucleus::GetN_asInt(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4ReactionProduct::Lorentz(), G4InuclParticleNames::pp, G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTOF(), and G4ReactionProduct::SetTotalEnergy().

Referenced by G4InelasticInteraction::CalculateMomenta().

  {
    // 
    // derived from original FORTRAN code TWOB by H. Fesefeldt (15-Sep-1987)
    //
    // Generation of momenta for elastic and quasi-elastic 2 body reactions
    //
    // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
    // The b values are parametrizations from experimental data.
    // Not available values are taken from those of similar reactions.
    //
    
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);    
    static const G4double expxu =  82.;           // upper bound for arg. of exp
    static const G4double expxl = -expxu;         // lower bound for arg. of exp
    
    const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
    const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
    const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
    const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
    G4double currentMass = currentParticle.GetMass()/GeV;
    G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;

    targetMass = targetParticle.GetMass()/GeV;
    const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
    
    G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
    G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;

    G4double cmEnergy = std::sqrt( currentMass*currentMass +
                              targetMass*targetMass +
                              2.0*targetMass*etCurrent );  // in GeV

    //if( (pOriginal < 0.1) ||
    //    (centerofmassEnergy < 0.01) ) // 2-body scattering not possible
    // Continue with original particle, but spend the nuclear evaporation energy
    //  targetParticle.SetMass( 0.0 );  // flag that the target doesn't exist
    //else                           // Two-body scattering is possible

    if( (pCurrent < 0.1) || (cmEnergy < 0.01) ) // 2-body scattering not possible
    {
      targetParticle.SetMass( 0.0 );  // flag that the target particle doesn't exist
    }
    else
    {
// moved this if-block to a later stage, i.e. to the assignment of the scattering angle
// @@@@@ double-check.
//      if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" || 
//          targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
//        if( G4UniformRand() < 0.5 )
//          targetParticle.SetDefinitionAndUpdateE( aNeutron );
//        else
//          targetParticle.SetDefinitionAndUpdateE( aProton );
//        targetHasChanged = true;
//        targetMass = targetParticle.GetMass()/GeV;
//      }
      //
      // Set masses and momenta for final state particles
      //
      G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
      pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
      
      if( pf < 0.001 )
      {
        for(G4int i=0; i<vecLen; i++) delete vec[i];
        vecLen = 0;
        throw G4HadronicException(__FILE__, __LINE__, "G4ReactionDynamics::TwoBody: pf is too small ");
      }
      
      pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
      //
      // Set beam and target in centre of mass system
      //
      G4ReactionProduct pseudoParticle[3];
      
      if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
          targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
        pseudoParticle[0].SetMass( targetMass*GeV );
        pseudoParticle[0].SetTotalEnergy( etOriginal*GeV );
        pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
      
        pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
        pseudoParticle[1].SetMass( mOriginal*GeV );
        pseudoParticle[1].SetKineticEnergy( 0.0 );

      } else {
        pseudoParticle[0].SetMass( currentMass*GeV );
        pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
        pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
      
        pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
        pseudoParticle[1].SetMass( targetMass*GeV );
        pseudoParticle[1].SetKineticEnergy( 0.0 );
      }
      //
      // Transform into centre of mass system
      //
      pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
      pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
      pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
      //
      // Set final state masses and energies in centre of mass system
      //
      currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
      targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
      //
      // Set |t| and |tmin|
      //
      const G4double cb = 0.01;
      const G4double b1 = 4.225;
      const G4double b2 = 1.795;
      //
      // Calculate slope b for elastic scattering on proton/neutron
      //
      G4double b = std::max( cb, b1+b2*std::log(pOriginal) );     
      G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
      
      G4double exindt = -1.0;
      exindt += std::exp(std::max(-btrang,expxl));
      //
      // Calculate sqr(std::sin(teta/2.) and std::cos(teta), set azimuth angle phi
      //
      G4double ctet = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) / btrang;
      if( std::fabs(ctet) > 1.0 )ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
      G4double stet = std::sqrt( (1.0-ctet)*(1.0+ctet) );
      G4double phi = twopi * G4UniformRand();
      //
      // Calculate final state momenta in centre of mass system
      //
      if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
          targetParticle.GetDefinition()->GetParticleSubType() == "pi") {

        currentParticle.SetMomentum( -pf*stet*std::sin(phi)*GeV,
                                     -pf*stet*std::cos(phi)*GeV,
                                     -pf*ctet*GeV );
      } else {

        currentParticle.SetMomentum( pf*stet*std::sin(phi)*GeV,
                                     pf*stet*std::cos(phi)*GeV,
                                     pf*ctet*GeV );
      }
      targetParticle.SetMomentum( currentParticle.GetMomentum() * (-1.0) );
      //
      // Transform into lab system
      //
      currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
      targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
      
      Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
      
      G4double pp, pp1, ekin;
      if( atomicWeight >= 1.5 )
      {
        const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
        pp1 = currentParticle.GetMomentum().mag()/MeV;
        if( pp1 >= 1.0 )
        {
          ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
          ekin = std::max( 0.0001*GeV, ekin );
          currentParticle.SetKineticEnergy( ekin*MeV );
          pp = currentParticle.GetTotalMomentum()/MeV;
          currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
        }
        pp1 = targetParticle.GetMomentum().mag()/MeV;
        if( pp1 >= 1.0 )
        {
          ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
          ekin = std::max( 0.0001*GeV, ekin );
          targetParticle.SetKineticEnergy( ekin*MeV );
          pp = targetParticle.GetTotalMomentum()/MeV;
          targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
        }
      }
    }

    // Get number of final state nucleons and nucleons remaining in
    // target nucleus
    
    std::pair<G4int, G4int> finalStateNucleons = 
      GetFinalStateNucleons(originalTarget, vec, vecLen);
    G4int protonsInFinalState = finalStateNucleons.first;
    G4int neutronsInFinalState = finalStateNucleons.second;

    G4int PinNucleus = std::max(0, 
      targetNucleus.GetZ_asInt() - protonsInFinalState);
    G4int NinNucleus = std::max(0,
      targetNucleus.GetN_asInt() - neutronsInFinalState);

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    if( atomicWeight >= 1.5 )
    {
      // Add black track particles
      // npnb is number of proton/neutron black track particles
      // ndta is the number of deuterons, tritons, and alphas produced
      // epnb is the kinetic energy available for proton/neutron black track particles
      // edta is the kinetic energy available for deuteron/triton/alpha particles
      //
      G4double epnb, edta;
      G4int npnb=0, ndta=0;
      
      epnb = targetNucleus.GetPNBlackTrackEnergy();            // was enp1 in fortran code
      edta = targetNucleus.GetDTABlackTrackEnergy();           // was enp3 in fortran code
      const G4double pnCutOff = 0.0001;       // GeV
      const G4double dtaCutOff = 0.0001;      // GeV
      const G4double kineticMinimum = 0.0001;
      const G4double kineticFactor = -0.010;
      G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
      if( epnb >= pnCutOff )
      {
        npnb = Poisson( epnb/0.02 );
        if( npnb > atomicWeight )npnb = G4int(atomicWeight);
        if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
        npnb = std::min( npnb, 127-vecLen );
      }
      if( edta >= dtaCutOff )
      {
        ndta = G4int(2.0 * std::log(atomicWeight));
        ndta = std::min( ndta, 127-vecLen );
      }

      if (npnb == 0 && ndta == 0) npnb = 1;

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);

      AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum, 
                             kineticFactor, modifiedOriginal, 
                             PinNucleus, NinNucleus, targetNucleus,
                             vec, vecLen);
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    //
    //  calculate time delay for nuclear reactions
    //
    if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
      currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
    else
      currentParticle.SetTOF( 1.0 );
    return;
  }

G4bool G4ReactionDynamics::TwoCluster	(	G4FastVector< G4ReactionProduct, GHADLISTSIZE > &	vec,
		G4int &	vecLen,
		G4ReactionProduct &	modifiedOriginal,
		const G4HadProjectile *	originalIncident,
		G4ReactionProduct &	currentParticle,
		G4ReactionProduct &	targetParticle,
		const G4DynamicParticle *	originalTarget,
		const G4Nucleus &	targetNucleus,
		G4bool &	incidentHasChanged,
		G4bool &	targetHasChanged,
		G4bool	leadFlag,
		G4ReactionProduct &	leadingStrangeParticle
	)

Definition at line 1434 of file G4ReactionDynamics.cc.

References G4UniformRand, GenerateNBodyEvent(), G4Nucleus::GetA_asInt(), G4Nucleus::GetAnnihilationDTABlackTrackEnergy(), G4Nucleus::GetAnnihilationPNBlackTrackEnergy(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4Nucleus::GetN_asInt(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetSide(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4FastVector< Type, N >::Initialize(), G4Lambda::Lambda(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4InuclParticleNames::pp, G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetSide(), and G4ReactionProduct::SetTotalEnergy().

Referenced by G4InelasticInteraction::CalculateMomenta().

  {
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    // derived from original FORTRAN code TWOCLU by H. Fesefeldt (11-Oct-1987)
    //
    //  Generation of X- and PT- values for incident, target, and all secondary particles 
    //  
    //  A simple two cluster model is used.
    //  This should be sufficient for low energy interactions.
    //

    // + debugging
    // raise(SIGSEGV);
    // - debugging

    G4int i;
    G4ParticleDefinition *aProton = G4Proton::Proton();
    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
    G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
    G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
    G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
    G4bool veryForward = false;

    const G4double protonMass = aProton->GetPDGMass()/MeV;
    const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
    const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
    const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
    const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
    G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
    G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
                                        targetMass*targetMass +
                                        2.0*targetMass*etOriginal );  // GeV
    G4double currentMass = currentParticle.GetMass()/GeV;
    targetMass = targetParticle.GetMass()/GeV;

    if( currentMass == 0.0 && targetMass == 0.0 )
    {
      G4double ek = currentParticle.GetKineticEnergy();
      G4ThreeVector mom = currentParticle.GetMomentum();
      currentParticle = *vec[0];
      targetParticle = *vec[1];
      for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
      if(vecLen<2) 
      {
        for(G4int j=0; j<vecLen; j++) delete vec[j];
        vecLen = 0;
        throw G4HadReentrentException(__FILE__, __LINE__,
        "G4ReactionDynamics::TwoCluster: Negative number of particles");
      }
      delete vec[vecLen-1];
      delete vec[vecLen-2];
      vecLen -= 2;
      currentMass = currentParticle.GetMass()/GeV;
      targetMass = targetParticle.GetMass()/GeV;
      incidentHasChanged = true;
      targetHasChanged = true;
      currentParticle.SetKineticEnergy( ek );
      currentParticle.SetMomentum( mom );
      veryForward = true;
    }

    const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
    const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
    //
    // particles have been distributed in forward and backward hemispheres
    // in center of mass system of the hadron nucleon interaction
    //
    // incident is always in forward hemisphere
    G4int forwardCount = 1;            // number of particles in forward hemisphere
    currentParticle.SetSide( 1 );
    G4double forwardMass = currentParticle.GetMass()/GeV;
    G4double cMass = forwardMass;
    
    // target is always in backward hemisphere
    G4int backwardCount = 1;           // number of particles in backward hemisphere
    G4int backwardNucleonCount = 1;    // number of nucleons in backward hemisphere
    targetParticle.SetSide( -1 );
    G4double backwardMass = targetParticle.GetMass()/GeV;
    G4double bMass = backwardMass;
    
    for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 ); // added by JLC, 2Jul97
      // to take care of the case where vec has been preprocessed by GenerateXandPt
      // and some of them have been set to -2 or -3
      if( vec[i]->GetSide() == -1 )
      {
        ++backwardCount;
        backwardMass += vec[i]->GetMass()/GeV;
      }
      else
      {
        ++forwardCount;
        forwardMass += vec[i]->GetMass()/GeV;
      }
    }
    //
    // nucleons and some pions from intranuclear cascade
    //
    G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
    if(term1 < 0) term1 = 0.0001; // making sure xtarg<0;
    const G4double afc = 0.312 + 0.2 * std::log(term1);
    G4double xtarg;
    if( centerofmassEnergy < 2.0+G4UniformRand() )        // added +2 below, JLC 4Jul97
      xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
    else
      xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
    if( xtarg <= 0.0 )xtarg = 0.01;
    G4int nuclearExcitationCount = Poisson( xtarg );
    if(atomicWeight<1.0001) nuclearExcitationCount = 0;
    G4int extraNucleonCount = 0;
    G4double extraMass = 0.0;
    G4double extraNucleonMass = 0.0;
    if( nuclearExcitationCount > 0 )
    {
      G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );     
      const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
      //
      //  NOTE: in TWOCLU, these new particles were given negative codes
      //        here we use  NewlyAdded = true  instead
      //
//      G4ReactionProduct *pVec = new G4ReactionProduct [nuclearExcitationCount];
      for( i=0; i<nuclearExcitationCount; ++i )
      {
        G4ReactionProduct* pVec = new G4ReactionProduct();
        if( G4UniformRand() < nucsup[momentumBin] )  // add proton or neutron
        {
          if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
            pVec->SetDefinition( aProton );
          else
            pVec->SetDefinition( aNeutron );
          ++backwardNucleonCount;
          ++extraNucleonCount;
          extraNucleonMass += pVec->GetMass()/GeV;
        }
        else
        {                                       // add a pion
          G4double ran = G4UniformRand();
          if( ran < 0.3181 )
            pVec->SetDefinition( aPiPlus );
          else if( ran < 0.6819 )
            pVec->SetDefinition( aPiZero );
          else
            pVec->SetDefinition( aPiMinus );
        }
        pVec->SetSide( -2 );    // backside particle
        extraMass += pVec->GetMass()/GeV;
        pVec->SetNewlyAdded( true );
        vec.SetElement( vecLen++, pVec );
        // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
    G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
    G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
    G4bool secondaryDeleted;
    G4double pMass;

    while( eAvailable <= 0.0 )   // must eliminate a particle
    {
      secondaryDeleted = false;
      for( i=(vecLen-1); i>=0; --i )
      {
        if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
        {
          pMass = vec[i]->GetMass()/GeV;
          for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];     // shift up
          --forwardCount;
          forwardEnergy += pMass;
          forwardMass -= pMass;
          secondaryDeleted = true;
          break;
        }
        else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
        {
          pMass = vec[i]->GetMass()/GeV;
          for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
          --backwardCount;
          backwardEnergy += pMass;
          backwardMass -= pMass;
          secondaryDeleted = true;
          break;
        }
      }           // breaks go down to here
      if( secondaryDeleted )
      {      
        delete vec[vecLen-1];
        --vecLen;
        // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      }
      else
      {
        if( vecLen == 0 )
        {
          return false;  // all secondaries have been eliminated
        }
        if( targetParticle.GetSide() == -1 )
        {
          pMass = targetParticle.GetMass()/GeV;
          targetParticle = *vec[0];
          for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
          --backwardCount;
          backwardEnergy += pMass;
          backwardMass -= pMass;
          secondaryDeleted = true;
        }
        else if( targetParticle.GetSide() == 1 )
        {
          pMass = targetParticle.GetMass()/GeV;
          targetParticle = *vec[0];
          for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
          --forwardCount;
          forwardEnergy += pMass;
          forwardMass -= pMass;
          secondaryDeleted = true;
        }
        if( secondaryDeleted )
        {
          delete vec[vecLen-1];
          --vecLen;
          // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
        else
        {
          if( currentParticle.GetSide() == -1 )
          {
            pMass = currentParticle.GetMass()/GeV;
            currentParticle = *vec[0];
            for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
            --backwardCount;
            backwardEnergy += pMass;
            backwardMass -= pMass;
            secondaryDeleted = true;
          }
          else if( currentParticle.GetSide() == 1 )
          {
            pMass = currentParticle.GetMass()/GeV;
            currentParticle = *vec[0];
            for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
            --forwardCount;
            forwardEnergy += pMass;
            forwardMass -= pMass;
            secondaryDeleted = true;
          }
          if( secondaryDeleted )
          {
            delete vec[vecLen-1];
            --vecLen;
            // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
          }
          else break;
        }
      }
      eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
    }
    //
    // This is the start of the TwoCluster function
    //  Choose masses for the 3 clusters:
    //   forward cluster
    //   backward meson cluster
    //   backward nucleon cluster
    //
    G4double rmc = 0.0, rmd = 0.0;
    const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
    const G4double gpar[] = { 2.6, 2.6, 1.8, 1.30, 1.20 };
    
    if (forwardCount <= 0 || backwardCount <= 0) return false;  // array bounds protection

    if (forwardCount == 1) rmc = forwardMass;
    else 
    {
      G4int ntc = std::max(1, std::min(5,forwardCount))-1; // check if offset by 1 @@
      rmc = forwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
    }

    if (backwardCount == 1) rmd = backwardMass;
    else
    {
      G4int ntc = std::max(1, std::min(5,backwardCount)); // check, if offfset by 1 @@
      rmd = backwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
    }

    while( rmc+rmd > centerofmassEnergy )
    {
      if( (rmc <= forwardMass) && (rmd <= backwardMass) )
      {
        G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
        rmc *= temp;
        rmd *= temp;
      }
      else
      {
        rmc = 0.1*forwardMass + 0.9*rmc;
        rmd = 0.1*backwardMass + 0.9*rmd;
      }
    }

    G4ReactionProduct pseudoParticle[8];
    for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
    
    pseudoParticle[1].SetMass( mOriginal*GeV );
    pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
    pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
    
    pseudoParticle[2].SetMass( protonMass*MeV );
    pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
    pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
    //
    //  transform into centre of mass system
    //
    pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
    pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
    pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
    
    const G4double pfMin = 0.0001;
    G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
    pf *= pf;
    pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
    pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
    //
    //  set final state masses and energies in centre of mass system
    //
    pseudoParticle[3].SetMass( rmc*GeV );
    pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
    
    pseudoParticle[4].SetMass( rmd*GeV );
    pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
    //
    // set |T| and |TMIN|
    //
    const G4double bMin = 0.01;
    const G4double b1 = 4.0;
    const G4double b2 = 1.6;

    G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
    G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
    G4double factor = 1.0 - std::exp(-dtb);
    G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
    costheta = std::max(-1.0, std::min(1.0, costheta));
    G4double sintheta = std::sqrt(1.0-costheta*costheta);
    G4double phi = G4UniformRand() * twopi;

    //
    // calculate final state momenta in centre of mass system
    //
    pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
                                   pf*sintheta*std::sin(phi)*GeV,
                                   pf*costheta*GeV );

    pseudoParticle[4].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
    //
    // simulate backward nucleon cluster in lab. system and transform in cms
    //
    G4double pp, pp1;
    if( nuclearExcitationCount > 0 )
    {
      const G4double ga = 1.2;
      G4double ekit1 = 0.04;
      G4double ekit2 = 0.6;
      if( ekOriginal <= 5.0 )
      {
        ekit1 *= ekOriginal*ekOriginal/25.0;
        ekit2 *= ekOriginal*ekOriginal/25.0;
      }
      const G4double a = (1.0-ga)/(std::pow(ekit2,(1.0-ga)) - std::pow(ekit1,(1.0-ga)));
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == -2 )
        {
          G4double kineticE =
            std::pow( (G4UniformRand()*(1.0-ga)/a+std::pow(ekit1,(1.0-ga))), (1.0/(1.0-ga)) );
          vec[i]->SetKineticEnergy( kineticE*GeV );
          G4double vMass = vec[i]->GetMass()/MeV;
          G4double totalE = kineticE*GeV + vMass;
          pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
          G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
          G4double sint = std::sqrt( std::max( 0.0, (1.0-cost*cost) ) );
          phi = twopi*G4UniformRand();
          vec[i]->SetMomentum( pp*sint*std::sin(phi)*MeV,
                               pp*sint*std::cos(phi)*MeV,
                               pp*cost*MeV );
          vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
        }
      }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    //
    // fragmentation of forward cluster and backward meson cluster
    //
    currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
    currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
    
    targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
    targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
    
    pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
    pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
    pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
    
    pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
    pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
    pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
    
    G4double wgt;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    if( forwardCount > 1 )     // tempV will contain the forward particles
    {
      G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
      tempV.Initialize( forwardCount );
      G4bool constantCrossSection = true;
      G4int tempLen = 0;
      if( currentParticle.GetSide() == 1 )
        tempV.SetElement( tempLen++, &currentParticle );
      if( targetParticle.GetSide() == 1 )
        tempV.SetElement( tempLen++, &targetParticle );
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == 1 )
        {
          if( tempLen < 18 )
            tempV.SetElement( tempLen++, vec[i] );
          else
          {
            vec[i]->SetSide( -1 );
            continue;
          }
        }
      }
      if( tempLen >= 2 )
      {
        wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
                                  constantCrossSection, tempV, tempLen );
        if( currentParticle.GetSide() == 1 )
          currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
        if( targetParticle.GetSide() == 1 )
          targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
        for( i=0; i<vecLen; ++i )
        {
          if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
        }
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    if( backwardCount > 1 )   //  tempV will contain the backward particles,
    {                         //  but not those created from the intranuclear cascade
      G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
      tempV.Initialize( backwardCount );
      G4bool constantCrossSection = true;
      G4int tempLen = 0;
      if( currentParticle.GetSide() == -1 )
        tempV.SetElement( tempLen++, &currentParticle );
      if( targetParticle.GetSide() == -1 )
        tempV.SetElement( tempLen++, &targetParticle );
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == -1 )
        {
          if( tempLen < 18 )
            tempV.SetElement( tempLen++, vec[i] );
          else
          {
            vec[i]->SetSide( -2 );
            vec[i]->SetKineticEnergy( 0.0 );
            vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
            continue;
          }
        }
      }
      if( tempLen >= 2 )
      {
        wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
                                  constantCrossSection, tempV, tempLen );
        if( currentParticle.GetSide() == -1 )
          currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
        if( targetParticle.GetSide() == -1 )
          targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
        for( i=0; i<vecLen; ++i )
        {
          if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
        }
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    //
    // Lorentz transformation in lab system
    //
    currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
    targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
    for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    //
    // sometimes the leading strange particle is lost, set it back
    //
    G4bool dum = true;
    if( leadFlag )
    {
      // leadFlag will be true
      //  iff original particle is at least as heavy as K+ and not a proton or 
      //  neutron AND if incident particle is at least as heavy as K+ and it is
      //  not a proton or neutron leadFlag is set to the incident particle
      //  or
      //  target particle is at least as heavy as K+ and it is not a proton or 
      //  neutron leadFlag is set to the target particle
      //
      if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
        dum = false;
      else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
        dum = false;
      else
      {
        for( i=0; i<vecLen; ++i )
        {
          if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
          {
            dum = false;
            break;
          }
        }
      }
      if( dum )
      {
        G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
        G4double ekin;
        if( ((leadMass <  protonMass) && (targetParticle.GetMass()/MeV <  protonMass)) ||
            ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
        {
          ekin = targetParticle.GetKineticEnergy()/GeV;
          pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
          targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
          targetParticle.SetKineticEnergy( ekin*GeV );
          pp = targetParticle.GetTotalMomentum()/MeV;   // new momentum
          if( pp1 < 1.0e-3 )
          {
            G4double costheta2 = 2.*G4UniformRand() - 1.;
            G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
            G4double phi2 = twopi*G4UniformRand();
            targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
                                        pp*sintheta2*std::sin(phi2)*MeV,
                                        pp*costheta2*MeV ) ;
          }
          else
            targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
          
          targetHasChanged = true;
        }
        else
        {
          ekin = currentParticle.GetKineticEnergy()/GeV;
          pp1 = currentParticle.GetMomentum().mag()/MeV;
          currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
          currentParticle.SetKineticEnergy( ekin*GeV );
          pp = currentParticle.GetTotalMomentum()/MeV;
          if( pp1 < 1.0e-3 )
          {
            G4double costheta2 = 2.*G4UniformRand() - 1.;
            G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
            G4double phi2 = twopi*G4UniformRand();
            currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
                                         pp*sintheta2*std::sin(phi2)*MeV,
                                         pp*costheta2*MeV ) ;
          }
          else
            currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
          
          incidentHasChanged = true;
        }
      }
    }    // end of if( leadFlag )

    // Get number of final state nucleons and nucleons remaining in
    // target nucleus
    
    std::pair<G4int, G4int> finalStateNucleons = 
      GetFinalStateNucleons(originalTarget, vec, vecLen);

    G4int protonsInFinalState = finalStateNucleons.first;
    G4int neutronsInFinalState = finalStateNucleons.second;

    G4int numberofFinalStateNucleons = 
      protonsInFinalState + neutronsInFinalState;

    if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
        targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
        originalIncident->GetDefinition()->GetPDGMass() < 
                                   G4Lambda::Lambda()->GetPDGMass())
      numberofFinalStateNucleons++;

    numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);

    G4int PinNucleus = std::max(0, 
      targetNucleus.GetZ_asInt() - protonsInFinalState);
    G4int NinNucleus = std::max(0,
      targetNucleus.GetN_asInt() - neutronsInFinalState);
    //
    //  for various reasons, the energy balance is not sufficient,
    //  check that,  energy balance, angle of final system, etc.
    //
    pseudoParticle[4].SetMass( mOriginal*GeV );
    pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
    pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
    
    G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
    G4int diff = 0;
    if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() )  diff = 1;
    if(numberofFinalStateNucleons == 1) diff = 0;
    pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
    pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
    pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
    
    G4double theoreticalKinetic =
      pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
    
    pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
    pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
    pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
     
    if( vecLen < 16 )
    {
      G4ReactionProduct tempR[130];
      tempR[0] = currentParticle;
      tempR[1] = targetParticle;
      for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];

      G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
      tempV.Initialize( vecLen+2 );
      G4bool constantCrossSection = true;
      G4int tempLen = 0;
      for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );

      if( tempLen >= 2 )
      {
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
                                  pseudoParticle[5].GetTotalEnergy()/MeV,
                                  constantCrossSection, tempV, tempLen );
        if (wgt == -1) {
          G4double Qvalue = 0;
          for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
          wgt = GenerateNBodyEvent( Qvalue/MeV,
                                    constantCrossSection, tempV, tempLen );
        }
        theoreticalKinetic = 0.0;
        for( i=0; i<vecLen+2; ++i )
        {
          pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
          pseudoParticle[7].SetMass( tempV[i]->GetMass() );
          pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
          pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
          theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
        }
      }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    else
    {
      theoreticalKinetic -=
        ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
      for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
    }
    G4double simulatedKinetic =
      currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
    for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
    //
    // make sure that kinetic energies are correct
    // the backward nucleon cluster is not produced within proper kinematics!!!
    //
    
    if( simulatedKinetic != 0.0 )
    {
      wgt = (theoreticalKinetic)/simulatedKinetic;
      currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
      pp = currentParticle.GetTotalMomentum()/MeV;
      pp1 = currentParticle.GetMomentum().mag()/MeV;
      if( pp1 < 0.001*MeV )
      {
        G4double costheta2 = 2.*G4UniformRand() - 1.;
        G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
        G4double phi2 = twopi*G4UniformRand();
        currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
                                     pp*sintheta2*std::sin(phi2)*MeV,
                                     pp*costheta2*MeV ) ;
      }
      else
        currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
      
      targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
      pp = targetParticle.GetTotalMomentum()/MeV;
      pp1 = targetParticle.GetMomentum().mag()/MeV;
      if( pp1 < 0.001*MeV )
      {
        G4double costheta2 = 2.*G4UniformRand() - 1.;
        G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
        G4double phi2 = twopi*G4UniformRand();
        targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
                                    pp*sintheta2*std::sin(phi2)*MeV,
                                    pp*costheta2*MeV ) ;
      }
      else
        targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
      
      for( i=0; i<vecLen; ++i )
      {
        vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
        pp = vec[i]->GetTotalMomentum()/MeV;
        pp1 = vec[i]->GetMomentum().mag()/MeV;
        if( pp1 < 0.001 )
        {
          G4double costheta2 = 2.*G4UniformRand() - 1.;
          G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
          G4double phi2 = twopi*G4UniformRand();
          vec[i]->SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
                               pp*sintheta2*std::sin(phi2)*MeV,
                               pp*costheta2*MeV ) ;
        }
        else
          vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);

    Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
            modifiedOriginal, originalIncident, targetNucleus,
            currentParticle, targetParticle, vec, vecLen );
    //
    //  add black track particles
    //  the total number of particles produced is restricted to 198
    //  this may have influence on very high energies
    //
    if( atomicWeight >= 1.5 )
    {
      // npnb is number of proton/neutron black track particles
      // ndta is the number of deuterons, tritons, and alphas produced
      // epnb is the kinetic energy available for proton/neutron black track 
      //   particles
      // edta is the kinetic energy available for deuteron/triton/alpha 
      //   particles

      G4int npnb = 0;
      G4int ndta = 0;

      G4double epnb, edta;
      if (veryForward) {
        epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
        edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
      } else {
        epnb = targetNucleus.GetPNBlackTrackEnergy();
        edta = targetNucleus.GetDTABlackTrackEnergy();
      }

      const G4double pnCutOff = 0.001;     // GeV
      const G4double dtaCutOff = 0.001;    // GeV
      const G4double kineticMinimum = 1.e-6;
      const G4double kineticFactor = -0.005;
      
      G4double sprob = 0.0;   // sprob = probability of self-absorption in 
                              // heavy molecules
      const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
      if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
      
      if( epnb >= pnCutOff )
      {
        npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
        if( numberofFinalStateNucleons + npnb > atomicWeight )
          npnb = G4int(atomicWeight - numberofFinalStateNucleons);
        npnb = std::min( npnb, 127-vecLen );
      }
      if( edta >= dtaCutOff )
      {
        ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
        ndta = std::min( ndta, 127-vecLen );
      }
      if (npnb == 0 && ndta == 0) npnb = 1;

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);

      AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum, 
                             kineticFactor, modifiedOriginal, 
                             PinNucleus, NinNucleus, targetNucleus,
                             vec, vecLen );
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
    //  MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
    //
    //  calculate time delay for nuclear reactions
    //
    if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
      currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
    else
      currentParticle.SetTOF( 1.0 );
    
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    return true;
  }

---

The documentation for this class was generated from the following files:

• G4ReactionDynamics.hh
• G4ReactionDynamics.cc

---