FullModelReactionDynamics.cxx

Go to the documentation of this file.

//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in:                      *
// *   http://cern.ch/geant4/license                                  *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.                                                             *
// *                                                                  *
// * This  code  implementation is the  intellectual property  of the *
// * GEANT4 collaboration.                                            *
// * By copying,  distributing  or modifying the Program (or any work *
// * based  on  the Program)  you indicate  your  acceptance of  this *
// * statement, and all its terms.                                    *
// ********************************************************************
//
//
//
// Hadronic Process: Reaction Dynamics
// original by H.P. Wellisch
// modified by J.L. Chuma, TRIUMF, 19-Nov-1996
// Last modified: 27-Mar-1997
// modified by H.P. Wellisch, 24-Apr-97
// H.P. Wellisch, 25.Apr-97: Side of current and target particle taken into account
// H.P. Wellisch, 29.Apr-97: Bug fix in NuclearReaction. (pseudo1 was without energy)
// J.L. Chuma, 30-Apr-97:  Changed return value for GenerateXandPt.  It seems possible
//                         that GenerateXandPt could eliminate some secondaries, but
//                         still finish its calculations, thus we would not want to
//                         then use TwoCluster to again calculate the momenta if vecLen
//                         was less than 6.
// J.L. Chuma, 10-Jun-97:  Modified NuclearReaction.  Was not creating ReactionProduct's
//                         with the new operator, thus they would be meaningless when
//                         going out of scope.
// J.L. Chuma, 20-Jun-97:  Modified GenerateXandPt and TwoCluster to fix units problems
// J.L. Chuma, 23-Jun-97:  Modified ProduceStrangeParticlePairs to fix units problems
// J.L. Chuma, 26-Jun-97:  Modified ProduceStrangeParticlePairs to fix array indices
//                         which were sometimes going out of bounds
// J.L. Chuma, 04-Jul-97:  Many minor modifications to GenerateXandPt and TwoCluster
// J.L. Chuma, 06-Aug-97:  Added original incident particle, before Fermi motion and
//                         evaporation effects are included, needed for self absorption
//                         and corrections for single particle spectra (shower particles)
// logging stopped 1997
// J. Allison, 17-Jun-99:  Replaced a min function to get correct behaviour on DEC.
 
#include "FullModelReactionDynamics.hh"
#include "G4Version.hh"
#include "G4AntiProton.hh"
#include "G4AntiNeutron.hh"
#include "Randomize.hh"
#include <iostream>
#if G4VERSION_NUMBER < 1100
#include "G4HadReentrentException.hh"
#else
#include "G4HadronicException.hh"
#endif
#include <signal.h>
//#include "G4ParticleTable.hh"
 
// #include "DumpFrame.hh"
 
/*         G4double GetQValue(G4ReactionProduct * aSec)
           {
           double QValue=0;
           if(aSec->GetDefinition()->GetParticleType() == "baryon")
           {
           if(aSec->GetDefinition()->GetBaryonNumber() < 0)
           {
           QValue = aSec->GetTotalEnergy();
           QValue += G4Neutron::Neutron()->GetPDGMass();
           }
           else
           {
           G4double ss = 0;
           ss +=aSec->GetDefinition()->GetPDGMass();
           if(aSec->GetDefinition() == G4Proton::Proton())
           {
           ss -=G4Proton::Proton()->GetPDGMass();
           }
           else
           {
           ss -=G4Neutron::Neutron()->GetPDGMass();
           }
           ss += aSec->GetKineticEnergy();
           QValue = ss;
           }
           }
           else if(aSec->GetDefinition()->GetPDGEncoding() == 0)
           {
           QValue = aSec->GetKineticEnergy();
           }
           else
           {
           QValue = aSec->GetTotalEnergy();
           }
           return QValue;
           }
*/
 
