// $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, UInt_t seed) 
  : TObject(),
    fStack(stack),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    frnd(0)
{
// Standard constructor
  fprimColl = new TClonesArray("PamVMCPrimary");
  frnd = new TRandom3(seed);

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

} 

// PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(UInt_t seed)
//   : TObject(),
//     fStack(0),
//     fevno(0),
//     fmass(0.),
//     fcharge(0.),
//     fprimColl(0),
//     frnd(0)
// {    
//   frnd = new TRandom3(seed);
// // 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()
  : 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)); // ritabrata
  //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); // RADIANTS
  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);

}


// Spherical distribution -- Test by Cecilia P march 2009 -----------------------------///
void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, 
					  Double_t zmin, Double_t zmax)
{
  Bool_t trkGood = kFALSE;
  Double_t theta = 999.;
  Double_t phi = 0.;
  Double_t x2,y2,x3,y3;

  //static const Double_t rad2deg = 57.2958;

  //static const Double_t s1_xmax=20.4, s1_ymax=16.5, s1_pz=102.8866;
  // calo cavity 8.07x6.57 with z2=71.6 z3=27.4 circa
  static const Double_t s2_xmax=7.8,  s2_ymax=6.0,  s2_z=73.489;
  static const Double_t s3_xmax=8.0,  s3_ymax=6.0,  s3_z=25.3159;

  // Generate random position unif. distr.
  // GenPosition(-25.,25., 25.,25., 108.0,108.0);
//   Double_t x0= fprim.fX0;
//   Double_t y0= fprim.fY0;
//   Double_t z0= fprim.fZ0;
  Double_t x0=  frnd->Uniform(xmin,xmax);
  Double_t y0=  frnd->Uniform(ymin,ymax);
  Double_t z0=  frnd->Uniform(zmin,zmax);
  Int_t posGood = 0;


  while (trkGood!=kTRUE)
    {
      // Generate phi theta
      theta=999.; // init
      while (theta>=0.82) //take only theta smaller than 30deg=0.52rad
	{
	  Double_t xcos = sqrt( frnd->Uniform(0.,1.) );
	  theta = acos(xcos);   //RAD
	}
      phi = frnd->Uniform(0.,2.*Pi());
 
      // Calculate xy at beginning/end of magnet cavity
      Double_t fact2 = (s2_z-z0)/cos(theta);
      x2 = x0 + fabs(fact2) *  sin(theta) * cos(phi);
      y2 = y0 + fabs(fact2) *  sin(theta) * sin(phi);
      Double_t fact3 = (s3_z-z0)/cos(theta);
      x3 = x0 + fabs(fact3) *  sin(theta) * cos(phi);
      y3 = y0 + fabs(fact3) *  sin(theta) * sin(phi);
	 
      // Test condition on the direction
      if ( Abs(x2) <= Abs(s2_xmax) && Abs(y2) <= Abs(s2_ymax) &&
	   Abs(x3) <= Abs(s3_xmax) && Abs(y3) <= Abs(s3_ymax) ) {
	trkGood = kTRUE;
	//cout<<" x/y0= "<<x0<<" "<<y0<<"  x/y2= "<<fact2*sin(theta)*cos(phi)<<" "<<x2<<"  xy3= "<<
	//  fact3*sin(theta)*cos(phi)<<" "<<x3<<endl;
      }
	 
      //if from current x0y0z0 the condition are not satysfied for any angle => def new start position
      posGood++;
      if ( posGood == 100 )
	{
	  x0=  frnd->Uniform(xmin,xmax);
	  y0=  frnd->Uniform(ymin,ymax);
	  z0=  frnd->Uniform(zmin,zmax);
	  posGood = 0;
	}
      
    } 

  // Set direction and position:
  SetDirection(theta, phi);
  SetPosition(x0, y0, z0);
  
  return;
}
//--- end Test by Cecilia ------------------------------------///




void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax)
{
  // all angles in RADIANTS
	Double_t theta, phi, y, f;
	phi = frnd->Uniform(Phmin,Phmax);
	
	do
	{	y = frnd->Uniform(0.,1.);
		theta = frnd->Uniform(Thmin,Thmax);
		f = Sin(theta);
	} while (y>f);
	
	//SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta))); // ritabrata
	SetPosition(1.,1.,130.); // cecilia


	//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)) ,  frnd->Uniform(0.,2*Pi()) );
	//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., 2*Pi() );
	
		Double_t th = fprim.fTHETA;
		Double_t ph = 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) - fprim.fX0;
		y1 = s1_pz*Tan(th)*Sin(ph) - fprim.fY0;
		x2 = s2_pz*Tan(th)*Cos(ph) - fprim.fX0;
		y2 = s2_pz*Tan(th)*Sin(ph) - fprim.fY0;
		x3 = s3_pz*Tan(th)*Cos(ph) - fprim.fX0;
		y3 = s3_pz*Tan(th)*Sin(ph) - fprim.fY0;

		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;

}