#include "PrimaryGeneratorMessenger.hh"
#include "PrimaryGeneratorAction.hh"
#include "DetectorConstruction.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "globals.hh"
#include "Randomize.hh"

PrimaryGeneratorAction::PrimaryGeneratorAction(
                            DetectorConstruction* DC)
  :Detector(DC),rndmFlag("off")			    
{
  G4int n_particle = 1;
  particleGun = new G4ParticleGun(n_particle);
  //create a messenger for this class
  gunMessenger = new PrimaryGeneratorMessenger(this);
// default particle

  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4String particleName;
  G4ParticleDefinition* particle 
    = particleTable->FindParticle(particleName="proton");
  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
  particleGun->SetParticleEnergy(20.*GeV);
  G4double position = -(Detector->GetWorld_Z())/2.;
  particleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, position));
}

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

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  if (rndmFlag == "on"){
   
    G4double z = Detector->GetWorld_XY()*(G4UniformRand()-0.5);
    G4double y = Detector->GetWorld_XY()*(G4UniformRand()-0.5);
    G4double x = Detector->GetWorld_Z()/2.;  
    particleGun->SetParticlePosition(G4ThreeVector(x,y,z));

  //k  G4int i=0;
  //k G4double theta=pi/10*(i+G4UniformRand());
    G4double phi=G4UniformRand()*twopi;
   



    G4double theta;
    do 
      theta = G4UniformRand() * halfpi;   
    while( G4UniformRand() > 2.*cos(theta)*sin(theta));
   
    G4double px=sin(theta)*sin(phi);
     
    G4double py=sin(theta)*cos(phi);
 
    G4double pz=cos(theta);

    particleGun->SetParticleMomentumDirection(G4ThreeVector(-fabs(pz),py,px));

 //   G4cout<<"y0 "<<y0<< " z0 "<<z0<<" Theta " << theta<<" Phi "<<phi<<G4endl;

   //    G4double E=pow(6.0255959E-6 -5.91594E-6*G4UniformRand(),-0.5747712644);
   //   G4double E=pow(0.000316228 -0.000310604*G4UniformRand(),-0.571429);
        
   //    particleGun->SetParticleEnergy(E*GeV);
   //    G4cout<<E<<"GeV"<<G4endl;
     
    //l  particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
  }
  particleGun->GeneratePrimaryVertex(anEvent);
}


// 2005 by G.I.Vasilyev