G4bool FullModelReactionDynamics::GenerateXandPt(
                                                 G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                                 G4int &vecLen,
                                                 G4ReactionProduct &modifiedOriginal,   // Fermi motion & evap. effects included
                                                 const G4HadProjectile *originalIncident,   // the original incident particle
                                                 G4ReactionProduct &currentParticle,
                                                 G4ReactionProduct &targetParticle,
                                                 const G4Nucleus &targetNucleus,
                                                 G4bool &incidentHasChanged,
                                                 G4bool &targetHasChanged,
                                                 G4bool leadFlag,
                                                 G4ReactionProduct &leadingStrangeParticle )
{
  //
  // 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 *anAntiProton = G4AntiProton::AntiProton();
  // G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
  G4ParticleDefinition *aProton = G4Proton::Proton();
  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
  G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
  G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
 
  G4int i, l;
  //G4double forVeryForward = 0.; // not used.
  G4bool veryForward = false;
 
  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/CLHEP::GeV;
  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/CLHEP::GeV;
  const G4double mOriginal = modifiedOriginal.GetMass()/CLHEP::GeV;
  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/CLHEP::GeV;
  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/CLHEP::GeV;
  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
                                           targetMass*targetMass +
                                           2.0*targetMass*etOriginal );  // GeV
  G4double currentMass = currentParticle.GetMass()/CLHEP::GeV;
  targetMass = targetParticle.GetMass()/CLHEP::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 )       // annihilation
    {
      // no kinetic energy in target .....
      G4double ek = currentParticle.GetKineticEnergy();
      G4ThreeVector m = 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()/CLHEP::GeV;
      targetMass = targetParticle.GetMass()/CLHEP::GeV;
      incidentHasChanged = true;
      targetHasChanged = true;
      currentParticle.SetKineticEnergy( ek );
      currentParticle.SetMomentum( m );
      //forVeryForward = aProton->GetPDGMass(); // not used.
      veryForward = true;
    }
  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());//targetNucleus.GetN();
  //    G4cout <<"Atomic weight is: "<<atomicWeight<<G4endl;
  const G4double atomicNumber = targetNucleus.GetZ_asInt();
  const G4double protonMass = aProton->GetPDGMass()/CLHEP::MeV;
  if( (originalIncident->GetDefinition() == aKaonMinus ||
       originalIncident->GetDefinition() == aKaonZeroL ||
       originalIncident->GetDefinition() == aKaonZeroS ||
       originalIncident->GetDefinition() == aKaonPlus) &&
      G4UniformRand() >= 0.7 )
    {
      G4ReactionProduct temp = currentParticle;
      currentParticle = targetParticle;
      targetParticle = temp;
      incidentHasChanged = true;
      targetHasChanged = true;
      currentMass = currentParticle.GetMass()/CLHEP::GeV;
      targetMass = targetParticle.GetMass()/CLHEP::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 );
 
  // PROBLEMET ER HER!!!
 
  G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
 
  if(freeEnergy<0)
    {
      G4cout<<"Free energy < 0!"<<G4endl;
      G4cout<<"E_CMS = "<<centerofmassEnergy<<" GeV"<<G4endl;
      G4cout<<"m_curr = "<<currentMass<<" GeV"<<G4endl;
      G4cout<<"m_orig = "<<mOriginal<<" GeV"<<G4endl;
      G4cout<<"m_targ = "<<targetMass<<" GeV"<<G4endl;
      G4cout<<"E_free = "<<freeEnergy<<" GeV"<<G4endl;
    }
 
  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()/CLHEP::GeV;
        } else {
        ++forwardCount;
        forwardEnergy -= vec[i]->GetMass()/CLHEP::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;
  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()/CLHEP::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
      //
      //      G4ReactionProduct *pVec = new G4ReactionProduct [nuclearExcitationCount];
      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()/CLHEP::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()/CLHEP::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()/CLHEP::GeV;
                  for( G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1];    // shift up
                  --forwardCount;
                  G4ReactionProduct *temp = vec[vecLen-1];
                  delete temp;
                  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 )
        {
          forwardEnergy += currentParticle.GetMass()/CLHEP::GeV;
          currentParticle.SetDefinitionAndUpdateE( targetParticle.GetDefinition() );
          targetParticle.SetDefinitionAndUpdateE( vec[0]->GetDefinition() );
          // above two lines modified 20-oct-97: were just simple equalities
          --forwardCount;
          for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
          G4ReactionProduct *temp = vec[vecLen-1];
          delete temp;
          if( --vecLen == 0 )return false;  // all the secondaries have been eliminated
          break;
        }
    }
  // 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;
                      backwardEnergy -= vec[i]->GetTotalEnergy()/CLHEP::GeV;
                    }
                  backwardEnergy += vec[i]->GetTotalEnergy()/CLHEP::GeV;
                  for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];   // shift up
                  --backwardCount;
                  G4ReactionProduct *temp = vec[vecLen-1];
                  delete temp;
                  if( --vecLen == 0 )return false;  // all the secondaries have been eliminated
                  break;
                }
            }
        }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      if( backwardParticlesLeft == 0 )
        {
          backwardEnergy += targetParticle.GetMass()/CLHEP::GeV;
          targetParticle = *vec[0];
          --backwardCount;
          for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
          G4ReactionProduct *temp = vec[vecLen-1];
          delete temp;
          if( --vecLen == 0 )return false;  // all the secondaries have been eliminated
          break;
        }
    }
  // 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*CLHEP::GeV );
  pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*CLHEP::GeV );
  pseudoParticle[0].SetTotalEnergy(
                                   std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*CLHEP::GeV );
 
  pseudoParticle[1].SetMass( protonMass*CLHEP::MeV );        // this could be targetMass
  pseudoParticle[1].SetTotalEnergy( protonMass*CLHEP::MeV );
 
  pseudoParticle[3].SetMass( protonMass*(1+extraNucleonCount)*CLHEP::MeV );
  pseudoParticle[3].SetTotalEnergy( protonMass*(1+extraNucleonCount)*CLHEP::MeV );
 
  pseudoParticle[8].SetMomentum( 1.0*CLHEP::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, rthnve, phinve, 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()*CLHEP::twopi;
      if( vec[i]->GetNewlyAdded() )           // added from intranuclear cascade
        {
          if( vec[i]->GetSide() == -2 )         //  is a nucleon
            {
              if( backwardNucleonCount < 18 )
                {
                  if( vec[i]->GetDefinition() == G4PionMinus::PionMinus() ||
                      vec[i]->GetDefinition() == G4PionPlus::PionPlus() ||
                      vec[i]->GetDefinition() == G4PionZero::PionZero() )
                    {
                      for(G4int i=0; i<vecLen; i++) delete vec[i];
                      vecLen = 0;
#if G4VERSION_NUMBER < 1100
                      throw G4HadReentrentException(__FILE__, __LINE__,
#else
                      throw G4HadronicException(__FILE__, __LINE__,
#endif
                                                    "FullModelReactionDynamics::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()/CLHEP::GeV;
      G4double ran = -std::log(1.0-G4UniformRand())/3.5;
      if( vec[i]->GetSide() == -2 )   // backward nucleon
        {
          if( vec[i]->GetDefinition() == aKaonMinus ||
              vec[i]->GetDefinition() == aKaonZeroL ||
              vec[i]->GetDefinition() == aKaonZeroS ||
              vec[i]->GetDefinition() == aKaonPlus ||
              vec[i]->GetDefinition() == aPiMinus ||
              vec[i]->GetDefinition() == aPiZero ||
              vec[i]->GetDefinition() == aPiPlus )
            {
              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() == aPiMinus ||
            vec[i]->GetDefinition() == aPiZero ||
            vec[i]->GetDefinition() == aPiPlus )
          {
            aspar = 0.75;
            pt = std::sqrt( std::pow( ran, 1.7 ) );
          } else if( vec[i]->GetDefinition() == aKaonMinus ||
                     vec[i]->GetDefinition() == aKaonZeroL ||
                     vec[i]->GetDefinition() == aKaonZeroS ||
                     vec[i]->GetDefinition() == aKaonPlus )
          {
            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)*CLHEP::GeV, pt*std::sin(phi)*CLHEP::GeV );
      for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
      if( vec[i]->GetSide() > 0 )
        et = pseudoParticle[0].GetTotalEnergy()/CLHEP::GeV;
      else
        et = pseudoParticle[1].GetTotalEnergy()/CLHEP::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)*CLHEP::GeV, pt*std::sin(phi)*CLHEP::GeV );
          //
          //   start of inner iteration loop
          //
          while( ++innerCounter < 7 )
            {
              ran = G4UniformRand()*dndl[19];
              l = 1;
              while(l<19 && ran>=dndl[l]) l++;
              x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1])/2.) );
              if( vec[i]->GetSide() < 0 )x *= -1.;
              vec[i]->SetMomentum( x*et*CLHEP::GeV );              // set the z-momentum
              totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
              vec[i]->SetTotalEnergy( totalEnergy*CLHEP::GeV );
              kineticEnergy = vec[i]->GetKineticEnergy()/CLHEP::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()/CLHEP::MeV < 0.0 )phi = CLHEP::twopi - phi;
                      phi += CLHEP::pi + normal()*CLHEP::pi/12.0;
                      if( phi > CLHEP::twopi )phi -= CLHEP::twopi;
                      if( phi < 0.0 )phi = CLHEP::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()/CLHEP::MeV < 0.0 )phi = CLHEP::twopi - phi;
                    phi += CLHEP::pi + normal() * CLHEP::pi / 12.0;
                    if( phi > CLHEP::twopi )phi -= CLHEP::twopi;
                    if( phi < 0.0 )phi = CLHEP::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() == aKaonMinus ||
                      vec[l]->GetDefinition() == aKaonZeroL ||
                      vec[l]->GetDefinition() == aKaonZeroS ||
                      vec[l]->GetDefinition() == aKaonPlus ||
                      vec[l]->GetDefinition() == aPiMinus ||
                      vec[l]->GetDefinition() == aPiZero ||
                      vec[l]->GetDefinition() == aPiPlus )
                    {
                      G4double tempMass = vec[l]->GetMass()/CLHEP::MeV;
                      totalEnergy = 0.95*vec[l]->GetTotalEnergy()/CLHEP::MeV + 0.05*tempMass;
                      totalEnergy = std::max( tempMass, totalEnergy );
                      vec[l]->SetTotalEnergy( totalEnergy*CLHEP::MeV );
                      pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
                      pp1 = vec[l]->GetMomentum().mag()/CLHEP::MeV;
                      if( pp1 < 1.0e-6*CLHEP::GeV )
                        {
                          G4double rthnve = CLHEP::pi*G4UniformRand();
                          G4double phinve = CLHEP::twopi*G4UniformRand();
                          G4double srth = std::sin(rthnve);
                          vec[l]->SetMomentum( pp*srth*std::cos(phinve)*CLHEP::MeV,
                                               pp*srth*std::sin(phinve)*CLHEP::MeV,
                                               pp*std::cos(rthnve)*CLHEP::MeV ) ;
                        } else {
                        vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
                      }
                      G4double px = vec[l]->GetMomentum().x()/CLHEP::MeV;
                      G4double py = vec[l]->GetMomentum().y()/CLHEP::MeV;
                      pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/CLHEP::GeV;
                      if( vec[l]->GetSide() > 0 )
                        {
                          forwardKinetic += vec[l]->GetKineticEnergy()/CLHEP::GeV;
                          pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
                        } else {
                        backwardKinetic += vec[l]->GetKineticEnergy()/CLHEP::GeV;
                        pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
                      }
                    }
                }
            }
        }    // 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;
                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()/CLHEP::MeV < 0.0 )phi = CLHEP::twopi - phi;
          phi += CLHEP::pi + normal() * CLHEP::pi / 12.0;
          if( phi > CLHEP::twopi )phi -= CLHEP::twopi;
          if( phi < 0.0 )phi = CLHEP::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()*CLHEP::twopi;
  G4double ran = -std::log(1.0-G4UniformRand());
  if( currentParticle.GetDefinition() == aPiMinus ||
      currentParticle.GetDefinition() == aPiZero ||
      currentParticle.GetDefinition() == aPiPlus )
    {
      aspar = 0.60;
      pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
    } else if( currentParticle.GetDefinition() == aKaonMinus ||
               currentParticle.GetDefinition() == aKaonZeroL ||
               currentParticle.GetDefinition() == aKaonZeroS ||
               currentParticle.GetDefinition() == aKaonPlus )
    {
      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)*CLHEP::GeV, pt*std::sin(phi)*CLHEP::GeV );
  et = pseudoParticle[0].GetTotalEnergy()/CLHEP::GeV;
  dndl[0] = 0.0;
  vecMass = currentParticle.GetMass()/CLHEP::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( (l<20) && (ran>dndl[l]) )l++;
  l = std::min( 19, l );
  x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1])/2.) );
  currentParticle.SetMomentum( x*et*CLHEP::GeV );                 // set the z-momentum
  if( forwardEnergy < forwardKinetic )
    totalEnergy = vecMass + 0.04*std::fabs(normal());
  else
    totalEnergy = vecMass + forwardEnergy - forwardKinetic;
  currentParticle.SetTotalEnergy( totalEnergy*CLHEP::GeV );
  pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*CLHEP::GeV;
  pp1 = currentParticle.GetMomentum().mag()/CLHEP::MeV;
  if( pp1 < 1.0e-6*CLHEP::GeV )
    {
      G4double rthnve = CLHEP::pi*G4UniformRand();
      G4double phinve = CLHEP::twopi*G4UniformRand();
      G4double srth = std::sin(rthnve);
      currentParticle.SetMomentum( pp*srth*std::cos(phinve)*CLHEP::MeV,
                                   pp*srth*std::sin(phinve)*CLHEP::MeV,
                                   pp*std::cos(rthnve)*CLHEP::MeV );
    } else {
    currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
  }
  pseudoParticle[4] = pseudoParticle[4] + currentParticle;
  //
  // this finishes the current particle
  // now for the 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()/CLHEP::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)*CLHEP::GeV, pt*std::sin(phi)*CLHEP::GeV );
      for( G4int j=0; j<20; ++j )binl[j] = (j-1.)/(19.*pt);
      et = pseudoParticle[1].GetTotalEnergy()/CLHEP::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( ( l < 20 ) && ( ran >= dndl[l] ) ) l++;
              l = std::min( 19, l );
              x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1])/2.) );
              if( targetParticle.GetSide() < 0 )x *= -1.;
              targetParticle.SetMomentum( x*et*CLHEP::GeV );                // set the z-momentum
              totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
              targetParticle.SetTotalEnergy( totalEnergy*CLHEP::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()/CLHEP::MeV < 0.0 )phi = CLHEP::twopi - phi;
                      phi += CLHEP::pi + normal() * CLHEP::pi / 12.0;
                      if( phi > CLHEP::twopi )phi -= CLHEP::twopi;
                      if( phi < 0.0 )phi = CLHEP::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*CLHEP::GeV );
                  pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*CLHEP::GeV;
                  pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV;
                  if( pp1 < 1.0e-6*CLHEP::GeV )
                    {
                      G4double rthnve = CLHEP::pi*G4UniformRand();
                      G4double phinve = CLHEP::twopi*G4UniformRand();
                      G4double srth = std::sin(rthnve);
                      targetParticle.SetMomentum( pp*srth*std::cos(phinve)*CLHEP::MeV,
                                                  pp*srth*std::sin(phinve)*CLHEP::MeV,
                                                  pp*std::cos(rthnve)*CLHEP::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() == aKaonMinus ||
                      vec[l]->GetDefinition() == aKaonZeroL ||
                      vec[l]->GetDefinition() == aKaonZeroS ||
                      vec[l]->GetDefinition() == aKaonPlus ||
                      vec[l]->GetDefinition() == aPiMinus ||
                      vec[l]->GetDefinition() == aPiZero ||
                      vec[l]->GetDefinition() == aPiPlus )
                    {
                      G4double tempMass = vec[l]->GetMass()/CLHEP::GeV;
                      totalEnergy =
                        std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/CLHEP::GeV + 0.05*tempMass );
                      vec[l]->SetTotalEnergy( totalEnergy*CLHEP::GeV );
                      pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*CLHEP::GeV;
                      pp1 = vec[l]->GetMomentum().mag()/CLHEP::MeV;
                      if( pp1 < 1.0e-6*CLHEP::GeV )
                        {
                          G4double rthnve = CLHEP::pi*G4UniformRand();
                          G4double phinve = CLHEP::twopi*G4UniformRand();
                          G4double srth = std::sin(rthnve);
                          vec[l]->SetMomentum( pp*srth*std::cos(phinve)*CLHEP::MeV,
                                               pp*srth*std::sin(phinve)*CLHEP::MeV,
                                               pp*std::cos(rthnve)*CLHEP::MeV );
                        }
                      else
                        vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
 
                      pt = std::max( 0.001*CLHEP::GeV, std::sqrt( sqr(vec[l]->GetMomentum().x()/CLHEP::MeV) +
                                                                  sqr(vec[l]->GetMomentum().y()/CLHEP::MeV) ) )/CLHEP::GeV;
                      if( vec[l]->GetSide() > 0)
                        {
                          forwardKinetic += vec[l]->GetKineticEnergy()/CLHEP::GeV;
                          pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
                        } else {
                        backwardKinetic += vec[l]->GetKineticEnergy()/CLHEP::GeV;
                        pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
                      }
                    }
                }
            }
        }    // closes outer loop
      //      if( eliminateThisParticle )  // not enough energy, eliminate target
      //      {
      //        G4cerr << "Warning: eliminating target particle" << G4endl;
      //        exit( EXIT_FAILURE );
      //      }
    }
  //
  //  this finishes the target particle
  // backward nucleons produced 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/CLHEP::GeV );
      if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
      vecMass = targetParticle.GetMass()/CLHEP::GeV;
      totalEnergy = ekin+vecMass;
      targetParticle.SetTotalEnergy( totalEnergy*CLHEP::GeV );
      pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*CLHEP::GeV;
      pp1 = pseudoParticle[6].GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 1.0e-6*CLHEP::GeV )
        {
          rthnve = CLHEP::pi*G4UniformRand();
          phinve = CLHEP::twopi*G4UniformRand();
          G4double srth = std::sin(rthnve);
          targetParticle.SetMomentum( pp*srth*std::cos(phinve)*CLHEP::MeV,
                                      pp*srth*std::sin(phinve)*CLHEP::MeV,
                                      pp*std::cos(rthnve)*CLHEP::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()/CLHEP::GeV;
      for( i=0; i<vecLen; ++i )
        {
          if( vec[i]->GetSide() == -3 )rmb0 += vec[i]->GetMass()/CLHEP::GeV;
        }
      rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()),cpar[tempCount]) / gpar[tempCount];
      totalEnergy = pseudoParticle[6].GetTotalEnergy()/CLHEP::GeV;
      vecMass = std::min( rmb, totalEnergy );
      pseudoParticle[6].SetMass( vecMass*CLHEP::GeV );
      pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*CLHEP::GeV;
      pp1 = pseudoParticle[6].GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 1.0e-6*CLHEP::GeV )
        {
          rthnve = CLHEP::pi * G4UniformRand();
          phinve = CLHEP::twopi * G4UniformRand();
          G4double srth = std::sin(rthnve);
          pseudoParticle[6].SetMomentum( -pp*srth*std::cos(phinve)*CLHEP::MeV,
                                         -pp*srth*std::sin(phinve)*CLHEP::MeV,
                                         -pp*std::cos(rthnve)*CLHEP::MeV );
        }
      else
        pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
 
      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> 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;
          throw std::runtime_error("FullModelReactionDynamics::GenerateXandPt: "
                                   "tempLen is not the same as backwardNucleonCount");
        }
      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);
 
  G4int numberofFinalStateNucleons = 0;
  if( currentParticle.GetDefinition() ==aProton ||
      currentParticle.GetDefinition() == aNeutron ||
      currentParticle.GetDefinition() == G4SigmaMinus::SigmaMinus()||
      currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus()||
      currentParticle.GetDefinition() == G4SigmaZero::SigmaZero()||
      currentParticle.GetDefinition() == G4XiZero::XiZero()||
      currentParticle.GetDefinition() == G4XiMinus::XiMinus()||
      currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus()||
      currentParticle.GetDefinition() == G4Lambda::Lambda()) ++numberofFinalStateNucleons;
  currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
 
  if( targetParticle.GetDefinition() ==aProton ||
      targetParticle.GetDefinition() == aNeutron ||
      targetParticle.GetDefinition() == G4Lambda::Lambda() ||
      targetParticle.GetDefinition() == G4XiZero::XiZero()||
      targetParticle.GetDefinition() == G4XiMinus::XiMinus()||
      targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus()||
      targetParticle.GetDefinition() == G4SigmaZero::SigmaZero()||
      targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus()||
      targetParticle.GetDefinition() == G4SigmaMinus::SigmaMinus()) ++numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiProton::AntiProton()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiNeutron::AntiNeutron()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiSigmaMinus::AntiSigmaMinus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiSigmaPlus::AntiSigmaPlus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiSigmaZero::AntiSigmaZero()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiXiZero::AntiXiZero()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiXiMinus::AntiXiMinus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiOmegaMinus::AntiOmegaMinus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiLambda::AntiLambda()) --numberofFinalStateNucleons;
  targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
 
  for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetDefinition() ==aProton ||
          vec[i]->GetDefinition() == aNeutron ||
          vec[i]->GetDefinition() == G4Lambda::Lambda() ||
          vec[i]->GetDefinition() == G4XiZero::XiZero() ||
          vec[i]->GetDefinition() == G4XiMinus::XiMinus() ||
          vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
          vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus()||
          vec[i]->GetDefinition() == G4SigmaZero::SigmaZero()||
          vec[i]->GetDefinition() == G4SigmaMinus::SigmaMinus()) ++numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiProton::AntiProton()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiNeutron::AntiNeutron()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiSigmaMinus::AntiSigmaMinus()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiSigmaPlus::AntiSigmaPlus()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiSigmaZero::AntiSigmaZero()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiLambda::AntiLambda()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiXiZero::AntiXiZero()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiXiMinus::AntiXiMinus()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiOmegaMinus::AntiOmegaMinus()) --numberofFinalStateNucleons;
      vec[i]->Lorentz( *vec[i], pseudoParticle[1] );
    }
  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  if(veryForward) numberofFinalStateNucleons++;
  numberofFinalStateNucleons = std::max( 1, numberofFinalStateNucleons );
  //
  // 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() == aKaonMinus ||
             leadingStrangeParticle.GetDefinition() == aKaonZeroL ||
             leadingStrangeParticle.GetDefinition() == aKaonZeroS ||
             leadingStrangeParticle.GetDefinition() == aKaonPlus ||
             leadingStrangeParticle.GetDefinition() == aPiMinus ||
             leadingStrangeParticle.GetDefinition() == aPiZero ||
             leadingStrangeParticle.GetDefinition() == aPiPlus);
          G4bool targetTest = false;
          /*          (targetParticle.GetDefinition() == aKaonMinus ||
                      targetParticle.GetDefinition() == aKaonZeroL ||
                      targetParticle.GetDefinition() == aKaonZeroS ||
                      targetParticle.GetDefinition() == aKaonPlus ||
                      targetParticle.GetDefinition() == aPiMinus ||
                      targetParticle.GetDefinition() == aPiZero ||
                      targetParticle.GetDefinition() == aPiPlus);
          */
          // 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 )
 
  pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*CLHEP::GeV );
  pseudoParticle[3].SetMass( mOriginal*CLHEP::GeV );
  pseudoParticle[3].SetTotalEnergy(
                                   std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*CLHEP::GeV );
 
  //    G4double ekin0 = pseudoParticle[3].GetKineticEnergy()/CLHEP::GeV;
 
  const 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)*CLHEP::MeV );
  pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*CLHEP::MeV );
 
  G4double theoreticalKinetic =
    pseudoParticle[3].GetTotalEnergy()/CLHEP::MeV +
    pseudoParticle[4].GetTotalEnergy()/CLHEP::MeV -
    currentParticle.GetMass()/CLHEP::MeV -
    targetParticle.GetMass()/CLHEP::MeV;
 
  G4double simulatedKinetic =
    currentParticle.GetKineticEnergy()/CLHEP::MeV + targetParticle.GetKineticEnergy()/CLHEP::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()/CLHEP::MeV;
      theoreticalKinetic -= vec[i]->GetMass()/CLHEP::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];
      //G4ReactionProduct *tempR = new G4ReactionProduct [vecLen+2];
      tempR[0] = currentParticle;
      tempR[1] = targetParticle;
      for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> 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()/CLHEP::MeV+
                                pseudoParticle[4].GetTotalEnergy()/CLHEP::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()/CLHEP::MeV;
            }
        }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      //delete [] tempR;
    }
  //
  //  Make sure, that the kinetic energies are correct
  //
  if( simulatedKinetic != 0.0 )
    {
      wgt = (theoreticalKinetic)/simulatedKinetic;
      theoreticalKinetic = currentParticle.GetKineticEnergy()/CLHEP::MeV * wgt;
      simulatedKinetic = theoreticalKinetic;
      currentParticle.SetKineticEnergy( theoreticalKinetic*CLHEP::MeV );
      pp = currentParticle.GetTotalMomentum()/CLHEP::MeV;
      pp1 = currentParticle.GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 1.0e-6*CLHEP::GeV )
        {
          rthnve = CLHEP::pi*G4UniformRand();
          phinve = CLHEP::twopi*G4UniformRand();
          currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                       pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                       pp*std::cos(rthnve)*CLHEP::MeV );
        }
      else
        {
          currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
        }
      theoreticalKinetic = targetParticle.GetKineticEnergy()/CLHEP::MeV * wgt;
      targetParticle.SetKineticEnergy( theoreticalKinetic*CLHEP::MeV );
      simulatedKinetic += theoreticalKinetic;
      pp = targetParticle.GetTotalMomentum()/CLHEP::MeV;
      pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      if( pp1 < 1.0e-6*CLHEP::GeV )
        {
          rthnve = CLHEP::pi*G4UniformRand();
          phinve = CLHEP::twopi*G4UniformRand();
          targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                      pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                      pp*std::cos(rthnve)*CLHEP::MeV );
        } else {
        targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
      }
      for( i=0; i<vecLen; ++i )
        {
          theoreticalKinetic = vec[i]->GetKineticEnergy()/CLHEP::MeV * wgt;
          simulatedKinetic += theoreticalKinetic;
          vec[i]->SetKineticEnergy( theoreticalKinetic*CLHEP::MeV );
          pp = vec[i]->GetTotalMomentum()/CLHEP::MeV;
          pp1 = vec[i]->GetMomentum().mag()/CLHEP::MeV;
          if( pp1 < 1.0e-6*CLHEP::GeV )
            {
              rthnve = CLHEP::pi*G4UniformRand();
              phinve = CLHEP::twopi*G4UniformRand();
              vec[i]->SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                   pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                   pp*std::cos(rthnve)*CLHEP::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
      //
      G4double epnb, edta;
      G4int npnb = 0;
      G4int ndta = 0;
 
      epnb = targetNucleus.GetPNBlackTrackEnergy();            // was enp1 in fortran code
      edta = targetNucleus.GetDTABlackTrackEnergy();           // was enp3 in fortran code
      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()/CLHEP::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 );
        }
      G4double spall = numberofFinalStateNucleons;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
 
      AddBlackTrackParticles( epnb, npnb, edta, ndta, sprob, kineticMinimum, kineticFactor,
                              modifiedOriginal, spall, targetNucleus,
                              vec, vecLen );
 
      //       G4double jpw=0;
      //       jpw+=GetQValue(&currentParticle);
      //       jpw+=GetQValue(&targetParticle);
      //       for( i=0; i<vecLen; ++i )jpw += GetQValue(vec[i]);
      //       G4cout << "JPW ### "<<jpw<<G4endl;
      // 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 FullModelReactionDynamics::SuppressChargedPions(
                                                     G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                                     G4int &vecLen,
                                                     const G4ReactionProduct &modifiedOriginal,
                                                     G4ReactionProduct &currentParticle,
                                                     //   G4ReactionProduct &targetParticle,
                                                     const G4Nucleus &targetNucleus,
                                                     G4bool &incidentHasChanged
                                                     //   G4bool &targetHasChanged
                                                     )
{
  // this code was originally in the fortran code TWOCLU
  //
  // suppress charged pions, for various reasons
  //
  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());//targetNucleus.GetN();
  const G4double atomicNumber = targetNucleus.GetZ_asInt();
  const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/CLHEP::GeV;
 
  // G4ParticleDefinition *aGamma = G4Gamma::Gamma();
  G4ParticleDefinition *aProton = G4Proton::Proton();
  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
 
  const G4bool antiTest =
    modifiedOriginal.GetDefinition() != anAntiProton &&
    modifiedOriginal.GetDefinition() != anAntiNeutron;
  if( antiTest && (
                   // currentParticle.GetDefinition() == aGamma ||
                   currentParticle.GetDefinition() == aPiPlus ||
                   currentParticle.GetDefinition() == aPiMinus ) &&
      ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
      ( G4UniformRand() <= atomicWeight/300.0 ) )
    {
      if( G4UniformRand() > atomicNumber/atomicWeight )
        currentParticle.SetDefinitionAndUpdateE( aNeutron );
      else
        currentParticle.SetDefinitionAndUpdateE( aProton );
      incidentHasChanged = true;
    }
  /*    if( antiTest && (
  // targetParticle.GetDefinition() == aGamma ||
  targetParticle.GetDefinition() == aPiPlus ||
  targetParticle.GetDefinition() == aPiMinus ) &&
  ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
  ( G4UniformRand() <= atomicWeight/300.0 ) )
  {
  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() == aGamma ||
                       vec[i]->GetDefinition() == aPiPlus ||
                       vec[i]->GetDefinition() == aPiMinus ) &&
          ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
          ( G4UniformRand() <= atomicWeight/300.0 ) )
        {
          if( G4UniformRand() > atomicNumber/atomicWeight )
            vec[i]->SetDefinitionAndUpdateE( aNeutron );
          else
            vec[i]->SetDefinitionAndUpdateE( aProton );
          // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
    }
  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
 
G4bool FullModelReactionDynamics::TwoCluster(
                                             G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                             G4int &vecLen,
                                             G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effects included
                                             const G4HadProjectile *originalIncident, // the original incident particle
                                             G4ReactionProduct &currentParticle,
                                             G4ReactionProduct &targetParticle,
                                             const G4Nucleus &targetNucleus,
                                             G4bool &incidentHasChanged,
                                             G4bool &targetHasChanged,
                                             G4bool leadFlag,
                                             G4ReactionProduct &leadingStrangeParticle )
{
  // 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 *aPiMinus = G4PionMinus::PionMinus();
  G4ParticleDefinition *aProton = G4Proton::Proton();
  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
  const G4double protonMass = aProton->GetPDGMass()/CLHEP::MeV;
  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/CLHEP::GeV;
  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/CLHEP::GeV;
  const G4double mOriginal = modifiedOriginal.GetMass()/CLHEP::GeV;
  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/CLHEP::GeV;
  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/CLHEP::GeV;
  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
                                           targetMass*targetMass +
                                           2.0*targetMass*etOriginal );  // GeV
  G4double currentMass = currentParticle.GetMass()/CLHEP::GeV;
  targetMass = targetParticle.GetMass()/CLHEP::GeV;
 
  if( currentMass == 0.0 && targetMass == 0.0 )
    {
      G4double ek = currentParticle.GetKineticEnergy();
      G4ThreeVector m = 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 i=0; i<vecLen; i++) delete vec[i];
          vecLen = 0;
#if G4VERSION_NUMBER < 1100
          throw G4HadReentrentException(__FILE__, __LINE__,
#else
          throw G4HadronicException(__FILE__, __LINE__,
#endif
                                        "FullModelReactionDynamics::TwoCluster: Negative number of particles");
        }
      delete vec[vecLen-1];
      delete vec[vecLen-2];
      vecLen -= 2;
      currentMass = currentParticle.GetMass()/CLHEP::GeV;
      targetMass = targetParticle.GetMass()/CLHEP::GeV;
      incidentHasChanged = true;
      targetHasChanged = true;
      currentParticle.SetKineticEnergy( ek );
      currentParticle.SetMomentum( m );
    }
  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());//targetNucleus.GetN();
  const G4double atomicNumber = 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()/CLHEP::GeV;
 
  // target is always in backward hemisphere
  G4int backwardCount = 1;           // number of particles in backward hemisphere
  targetParticle.SetSide( -1 );
  G4double backwardMass = targetParticle.GetMass()/CLHEP::GeV;
 
  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()/CLHEP::GeV;
        }
      else
        {
          ++forwardCount;
          forwardMass += vec[i]->GetMass()/CLHEP::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;
  //G4double extraMass = 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 )
                // HPW: looks like a gheisha bug
                pVec->SetDefinition( aProton );
              else
                pVec->SetDefinition( aNeutron );
            }
          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()/CLHEP::GeV;
          pVec->SetNewlyAdded( true );
          vec.SetElement( vecLen++, pVec );
          // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
    }
  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  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()/CLHEP::GeV;
              for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];     // shift up
              --forwardCount;
              forwardMass -= pMass;
              secondaryDeleted = true;
              break;
            }
          else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
            {
              pMass = vec[i]->GetMass()/CLHEP::GeV;
              for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
              --backwardCount;
              backwardMass -= pMass;
              secondaryDeleted = true;
              break;
            }
        }           // breaks go down to here
      if( secondaryDeleted )
        {
          G4ReactionProduct *temp = vec[vecLen-1];
          delete temp;
          --vecLen;
          // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
      else
        {
          if( vecLen == 0 )
            {
              return false;  // all secondaries have been eliminated
            }
          if( targetParticle.GetSide() == -1 )
            {
              pMass = targetParticle.GetMass()/CLHEP::GeV;
              targetParticle = *vec[0];
              for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
              --backwardCount;
              backwardMass -= pMass;
              secondaryDeleted = true;
            }
          else if( targetParticle.GetSide() == 1 )
            {
              pMass = targetParticle.GetMass()/CLHEP::GeV;
              targetParticle = *vec[0];
              for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
              --forwardCount;
              forwardMass -= pMass;
              secondaryDeleted = true;
            }
          if( secondaryDeleted )
            {
              G4ReactionProduct *temp = vec[vecLen-1];
              delete temp;
              --vecLen;
              // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
            }
          else
            {
              if( currentParticle.GetSide() == -1 )
                {
                  pMass = currentParticle.GetMass()/CLHEP::GeV;
                  currentParticle = *vec[0];
                  for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
                  --backwardCount;
                  backwardMass -= pMass;
                  secondaryDeleted = true;
                }
              else if( currentParticle.GetSide() == 1 )
                {
                  pMass = currentParticle.GetMass()/CLHEP::GeV;
                  currentParticle = *vec[0];
                  for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
                  --forwardCount;//This line can cause infinite loop
                  forwardMass -= pMass;
                  secondaryDeleted = true;
                }
              if( secondaryDeleted )
                {
                  G4ReactionProduct *temp = vec[vecLen-1];
                  delete temp;
                  --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; // rme = 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) return false;
 
  if( forwardCount == 1 )rmc = forwardMass;
  else
    {
      //      G4int ntc = std::min(5,forwardCount); // check if offset by 1 @@
      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::min(5,backwardCount); // check, if offfset by 1 @@
      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;
        }
    }
  // note that rme is never used below this section
  //
  //if( nuclearExcitationCount == 0 )rme = 0.0;
  //else if( nuclearExcitationCount == 1 )rme = extraMass;
  //else
  //{
  //  G4int ntc = std::min(5,nuclearExcitationCount)-1;
  //  rme = extraMass + std::pow(-std::log(1.-G4UniformRand()),cpar[ntc])/gpar[ntc];
  //}
  //
  //  Set beam, target of first interaction in centre of mass system
  //
  G4ReactionProduct pseudoParticle[8];
  for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
 
  pseudoParticle[1].SetMass( mOriginal*CLHEP::GeV );
  pseudoParticle[1].SetTotalEnergy( etOriginal*CLHEP::GeV );
  pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*CLHEP::GeV );
 
  pseudoParticle[2].SetMass( protonMass*CLHEP::MeV );
  pseudoParticle[2].SetTotalEnergy( protonMass*CLHEP::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*CLHEP::GeV );
  pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*CLHEP::GeV );
 
  pseudoParticle[4].SetMass( rmd*CLHEP::GeV );
  pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*CLHEP::GeV );
  //
  // set |T| and |TMIN|
  //
  const G4double bMin = 0.01;
  const G4double b1 = 4.0;
  const G4double b2 = 1.6;
  G4double t = std::log( 1.0-G4UniformRand() ) / std::max( bMin, b1+b2*std::log(pOriginal) );
  G4double t1 =
    pseudoParticle[1].GetTotalEnergy()/CLHEP::GeV - pseudoParticle[3].GetTotalEnergy()/CLHEP::GeV;
  G4double pin = pseudoParticle[1].GetMomentum().mag()/CLHEP::GeV;
  G4double tacmin = t1*t1 - (pin-pf)*(pin-pf);
  //
  // calculate (std::sin(teta/2.)^2 and std::cos(teta), set azimuth angle phi
  //
  const G4double smallValue = 1.0e-10;
  G4double dumnve = 4.0*pin*pf;
  if( dumnve == 0.0 )dumnve = smallValue;
  G4double ctet = std::max( -1.0, std::min( 1.0, 1.0+2.0*(t-tacmin)/dumnve ) );
  dumnve = std::max( 0.0, 1.0-ctet*ctet );
  G4double stet = std::sqrt(dumnve);
  G4double phi = G4UniformRand() * CLHEP::twopi;
  //
  // calculate final state momenta in centre of mass system
  //
  pseudoParticle[3].SetMomentum( pf*stet*std::sin(phi)*CLHEP::GeV,
                                 pf*stet*std::cos(phi)*CLHEP::GeV,
                                 pf*ctet*CLHEP::GeV );
  pseudoParticle[4].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
  //
  // simulate backward nucleon cluster in lab. system and transform in cms
  //
  G4double pp, pp1, rthnve, phinve;
  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*CLHEP::GeV );
              G4double vMass = vec[i]->GetMass()/CLHEP::MeV;
              G4double totalE = kineticE + 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 = CLHEP::twopi*G4UniformRand();
              vec[i]->SetMomentum( pp*sint*std::sin(phi)*CLHEP::MeV,
                                   pp*sint*std::cos(phi)*CLHEP::MeV,
                                   pp*cost*CLHEP::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,MYGHADLISTSIZE> 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()/CLHEP::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,MYGHADLISTSIZE> 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()/CLHEP::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
  //
  G4int numberofFinalStateNucleons = 0;
  if( currentParticle.GetDefinition() ==aProton ||
      currentParticle.GetDefinition() == aNeutron ||
      currentParticle.GetDefinition() == aSigmaMinus||
      currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus()||
      currentParticle.GetDefinition() == G4SigmaZero::SigmaZero()||
      currentParticle.GetDefinition() == G4XiZero::XiZero()||
      currentParticle.GetDefinition() == G4XiMinus::XiMinus()||
      currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus()||
      currentParticle.GetDefinition() == G4Lambda::Lambda()) ++numberofFinalStateNucleons;
  currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
  if( targetParticle.GetDefinition() ==aProton ||
      targetParticle.GetDefinition() == aNeutron ||
      targetParticle.GetDefinition() == G4Lambda::Lambda() ||
      targetParticle.GetDefinition() == G4XiZero::XiZero()||
      targetParticle.GetDefinition() == G4XiMinus::XiMinus()||
      targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus()||
      targetParticle.GetDefinition() == G4SigmaZero::SigmaZero()||
      targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus()||
      targetParticle.GetDefinition() == aSigmaMinus) ++numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiProton::AntiProton()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiNeutron::AntiNeutron()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiSigmaMinus::AntiSigmaMinus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiSigmaPlus::AntiSigmaPlus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiSigmaZero::AntiSigmaZero()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiXiZero::AntiXiZero()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiXiMinus::AntiXiMinus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiOmegaMinus::AntiOmegaMinus()) --numberofFinalStateNucleons;
  if( targetParticle.GetDefinition() ==G4AntiLambda::AntiLambda()) --numberofFinalStateNucleons;
  targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
  for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetDefinition() ==aProton ||
          vec[i]->GetDefinition() == aNeutron ||
          vec[i]->GetDefinition() == G4Lambda::Lambda() ||
          vec[i]->GetDefinition() == G4XiZero::XiZero() ||
          vec[i]->GetDefinition() == G4XiMinus::XiMinus() ||
          vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
          vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus()||
          vec[i]->GetDefinition() == G4SigmaZero::SigmaZero()||
          vec[i]->GetDefinition() == aSigmaMinus) ++numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiProton::AntiProton()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiNeutron::AntiNeutron()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiSigmaMinus::AntiSigmaMinus()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiSigmaPlus::AntiSigmaPlus()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiSigmaZero::AntiSigmaZero()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiLambda::AntiLambda()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiXiZero::AntiXiZero()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiXiMinus::AntiXiMinus()) --numberofFinalStateNucleons;
      if( vec[i]->GetDefinition() ==G4AntiOmegaMinus::AntiOmegaMinus()) --numberofFinalStateNucleons;
      vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
    }
  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  numberofFinalStateNucleons = std::max( 1, numberofFinalStateNucleons );
  //
  // 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()/CLHEP::MeV;
          G4double ekin;
          if( ((leadMass <  protonMass) && (targetParticle.GetMass()/CLHEP::MeV <  protonMass)) ||
              ((leadMass >= protonMass) && (targetParticle.GetMass()/CLHEP::MeV >= protonMass)) )
            {
              ekin = targetParticle.GetKineticEnergy()/CLHEP::GeV;
              pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV; // old momentum
              targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
              targetParticle.SetKineticEnergy( ekin*CLHEP::GeV );
              pp = targetParticle.GetTotalMomentum()/CLHEP::MeV;   // new momentum
              if( pp1 < 1.0e-3 )
                {
                  rthnve = CLHEP::pi*G4UniformRand();
                  phinve = CLHEP::twopi*G4UniformRand();
                  targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(rthnve)*CLHEP::MeV,
                                              pp*std::sin(rthnve)*std::sin(rthnve)*CLHEP::MeV,
                                              pp*std::cos(rthnve)*CLHEP::MeV );
                }
              else
                targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
 
              targetHasChanged = true;
            }
          else
            {
              ekin = currentParticle.GetKineticEnergy()/CLHEP::GeV;
              pp1 = currentParticle.GetMomentum().mag()/CLHEP::MeV;
              currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
              currentParticle.SetKineticEnergy( ekin*CLHEP::GeV );
              pp = currentParticle.GetTotalMomentum()/CLHEP::MeV;
              if( pp1 < 1.0e-3 )
                {
                  rthnve = CLHEP::pi*G4UniformRand();
                  phinve = CLHEP::twopi*G4UniformRand();
                  currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(rthnve)*CLHEP::MeV,
                                               pp*std::sin(rthnve)*std::sin(rthnve)*CLHEP::MeV,
                                               pp*std::cos(rthnve)*CLHEP::MeV );
                }
              else
                currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
 
              incidentHasChanged = true;
            }
        }
    }    // end of if( leadFlag )
  //
  //  for various reasons, the energy balance is not sufficient,
  //  check that,  energy balance, angle of final system, etc.
  //
  pseudoParticle[4].SetMass( mOriginal*CLHEP::GeV );
  pseudoParticle[4].SetTotalEnergy( etOriginal*CLHEP::GeV );
  pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*CLHEP::GeV );
 
  const 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)*CLHEP::MeV );
  pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*CLHEP::MeV );
 
  //    G4double ekin0 = pseudoParticle[4].GetKineticEnergy()/CLHEP::GeV;
  G4double theoreticalKinetic =
    pseudoParticle[4].GetTotalEnergy()/CLHEP::GeV + pseudoParticle[5].GetTotalEnergy()/CLHEP::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];
      //G4ReactionProduct *tempR = new G4ReactionProduct [vecLen+2];
      tempR[0] = currentParticle;
      tempR[1] = targetParticle;
      for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
 
      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> 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()/CLHEP::MeV+pseudoParticle[5].GetTotalEnergy()/CLHEP::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()/CLHEP::GeV;
            }
        }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      //delete [] tempR;
    }
  else
    {
      theoreticalKinetic -=
        ( currentParticle.GetMass()/CLHEP::GeV + targetParticle.GetMass()/CLHEP::GeV );
      for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/CLHEP::GeV;
    }
  G4double simulatedKinetic =
    currentParticle.GetKineticEnergy()/CLHEP::GeV + targetParticle.GetKineticEnergy()/CLHEP::GeV;
  for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/CLHEP::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()/CLHEP::MeV;
      pp1 = currentParticle.GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 0.001*CLHEP::MeV )
        {
          rthnve = CLHEP::pi * G4UniformRand();
          phinve = CLHEP::twopi * G4UniformRand();
          currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                       pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                       pp*std::cos(rthnve)*CLHEP::MeV );
        }
      else
        currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
 
      targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
      pp = targetParticle.GetTotalMomentum()/CLHEP::MeV;
      pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 0.001*CLHEP::MeV )
        {
          rthnve = CLHEP::pi * G4UniformRand();
          phinve = CLHEP::twopi * G4UniformRand();
          targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                      pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                      pp*std::cos(rthnve)*CLHEP::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()/CLHEP::MeV;
          pp1 = vec[i]->GetMomentum().mag()/CLHEP::MeV;
          if( pp1 < 0.001 )
            {
              rthnve = CLHEP::pi * G4UniformRand();
              phinve = CLHEP::twopi * G4UniformRand();
              vec[i]->SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                   pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                   pp*std::cos(rthnve)*CLHEP::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
      //
      G4double epnb, edta;
      G4int npnb = 0;
      G4int ndta = 0;
 
      epnb = targetNucleus.GetPNBlackTrackEnergy();            // was enp1 in fortran code
      edta = targetNucleus.GetDTABlackTrackEnergy();           // was enp3 in fortran code
      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()/CLHEP::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 );
        }
      G4double spall = numberofFinalStateNucleons;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
 
      AddBlackTrackParticles( epnb, npnb, edta, ndta, sprob, kineticMinimum, kineticFactor,
                              modifiedOriginal, spall, 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;
}
 
void FullModelReactionDynamics::TwoBody(
                                        G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                        G4int &vecLen,
                                        G4ReactionProduct &modifiedOriginal,
                                        const G4DynamicParticle */*originalTarget*/,
                                        G4ReactionProduct &currentParticle,
                                        G4ReactionProduct &targetParticle,
                                        const G4Nucleus &targetNucleus,
                                        G4bool &/* targetHasChanged*/ )
{
  //    G4cout<<"TwoBody called"<<G4endl;
  //
  // 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.
  //
  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
  G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
  G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
 
  // 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()/CLHEP::GeV;
  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/CLHEP::GeV;
  G4double currentMass = currentParticle.GetMass()/CLHEP::GeV;
  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/CLHEP::GeV;
 
  targetMass = targetParticle.GetMass()/CLHEP::GeV;
  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());//targetNucleus.GetN();
  //    G4cout<<"Atomic weight is found to be: "<<atomicWeight<<G4endl;
  G4double etCurrent = currentParticle.GetTotalEnergy()/CLHEP::GeV;
  G4double pCurrent = currentParticle.GetTotalMomentum()/CLHEP::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() == aKaonMinus ||
      //          targetParticle.GetDefinition() == aKaonZeroL ||
      //          targetParticle.GetDefinition() == aKaonZeroS ||
      //          targetParticle.GetDefinition() == aKaonPlus ||
      //          targetParticle.GetDefinition() == aPiMinus ||
      //          targetParticle.GetDefinition() == aPiZero ||
      //          targetParticle.GetDefinition() == aPiPlus )
      //      {
      //        if( G4UniformRand() < 0.5 )
      //          targetParticle.SetDefinitionAndUpdateE( aNeutron );
      //        else
      //          targetParticle.SetDefinitionAndUpdateE( aProton );
      //        targetHasChanged = true;
      //        targetMass = targetParticle.GetMass()/CLHEP::GeV;
      //      }
      //
      // Set masses and momenta for final state particles
      //
      /*
        G4cout<<"Check 0"<<G4endl;
        G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()/CLHEP::GeV<<G4endl;
        G4cout<<"target mass: "<<targetParticle.GetMass()/CLHEP::GeV<<G4endl;
        G4cout<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
        G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()/CLHEP::GeV<<G4endl;
        G4cout<<"current mass: "<<currentParticle.GetMass()/CLHEP::GeV<<G4endl;
        G4cout<<currentParticle.GetDefinition()->GetParticleName()<<G4endl;
      */
      G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
      pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
      //      G4cout << "pf: " << pf<< G4endl;
 
      if( pf <= 0.)// 0.001 )
        {
          for(G4int i=0; i<vecLen; i++) delete vec[i];
          vecLen = 0;
          throw G4HadronicException(__FILE__, __LINE__, "FullModelReactionDynamics::TwoBody: pf is too small ");
        }
 
      pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
      //
      // Set beam and target in centre of mass system
      //
      G4ReactionProduct pseudoParticle[3];
      /* //Marker
         if(//targetParticle.GetDefinition()->GetParticleType()=="rhadron")
         targetParticle.GetDefinition() == aKaonMinus ||
         targetParticle.GetDefinition() == aKaonZeroL ||
         targetParticle.GetDefinition() == aKaonZeroS ||
         targetParticle.GetDefinition() == aKaonPlus ||
         targetParticle.GetDefinition() == aPiMinus ||
         targetParticle.GetDefinition() == aPiZero ||
         targetParticle.GetDefinition() == aPiPlus )
         {
         G4cout<<"Particlecheck"<<G4endl;
         pseudoParticle[0].SetMass( targetMass*CLHEP::GeV );
         G4cout<<pseudoParticle[0].GetMass()<<G4endl;
         pseudoParticle[0].SetTotalEnergy( etOriginal*CLHEP::GeV );
         G4cout<<pseudoParticle[0].GetTotalEnergy()<<G4endl;
         pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*CLHEP::GeV );
 
         pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
         pseudoParticle[1].SetMass( mOriginal*CLHEP::GeV );
         G4cout<<pseudoParticle[1].GetMass()<<G4endl;
         pseudoParticle[1].SetKineticEnergy( 0.0 );
         G4cout<<pseudoParticle[1].GetTotalEnergy()<<G4endl;
         }
         else
         {*/
      pseudoParticle[0].SetMass( currentMass*CLHEP::GeV );
      pseudoParticle[0].SetTotalEnergy( etCurrent*CLHEP::GeV );
      pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*CLHEP::GeV );
 
      pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
      pseudoParticle[1].SetMass( targetMass*CLHEP::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)*CLHEP::GeV );
      targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*CLHEP::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 b = std::max( cb, b1+b2*std::log(ptemp) );
      G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/CLHEP::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 = CLHEP::twopi * G4UniformRand();
      //
      // Calculate final state momenta in centre of mass system
      //
      if(//targetParticle.GetDefinition()->GetParticleType()=="rhadron")
         targetParticle.GetDefinition() == aKaonMinus ||
         targetParticle.GetDefinition() == aKaonZeroL ||
         targetParticle.GetDefinition() == aKaonZeroS ||
         targetParticle.GetDefinition() == aKaonPlus ||
         targetParticle.GetDefinition() == aPiMinus ||
         targetParticle.GetDefinition() == aPiZero ||
         targetParticle.GetDefinition() == aPiPlus )
        {
          currentParticle.SetMomentum( -pf*stet*std::sin(phi)*CLHEP::GeV,
                                       -pf*stet*std::cos(phi)*CLHEP::GeV,
                                       -pf*ctet*CLHEP::GeV );
        }
      else
        {
          currentParticle.SetMomentum( pf*stet*std::sin(phi)*CLHEP::GeV,
                                       pf*stet*std::cos(phi)*CLHEP::GeV,
                                       pf*ctet*CLHEP::GeV );
        }
      targetParticle.SetMomentum( currentParticle.GetMomentum() * (-1.0) );
      //
      // Transform into lab system
      //
      currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
      targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
 
      /*
        G4cout<<"Check 1"<<G4endl;
        G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
        G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
        G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
        G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
      */
 
      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()/CLHEP::MeV;
          if( pp1 >= 1.0 )
            {
              ekin = currentParticle.GetKineticEnergy()/CLHEP::MeV - cfa*(1.0+0.5*normal())*CLHEP::GeV;
              ekin = std::max( 0.0001*CLHEP::GeV, ekin );
              currentParticle.SetKineticEnergy( ekin*CLHEP::MeV );
              pp = currentParticle.GetTotalMomentum()/CLHEP::MeV;
              currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
            }
          pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV;
          if( pp1 >= 1.0 )
            {
              ekin = targetParticle.GetKineticEnergy()/CLHEP::MeV - cfa*(1.0+normal()/2.)*CLHEP::GeV;
              ekin = std::max( 0.0001*CLHEP::GeV, ekin );
              targetParticle.SetKineticEnergy( ekin*CLHEP::MeV );
              pp = targetParticle.GetTotalMomentum()/CLHEP::MeV;
              targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
            }
        }
    }
  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  if( atomicWeight >= 1.5 )
    {
      // Add black track particles
      //  the procedure is somewhat different than in TwoCluster and GenerateXandPt.
      //  The reason is that we have to also simulate the nuclear reactions
      //  at low energies like a(h,p)b, a(h,p p)b, a(h,n)b etc.
      //
      // 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 );
          /*
            G4cout<<"A couple of Poisson numbers:"<<G4endl;
            for (int n=0;n!=10;n++) G4cout<<Poisson(epnb/0.02)<<", ";
            G4cout<<G4endl;
          */
          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 );
        }
      G4double spall = 0.0;
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
 
      /*
        G4cout<<"Check 2"<<G4endl;
        G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
        G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
        G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
        G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
 
        G4cout<<"------------------------------------------------------------------------"<<G4endl;
        G4cout<<"Atomic weight: "<<atomicWeight<<G4endl;
        G4cout<<"number of proton/neutron black track particles: "<<npnb<<G4endl;
        G4cout<<"number of deuterons, tritons, and alphas produced: "<<ndta <<G4endl;
        G4cout<<"kinetic energy available for proton/neutron black track particles: "<<epnb/CLHEP::GeV<<" GeV" <<G4endl;
        G4cout<<"kinetic energy available for deuteron/triton/alpha particles: "<<edta/CLHEP::GeV <<" GeV"<<G4endl;
        G4cout<<"------------------------------------------------------------------------"<<G4endl;
      */
 
      AddBlackTrackParticles( epnb, npnb, edta, ndta, sprob, kineticMinimum, kineticFactor,
                              modifiedOriginal, spall, 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;
}
 
