#include "IElectronCloud.hxx"
#include "IFieldManager.hxx"
#include <TRandom.h>
#include <iostream>
#include <cmath>
#include <HEPUnits.hxx>

using namespace std;

//***********************************************************
COMET::IElectronCloud::IElectronCloud(int nE, TVector3 pos, TVector3 vel, double t, int q) {
//***********************************************************

  //Save the start values
  fNbrEs = nE;
  fPos = pos;
  fStartPos = pos;
  fVel = vel;
  fT = t;
  fStartT = t;
  fCharge = q;

  Initialize();

}

//***********************************************************
COMET::IElectronCloud::~IElectronCloud() {
//***********************************************************
}

//***********************************************************
void COMET:: IElectronCloud::Initialize() {
//***********************************************************
  fDrifted = false;
  fMeanDriftXVelocity = 0;

  //Sigma of dot is sqrt of (radius squared divided by 4)  (in mm) 
  fYInitialSigma = 2;
  fZInitialSigma = 2;
  fXInitialSigma = 0;
  fRDisk = 0;
  fCornersAreSet = false;
  //Do not initialize field manager here! That is done in IDriftVolume
}

//***********************************************************
int COMET::IElectronCloud::GetNumberElectrons() {
//***********************************************************
  return fNbrEs;
}
//***********************************************************
int COMET::IElectronCloud::GetCharge() {
//***********************************************************
  return fCharge;
}
//***********************************************************
void COMET::IElectronCloud::SetNumberElectrons(int nE) {
//***********************************************************
  fNbrEs = nE;
}
//***********************************************************
TVector3 COMET::IElectronCloud::GetPosition() {
//***********************************************************
  return fPos;
}
//***********************************************************
void COMET::IElectronCloud::SetPosition(TVector3 pos) {
//***********************************************************
  fPos = pos;
}
//***********************************************************
TVector3 COMET::IElectronCloud::GetVelocity() {
//***********************************************************
  return fVel;
}
//***********************************************************
void COMET::IElectronCloud::SetVelocity(TVector3 vel) {
//***********************************************************
  fVel = vel;
}
//***********************************************************
double COMET::IElectronCloud::GetMeanDriftXVelocity() {
//***********************************************************
  return fMeanDriftXVelocity;
}
//***********************************************************
void COMET::IElectronCloud::SetMeanDriftXVelocity(double vel) {
//***********************************************************
  fMeanDriftXVelocity = vel;
}
//***********************************************************
double COMET::IElectronCloud::GetTime() {
//***********************************************************
  return fT;
}
//***********************************************************
void COMET::IElectronCloud::SetTime(double time){ 
//***********************************************************
  fT = time;
}
//***********************************************************
TVector3 COMET::IElectronCloud::GetStartPosition() {
//***********************************************************
  return fStartPos;
}
//***********************************************************
bool COMET::IElectronCloud::IsDrifted() {
//***********************************************************
  return fDrifted;
}
//***********************************************************
void COMET::IElectronCloud::SetIsDrifted(bool ok) {
//***********************************************************
  fDrifted = ok;;
}
//***********************************************************
int COMET::IElectronCloud::NumberDriftSteps() {
//***********************************************************
  return fNbrSteps;
}
//***********************************************************
void COMET::IElectronCloud::SetNumberDriftSteps(int steps) {
//***********************************************************
  fNbrSteps = steps;
}
//***********************************************************
double COMET::IElectronCloud::GetXDriftLength() {
//***********************************************************
  return fabs((fStartPos - fPos).X()); //returns 0.0 if not drifted
}

//***********************************************************
double COMET::IElectronCloud::GetTotDriftLength() {
//***********************************************************
  return (fStartPos - fPos).Mag();  //returns 0.0 if not drifted
}


//***********************************************************
double COMET::IElectronCloud::GetRandomArrivalTime(double longDiffCoef) {
//***********************************************************
  if (IsDrifted()) {
    //Get the x-displacement due to E-field
    double totalDx_E = GetXDriftLength();
    if (totalDx_E > 0.000001) {
      //Standard deviation 
      double sigma = longDiffCoef*sqrt( totalDx_E );
      //Get dx displacemenent due to longitudinal diffusion
      double dx = gRandom->Gaus(0, sigma + fXInitialSigma);
      //cout<<"IElectronCloud:  dx = "<<dx<<"   sigma = "<<sigma<<"    fXInitialSigma = "<<fXInitialSigma<<"  mead drift velocity = "
      //	  <<fMeanDriftXVelocity<<"    time "<<GetTime()<<endl;
      //Return diffused arrival time around the clouds arrival time
      if (fabs(fMeanDriftXVelocity) > 0.000000001) {
	return GetTime() + dx/fabs(fMeanDriftXVelocity);
      }
    }
  }
  cout<<"IElectronCloud::GetRandomArrivalTime; Returning -1"<<endl;
  return -1.;
}

