#include "G4SystemOfUnits.hh"
#include "G4ThreeVector.hh"
#include "G4Event.hh"
#include "G4ParticleTable.hh"
#include "Randomize.hh"

#include "PrimaryGeneratorAction.hh"

PrimaryGeneratorAction::PrimaryGeneratorAction()
{
  fParticleGun  = new G4ParticleGun(1);
}

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  delete fParticleGun;
}

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  G4int i=0;
  G4double x0,y0,z0,theta,phi;
  G4double test;
  
  x0=y0=z0=0.0;

  phi = CLHEP::twopi*CLHEP::RandFlat::shoot(0.0,1.0)*rad; 

  while(i==0){
     theta = CLHEP::RandFlat::shoot(0.0,CLHEP::pi)*rad;
     test = CLHEP::RandFlat::shoot(0.0,1.0);
     if(sin(theta)>test)i++;
  }

  //Particle direction
  G4double kx, ky, kz; //components
  kx=cos(phi)*sin(theta);
  ky=sin(phi)*sin(theta);
  kz=cos(theta);
  G4ThreeVector dir_vec (kx,ky,kz);

  //Polarization 
  G4ThreeVector polar = Polarisation(dir_vec);
  
  fParticleGun->SetParticleEnergy(4*GeV);
  fParticleGun->SetParticleMomentumDirection(dir_vec);
  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
  fParticleGun->SetParticlePolarization(polar); 
  
  G4ParticleDefinition* particle=
    G4ParticleTable::GetParticleTable()->FindParticle("mu-");
    //G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
  
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->GeneratePrimaryVertex(anEvent);

}

G4ThreeVector PrimaryGeneratorAction::Polarisation(G4ThreeVector d){

  if(d.mag()!=1.0) d = d.unit();

  G4ThreeVector a = d.orthogonal();
  G4ThreeVector b = a.cross(d); 
  a = a.unit();
  b = b.unit();
  G4double theta = CLHEP::twopi*G4UniformRand();
  return cos(theta)*a + sin(theta)*b;

}
