#include <Randomize.hh>
#include <G4OpProcessSubType.hh>

#include <RAT/OpRayleigh.hh>
#include <RAT/UserTrackInformation.hh>

using CLHEP::twopi;

using namespace RAT;

OpRayleigh::OpRayleigh( const G4String& processName,
                          G4ProcessType type)
  : G4VDiscreteProcess( processName, type )
{
  SetProcessSubType( fOpRayleigh );

  fPhysicsTable = NULL;
  BuildThePhysicsTable();
}

OpRayleigh::~OpRayleigh()
{
  if( fPhysicsTable != NULL )
    {
      fPhysicsTable->clearAndDestroy();
      delete fPhysicsTable;
    }
}

G4VParticleChange*
OpRayleigh::PostStepDoIt( const G4Track& track,
                          const G4Step& step )
{
  aParticleChange.Initialize( track );
  const G4DynamicParticle* aParticle = track.GetDynamicParticle();

  double psi = 0.0;
  do
    {
       psi = twopi * G4UniformRand() / 2.0;
    }
  while( G4UniformRand() > pow( sin( psi ) , 3 ) );
  const double phi = G4UniformRand() * twopi;

  const G4ThreeVector oldMomentum = aParticle->GetMomentumDirection().unit();
  const G4ThreeVector oldPolarisation = aParticle->GetPolarization().unit();
  const G4ThreeVector e2 = oldMomentum.cross( oldPolarisation );

  // oldPolarisation is Z, oldMomentum is x and e2 y
  G4ThreeVector newMomentum = cos( psi ) * oldPolarisation + sin( psi ) * cos( phi ) * oldMomentum + sin( psi ) * sin( phi ) * e2;
  newMomentum.unit();

  G4ThreeVector newPolarisation = oldPolarisation - cos( psi ) * newMomentum;
  newPolarisation.unit();

  aParticleChange.ProposePolarization( newPolarisation );
  aParticleChange.ProposeMomentumDirection( newMomentum );


  RAT::UserTrackInformation* trackHistory = dynamic_cast<RAT::UserTrackInformation*>(track.GetUserInformation());
  if (trackHistory) {
    trackHistory->AddToHistory(RAT::DS::MCPE::hRayleighScatt);
  }

  return G4VDiscreteProcess::PostStepDoIt( track, step );
}

G4double
OpRayleigh::GetMeanFreePath( const G4Track& track,
                             G4double,
                             G4ForceCondition* )
{
  const G4DynamicParticle* particle = track.GetDynamicParticle();
  const G4double photonMomentum = particle->GetTotalMomentum();
  const G4Material* material = track.GetMaterial();
  G4PhysicsOrderedFreeVector* rayleigh = static_cast<G4PhysicsOrderedFreeVector*>( (*fPhysicsTable)( material->GetIndex() ) );

  G4double rsLength = DBL_MAX;
  if( rayleigh != NULL )
    rsLength = rayleigh->Value( photonMomentum );
  return rsLength;
}

void
OpRayleigh::BuildThePhysicsTable()
{
  if( fPhysicsTable )
    return;
  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
  const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
  fPhysicsTable = new G4PhysicsTable( numOfMaterials );

  for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
    {
      G4Material* material = (*theMaterialTable)[iMaterial];
      G4MaterialPropertiesTable* materialProperties = material->GetMaterialPropertiesTable();
      G4PhysicsOrderedFreeVector* rayleigh = NULL;
      if( materialProperties != NULL )
        {
          rayleigh = materialProperties->GetProperty( "RSLENGTH" );
        }
      fPhysicsTable->insertAt( iMaterial, rayleigh );
    }
}