//
// ********************************************************************
// * 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: G4DNAProcess.icc,v 1.12 2009/01/20 07:50:28 sincerti Exp $
// GEANT4 tag $Name: geant4-09-04-ref-00 $
// 
// Contact Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
//
// Reference: TNS Geant4-DNA paper
//

// History:
// -----------
// Date         Name              Modification
// 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
//
// -------------------------------------------------------------------


template <class TCrossSection,class TFinalState>
G4double G4DNAProcess<TCrossSection,TFinalState>::GetMeanFreePath(const G4Track& track, 
								  G4double /* previousStepSize */, 
								  G4ForceCondition* /* condition */)
{ 
  G4double meanFreePath = DBL_MAX;

  // Assume the interacting medium to be water; one of the elements must be oxygen
  G4Material* material(track.GetMaterial());
  size_t i = material->GetNumberOfElements();
  while (i>0)
    {
      i--;
      const G4Element* element(material->GetElement(i));
      if (element->GetZ() == 8.)
	{
	  // Number of oxygen atoms per volume = number of water molecules per volume
	  G4double density = material->GetAtomicNumDensityVector()[i];
	  // G4cout << "density = " << density << G4endl;
	  if (density > 0.)
	    {
	      G4double cross = crossSection.CrossSection(track);
	      if (cross > 0.0) meanFreePath = 1. / (density*cross);
              if (meanFreePath == 0.) meanFreePath = DBL_MIN;
	      return meanFreePath;
	    }
	}
    } // end while
  
  // If it ends up here, it means that the material is not water
  G4Exception("G4DNAProcess::GetMeanFreePath - material is not water");
  // One does not really need a return statement here
  return DBL_MAX;
}

								  
template <class TCrossSection,class TFinalState>
G4VParticleChange* G4DNAProcess<TCrossSection,TFinalState>::PostStepDoIt(const G4Track& track, const G4Step& step)
{
  aParticleChange.Initialize(track);
  
  // G4cout << "Track initialized" << G4endl;

  // Interaction product
  const G4FinalStateProduct& product = finalState.GenerateFinalState(track,step);

  // Number of secondary products to be generated
  G4int nSecondaries  = product.NumberOfSecondaries();
  aParticleChange.SetNumberOfSecondaries(nSecondaries);
  
  // Secondaries
  for (G4int l = 0; l<nSecondaries; l++ )
    {
      G4DynamicParticle* particle = product.GetSecondaries()[l];
      if (particle != 0) 
	{
//	  aParticleChange.SetNumberOfSecondaries(nSecondaries);
	  aParticleChange.AddSecondary(particle);
	}
    }
  
  // Take care of incident particle to be killed, if necessary; dump its energy deposit locally
  G4double deposit = product.GetEnergyDeposit();
  if (deposit > 0.0) aParticleChange.ProposeLocalEnergyDeposit(deposit);
  
  if (product.PrimaryParticleIsKilled())		
    {
      aParticleChange.ProposeTrackStatus(fStopAndKill);
      aParticleChange.ProposeEnergy(0.);
      aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
      
      if (product.PrimaryParticleIsKilledAndDoNotDepositEnergy())
      {
        aParticleChange.ProposeLocalEnergyDeposit(deposit);
      }
      else
      {
        aParticleChange.ProposeLocalEnergyDeposit(track.GetKineticEnergy() + deposit);
      }
      
    }
  else
    {
      // Modify incident particle kinematics taking into account the generated products
      
      // ---- MGP ---- Temporary: assume at most one secondary product
      //               Sebastien: please check if consistent with current models or generalize
      
      // Primary particle momentum and kinetic energy
      G4ThreeVector primaryMomentum = track.GetMomentum();
      G4double primaryKineticEnergy = track.GetKineticEnergy();
      
      // Secondary product momentum and energy
      
      G4double secondaryKineticEnergy = 0.;
      if (nSecondaries >0 )
	{
	  G4DynamicParticle* secondary = product.GetSecondaries()[0];
	  secondaryKineticEnergy = secondary->GetKineticEnergy();
	  
	  // Calculate new primary particle kinetic energy
	  G4double finalKineticEnergy = primaryKineticEnergy - secondaryKineticEnergy - deposit;
	  
	  if (finalKineticEnergy <= 0.0) 
	    {
	      // Primary particle is stopped; kill it
	      aParticleChange.ProposeTrackStatus(fStopAndKill);
	      aParticleChange.ProposeEnergy(0.);
	      aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
	    } 
	  else 
	    { 
	      // Calculate new primary particle momentum: difference between original primary one and secondary 
	      G4ThreeVector secondaryMomentum = secondary->GetMomentum();
	      G4ThreeVector finalMomentum = primaryMomentum - secondaryMomentum;
	      G4ThreeVector finalDirection = finalMomentum.unit();
	      aParticleChange.ProposeMomentumDirection(finalDirection);
	      aParticleChange.ProposeEnergy(finalKineticEnergy);
	    }
	}
      else
	{
	  // Check whether primary particle is modified
	  if (product.PrimaryParticleIsModified())
	    {
	      G4ThreeVector finalDirection = product.GetModifiedDirection();
	      aParticleChange.ProposeMomentumDirection(finalDirection);
	      G4double finalKineticEnergy = product.GetModifiedEnergy();
	      aParticleChange.ProposeEnergy(finalKineticEnergy);	      
	    }
	}  
	
    }
  
  return G4VDiscreteProcess::PostStepDoIt(track,step );
}