G4double FullModelReactionDynamics::GenerateNBodyEvent(
                                                       const G4double totalEnergy,                // MeV
                                                       const G4bool constantCrossSection,
                                                       G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                                       G4int &vecLen )
{
  //      // 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 FullModelReactionDynamics::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 *mass = new G4double [vecLen];    // mass of each particle
  //G4double *energy = new G4double [vecLen];  // total energy of each particle
  //G4double **pcm;           // pcm is an array with 3 rows and vecLen columns
  //pcm = new G4double * [3];
  //for( i=0; i<3; ++i )pcm[i] = new G4double [vecLen];
 
  G4double totalMass = 0.0;
  G4double sm[18];
  //G4double *sm = new G4double [vecLen];
 
  for( i=0; i<vecLen; ++i )
    {
      mass[i] = vec[i]->GetMass()/CLHEP::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/CLHEP::GeV;
  if( totalMass > totalE )
    {
      //G4cerr << "*** Error in FullModelReactionDynamics::GenerateNBodyEvent" << G4endl;
      //G4cerr << "    total mass (" << totalMass*GeV << "MeV) > total energy ("
      //     << totalEnergy << "MeV)" << G4endl;
      totalE = totalMass;
      //delete [] mass;
      //delete [] energy;
      //for( i=0; i<3; ++i )delete [] pcm[i];
      //delete [] pcm;
      //delete [] sm;
      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) = CLHEP::pi*(2*CLHEP::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, s, 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 = CLHEP::twopi*G4UniformRand();
      cb = std::cos(bang);
      sb = std::sin(bang);
      c = 2.0*G4UniformRand() - 1.0;
      s = 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*s;                           //  rotation
              pcm[1][j] = s0*s + 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*s;                           //  rotation
              pcm[1][j] = s0*s + 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]*CLHEP::GeV, pcm[1][i]*CLHEP::GeV, pcm[2][i]*CLHEP::GeV );
      vec[i]->SetTotalEnergy( energy[i]*CLHEP::GeV );
    }
  //delete [] mass;
  //delete [] energy;
  //for( i=0; i<3; ++i )delete [] pcm[i];
  //delete [] pcm;
  //delete [] emm;
  //delete [] sm;
  //delete [] pd;
  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  return weight;
}
 