//***********************************************************
TVector3 COMET::IElectronCloud::GetAnElectronsPosition(double tranDiffCoef, double longDiffCoef, bool calibtarget) {
//***********************************************************
  double dXDrift = GetXDriftLength();
  //Standard deviations 

  
  double tranSigma = tranDiffCoef*sqrt( dXDrift );
  double longSigma = longDiffCoef*sqrt( dXDrift );
 
  

  double ySigma;
  double zSigma;
 
  if(!calibtarget){
     ySigma = tranSigma ;
     zSigma = tranSigma;
  } else{
    //Add the transdiffusion and dot radius in quadrature. 
    ySigma = TMath::Sqrt(tranSigma*tranSigma + fYInitialSigma*fYInitialSigma);
    zSigma = TMath::Sqrt(tranSigma*tranSigma + fZInitialSigma*fZInitialSigma);
  }
  
  double xSigma = longSigma + fXInitialSigma;

    
  //Original positions (if cloud did not origin from a perfect point)
  //TVector2 localPosOnDisk(GetRandomPositionOnDisk());
  //TVector2 localPosOnParallelogram(GetRandomPositionOnParallelogram());
  //TVector2 localPos = localPosOnDisk; // + localPosOnParallelogram; //Will probably never be both, but just in case

  //TEMP
  //TEMP assume  cloud from origin from a perfect point (dot center)
  //localPos = TVector2(0,0);  //Take away this if we dont want to assume that all electrons
                             //in the cloud originated from center of the dot
  //cout<<"tranDiffCoef = "<<tranDiffCoef<<endl;
  //cout<<"dXDrift = "<<dXDrift<<endl;
  //cout<<"tranSigma = "<<tranSigma<<endl;
  //cout<<"ySigma = "<<ySigma<<endl;
  //cout<<"zSigma = "<<zSigma<<endl;
  //END TEMP


  //Return diffused arrival time around the clouds arrival position + orgin local position (if any)
  double x = 0;
  if (fDrifted) {  //If drifted all electrons of the cloud is on the pad plane, then use 
    x = fPos.X();  //GetRandomArrivalTime to get the diffusion in x
  } else {
    x = fPos.X() + gRandom->Gaus(0, xSigma);
  }
  return TVector3 (x, 
		   fPos.Y() + gRandom->Gaus(0, ySigma), 
		   fPos.Z() + gRandom->Gaus(0, zSigma));
}

//***********************************************************
TVector2 COMET::IElectronCloud::GetRandomPositionOnDisk() {
//***********************************************************
  //For now assume homogenous distribution of the photo electrons on the dot
  TVector2 localPos(0,0);
  if (fRDisk > 0.0) {
    localPos.SetMagPhi(fRDisk*gRandom->Uniform(), 2*TMath::Pi()*gRandom->Uniform());
  }
  return localPos;
}

//***********************************************************  
TVector2 COMET::IElectronCloud::GetRandomPositionOnParallelogram() {
//***********************************************************

  if (!fCornersAreSet) return TVector2(0,0);
  //For now this only works if it is a parallelogram, i.e. two pairs of parallel sides, and if the corners
  //are put consecutively in the vector, so check this  (this could easily be changed if needed)
  if (fCorners.size() != 4) cout<<"ERROR in IElectronCloud::GetRandomPositionOnParallelogram, FIXME: allow for more corners than 4"<<endl; 
  double side1 = (fCorners[1] - fCorners[0]).Mag();
  double side2 = (fCorners[2] - fCorners[1]).Mag();
  double side3 = (fCorners[3] - fCorners[2]).Mag();
  double side4 = (fCorners[0] - fCorners[3]).Mag();
  if ((fabs(side1 - side3) > 0.001)||(fabs(side2 - side4) > 0.001)) {
    cout<<"ERROR in IElectronCloud::GetRandomPositionOnFlatShape, FIXME: allow for arbitrary strips"<<endl;
    exit(1);
  } 
  if (TVector3(fCorners[1]-fCorners[0]).Angle(TVector3(fCorners[2]-fCorners[0])) >
      TVector3(fCorners[1]-fCorners[0]).Angle(TVector3(fCorners[3]-fCorners[0]))) {
    cout<<"ERROR in IElectronCloud::GetRandomPositionOnParallelogram, FIXME: allow for non consecutive corners"<<endl;
    exit(1);
  }
  //For now assume homogenous distribution of the photo electrons on the strip-parts
  //Start as if corner 0 was origin 
  //Move along the side between corner 0 and 1
  TVector3 pos = gRandom->Uniform()*(fCorners[1] - fCorners[0]).Mag()*(fCorners[1] - fCorners[0]).Unit();
  //Then move along the side between corner 0 and 3
  pos = pos + gRandom->Uniform()*(fCorners[3] - fCorners[0]).Mag()*(fCorners[3] - fCorners[0]).Unit();
  //Now get the vector from corner 0 to the center in similar way
  TVector3 center = 0.5*(fCorners[1] - fCorners[0]).Mag()*(fCorners[1] - fCorners[0]).Unit();
  center = center + 0.5*(fCorners[3] - fCorners[0]).Mag()*(fCorners[3] - fCorners[0]).Unit();
  //Get the vector pointing from center (that we will now use as origo) to the position
  pos = (pos - center);
  //Return a local postion where center is origo
  return TVector2(pos.Z(), pos.Y());
}

//***********************************************************
void COMET::IElectronCloud::UnDrift() {
//***********************************************************
  fPos = fStartPos;
  fT = fStartT;
  fVel = TVector3(0,0,0);
  fMeanDriftXVelocity = 0;
  fDrifted = false;

}