// -------------------------------------------------------------------------
// -----            FairBoxGeneratorMOD source file                        -----
// -----          Created 09/09/04  by Yu.Kharlov
// -------------------------------------------------------------------------

/* $Id: FairBoxGeneratorMOD.cxx,v 1.4 2006/07/18 09:28:06 prokudin Exp $ */

/* History of cvs commits:
 *
 * $Log: FairBoxGeneratorMOD.cxx,v $
 * Revision 1.4  2006/07/18 09:28:06  prokudin
 * Should be * instead /
 *
 * Revision 1.3  2006/07/14 11:23:57  kharlov
 * Add protection for simultaneously set ranges; split vertex and kinematics ranges
 *
 * Revision 1.2  2006/03/29 16:25:50  kharlov
 * New functionality added
 *
 */

#include "FairBoxGeneratorMOD.h"

#include "FairPrimaryGenerator.h"

#include "TRandom.h"
#include "TParticlePDG.h"
#include "TDatabasePDG.h"
#include "TMath.h"

// ------------------------------------------------------------------------
FairBoxGeneratorMOD::FairBoxGeneratorMOD() :
  FairGenerator(),
  fPDGType(0),fMult(0),fPDGMass(0),fPtMin(0),fPtMax(0),
  fPhiMin(0),fPhiMax(0),fEtaMin(0),fEtaMax(0),fYMin(0),fYMax(0),
  fPMin(0),fPMax(0),fThetaMin(0),fThetaMax(0),fX(0),fY(0),fZ(0),
  fX1(0),fX2(0),fY1(0),fY2(0),fZ1(0),fZ2(0),fCX(0),fCY(0),
  fEtaRangeIsSet(0),fYRangeIsSet(0),fPhiRangeIsSet(0),fThetaRangeIsSet(0),
  fCosThetaIsSet(0),fPtRangeIsSet(0),fPRangeIsSet(0),fCosmicsIsSet(0),fSphereIsSet(0),
  fPointVtxIsSet(0),fBoxVtxIsSet(0),fCylVtxIsSet(0),fDebug(0),fPlane1VtxIsSet(0),fDonutIsSet(0)
{
  // Default constructor
}

// ------------------------------------------------------------------------
FairBoxGeneratorMOD::FairBoxGeneratorMOD(Int_t pdgid, Int_t mult) :
  fPDGType(pdgid),fMult(mult),fPDGMass(0),fPtMin(0),fPtMax(0),
  fPhiMin(0),fPhiMax(0),fEtaMin(0),fEtaMax(0),fYMin(0),fYMax(0),
  fPMin(0),fPMax(0),fThetaMin(0),fThetaMax(0),fX(0),fY(0),fZ(0),
  fX1(0),fX2(0),fY1(0),fY2(0),fZ1(0),fZ2(0),
  fEtaRangeIsSet(0), fYRangeIsSet(0),fPhiRangeIsSet(0),fThetaRangeIsSet(0),
  fCosThetaIsSet(0), fPtRangeIsSet(0),fPRangeIsSet(0),fCosmicsIsSet(0),fSphereIsSet(0),
  fPointVtxIsSet(0),fBoxVtxIsSet(0),fCylVtxIsSet(0),fDebug(0),fPlane1VtxIsSet(0),fDonutIsSet(0)
{
  // Constructor. Set default kinematics limits
  SetPhiRange  ();
}
// ------------------------------------------------------------------------
Bool_t  FairBoxGeneratorMOD::Init()
{
  // Initialize generator


  if (fPhiMax-fPhiMin>360)
    Fatal("Init()","FairBoxGeneratorMOD: phi range is too wide: %f<phi<%f",
          fPhiMin,fPhiMax);
  if (fPRangeIsSet && fPtRangeIsSet) {
    Fatal("Init()","FairBoxGeneratorMOD: Cannot set P and Pt ranges simultaneously");
  }
  if (fPRangeIsSet && fYRangeIsSet) {
    Fatal("Init()","FairBoxGeneratorMOD: Cannot set P and Y ranges simultaneously");
  }
  if ( (fThetaRangeIsSet && fYRangeIsSet) ||
       (fThetaRangeIsSet && fEtaRangeIsSet) ||
       (fYRangeIsSet     && fEtaRangeIsSet) ) {
    Fatal("Init()","FairBoxGeneratorMOD: Cannot set Y, Theta or Eta ranges simultaneously");
  }
  if (fPointVtxIsSet && fBoxVtxIsSet) {
    Fatal("Init()","FairBoxGeneratorMOD: Cannot set point and box vertices simultaneously");
  }

  // Check for particle type
  TDatabasePDG* pdgBase = TDatabasePDG::Instance();
  TParticlePDG* particle = pdgBase->GetParticle(fPDGType);

  if (! particle) { Fatal("FairBoxGeneratorMOD","PDG code %d not defined.",fPDGType); }

  if (fCylVtxIsSet)
    fFunc=new TF1("fFunc","sin(x)^2",fCPhi1*TMath::DegToRad(),fCPhi2*TMath::DegToRad());
  if (fCosmicsIsSet)
    fFunc=new TF1("fFunc","sin(x)*cos(x)^2",fCPhi1*TMath::DegToRad(),fCPhi2*TMath::DegToRad());
  else
    fFunc=new TF1("fFunc","x",fCPhi1*TMath::DegToRad(),fCPhi2*TMath::DegToRad());

  fPDGMass = particle->Mass();
  // printf("particle->Mass() = %f \n", fPDGMass);
  return kTRUE;

}