G4double
FullModelReactionDynamics::normal()
{
  G4double ran = -6.0;
  for( G4int i=0; i<12; ++i )ran += G4UniformRand();
  return ran;
}
 
G4int
FullModelReactionDynamics::Poisson( G4double x )  // generation of poisson distribution
{
  G4int iran;
  G4double ran;
 
  if( x > 9.9 )    // use normal distribution with sigma^2 = <x>
    iran = static_cast<G4int>(std::max( 0.0, x+normal()*std::sqrt(x) ) );
  else {
    G4int mm = G4int(5.0*x);
    if( mm <= 0 )   // for very small x try iran=1,2,3
      {
        G4double p1 = x*std::exp(-x);
        G4double p2 = x*p1/2.0;
        G4double p3 = x*p2/3.0;
        ran = G4UniformRand();
        if( ran < p3 )
          iran = 3;
        else if( ran < p2 )   // this is original Geisha, it should be ran < p2+p3
          iran = 2;
        else if( ran < p1 )   // should be ran < p1+p2+p3
          iran = 1;
        else
          iran = 0;
      }
    else
      {
        iran = 0;
        G4double r = std::exp(-x);
        ran = G4UniformRand();
        if( ran > r )
          {
            G4double rrr;
            G4double rr = r;
            for( G4int i=1; i<=mm; ++i )
              {
                iran++;
                if( i > 5 )   // Stirling's formula for large numbers
                  rrr = std::exp(i*std::log(x)-(i+0.5)*std::log((G4double)i)+i-0.9189385);
                else
                  rrr = std::pow(x,i)/Factorial(i);
                rr += r*rrr;
                if( ran <= rr )break;
              }
          }
      }
  }
  return iran;
}
 
