// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 
//

#include <TVirtualMC.h>
#include <TVirtualMCStack.h>
#include <TPDGCode.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TVector3.h>
#include <TMath.h>

#include "PamVMCPrimaryGenerator.h"

ClassImp(PamVMCPrimaryGenerator)

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) 
  : TObject(),
    fStack(stack),
    fPdg(kProton),
    fKinEnergy(0.05),//100 GeV
    fDirX(0.),
    fDirY(0.),
    fDirZ(-1.),
    fPolAngle(0.),
    fNofPrimaries(1)
{
// Standard constructor

} 

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
  : TObject(),
    fStack(0),
    fPdg(0),
    fKinEnergy(0.),
    fDirX(0.),
    fDirY(0.),
    fDirZ(0.),
    fPolAngle(0.),
    fNofPrimaries(0)
{    
// Default constructor
}

PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() 
{
// Destructor  
}

// private methods

#include <Riostream.h>

void PamVMCPrimaryGenerator::GeneratePrimary()
{    
// Add one primary particle to the user stack (derived from TVirtualMCStack).
  
  // Track ID (filled by stack)
  Int_t ntr;
 
  // Option: to be tracked
  Int_t toBeDone = 1; 
 
  // Particle type
  Int_t pdg  = fPdg;
  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fPdg);
 
  // Position
  Double_t vx  = 1.; 
  Double_t vy  = 1.; 
  Double_t vz =  130.;
  Double_t tof = 0.;

  // Energy (in GeV)
  Double_t kinEnergy = fKinEnergy; 
  Double_t mass = particlePDG->Mass();      
  Double_t e  = mass + kinEnergy;
 
  // Particle momentum
  Double_t p0, px, py, pz;
  p0 = sqrt(e*e - mass*mass); 
  px = p0 * fDirX;
  py = p0 * fDirY; 
  pz = p0 * fDirZ; 
 
  // Polarization
  TVector3 polar;
//  if ( fPdg == 50000050 ) {
//    TVector3 normal (1., 0., 0.);
//    TVector3 kphoton = TVector3(fDirX, fDirY, fDirZ);
//    TVector3 product = normal.Cross(kphoton); 
//    Double_t modul2  = product*product;
 
//    TVector3 e_perpend (0., 0., 1.);
//    if (modul2 > 0.) e_perpend = (1./sqrt(modul2))*product; 
//    TVector3 e_paralle = e_perpend.Cross(kphoton);
 
///    polar =   TMath::Cos(fPolAngle*TMath::DegToRad())*e_paralle 
  //           + TMath::Sin(fPolAngle*TMath::DegToRad())*e_perpend;
  //  }	
  //else 
  //  Warning("GeneratePrimary",
  //          "The primary particle is not an opticalphoton");

  // Add particle to stack 
  fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, vx, vy, vz, tof, 
                   polar.X(), polar.Y(), polar.Z(), 
                   kPPrimary, ntr, 1., 0);
}

// public methods

void PamVMCPrimaryGenerator::GeneratePrimaries()
{    
// Fill the user stack (derived from TVirtualMCStack) with primary particle

  for (Int_t i=0; i<fNofPrimaries; i++) GeneratePrimary();  
}

void PamVMCPrimaryGenerator::SetDirection(
                              Double_t dirX, Double_t dirY, Double_t dirZ) 
{ 
// Set normalized direction

  Double_t norm = TMath::Sqrt(dirX*dirX + dirY*dirY + dirZ*dirZ);
  
  fDirX = dirX/norm; 
  fDirY = dirY/norm;  
  fDirZ = dirZ/norm;
}