// $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"

using TMath::Sqrt;
using TMath::Sin;
using TMath::Cos;
using TMath::ATan;
using TMath::Log;
using TMath::Power;
using TMath::Exp;

ClassImp(PamVMCPrimary)

PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
{
  a.fPDG=b.fPDG;
  a.fX0=b.fX0;
  a.fY0=b.fY0;
  a.fZ0=b.fZ0;
  a.fTHETA=b.fTHETA;
  a.fPHI=b.fPHI;
  a.fP0=b.fP0;
  a.fGOOD=b.fGOOD;

  return a;
}


ClassImp(PamVMCPrimaryGenerator)

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) 
  : TObject(),
    fStack(stack),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    frnd(0)
{
// Standard constructor
  fprimColl = new TClonesArray("PamVMCPrimary");
  frnd = new TRandom3(0);

  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV

} 

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
  : TObject(),
    fStack(0),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    fprimColl(0),
    frnd(0)
{    
  frnd = new TRandom3(0);
// Default constructor
  //Default primary proton
  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV
}

PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() 
{
// Destructor
  delete frnd;
  delete fprimColl;
}

// 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  = fprim.fPDG;
  
  Double_t fvx, fvy, fvz;
  fvx=fprim.fX0;
  fvy=fprim.fY0;
  fvz=fprim.fZ0;

  // Position
  
  Double_t tof = 0.;

  // Energy (in GeV)
  Double_t kinEnergy = MomentumToKinE(fprim.fP0);     
  Double_t e  = fmass + kinEnergy;
 
  // Particle momentum
  Double_t  px, py, pz;
   
  px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI);
  py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
  pz = -fprim.fP0*Cos(fprim.fTHETA);
 
  // Polarization
  TVector3 polar;

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

  PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
 
  *pc = fprim;
}


void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
  fprim.fPDG=pdg;
  //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
  fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
  fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
}

void PamVMCPrimaryGenerator::SetMomentum(
                              Double_t px, Double_t py, Double_t pz) 
{
  fprim.fP0= Sqrt(px*px+py*py+pz*pz);
  fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz);
  fprim.fPHI=ATan(py/px);
}
   				     
void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
{
  if(isEnergy) {
    fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
  } else{
    fprim.fP0=frnd->Uniform(PEmin,PEmax);
  }

}

void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
{
  Double_t alpha = 1.+gamma; //integral spectral index
  if(alpha==0.){
    fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
  } else {
    if(PEmin==0.) PEmin=1.E-10;
    fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
  }

  if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);

}