#ifdef USE_CDFEDM2
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Component: FourVectorCalculator.cc
//
// Purpose:   This class inherits from TowerCalculator. It sums the EM and hadronic compartments as if they were four-vectors.
//
// Created:   15/09/1999  Pierre Savard
//
// History:
//            29/08/2000  Matthias Toennesmann - Wrote makePhysicsTower function. Towers are only accepted if their Pt (not their
//                                               Et!) is greater than _etThreshold (as described in CDF Note 5293, Section 4.3.1)!
//            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 -- add iPhi index to _locs array, to 
//                        work with new version of Locations.cc.
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

//namespace calor{

FourVectorCalculator::FourVectorCalculator()
{
  _etThreshold = 0;
  _calcName = "FourVector";
}

FourVectorCalculator::FourVectorCalculator(float etThreshold)
{
  _etThreshold = etThreshold;
  _calcName = "FourVector";
}
PhysicsTower* FourVectorCalculator::makePhysicsTower(EnergyData* eData, Locations* locs) const
{
  return makeCorrPhysicsTower(eData,locs,1.0,1.0);
}

PhysicsTower* FourVectorCalculator::makeCorrPhysicsTower(EnergyData* eData, Locations* locs,
					      float correm,float corrhad) const
{
  float emEnergy = eData->emEnergy();
  float hadEnergy = eData->hadEnergy();

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

  size_t iEta = eData->iEta();
  size_t iPhi = eData->iPhi();

  float emPhi = eData->emPhi();
  float hadPhi = eData->hadPhi();

  // Take sinTheta from iEta (CalData) 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;
  }

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

  float pt = sqrt(px*px + py*py);

  if(pt < _etThreshold) return 0;   // do not waste time making this tower

  float pz;
  if(eData->fourMomentum().rho() == 0.0)
    pz = emEnergy*cos(locs->emTheta(iEta,iPhi)) + hadEnergy*cos(locs->hadTheta(iEta,iPhi));
  else
    pz = eData->fourMomentum().z();
  float totEnergy = emEnergy + hadEnergy;

  HepLorentzVector fourMomentum(px,py,pz,totEnergy);

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

  float totP = fourMomentum.rho();
  float totEt = totEnergy/totP*pt;

  float emEta, hadEta;
  if(eData->fourMomentum().rho() == 0.0){
    emEta = locs->emEta(iEta,iPhi);
    hadEta = locs->hadEta(iEta,iPhi);
  }
  else{
    emEta = eData->fourMomentum().pseudoRapidity();
    hadEta = emEta;
  }
  float totEta = fourMomentum.pseudoRapidity();

  float totPhi = fourMomentum.phi();
  if(totPhi < 0) totPhi += TWOPI;

  return new PhysicsTower(totEta,emEta,hadEta,totPhi,emPhi,hadPhi,emEnergy,hadEnergy,fourMomentum,emEt,hadEt,totEt,iEta,iPhi);
}

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