G4int
FullModelReactionDynamics::Factorial( G4int n )
{   // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
  G4int m = std::min(n,10);
  G4int result = 1;
  if( m <= 1 )return result;
  for( G4int i=2; i<=m; ++i )result *= i;
  return result;
}
 
void FullModelReactionDynamics::Defs1(
                                      const G4ReactionProduct &modifiedOriginal,
                                      G4ReactionProduct &currentParticle,
                                      G4ReactionProduct &targetParticle,
                                      G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                      G4int &vecLen )
{
  const G4double pjx = modifiedOriginal.GetMomentum().x()/CLHEP::MeV;
  const G4double pjy = modifiedOriginal.GetMomentum().y()/CLHEP::MeV;
  const G4double pjz = modifiedOriginal.GetMomentum().z()/CLHEP::MeV;
  const G4double p = modifiedOriginal.GetMomentum().mag()/CLHEP::MeV;
  if( pjx*pjx+pjy*pjy > 0.0 )
    {
      G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
      cost = pjz/p;
      sint = 0.5 * ( std::sqrt(std::abs((1.0-cost)*(1.0+cost))) + std::sqrt(pjx*pjx+pjy*pjy)/p );
      if( pjy < 0.0 )
        ph = 3*CLHEP::halfpi;
      else
        ph = CLHEP::halfpi;
      if( std::abs( pjx ) > 0.001*CLHEP::MeV )ph = std::atan2(pjy,pjx);
      cosp = std::cos(ph);
      sinp = std::sin(ph);
      pix = currentParticle.GetMomentum().x()/CLHEP::MeV;
      piy = currentParticle.GetMomentum().y()/CLHEP::MeV;
      piz = currentParticle.GetMomentum().z()/CLHEP::MeV;
      currentParticle.SetMomentum( cost*cosp*pix*CLHEP::MeV - sinp*piy+sint*cosp*piz*CLHEP::MeV,
                                   cost*sinp*pix*CLHEP::MeV + cosp*piy+sint*sinp*piz*CLHEP::MeV,
                                   -sint*pix*CLHEP::MeV     + cost*piz*CLHEP::MeV );
      pix = targetParticle.GetMomentum().x()/CLHEP::MeV;
      piy = targetParticle.GetMomentum().y()/CLHEP::MeV;
      piz = targetParticle.GetMomentum().z()/CLHEP::MeV;
      targetParticle.SetMomentum( cost*cosp*pix*CLHEP::MeV - sinp*piy+sint*cosp*piz*CLHEP::MeV,
                                  cost*sinp*pix*CLHEP::MeV + cosp*piy+sint*sinp*piz*CLHEP::MeV,
                                  -sint*pix*CLHEP::MeV     + cost*piz*CLHEP::MeV );
      for( G4int i=0; i<vecLen; ++i )
        {
          pix = vec[i]->GetMomentum().x()/CLHEP::MeV;
          piy = vec[i]->GetMomentum().y()/CLHEP::MeV;
          piz = vec[i]->GetMomentum().z()/CLHEP::MeV;
          vec[i]->SetMomentum( cost*cosp*pix*CLHEP::MeV - sinp*piy+sint*cosp*piz*CLHEP::MeV,
                               cost*sinp*pix*CLHEP::MeV + cosp*piy+sint*sinp*piz*CLHEP::MeV,
                               -sint*pix*CLHEP::MeV     + cost*piz*CLHEP::MeV );
        }
    }
  else
    {
      if( pjz < 0.0 )
        {
          currentParticle.SetMomentum( -currentParticle.GetMomentum().z() );
          targetParticle.SetMomentum( -targetParticle.GetMomentum().z() );
          for( G4int i=0; i<vecLen; ++i )
            vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
        }
    }
}
 
