//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * 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.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// $Id: G4RPGTwoCluster.cc,v 1.5 2008/06/09 18:13:35 dennis Exp $
// GEANT4 tag $Name: geant4-09-02 $
//

#include "G4RPGTwoCluster.hh"
#include "Randomize.hh"
#include "G4Poisson.hh"
#include <iostream>
#include "G4HadReentrentException.hh"
#include <signal.h>


G4RPGTwoCluster::G4RPGTwoCluster()
  : G4RPGReaction() {}


G4bool G4RPGTwoCluster::
ReactionStage(const G4HadProjectile* originalIncident,
              G4ReactionProduct& modifiedOriginal,
              G4bool& incidentHasChanged,
              const G4DynamicParticle* originalTarget,
              G4ReactionProduct& targetParticle,
              G4bool& targetHasChanged,
              const G4Nucleus& targetNucleus,
              G4ReactionProduct& currentParticle,
              G4FastVector<G4ReactionProduct,256>& vec,
              G4int& vecLen,
              G4bool leadFlag,
              G4ReactionProduct& leadingStrangeParticle)
{
  // Derived from H. Fesefeldt's FORTRAN code TWOCLU
  //
  // A simple two cluster model is used to generate x- and pt- values for 
  // incident, target, and all secondary particles. 
  // This should be sufficient for low energy interactions.
  //

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

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

  if( currentMass == 0.0 && targetMass == 0.0 )
  {
    G4double ek = currentParticle.GetKineticEnergy();
    G4ThreeVector 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;
      throw G4HadReentrentException(__FILE__, __LINE__,
      "G4RPGTwoCluster::ReactionStage : Negative number of particles");
    }
    delete vec[vecLen-1];
    delete vec[vecLen-2];
    vecLen -= 2;
    currentMass = currentParticle.GetMass()/GeV;
    targetMass = targetParticle.GetMass()/GeV;
    incidentHasChanged = true;
    targetHasChanged = true;
    currentParticle.SetKineticEnergy( ek );
    currentParticle.SetMomentum( m );
    veryForward = true;
  }

  const G4double atomicWeight = targetNucleus.GetN();
  const G4double atomicNumber = targetNucleus.GetZ();
  //
  // particles have been distributed in forward and backward hemispheres
  // in center of mass system of the hadron nucleon interaction
  //

  // Incident particle always in forward hemisphere

  G4int forwardCount = 1;        // number of particles in forward hemisphere
  currentParticle.SetSide( 1 );
  G4double forwardMass = currentParticle.GetMass()/GeV;
  G4double cMass = forwardMass;
    
  // Target particle always in backward hemisphere

  G4int backwardCount = 1;       // number of particles in backward hemisphere
  targetParticle.SetSide( -1 );
  G4double backwardMass = targetParticle.GetMass()/GeV;
  G4double bMass = backwardMass;

  //  G4int backwardNucleonCount = 1;  // number of nucleons in backward hemisphere
    
  for( i=0; i<vecLen; ++i )
  {
    if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 );   // to take care of 
    // case where vec has been preprocessed by GenerateXandPt
    // and some of them have been set to -2 or -3
    if( vec[i]->GetSide() == -1 )
    {
      ++backwardCount;
      backwardMass += vec[i]->GetMass()/GeV;
    }
    else
    {
      ++forwardCount;
      forwardMass += vec[i]->GetMass()/GeV;
    }
  }

  //
  // Add nucleons and some pions from intra-nuclear 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 = G4Poisson( xtarg );

  if(atomicWeight<1.0001) nuclearExcitationCount = 0;
  //  G4int extraNucleonCount = 0;
  //  G4double extraMass = 0.0;
  //  G4double extraNucleonMass = 0.0;
  if( nuclearExcitationCount > 0 )
  {
    G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );     
    const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
    //
    //  NOTE: in TWOCLU, these new particles were given negative codes
    //        here we use  NewlyAdded = true  instead
    //
    for( i=0; i<nuclearExcitationCount; ++i )
    {
      G4ReactionProduct* pVec = new G4ReactionProduct();
      if( G4UniformRand() < nucsup[momentumBin] )  // add proton or neutron
      {
        if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
          pVec->SetDefinition( aProton );
        else
          pVec->SetDefinition( aNeutron );
	// Not used        ++backwardNucleonCount;
	// Not used        ++extraNucleonCount;
	// Not used        extraNucleonMass += pVec->GetMass()/GeV;
      }
      else
      {                                       // add a pion
        G4double ran = G4UniformRand();
        if( ran < 0.3181 )
          pVec->SetDefinition( aPiPlus );
        else if( ran < 0.6819 )
          pVec->SetDefinition( aPiZero );
        else
          pVec->SetDefinition( aPiMinus );

	// DHW: add following two lines to correct energy balance
	//        ++backwardCount;
	//        backwardMass += pVec->GetMass()/GeV;
      }
      pVec->SetSide( -2 );    // backside particle
      // Not used     extraMass += pVec->GetMass()/GeV;
      pVec->SetNewlyAdded( true );
      vec.SetElement( vecLen++, pVec );
    }
  }

  // Masses of particles added from cascade not included in energy balance.
  // That's correct for nucleons from the intra-nuclear cascade but not for 
  // pions from the cascade.
 
  G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
  G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
  G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
  G4bool secondaryDeleted;
  G4double pMass;

  while( eAvailable <= 0.0 )   // must eliminate a particle
  {
    secondaryDeleted = false;
    for( i=(vecLen-1); i>=0; --i )
    {
      if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
      {
        pMass = vec[i]->GetMass()/GeV;
        for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];     // shift up
        --forwardCount;
        forwardEnergy += pMass;
        forwardMass -= pMass;
        secondaryDeleted = true;
        break;
      }
      else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
      {
        pMass = vec[i]->GetMass()/GeV;
        for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
        --backwardCount;
        backwardEnergy += pMass;
        backwardMass -= pMass;
        secondaryDeleted = true;
        break;
      }
    }  // breaks go down to here

    if( secondaryDeleted )
    {
      delete vec[vecLen-1];
      --vecLen;
        // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
    }
    else
    {
      if( vecLen == 0 ) return false;  // all secondaries have been eliminated
      if( targetParticle.GetSide() == -1 )
      {
        pMass = targetParticle.GetMass()/GeV;
        targetParticle = *vec[0];
        for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
        --backwardCount;
        backwardEnergy += pMass;
        backwardMass -= pMass;
        secondaryDeleted = true;
      }
      else if( targetParticle.GetSide() == 1 )
      {
        pMass = targetParticle.GetMass()/GeV;
        targetParticle = *vec[0];
        for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
        --forwardCount;
        forwardEnergy += pMass;
        forwardMass -= pMass;
        secondaryDeleted = true;
      }

      if( secondaryDeleted )
      {
        delete vec[vecLen-1];
        --vecLen;
      }
      else
      {
        if( currentParticle.GetSide() == -1 )
        {
          pMass = currentParticle.GetMass()/GeV;
          currentParticle = *vec[0];
          for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
          --backwardCount;
          backwardEnergy += pMass;
          backwardMass -= pMass;
          secondaryDeleted = true;
        }
        else if( currentParticle.GetSide() == 1 )
        {
          pMass = currentParticle.GetMass()/GeV;
          currentParticle = *vec[0];
          for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
          --forwardCount;
          forwardEnergy += pMass;
          forwardMass -= pMass;
          secondaryDeleted = true;
        }
        if( secondaryDeleted )
        {
          delete vec[vecLen-1];
          --vecLen;
        }
        else break;

      }  // secondary not deleted 
    }  // secondary not deleted

    eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
  } // while

  //
  // This is the start of the TwoCluster function
  // Choose multi-particle resonance masses by sampling 
  //    P(M) = gc[g(M-M0)]**(c-1) *exp[-(g(M-M0))**c] 
  // for M > M0
  //
  // Use for the forward and backward clusters, but not
  // the cascade cluster

  const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
  const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
  G4int ntc = 0;

  if (forwardCount < 1 || backwardCount < 1) return false;  // array bounds protection

  G4double rmc = forwardMass;
  if (forwardCount > 1) {
    ntc = std::min(3,forwardCount-2);
    rmc += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
  }
  G4double rmd = backwardMass;
  if( backwardCount > 1 ) {
    ntc = std::min(3,backwardCount-2);
    rmd += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
  }

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

  G4ReactionProduct pseudoParticle[8];
  for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
    
  pseudoParticle[1].SetMass( mOriginal*GeV );
  pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
  pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
    
  pseudoParticle[2].SetMass( protonMass*MeV );
  pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
  pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
  //
  //  transform into center of mass system
  //
  pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
  pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
  pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );

  // Calculate cm momentum for forward and backward masses
  // W = sqrt(pf*pf + rmc*rmc) + sqrt(pf*pf + rmd*rmd)
  // Solve for pf

  const G4double pfMin = 0.0001;
  G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
  pf *= pf;
  pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
  pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
  //
  //  set final state masses and energies in centre of mass system
  //
  pseudoParticle[3].SetMass( rmc*GeV );
  pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
    
  pseudoParticle[4].SetMass( rmd*GeV );
  pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );

  //
  // Get cm scattering angle by sampling t from tmin to tmax
  //
  const G4double bMin = 0.01;
  const G4double b1 = 4.0;
  const G4double b2 = 1.6;
  G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
  G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
  G4double factor = 1.0 - std::exp(-dtb);
  G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;

  costheta = std::max(-1.0, std::min(1.0, costheta));
  G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
  G4double phi = G4UniformRand() * twopi;
  //
  // calculate final state momenta in centre of mass system
  //
  pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
                                 pf*sintheta*std::sin(phi)*GeV,
                                 pf*costheta*GeV );
  pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());

  // Backward cluster of nucleons and pions from intra-nuclear cascade
  // Set up in lab system and transform to cms

  G4double pp, pp1;
  if( nuclearExcitationCount > 0 )
  {
    const G4double ga = 1.2;
    G4double ekit1 = 0.04;
    G4double ekit2 = 0.6;   // Max KE of cascade particle
    if( ekOriginal <= 5.0 )
    {
      ekit1 *= ekOriginal*ekOriginal/25.0;
      ekit2 *= ekOriginal*ekOriginal/25.0;
    }
    G4double scale = std::pow(ekit2/ekit1, 1.0-ga) - 1.0;
    for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetSide() == -2 )
      {
        G4double kineticE = ekit1*std::pow((1.0 + G4UniformRand()*scale), 1.0/(1.0-ga) );
        vec[i]->SetKineticEnergy( kineticE*GeV );
        G4double vMass = vec[i]->GetMass()/MeV;
        G4double totalE = kineticE*GeV + vMass;
        pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
        G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
        G4double sint = std::sqrt(1.0-cost*cost);
        phi = twopi*G4UniformRand();
        vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
                             pp*sint*std::sin(phi)*MeV,
                             pp*cost*MeV );
        vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
      }
    }
  }

  //
  // Fragmentation of forward and backward clusters
  //

  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,256> tempV;
    tempV.Initialize( forwardCount );
    G4bool constantCrossSection = true;
    G4int tempLen = 0;
    if( currentParticle.GetSide() == 1 )
      tempV.SetElement( tempLen++, &currentParticle );
    if( targetParticle.GetSide() == 1 )
      tempV.SetElement( tempLen++, &targetParticle );
    for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetSide() == 1 )
      {
        if( tempLen < 18 )
          tempV.SetElement( tempLen++, vec[i] );
        else
        {
          vec[i]->SetSide( -1 );
          continue;
        }
      }
    }
    if( tempLen >= 2 )
    {
      wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
                                constantCrossSection, tempV, tempLen );
      if( currentParticle.GetSide() == 1 )
        currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
      if( targetParticle.GetSide() == 1 )
        targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
      }
    }
  }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  if( backwardCount > 1 )   //  tempV will contain the backward particles,
  {                         //  but not those created from the intranuclear cascade
    G4FastVector<G4ReactionProduct,256> tempV;
    tempV.Initialize( backwardCount );
    G4bool constantCrossSection = true;
    G4int tempLen = 0;
    if( currentParticle.GetSide() == -1 )
      tempV.SetElement( tempLen++, &currentParticle );
    if( targetParticle.GetSide() == -1 )
      tempV.SetElement( tempLen++, &targetParticle );
    for( i=0; i<vecLen; ++i )
    {
      if( vec[i]->GetSide() == -1 )
      {
        if( tempLen < 18 )
          tempV.SetElement( tempLen++, vec[i] );
        else
        {
          vec[i]->SetSide( -2 );
          vec[i]->SetKineticEnergy( 0.0 );
          vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
          continue;
        }
      }
    }
    if( tempLen >= 2 )
    {
      wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
                                constantCrossSection, tempV, tempLen );
      if( currentParticle.GetSide() == -1 )
        currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
      if( targetParticle.GetSide() == -1 )
        targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
      }
    }
  }

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  //
  // Lorentz transformation in lab system
  //
  currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
  targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
  for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );

      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  //
  // sometimes the leading strange particle is lost, set it back
  //
  G4bool dum = true;
  if( leadFlag )
  {
    // leadFlag will be true
    //  iff original particle is strange AND if incident particle is strange
    //  leadFlag is set to the incident particle
    //  or
    //  target particle is strange leadFlag is set to the target particle

    if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
      dum = false;
    else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
      dum = false;
    else
    {
      for( i=0; i<vecLen; ++i )
      {
        if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
        {
          dum = false;
          break;
        }
      }
    }
    if( dum )
    {
      G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
      G4double ekin;
      if( ((leadMass <  protonMass) && (targetParticle.GetMass()/MeV <  protonMass)) ||
          ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
      {
        ekin = targetParticle.GetKineticEnergy()/GeV;
        pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
        targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
        targetParticle.SetKineticEnergy( ekin*GeV );
        pp = targetParticle.GetTotalMomentum()/MeV;   // new momentum
        if( pp1 < 1.0e-3 ) {
          G4ThreeVector iso = Isotropic(pp);
          targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
        } else {
          targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
	}
        targetHasChanged = true;
      }
      else
      {
        ekin = currentParticle.GetKineticEnergy()/GeV;
        pp1 = currentParticle.GetMomentum().mag()/MeV;
        currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
        currentParticle.SetKineticEnergy( ekin*GeV );
        pp = currentParticle.GetTotalMomentum()/MeV;
        if( pp1 < 1.0e-3 ) {
          G4ThreeVector iso = Isotropic(pp);
          currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
        } else {
          currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
	}
        incidentHasChanged = true;
      }
    }
  }    // end of if( leadFlag )

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

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

  G4int numberofFinalStateNucleons = 
    protonsInFinalState + neutronsInFinalState;

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

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

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

  if( vecLen < 16 )
  {
    G4ReactionProduct tempR[130];
    tempR[0] = currentParticle;
    tempR[1] = targetParticle;
    for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];

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

    if( tempLen >= 2 )
    {
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
      wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
                                pseudoParticle[5].GetTotalEnergy()/MeV,
                                constantCrossSection, tempV, tempLen );
      if (wgt == -1) {
        G4double Qvalue = 0;
        for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
        wgt = GenerateNBodyEvent( Qvalue/MeV,
                                  constantCrossSection, tempV, tempLen );
      }
      theoreticalKinetic = 0.0;
      for( i=0; i<vecLen+2; ++i )
      {
        pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
        pseudoParticle[7].SetMass( tempV[i]->GetMass() );
        pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
        pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
        theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
      }
    }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  }
  else
  {
    theoreticalKinetic -=
      ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
    for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
  }
  G4double simulatedKinetic =
    currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
  for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;

  // make sure that kinetic energies are correct
  // the backward nucleon cluster is not produced within proper kinematics!!!
    
  if( simulatedKinetic != 0.0 )
  {
    wgt = (theoreticalKinetic)/simulatedKinetic;
    currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
    pp = currentParticle.GetTotalMomentum()/MeV;
    pp1 = currentParticle.GetMomentum().mag()/MeV;
    if( pp1 < 0.001*MeV ) {
      G4ThreeVector iso = Isotropic(pp);
      currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
    } else {
      currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
    }

    targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
    pp = targetParticle.GetTotalMomentum()/MeV;
    pp1 = targetParticle.GetMomentum().mag()/MeV;
    if( pp1 < 0.001*MeV ) {
      G4ThreeVector iso = Isotropic(pp);
      targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
    } else {
      targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
    }

    for( i=0; i<vecLen; ++i )
    {
      vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
      pp = vec[i]->GetTotalMomentum()/MeV;
      pp1 = vec[i]->GetMomentum().mag()/MeV;
      if( pp1 < 0.001 ) {
        G4ThreeVector iso = Isotropic(pp);
        vec[i]->SetMomentum( iso.x(), iso.y(), iso.z() );
      } else {
        vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
      }
    }
  }
      // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);

  Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
          modifiedOriginal, originalIncident, targetNucleus,
          currentParticle, targetParticle, vec, vecLen );

  //  Add black track particles
  //  the total number of particles produced is restricted to 198
  //  this may have influence on very high energies

  if( atomicWeight >= 1.5 )
  {
    // npnb is number of proton/neutron black track particles
    // ndta is the number of deuterons, tritons, and alphas produced
    // epnb is the kinetic energy available for proton/neutron black track 
    //   particles
    // edta is the kinetic energy available for deuteron/triton/alpha 
    //   particles

    G4int npnb = 0;
    G4int ndta = 0;

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

    const G4double pnCutOff = 0.001;     // GeV
    const G4double dtaCutOff = 0.001;    // GeV
    //    const G4double kineticMinimum = 1.e-6;
    //    const G4double kineticFactor = -0.005;
      
    //    G4double sprob = 0.0;   // sprob = probability of self-absorption in 
                            // heavy molecules
    // Not currently used (DHW 9 June 2008)  const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
    //    if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
      
    if( epnb >= pnCutOff )
    {
      npnb = G4Poisson((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 = G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
      ndta = std::min( ndta, 127-vecLen );
    }
    if (npnb == 0 && ndta == 0) npnb = 1;

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

    AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal, 
                           PinNucleus, NinNucleus, targetNucleus,
                           vec, vecLen );
    // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
  }

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

  return true;
}
 
 /* end of file */