// ------------------------------------------------------------------------
Bool_t FairBoxGeneratorMOD::ReadEvent(FairPrimaryGenerator* primGen)
{
  // Generate one event: produce primary particles emitted from one vertex.
  // Primary particles are distributed uniformly along
  // those kinematics variables which were limitted by setters.
  // if SetCosTheta() function is used, the distribution will be uniform in
  // cos(theta)

  Double32_t pabs=0, phi=0, pt=0, theta=0, eta, y, mt, px=0, py=0, pz=0, vtxphi=0, vtxtheta=0;

  // Generate particles
  for (Int_t k = 0; k < fMult; k++) 
    {
  
      phi = gRandom->Uniform(fPhiMin,fPhiMax) * TMath::DegToRad();

      if      (fPRangeIsSet ) { pabs = gRandom->Uniform(fPMin,fPMax); }
      else if (fPtRangeIsSet) { pt   = gRandom->Uniform(fPtMin,fPtMax); }

      if      (fThetaRangeIsSet) 
	{
	  if (fCosThetaIsSet)
	    theta = acos(gRandom->Uniform(cos(fThetaMin* TMath::DegToRad()),
					  cos(fThetaMax* TMath::DegToRad())));
	  else 
	    {
	      theta = gRandom->Uniform(fThetaMin,fThetaMax) * TMath::DegToRad();
	    }
	} 
      else if (fEtaRangeIsSet) 
	{
	  eta   = gRandom->Uniform(fEtaMin,fEtaMax);
	  theta = 2*TMath::ATan(TMath::Exp(-eta));
	} 
      else if (fYRangeIsSet) 
	{
	  y     = gRandom->Uniform(fYMin,fYMax);
	  mt = TMath::Sqrt(fPDGMass*fPDGMass + pt*pt);
	  pz = mt * TMath::SinH(y);
	}

      if (fThetaRangeIsSet || fEtaRangeIsSet) 
	{
	  if(fPRangeIsSet ) 
	    {
	      pz = pabs*TMath::Cos(theta);
	      pt = pabs*TMath::Sin(theta);
	    } 
	  else if (fPtRangeIsSet) 
	    {
	      pz = pt/TMath::Tan(theta);
	    }
	}

      if(fCylVtxIsSet)
	{
	  if (fDebug)
	    std::cout<<"using cylindrical generator\n";
	  vtxphi=gRandom->Uniform(fCPhi1,fCPhi2)*TMath::DegToRad();
	  fX=fCR*TMath::Cos(vtxphi);
	  fY=fCR*TMath::Sin(vtxphi);
	  fZ=gRandom->Uniform(fCZ,fCZ+fCLen);
	  phi=180*TMath::DegToRad()+vtxphi;
	  fX+=gRandom->Uniform(fX1,fX2);
	  fY+=gRandom->Uniform(fY1,fY2);
	}


      px = pt*TMath::Cos(phi);
      py = pt*TMath::Sin(phi);

      if (fBoxVtxIsSet) 
	{
	  if (fDebug)
	    std::cout<<"Using box generator\n";
	  fX = gRandom->Uniform(fX1,fX2);
	  fY = gRandom->Uniform(fY1,fY2);
	  fZ = gRandom->Uniform(fZ1,fZ2);
	}

      if (fDonutIsSet)
	{
	  double r=gRandom->Uniform(fX1,fX2);
	  double phi=gRandom->Uniform(fY1*TMath::DegToRad(),fY2*TMath::DegToRad());
	  fX = r*TMath::Cos(phi);
	  fY = r*TMath::Sin(phi);
	  fZ = gRandom->Uniform(fZ1,fZ2);
	}

      if (fPlane1VtxIsSet)
	{
	  if (fDebug)
	    std::cout<<"using plane generator\n";
	  fX=gRandom->Uniform(fX1,fX2);
	  if ((fX1-fX2)<0.000001)
	    fY=gRandom->Uniform(fY1,fY2);
	  else
	    fY=( ( (fY1-fY2)/(fX2-fX1) )*(fX2-fX) )+fY2;
	  std::cout<<"fX="<<fX<<" fY="<<fY<<"\n";
	  std::cout<<"phi:"<<TMath::ATan( ( (fY1-fY2)/(fX2-fX1) ) )*TMath::RadToDeg()<<"\n";
	  fZ=gRandom->Uniform(fZ1,fZ2);
	}

      if(fCosmicsIsSet)
	{
	  if (fSphereIsSet)
	    {
	      if (fThetaRangeIsSet)
		vtxtheta=acos(gRandom->Uniform(cos(fThetaMin*TMath::DegToRad()),cos(fThetaMax*TMath::DegToRad())));
	      else
		vtxtheta=gRandom->Uniform(0,180)*TMath::DegToRad();
	      if (fPhiRangeIsSet)
		vtxphi=gRandom->Uniform(fPhiMin,fPhiMax) * TMath::DegToRad();
	      else
		vtxphi=gRandom->Uniform(0,360) * TMath::DegToRad();
	    }
	  else
	    {
	      vtxtheta=fFunc->GetRandom();
	      vtxphi=gRandom->Uniform(0,360)*TMath::DegToRad();
	    }
	  //vtxtheta=90*TMath::DegToRad();
	  //vtxphi=0;
	  fX=fCX+fCR*TMath::Sin(vtxtheta)*TMath::Cos(vtxphi)+gRandom->Uniform(fX1,fX2);;
	  fZ=fCZ+fCR*TMath::Sin(vtxtheta)*TMath::Sin(vtxphi)+gRandom->Uniform(fZ1,fZ2);
	  fY=fCY+fCR*TMath::Cos(vtxtheta)+gRandom->Uniform(fY1,fY2);
	  //theta=vtxtheta;
	  //phi=180*TMath::DegToRad()+vtxphi;
	  px=-pabs*TMath::Sin(vtxtheta)*TMath::Cos(vtxphi);
	  pz=-pabs*TMath::Sin(vtxtheta)*TMath::Sin(vtxphi);
	  py=-pabs*TMath::Cos(vtxtheta);
	  theta=TMath::ACos(pz/pabs);
	  pt=pz*TMath::Tan(theta);
	  phi=TMath::ASin(px/pt);
	}


    
      if (fPDGvarAnti)
	{
	  if(gRandom->Uniform(-1,1)<0)
	    fPDGType = -1*fPDGType;
	}

      if (fDebug)
	printf("BoxGen: kf=%d, p=(%.2f, %.2f, %.2f) GeV, x=(%.1f, %.1f, %.1f) cm, phi=%.1f, theta=%.1f \n",
	       fPDGType, px, py, pz, fX, fY, fZ,phi*TMath::RadToDeg(),theta*TMath::RadToDeg());

      primGen->AddTrack(fPDGType, px, py, pz, fX, fY, fZ);
    }
  return kTRUE;

}
// ------------------------------------------------------------------------


ClassImp(FairBoxGeneratorMOD)