void FullModelReactionDynamics::Rotate(
                                       const G4double numberofFinalStateNucleons,
                                       const G4ThreeVector &temp,
                                       const G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effect included
                                       const G4HadProjectile *originalIncident, // original incident particle
                                       const G4Nucleus &targetNucleus,
                                       G4ReactionProduct &currentParticle,
                                       G4ReactionProduct &targetParticle,
                                       G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                       G4int &vecLen )
{
  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
  //
  //   Rotate in direction of z-axis, this does disturb in some way our
  //    inclusive distributions, but it is necessary for momentum conservation
  //
  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());//targetNucleus.GetN();
  const G4double logWeight = std::log(atomicWeight);
 
  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
 
  G4int i;
  G4ThreeVector pseudoParticle[4];
  for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
  pseudoParticle[0] = currentParticle.GetMomentum()
    + targetParticle.GetMomentum();
  for( i=0; i<vecLen; ++i )
    pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
  //
  //  Some smearing in transverse direction from Fermi motion
  //
  G4double pp, pp1;
  G4double alekw, p, rthnve, phinve;
  G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
 
  r1 = CLHEP::twopi*G4UniformRand();
  r2 = G4UniformRand();
  a1 = std::sqrt(-2.0*std::log(r2));
  ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*CLHEP::GeV;
  ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*CLHEP::GeV;
  G4ThreeVector Fermi(ran1, ran2, 0);
 
  pseudoParticle[0] = pseudoParticle[0]+Fermi; // all particles + Fermi
  pseudoParticle[2] = temp; // original in cms system
  pseudoParticle[3] = pseudoParticle[0];
 
  pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
  G4double rotation = 2.*CLHEP::pi*G4UniformRand();
  pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
  pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
  for(G4int ii=1; ii<=3; ii++)
    {
      p = pseudoParticle[ii].mag();
      if( p == 0.0 )
        pseudoParticle[ii]= G4ThreeVector( 0.0, 0.0, 0.0 );
      else
        pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
    }
 
  pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
  pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
  pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
  currentParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
 
  pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
  pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
  pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
  targetParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
 
  for( i=0; i<vecLen; ++i )
    {
      pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
      pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
      pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
      vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
    }
  //
  //  Rotate in direction of primary particle, subtract binding energies
  //   and make some further corrections if required
  //
  Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
  G4double ekin;
  G4double dekin = 0.0;
  G4double ek1 = 0.0;
  G4int npions = 0;
  if( atomicWeight >= 1.5 )            // self-absorption in heavy molecules
    {
      // corrections for single particle spectra (shower particles)
      //
      const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
      const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
      alekw = std::log( originalIncident->GetKineticEnergy()/CLHEP::GeV );
      exh = 1.0;
      if( alekw > alem[0] )   //   get energy bin
        {
          exh = val0[6];
          for( G4int j=1; j<7; ++j )
            {
              if( alekw < alem[j] ) // use linear interpolation/extrapolation
                {
                  G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
                  exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
                  break;
                }
            }
          exh = 1.0 - exh;
        }
      const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
      ekin = currentParticle.GetKineticEnergy()/CLHEP::GeV - cfa*(1+normal()/2.0);
      ekin = std::max( 1.0e-6, ekin );
      xxh = 1.0;
      /*      if( ( (modifiedOriginal.GetDefinition() == aPiPlus) ||
              (modifiedOriginal.GetDefinition() == aPiMinus) ) &&
              (currentParticle.GetDefinition() == aPiZero) &&
              (G4UniformRand() <= logWeight) )xxh = exh;*/
      dekin += ekin*(1.0-xxh);
      ekin *= xxh;
      /*      if( (currentParticle.GetDefinition() == aPiPlus) ||
              (currentParticle.GetDefinition() == aPiZero) ||
              (currentParticle.GetDefinition() == aPiMinus) )
              {
              ++npions;
              ek1 += ekin;
              }*/
      currentParticle.SetKineticEnergy( ekin*CLHEP::GeV );
      pp = currentParticle.GetTotalMomentum()/CLHEP::MeV;
      pp1 = currentParticle.GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 0.001*CLHEP::MeV )
        {
          rthnve = CLHEP::pi*G4UniformRand();
          phinve = CLHEP::twopi*G4UniformRand();
          currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                       pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                       pp*std::cos(rthnve)*CLHEP::MeV );
        }
      else
        currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
      ekin = targetParticle.GetKineticEnergy()/CLHEP::GeV - cfa*(1+normal()/2.0);
      ekin = std::max( 1.0e-6, ekin );
      xxh = 1.0;
      if( ( (modifiedOriginal.GetDefinition() == aPiPlus) ||
            (modifiedOriginal.GetDefinition() == aPiMinus) ) &&
          (targetParticle.GetDefinition() == aPiZero) &&
          (G4UniformRand() < logWeight) )xxh = exh;
      dekin += ekin*(1.0-xxh);
      ekin *= xxh;
      if( (targetParticle.GetDefinition() == aPiPlus) ||
          (targetParticle.GetDefinition() == aPiZero) ||
          (targetParticle.GetDefinition() == aPiMinus) )
        {
          ++npions;
          ek1 += ekin;
        }
      targetParticle.SetKineticEnergy( ekin*CLHEP::GeV );
      pp = targetParticle.GetTotalMomentum()/CLHEP::MeV;
      pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV;
      if( pp1 < 0.001*CLHEP::MeV )
        {
          rthnve = CLHEP::pi*G4UniformRand();
          phinve = CLHEP::twopi*G4UniformRand();
          targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                      pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                      pp*std::cos(rthnve)*CLHEP::MeV );
        }
      else
        targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
      for( i=0; i<vecLen; ++i )
        {
          ekin = vec[i]->GetKineticEnergy()/CLHEP::GeV - cfa*(1+normal()/2.0);
          ekin = std::max( 1.0e-6, ekin );
          xxh = 1.0;
          if( ( (modifiedOriginal.GetDefinition() == aPiPlus) ||
                (modifiedOriginal.GetDefinition() == aPiMinus) ) &&
              (vec[i]->GetDefinition() == aPiZero) &&
              (G4UniformRand() < logWeight) )xxh = exh;
          dekin += ekin*(1.0-xxh);
          ekin *= xxh;
          if( (vec[i]->GetDefinition() == aPiPlus) ||
              (vec[i]->GetDefinition() == aPiZero) ||
              (vec[i]->GetDefinition() == aPiMinus) )
            {
              ++npions;
              ek1 += ekin;
            }
          vec[i]->SetKineticEnergy( ekin*CLHEP::GeV );
          pp = vec[i]->GetTotalMomentum()/CLHEP::MeV;
          pp1 = vec[i]->GetMomentum().mag()/CLHEP::MeV;
          if( pp1 < 0.001*CLHEP::MeV )
            {
              rthnve = CLHEP::pi*G4UniformRand();
              phinve = CLHEP::twopi*G4UniformRand();
              vec[i]->SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                   pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                   pp*std::cos(rthnve)*CLHEP::MeV );
            }
          else
            vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
        }
    }
  if( (ek1 != 0.0) && (npions > 0) )
    {
      dekin = 1.0 + dekin/ek1;
      //
      //  first do the incident particle
      //
      if( (currentParticle.GetDefinition() == aPiPlus) ||
          (currentParticle.GetDefinition() == aPiZero) ||
          (currentParticle.GetDefinition() == aPiMinus) )
        {
          currentParticle.SetKineticEnergy(
                                           std::max( 0.001*CLHEP::MeV, dekin*currentParticle.GetKineticEnergy() ) );
          pp = currentParticle.GetTotalMomentum()/CLHEP::MeV;
          pp1 = currentParticle.GetMomentum().mag()/CLHEP::MeV;
          if( pp1 < 0.001 )
            {
              rthnve = CLHEP::pi*G4UniformRand();
              phinve = CLHEP::twopi*G4UniformRand();
              currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                           pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                           pp*std::cos(rthnve)*CLHEP::MeV );
            }
          else
            currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
        }
      /*      if( (targetParticle.GetDefinition() == aPiPlus) ||
              (targetParticle.GetDefinition() == aPiZero) ||
              (targetParticle.GetDefinition() == aPiMinus) )
              {
              targetParticle.SetKineticEnergy(
              std::max( 0.001*CLHEP::MeV, dekin*targetParticle.GetKineticEnergy() ) );
              pp = targetParticle.GetTotalMomentum()/CLHEP::MeV;
              pp1 = targetParticle.GetMomentum().mag()/CLHEP::MeV;
              if( pp1 < 0.001 )
              {
              rthnve = CLHEP::pi*G4UniformRand();
              phinve = CLHEP::twopi*G4UniformRand();
              targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
              pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
              pp*std::cos(rthnve)*CLHEP::MeV );
              }
              else
              targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
              }*/
      for( i=0; i<vecLen; ++i )
        {
          if( (vec[i]->GetDefinition() == aPiPlus) ||
              (vec[i]->GetDefinition() == aPiZero) ||
              (vec[i]->GetDefinition() == aPiMinus) )
            {
              vec[i]->SetKineticEnergy( std::max( 0.001*CLHEP::MeV, dekin*vec[i]->GetKineticEnergy() ) );
              pp = vec[i]->GetTotalMomentum()/CLHEP::MeV;
              pp1 = vec[i]->GetMomentum().mag()/CLHEP::MeV;
              if( pp1 < 0.001 )
                {
                  rthnve = CLHEP::pi*G4UniformRand();
                  phinve = CLHEP::twopi*G4UniformRand();
                  vec[i]->SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*CLHEP::MeV,
                                       pp*std::sin(rthnve)*std::sin(phinve)*CLHEP::MeV,
                                       pp*std::cos(rthnve)*CLHEP::MeV );
                }
              else
                vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
            }
        }
    }
}
 
