// $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::ACos;
using TMath::Tan;
using TMath::ATan;
using TMath::ATan2;
using TMath::Log;
using TMath::Power;
using TMath::Exp;
using TMath::Pi;
using TMath::Abs;

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((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI)); 
  py = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Sin((Pi()/180.)*(fprim.fPHI));
  pz = -fprim.fP0*Cos((Pi()/180.)*(fprim.fTHETA)); // converting in radian
  //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);

}

void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax)
{
	Double_t theta, phi, y, f;
	phi = (Pi()/180.)*frnd->Uniform(Phmin,Phmax);
	
	do
	{	y = frnd->Uniform(0.,1.);
		theta = (Pi()/180.)*frnd->Uniform(Thmin,Thmax);
		f = Sin(theta);
	} while (y>f);
	
	//theta = phi = (Pi()/180.)*45.;
	SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta)));

	//random distribution of theta phi in the angle at the vertex at (0,0,r)
	//by the S3 max distant corners.
	static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159;
	Double_t ang = ATan((Sqrt(s3_x*s3_x+s3_y*s3_y))/(r-s3_pz));
	
	//SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform((phi-ang),(phi+ang)))/(Pi()/180.));
	SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform(0.,2*Pi()))/(Pi()/180.));
	//SetDirection(0., 0.);
}

void PamVMCPrimaryGenerator::GenSphDist(Double_t r)
{
	//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
	static const Double_t s1_x=20.4, s1_y=16.5, s1_pz=102.8866;
	static const Double_t s2_x=9.0, s2_y=7.5, s2_pz=73.489;
	static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159;

	//calculate max theta and phi angles to be allowed (calculating wrt one 
	//corner of s3 to opposite corner of s1)
	Double_t rx = s1_x+s3_x;
	Double_t ry = s1_y+s3_y;
	Double_t rz = s1_pz-s3_pz;

	Double_t thmax = (180./Pi())*(ACos(rz/Sqrt(rx*rx+ry*ry+rz*rz)));
	Double_t phmax = (180./Pi())*(ATan2(ry,rx));
	//cout << "~~~~~~Theta max Phi max : " << thmax <<", "<< phmax << endl;

	//generate a track check and let it go only if it is passing through all
	//the 3 TOF
	Bool_t trkGood = kFALSE;

	do
	{
	//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
		GenSphDist(r, 0., thmax, 0., phmax);
		
		Double_t th = (Pi()/180.)*fprim.fTHETA;
		Double_t ph = (Pi()/180.)*fprim.fPHI;
	//cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl;

		Double_t x1, y1, x2, y2, x3, y3;
		
		x1 = s1_pz*Tan(th)*Cos(ph);
		y1 = s1_pz*Tan(th)*Sin(ph);
		x2 = s2_pz*Tan(th)*Cos(ph);
		y2 = s2_pz*Tan(th)*Sin(ph);
		x3 = s3_pz*Tan(th)*Cos(ph);
		y3 = s3_pz*Tan(th)*Sin(ph);

		if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) &&
		     Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) &&
		     Abs(x3) <= Abs(s3_x) && Abs(y3) <= Abs(s3_y) ) trkGood = kTRUE;
	//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
	}while (!trkGood);
	//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
}
