/*  Detector.h

    Detector Class  (Sep 2008)

    To simulate detector effects
*/

class Detector{
   public:
     Detector();
     void ecalFast(Particle&, double =0.0, double =0.0, double =0.0, double =0.0);
};

//------------------------------------------------------------------------------
// The constructor
//
Detector::Detector(){  }


//------------------------------------------------------------------------------
// ecalFast
// Simulates the effects of ECAL by smearing four-momenta of the photon
// according to: following energy and spatial resolutions:
//              dE = R*SQRT(E) + k*E    (energy  resolution)
//       dThetaPhi = A/sqrt(E) + B      (angular resolution)
// where E is energy(in GeV).
//
// For ALEPH: R=0.18 and k=0.009; A = 2.5 mrad.GeV^0.5 and B = 0.25 mrad
// For ATLAS: R=0.10 and k=0.007
// For CMS  : R=0.04 and k=0.005
//------------------------------------------------------------------------------
void Detector::ecalFast(Particle& par, double R, double k, double A, double B){

   // Get the four-momenta of the particle
   double pX   = par.px;
   double pY   = par.py;
   double pZ   = par.pz;
   double E    = par.E;
   double pAbs = par.pAbs;
   double pt   = par.pt;

   // Determine the direction of the particle
   double cosPhi   = pX/pt;
   double sinPhi   = pY/pt;
   double cosTheta = pZ/pAbs;
   double sinTheta = sqrt(1.0-cosTheta*cosTheta);

   // Energy resolution in GeV
   double dE = R*sqrt(E) + k*E;

   // Modify the energy of the particle by applying a gaussian smearing
   // and evaluate the components of the space-momentum. (Do not modify the direction)
   E    = Random::Gauss(E, dE);
   pAbs = sqrt(E*E-par.m*par.m);
   pX   = pAbs*sinTheta*cosPhi;
   pY   = pAbs*sinTheta*sinPhi;
   pZ   = pAbs*cosTheta;

   // Now modify the direction of the momentum vector to simulate angular resolution
   if(fabs(A)>0.0 && fabs(B)>0.0){

      // Spatial (angular) resolution in radians
      double dThetaPhi = (A/sqrt(E) + B)*1.0e-3;

      // transform momentum vector to a new coordinate system 
      // whose y axis is directed in to vector p
      double theta = acos(pY/pAbs);
      double phi   = atan(pZ/pX);
      if(pZ<0.0 && pX<0.0) phi += M_PI;
      if(pZ>0.0 && pX<0.0) phi += M_PI;
      double Spx   = cos(theta)*cos(phi)*pX - sin(theta)*pY + cos(theta)*sin(phi)*pZ;
      double Spy   = sin(theta)*cos(phi)*pX + cos(theta)*pY + sin(theta)*sin(phi)*pZ;
      double Spz   =-sin(phi)*pX + cos(phi)*pZ;
      double Sp    = sqrt(Spx*Spx + Spy*Spy + Spz*Spz);
 
      // smear (modify) the momentum vector in the new coordinate
      double alpha = 2.0*M_PI*Random::Flat();
      double beta  = Random::Gauss(0.0, dThetaPhi); // RMS value of smearing angle
      Spx = Sp*sin(beta)*cos(alpha);
      Spy = Sp*cos(beta);
      Spz = Sp*sin(beta)*sin(alpha);
 
      // inverse transform (back to original coordinate)
      pX = sin(theta)*cos(phi)*Spy + cos(theta)*cos(phi)*Spx - sin(phi)*Spz;
      pY = cos(theta)*Spy - sin(theta)*Spx;
      pZ = sin(theta)*sin(phi)*Spy + cos(theta)*sin(phi)*Spx + cos(phi)*Spz;
   }

   // Reset the four-momentum
   par.p4( pX, pY, pZ, E );
}