void FullModelReactionDynamics::AddBlackTrackParticles(
                                                       const G4double epnb,            // GeV
                                                       const G4int npnb,
                                                       const G4double edta,            // GeV
                                                       const G4int ndta,
                                                       const G4double sprob,
                                                       const G4double kineticMinimum,  // GeV
                                                       const G4double kineticFactor,   // GeV
                                                       const G4ReactionProduct &modifiedOriginal,
                                                       G4double spall,
                                                       const G4Nucleus &targetNucleus,
                                                       G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                                       G4int &vecLen )
{
  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
  //
  // 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
  //
 
  G4ParticleDefinition *aProton = G4Proton::Proton();
  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
  G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
  G4ParticleDefinition *aTriton = G4Triton::Triton();
  G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
 
  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/CLHEP::MeV;
  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());
  //std::cout<<"atomicWeight "<<atomicWeight<<" and GetN_asInt "<<targetNucleus.GetN_asInt()<<" and GetA_asInt "<<targetNucleus.GetA_asInt()<<" for GetN "<<targetNucleus.GetN()<<std::endl;
  const G4double atomicNumber = targetNucleus.GetZ_asInt();
  //std::cout<<"atomicNumber "<<atomicNumber<<" for GetZ "<<targetNucleus.GetZ()<<std::endl;
 
  const G4double ika1 = 3.6;
  const G4double ika2 = 35.56;
  const G4double ika3 = 6.45;
  const G4double sp1 = 1.066;
 
  G4int i;
  G4double pp;
  // G4double totalQ = 0;
  // G4double kinCreated = 0;
  G4double cfa = 0.025*((atomicWeight-1.0)/120.0) * std::exp(-(atomicWeight-1.0)/120.0);
  if( npnb > 0)  // first add protons and neutrons
    {
      G4double backwardKinetic = 0.0;
      G4int local_npnb = npnb;
      for( i=0; i<npnb; ++i ) if( G4UniformRand() < sprob ) local_npnb--;
      G4double ekin = epnb/local_npnb;
 
      for( i=0; i<local_npnb; ++i )
        {
          G4ReactionProduct * p1 = new G4ReactionProduct();
          if( backwardKinetic > epnb )
            {
              delete p1;
              break;
            }
          G4double ran = G4UniformRand();
          G4double kinetic = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
          if( kinetic < 0.0 )kinetic = -0.010*std::log(ran);
          backwardKinetic += kinetic;
          if( backwardKinetic > epnb )
            kinetic = std::max( kineticMinimum, epnb-(backwardKinetic-kinetic) );
          if( G4UniformRand() > (1.0-atomicNumber/atomicWeight) )
            p1->SetDefinition( aProton );
          else
            p1->SetDefinition( aNeutron );
          vec.SetElement( vecLen, p1 );
          ++spall;
          G4double cost = G4UniformRand() * 2.0 - 1.0;
          G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
          G4double phi = CLHEP::twopi * G4UniformRand();
          vec[vecLen]->SetNewlyAdded( true );
          vec[vecLen]->SetKineticEnergy( kinetic*CLHEP::GeV );
          // kinCreated+=kinetic;
          pp = vec[vecLen]->GetTotalMomentum()/CLHEP::MeV;
          vec[vecLen]->SetMomentum( pp*sint*std::sin(phi)*CLHEP::MeV,
                                    pp*sint*std::cos(phi)*CLHEP::MeV,
                                    pp*cost*CLHEP::MeV );
          vecLen++;
          // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
      if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*CLHEP::GeV) )
        {
          G4double ekw = ekOriginal/CLHEP::GeV;
          G4int ika, kk = 0;
          if( ekw > 1.0 )ekw *= ekw;
          ekw = std::max( 0.1, ekw );
          ika = G4int(ika1*std::exp((atomicNumber*atomicNumber/atomicWeight-ika2)/ika3)/ekw);
          if( ika > 0 )
            {
              for( i=(vecLen-1); i>=0; --i )
                {
                  if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
                    {
                      vec[i]->SetDefinitionAndUpdateE( aNeutron );  // modified 22-Oct-97
                      if( ++kk > ika )break;
                    }
                }
            }
        }
    }
  if( ndta > 0 )    //  now, try to add deuterons, tritons and alphas
    {
      G4double backwardKinetic = 0.0;
      G4int local_ndta=ndta;
      for( i=0; i<ndta; ++i )if( G4UniformRand() < sprob )local_ndta--;
      G4double ekin = edta/local_ndta;
 
      for( i=0; i<local_ndta; ++i )
        {
          G4ReactionProduct *p2 = new G4ReactionProduct();
          if( backwardKinetic > edta )
            {
              delete p2;
              break;
            }
          G4double ran = G4UniformRand();
          G4double kinetic = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
          if( kinetic < 0.0 )kinetic = kineticFactor*std::log(ran);
          backwardKinetic += kinetic;
          if( backwardKinetic > edta )kinetic = edta-(backwardKinetic-kinetic);
          if( kinetic < 0.0 )kinetic = kineticMinimum;
          G4double cost = 2.0*G4UniformRand() - 1.0;
          G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
          G4double phi = CLHEP::twopi*G4UniformRand();
          ran = G4UniformRand();
          if( ran <= 0.60 )
            p2->SetDefinition( aDeuteron );
          else if( ran <= 0.90 )
            p2->SetDefinition( aTriton );
          else
            p2->SetDefinition( anAlpha );
          spall += p2->GetMass()/CLHEP::GeV * sp1;
          if( spall > atomicWeight )
            {
              delete p2;
              break;
            }
          vec.SetElement( vecLen, p2 );
          vec[vecLen]->SetNewlyAdded( true );
          vec[vecLen]->SetKineticEnergy( kinetic*CLHEP::GeV );
          // kinCreated+=kinetic;
          pp = vec[vecLen]->GetTotalMomentum()/CLHEP::MeV;
          vec[vecLen++]->SetMomentum( pp*sint*std::sin(phi)*CLHEP::MeV,
                                      pp*sint*std::cos(phi)*CLHEP::MeV,
                                      pp*cost*CLHEP::MeV );
          // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
        }
    }
   // G4double delta = epnb+edta - kinCreated;
}
 
