#ifdef USE_CDFEDM2
//////////////////////////////////////////////////////////////////////////
//
// Component: MetCalculator.cc
// Purpose: This class inherits from TowerCalculator.  It sums the
//          EM and Hadronic compartements according to what was done
//          in Run 1.
//          
// Created: 15/09/99  Pierre Savard
// History: Aug 15, 2002 Jean-Francois Arguin: Change makePhysicsTower
//                       to accommodate PhysicsTower made from CdfTracks
//                       or HEPG bank, i.e.: if EnergyData has a non-null
//                       eta member, use it for PhysicsTower eta instead 
//                       of using ieta.
//          25/11/2002  Matthias Toennesmann - Changes triggered by the changes in 
//                      the EnergyData class: The EnergyData object now contains 
//                      a 4-vector instead of its pseudorapidity.
//                      The PhysicsTower now gets the 4-vector of the EnergyData 
//                      object. It is not made massless anymore (important if 
//                      HEPG is the source).
//          01/05/2003  Dan MacQueen --  changed _locs(ieta) to _locs(ieta,iphi),
//                      added vertex corrections.
//
//////////////////////////////////////////////////////////////////////////

#include <iostream>
#include "Calor/MetCalculator.hh"


// namespace calor {
  
MetCalculator::MetCalculator()
{
  _energyThreshold = 0.;
  _calcName = "Met";
}
MetCalculator::MetCalculator(float energyThreshold)
{
  _energyThreshold = energyThreshold;
  _calcName = "Met";
}