void FullModelReactionDynamics::MomentumCheck(
                                              const G4ReactionProduct &modifiedOriginal,
                                              G4ReactionProduct &currentParticle,
                                              G4ReactionProduct &targetParticle,
                                              G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                              G4int &vecLen )
{
  const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/CLHEP::MeV;
  G4double testMomentum = currentParticle.GetMomentum().mag()/CLHEP::MeV;
  G4double pMass;
  if( testMomentum >= pOriginal )
    {
      pMass = currentParticle.GetMass()/CLHEP::MeV;
      currentParticle.SetTotalEnergy(
                                     std::sqrt( pMass*pMass + pOriginal*pOriginal )*CLHEP::MeV );
      currentParticle.SetMomentum(
                                  currentParticle.GetMomentum() * (pOriginal/testMomentum) );
    }
  testMomentum = targetParticle.GetMomentum().mag()/CLHEP::MeV;
  if( testMomentum >= pOriginal )
    {
      pMass = targetParticle.GetMass()/CLHEP::MeV;
      targetParticle.SetTotalEnergy(
                                    std::sqrt( pMass*pMass + pOriginal*pOriginal )*CLHEP::MeV );
      targetParticle.SetMomentum(
                                 targetParticle.GetMomentum() * (pOriginal/testMomentum) );
    }
  for( G4int i=0; i<vecLen; ++i )
    {
      testMomentum = vec[i]->GetMomentum().mag()/CLHEP::MeV;
      if( testMomentum >= pOriginal )
        {
          pMass = vec[i]->GetMass()/CLHEP::MeV;
          vec[i]->SetTotalEnergy(
                                 std::sqrt( pMass*pMass + pOriginal*pOriginal )*CLHEP::MeV );
          vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
        }
    }
}
 
void FullModelReactionDynamics::ProduceStrangeParticlePairs(
                                                            G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> &vec,
                                                            G4int &vecLen,
                                                            const G4ReactionProduct &modifiedOriginal,
                                                            const G4DynamicParticle *originalTarget,
                                                            G4ReactionProduct &currentParticle,
                                                            G4ReactionProduct &targetParticle,
                                                            G4bool &incidentHasChanged,
                                                            G4bool &targetHasChanged )
{
  // 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()/CLHEP::GeV;
  const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/CLHEP::GeV;
  G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/CLHEP::GeV;
  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
                                           targetMass*targetMass +
                                           2.0*targetMass*etOriginal );  // GeV
  G4double currentMass = currentParticle.GetMass()/CLHEP::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()/CLHEP::GeV;
  const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/CLHEP::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
    {
      G4double 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()/CLHEP::GeV;
  targetMass = targetParticle.GetMass()/CLHEP::GeV;
 
  G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
  for( i=0; i<vecLen; ++i )
    {
      energyCheck -= vec[i]->GetMass()/CLHEP::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
FullModelReactionDynamics::NuclearReaction(
                                           G4FastVector<G4ReactionProduct,4> &vec,
                                           G4int &vecLen,
                                           const G4HadProjectile *originalIncident,
                                           const G4Nucleus &targetNucleus,
                                           const G4double theAtomicMass,
                                           const G4double *mass )
{
  // 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()/CLHEP::MeV;
  const G4double aNeutronMass = aNeutron->GetPDGMass()/CLHEP::MeV;
  const G4double aDeuteronMass = aDeuteron->GetPDGMass()/CLHEP::MeV;
  const G4double aTritonMass = aTriton->GetPDGMass()/CLHEP::MeV;
  const G4double anAlphaMass = anAlpha->GetPDGMass()/CLHEP::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*CLHEP::MeV )
    {
      G4double phinve = CLHEP::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()/CLHEP::MeV;
  G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/CLHEP::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 )
    {
      // atomic weight
      const G4double A = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());//targetNucleus.GetN();
      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*CLHEP::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__,
                                "FullModelReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
    }
  G4double ke = currentParticle.GetKineticEnergy()/CLHEP::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]*CLHEP::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*CLHEP::MeV );
  pseudo1.SetTotalEnergy( theAtomicMass*CLHEP::MeV );
  G4ReactionProduct pseudo2 = currentParticle + pseudo1;
  pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
  //
  // use phase space routine in centre of mass system
  //
  G4FastVector<G4ReactionProduct,MYGHADLISTSIZE> 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()/CLHEP::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()/CLHEP::MeV - aProtonMass < 0.1 )
    {
      v[0]->SetDefinition( aProton );
      particleIsDefined = true;
    }
  else if( v[0]->GetMass()/CLHEP::MeV - aNeutronMass < 0.1 )
    {
      v[0]->SetDefinition( aNeutron );
      particleIsDefined = true;
    }
  else if( v[0]->GetMass()/CLHEP::MeV - aDeuteronMass < 0.1 )
    {
      v[0]->SetDefinition( aDeuteron );
      particleIsDefined = true;
    }
  else if( v[0]->GetMass()/CLHEP::MeV - aTritonMass < 0.1 )
    {
      v[0]->SetDefinition( aTriton );
      particleIsDefined = true;
    }
  else if( v[0]->GetMass()/CLHEP::MeV - anAlphaMass < 0.1 )
    {
      v[0]->SetDefinition( anAlpha );
      particleIsDefined = true;
    }
  currentParticle.SetKineticEnergy(
                                   std::max( 0.001, currentParticle.GetKineticEnergy()/CLHEP::MeV ) );
  p = currentParticle.GetTotalMomentum();
  pp = currentParticle.GetMomentum().mag();
  if( pp <= 0.001*CLHEP::MeV )
    {
      G4double phinve = CLHEP::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()/CLHEP::MeV ) );
      p = v[0]->GetTotalMomentum();
      pp = v[0]->GetMomentum().mag();
      if( pp <= 0.001*CLHEP::MeV )
        {
          G4double phinve = CLHEP::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()/CLHEP::MeV ) );
  else
    v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/CLHEP::MeV ) );
 
  p = v[1]->GetTotalMomentum();
  pp = v[1]->GetMomentum().mag();
  if( pp <= 0.001*CLHEP::MeV )
    {
      G4double phinve = CLHEP::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()/CLHEP::MeV ) );
      else
        v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/CLHEP::MeV ) );
 
      p = v[2]->GetTotalMomentum();
      pp = v[2]->GetMomentum().mag();
      if( pp <= 0.001*CLHEP::MeV )
        {
          G4double phinve = CLHEP::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;
}
 
/* end of file */