PhysicsTower* MetCalculator::makePhysicsTower(EnergyData* eData, Locations* locs) const
{
  return makeCorrPhysicsTower(eData,locs,1.0,1.0);
}
PhysicsTower* MetCalculator::makeCorrPhysicsTower(EnergyData* eData, Locations* locs,
					      float correm,float corrhad) const
{
  size_t iEta = eData->iEta();
  size_t iPhi = eData->iPhi();
  float emEnergy = eData->emEnergy();
  float hadEnergy = eData->hadEnergy();
  if(emEnergy < _energyThreshold) emEnergy = 0.;  
  if(hadEnergy < _energyThreshold) hadEnergy = 0.;  

  if(emEnergy+hadEnergy == 0.) return 0;  // do not waste time making this tower
  
  // Take sinTheta from iEta & iPhi (locs) or eta (CdfTrack or HEPG)
  float emSinTheta, hadSinTheta;
  if(eData->fourMomentum().rho() == 0.0){
    emSinTheta = locs->emSinTheta(iEta,iPhi);
    hadSinTheta = locs->hadSinTheta(iEta,iPhi);
  }
  else{
    emSinTheta = sin(eData->fourMomentum().theta());
    hadSinTheta = emSinTheta;
  }

  // correction (offline LERs)
  emEnergy =emEnergy * correm;
  hadEnergy =hadEnergy * corrhad;

  float emEt = emEnergy * emSinTheta;
  float hadEt = hadEnergy * hadSinTheta;
  float totEt = emEt + hadEt;

  float emEta, hadEta, emTheta, hadTheta;
  if(eData->fourMomentum().rho() == 0.0){
   emEta = locs->emEta(iEta,iPhi);
   hadEta = locs->hadEta(iEta,iPhi);
   emTheta = locs->emTheta(iEta,iPhi);
   hadTheta = locs->hadTheta(iEta,iPhi);
  }
  else{
    emEta = eData->fourMomentum().pseudoRapidity();
    hadEta = emEta;
    emTheta = eData->fourMomentum().theta();
    hadTheta = emTheta;
  }
  float totEta = (emEta *emEt + hadEta *hadEt)/totEt;
  
  float emPhi = eData->emPhi();
  float hadPhi = eData->hadPhi();

  // Phis from locs are corrected for x,y,z vertex.
  // eData's phis assume x & y vertex of 0.0, but are corrected for the
  // tower centre of gravity.
  // Thus, we use eData's phis in conjunction with corrections from locs
  // to find the corrected phis.

  float emRxy = locs->emRxy(iEta,iPhi);
  float hadRxy = locs->hadRxy(iEta,iPhi);
  float xVertex = locs->threeVertex(0);
  float yVertex = locs->threeVertex(1);
  float emX;
  float emY;
  float hadX;
  float hadY;

  emY = emRxy*sin(emPhi) - yVertex;
  emX = emRxy*cos(emPhi) - xVertex;

  if (emX > 0.0 && emY >= 0.0){
    emPhi = atan(emY/emX);
  } else if (emX < 0.0 && emY > 0.0) {
    emPhi = PI - atan(emY/(-1.0*emX));
  } else if (emX < 0.0 && emY <= 0.0) {
    emPhi = PI + atan (emY/emX);
  } else if (emX > 0.0 && emY < 0.0) {
    emPhi = 2.0*PI - atan((-1.0*emY)/emX);
  } else if (emX = 0.0 && emY >= 0.0) {
    emPhi = PI/2.0;
  } else if (emX = 0.0 && emY < 0.0) {
    emPhi = 3.0*PI/2.0;
  }
    
  hadY = hadRxy*sin(hadPhi) - yVertex;
  hadX = hadRxy*cos(hadPhi) - xVertex;

  if (hadX > 0.0 && hadY >= 0.0){
    hadPhi = atan(hadY/hadX);
  } else if (hadX < 0.0 && hadY > 0.0) {
    hadPhi = PI - atan(hadY/(-1.0*hadX));
  } else if (hadX < 0.0 && hadY <= 0.0) {
    hadPhi = PI + atan (hadY/hadX);
  } else if (hadX > 0.0 && hadY < 0.0) {
    hadPhi = 2.0*PI - atan((-1.0*hadY)/hadX);
  } else if (hadX = 0.0 && hadY >= 0.0) {
    hadPhi = PI/2.0;
  } else if (hadX = 0.0 && hadY < 0.0) {
    hadPhi = 3.0*PI/2.0;
  }

  // float totPhi = (emPhi *emEt + hadPhi *hadEt)/totEt;
  // Calculation of totPhi and totEt will crash if totEt turns out to be zero,
  // but that shouldn't happen.
  // This is where we run into problems: if emPhi or hadPhi were initially
  // close to 2.0*PI, but the above vertex corrections introduced round-off
  // errors which moved one of them to a little above zero, totPhi will
  // end up close to PI.  This is bad.
  // What we do here is make it so if hadPhi is a little under 2*PI,
  // while emPhi is a little over zero, we use emPhi+2*PI in the weighted 
  // average.  Same if situations are reversed.

  float totPhi;
  if (hadPhi > 6.0 && emPhi < 0.28) { 
    totPhi = ((emPhi+2.0*PI) *emEt + hadPhi *hadEt)/totEt;
  } else if (emPhi > 6.0 && hadPhi < 0.28) {
    totPhi = (emPhi *emEt + (hadPhi+2.0*PI) *hadEt)/totEt;
  } else {
    totPhi = (emPhi *emEt + hadPhi *hadEt)/totEt;
  }
  if (totPhi > 2.0*PI) {
    totPhi = totPhi - 2.0*PI;
  }

  float px,py,pz;
  if(eData->fourMomentum().rho() == 0.0){
    px = emEnergy * cos(emPhi) * emSinTheta +
      hadEnergy * cos(hadPhi) * hadSinTheta;
    py = emEnergy * sin(emPhi) * emSinTheta +
      hadEnergy * sin(hadPhi) * hadSinTheta;
    pz = emEnergy * cos(emTheta) +
      hadEnergy * cos(hadTheta) ;
  }
  else{
    px = eData->fourMomentum().x();
    py = eData->fourMomentum().y();
    pz = eData->fourMomentum().z();
  }

  HepLorentzVector fourMomentum(px,py,pz,(emEnergy+hadEnergy));
  
  return new PhysicsTower(totEta, emEta, hadEta, totPhi, emPhi, 
                           hadPhi, emEnergy, hadEnergy, fourMomentum,
                           emEt, hadEt, totEt, iEta, iPhi);

}

// } // namespace calor
#endif  // USE_CDFEDM2