////////////////////////////////////////////////////////////////////
//
//     File: PhotonVariables.cc
//     Authors: Mike Kirby, Ray Culbertson
//     Description: High level analysis photon variables defined as functions
//     Created: 10/6/2000
//
////////////////////////////////////////////////////////////////////

//#include <cstdio>
//#include <math.h>
#include "Rtypes.h"
#include "HighLevelObjects/PhotonVariables.hh"
#include "HighLevelObjects/TrackingVariables.hh"
#include "CalorGeometry/CalParameters.hh"
#include "ElectronObjects/CdfEmObject.hh"
#include "ElectronObjects/EmCluster.hh"
#include "Electron/CemCorrAlg.hh"
#include "Electron/PemCorrAlg.hh"
#include "CalorObjects/CesCluster.hh"
#include "CalorObjects/Pes2dCluster.hh"
#include "Electron/emobj_alg.hh"
#include "CLHEP/Geometry/Point3D.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Tracks/track_cut.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Tracks/SimpleReconstructedTrack.hh"
#include "Trajectory/Helix.hh"
#include "TrackingSI/Utils/SiExpected.hh"

namespace PhotonVariables {


    // set lastVersion to the last version you know of to print changes
  int phoVariablesVersion(int lastVersion) {
    int currentVersion = 14;
    if(lastVersion<currentVersion){
      std::cout << "PhotonVariables verison history:" << std::endl;
    }
    if(lastVersion<1){
      std::cout << "V1: Added PES status variables, track Iso cleaned up." << std::endl;
      std::cout << "V1: Altered ShrMax position, phi, etc to include PES " << std::endl;
    }
    if(lastVersion<2){
      std::cout << "V2: fix for possible offset to x in phoCprCesSeeded" << std::endl;
    }
    if(lastVersion<3){
      std::cout << "V3: add tower seed phi and eta" << std::endl;
    }
    if(lastVersion<4){
      std::cout << "V4: add hadron timing" << std::endl;
    }
    if(lastVersion<5){
      std::cout << "V5: add correction to CES chiSq - significant change" << std::endl;
    }
    if(lastVersion<6){
      std::cout << "V6: add more CES variables" << std::endl;
    }
    if(lastVersion<7){
      std::cout << "V7: fix PES energy, add PES, PPR isolation variables" << std::endl;
    }
    if(lastVersion<8){
      std::cout << "V8: add Delta Z cut to tracking isolation" << std::endl;
    }
    if(lastVersion<9){
      std::cout << "V9: add photon-based Lshr and PEM 3x3 ChiSq" << std::endl;
    }
    if(lastVersion<10){
      std::cout << "V10: remove ces-seeded cpr from stat word" << std::endl;
    }
    if(lastVersion<11){
      std::cout << "V11: fix trigger tower Et for z vertex" << std::endl;
      std::cout << "V11: add 2-tower had/em" << std::endl;
    }
    if(lastVersion<12){
      std::cout << "V12: fix ces and cpr weights" << std::endl;
    }
    if(lastVersion<13){
      std::cout << "V13: switched to two-tower et, redefine stat word" << std::endl;
    }
    if(lastVersion<14){
      std::cout << "V14: added plug E correction" << std::endl;
    }
    if(lastVersion<15){
      std::cout << "V15: tweak IsoLcor, ad PES variables" << std::endl;
    }
    return currentVersion;
  }

    // bits set for each cut failed, by Run I convention
  int phoCutWord(const CdfEmObject& emObject, const double zvertex) {
    int word=0;
    if(phoDetector(emObject)!=CEM) word |= (1<<0);
    if(phoCorrectedEt(emObject,zvertex)<5.0) word |= (1<<1);
    double chiSq = phoCesChiSq(emObject);
    if(chiSq> 4.0) word |= (1<<2);
    if(chiSq>10.0) word |= (1<<3);
    if(chiSq>20.0) word |= (1<<4);
    if(phoLshr(emObject,zvertex)>0.2) word |= (1<<5);
    if(phoHadEm(emObject)>0.125) word |= (1<<6);
    double hademslide = 0.055+.00045*phoCorrectedEnergy(emObject);
    if(phoHadEm(emObject)>hademslide) word |= (1<<7);
    int n3d = phoN3D(emObject);
    if(n3d>0) word |= (1<<8);
    if(n3d>1) word |= (1<<9);
    double etcor = phoCorrectedEt(emObject,zvertex);
    double etiso = phoCo4PJW(emObject,zvertex);
    double eratio = (etcor>0.0? etiso/etcor : 0.0);
    if(eratio>0.10 && etcor>2.0) word |= (1<<10);
    if(etiso>2.0) word |= (1<<11);
    if(etcor<20.0) {
      if(eratio>0.1) word |= (1<<12);
    } else {
      if(etiso>2.0+0.005*etcor) word |= (1<<12);
    }
    double tiso = phoTrackIso(emObject,0.4,zvertex);
    if(tiso>2.0) word |= (1<<13);
    if(tiso>5.0+0.005*etcor) word |= (1<<14);
    double cesslide = phoCesSlide(emObject);
    double sinTh = phoSinTheta(emObject,zvertex);
    bool qces = false;
    //if( phoCesStripE2(emObject)<cesslide && 
    //	(phoCesWireE2(emObject)<cesslide || 
    // fabs( phoCesX(emObject) - phoCesX2(emObject))<7.0) )
    //qces = true;
    if( phoCesStripE2(emObject)*sinTh<cesslide && 
	phoCesWireE2(emObject)*sinTh<cesslide )
	qces = true;
    if(!qces) word |= (1<<15);
    double az = fabs(phoCesZ(emObject));
    if(az<9.0 || az>230.0) word |= (1<<16);
    if(phoPt(emObject)>1.0 + 0.005*etcor) word |= (1<<17);
    if(fabs(phoCesX(emObject))>21.0) word |= (1<<18);
    //double CprX,CprQ;
    // this has to be temorarily moved to StntupleMakerModule
    //int stat = phoCprCesSeeded(emObject, CprX, CprQ, 5, zvertex);
    //if(stat==0 && CprQ>500.0) word |= (1<<19);
    if(!phoHasShowerMax(emObject)) word |= (1<<20);
    qces = false;
    if( phoCesStripE2(emObject)*sinTh<1.0 && 
	phoCesWireE2(emObject)*sinTh<1.0 )
	qces = true;
    if(!qces) word |= (1<<21);

    //skip a few for central

    if(phoPesE(emObject)<=0.0) word |= (1<<25);

    return word;
  }

    // status word = 1000, 2000... + photon Et
  double phoStatus(const CdfEmObject& emObject, const double zvertex) {
    int loose;
    int looseiso;
    int tight;
    int tightiso;
    if(phoDetector(emObject)==CEM) {
      loose = (1<<20) + (1<<17) + (1<<9) + (1<<4) + (1<<7);
      tight = (1<<20) + (1<<18) + (1<<16) + (1<<17) + (1<<9) + (1<<4) + (1<<7);
      looseiso = (1<<10);
      tightiso = (1<<11) + (1<<14);
    } else {
      loose = (1<<20) + (1<<7)+ (1<<25);
      tight = (1<<20) + (1<<7)+ (1<<25);
      looseiso = (1<<10);
      tightiso = (1<<11);
    }

    int word = phoCutWord(emObject, zvertex);
    double status = 1000;  // 1000 for existing
    if((word & loose)==0) {
      status = 2000;
      if((word & looseiso)==0) {
	status = 3000;
      }
      if((word & tight)==0) {
	status = 4000;
	if((word & tightiso)==0) {
	  status = 5000;
	}
      }
    }
    status += phoCorrectedEt(emObject,zvertex);
    return status;
  }


    // Detector region (CEM=0, PEM=1 currently)"CalorGeometry/CalConstants.hh"
  Detector phoDetector(const CdfEmObject& emObject) {
    return emObject.getEmCluster()->region();
  }

    //the Et directly from z=0 and EmCluster centroid (for trigger)
  double phoTowerEt(const CdfEmObject& emObject) {
    const EmCluster* clus = emObject.getEmCluster();
    clus->setVertexZ(0.0);
    return clus->emEt();
  }

    //the Et directly from the EmCluster method
  double phoRawEt(const CdfEmObject& emObject, double zvertex ) {
    return phoRawEnergy(emObject)*phoSinTheta(emObject,zvertex);
  }

    //the Et corected for final offline calibrations
  double phoCorrectedEt(const CdfEmObject& emObject, double zvertex ) {
    return phoRawEt(emObject, zvertex )/phoEnergyCorrection(emObject);
  }

    //four vector using corrected energy
  HepLorentzVector phoFourVector(const CdfEmObject& emObject, 
			       double zvertex ){
    HepPoint3D p;
    if(phoHasShowerMax(emObject)) {   // then there was a CES cluster
      p = phoShowerMaxPosition(emObject);
    } else {
      p = phoCentroidPosition(emObject);
    }
    p = p - HepPoint3D(0.,0.,zvertex);
    double ce = phoCorrectedEnergy(emObject);
    p.setMag(ce);
    return HepLorentzVector(p.x(),p.y(),p.z(),ce);
  }

    // raw Energy
  double phoRawEnergy(const CdfEmObject& emObject) {
    double e = emobj_alg::CEMFourMomentum::TwoTower(emObject,false).e();
    if(e>0.001) return e;
    return (emObject.getEmCluster()->emEnergy());
  }

    // corrected Energy
  double phoCorrectedEnergy(const CdfEmObject& emObject) {
    return phoRawEnergy(emObject)/phoEnergyCorrection(emObject);
  }

    // the phi index of the seed tower
  int phoSeedPhi(const CdfEmObject& emObject) {
    return emObject.getEmCluster()->seedPhiAddr();
  }

    // the eta index of the seed tower
  int phoSeedEta(const CdfEmObject& emObject) {
    return emObject.getEmCluster()->seedEtaAddr();
  }

    // the phi using only towers
  double phoTowerPhi(const CdfEmObject& emObject) {
    return emObject.getEmCluster()->emPhi();
  }

    // the eta using only towers
  double phoTowerEta(const CdfEmObject& emObject) {
    const EmCluster* clus = emObject.getEmCluster();
    clus->setVertexZ(0.0);
    return clus->emEtaDetector();
  }

    //  the tower usage pattern
  int phoTowerPattern(const CdfEmObject& emObject) {
    int seedPhi = phoSeedPhi(emObject);
    int seedEta = phoSeedEta(emObject);
    bool qCen = phoDetector(emObject)==CEM;
    const EmCluster& clust = *(emObject.getEmCluster());
    const EmCluster::TowerIdList&  tlist = clust.towerIds();
    int pat = 0;
    for(EmCluster::TowerIdList::const_iterator it=tlist.begin(); 
	                                      it<tlist.end(); ++it) { 
      TowerKey tk = clust.towerKey(*it);
      if(qCen) {
	// this assumes one phi, up to 3 towers in eta
	pat |= ( 1<<( tk.iEta() - seedEta + 1) );
      } else {
	// the plug boundaries can't be summarized like this
	pat = 0;
      }
    } // end
    return pat;
  }

    // phi
  double phoPhi(const CdfEmObject& emObject) {
    HepPoint3D p;
    if(phoHasShowerMax(emObject)) {   // then there was a CES cluster
      p = phoShowerMaxPosition(emObject);
    } else {
      p = phoCentroidPosition(emObject);
    }
    double phi = p.phi();
    if(phi<0.0) phi += 2*M_PI;
    return phi;
  }

    // event eta
  double phoEta(const CdfEmObject& emObject, double zvertex) {

    HepPoint3D p;
    if(phoHasShowerMax(emObject)) {   // then there was a CES cluster
      p = phoShowerMaxPosition(emObject);
    } else {
      p = phoCentroidPosition(emObject);
    }
    p = p - HepPoint3D(0.,0.,zvertex);
    return p.pseudoRapidity();
  }

    // detector eta
  double phoDteta(const CdfEmObject& emObject) {
    return phoEta(emObject, 0.0);
  }

    // had/em ("3-tower")
  double phoHadEm(const CdfEmObject& emObject) {
    return (emObject.getEmCluster()->hadEm());
  }

  // had/em from trigger tower ("2-tower")
  double phoHadEmT(const CdfEmObject& emObject) {
    return emObject.getEmCluster()->hadEmTriggerTower();
  }

    // lshr this is the electron quantity, needs to be updated
  double phoLshr(const CdfEmObject& emObject, double zvertex, int nTowers) {
    if(phoDetector(emObject)!=CEM) return 99.0;
    double cesZ = phoCesZ(emObject);
    if(fabs(cesZ) > 230.0) return 99.0;
    double lshr = emobj_alg::CalculateLshr::lshr( *emObject.getEmCluster(), 
						  cesZ, zvertex, nTowers);
    return lshr;
  }


  double phoCo4(const CdfEmObject& emObject, const double zvertex) {
    const EmCluster* emc = emObject.getEmCluster();
    double eleSinTheta = emc->emEt()/emc->emEnergy();
    double eleIso = emc->totalIsolationEt4()/eleSinTheta;
    return eleIso*phoSinTheta(emObject,zvertex);
  }

  double phoCo7(const CdfEmObject& emObject, const double zvertex) {
    const EmCluster* emc = emObject.getEmCluster();
    double eleSinTheta = emc->emEt()/emc->emEnergy();
    double eleIso = emc->totalIsolationEt7()/eleSinTheta;
    return eleIso*phoSinTheta(emObject,zvertex);
  }

    // pt of highest pt track (order=1 is highest, 2 is next highest)
  double phoPt(const CdfEmObject& emObject, int order){
    double pt_max1 = -99.0;
    double pt_max2 = -99.0;
    //    const CdfTrack* best_track=0;
    //    const CdfTrack* best_track2=0;
    const CdfTrackView& ctv = emObject.matchingTracks();
    //vector of  CdfTrack_lnk 
    const CdfTrackView::CollType & vlc = ctv.contents();

    for (CdfTrackView::const_iterator it  = vlc.begin();
	                              it != vlc.end(); ++it) {
      double track_pt = (*it)->pt();

      if (track_pt >= pt_max1) {
	//	best_track2 = best_track;
	pt_max2     = pt_max1;
	//	best_track  = it->operator->();
	pt_max1     = track_pt;
      } else if (track_pt >= pt_max2) {
	//	best_track2 = it->operator->();
	pt_max2     = track_pt;
      }
	
    } // endif this track pt> pt_max2
    return (order == 1 ? pt_max1 : pt_max2);
  }


    // number of tracks pointing at cluster
    // this needs to be reviewed, not sure what the quantity is
  int phoN3D(const CdfEmObject& emObject){
    return emObject.numberTracks();
  }


    // Track isolation in cone, needs to be reviewed
  double phoTrackIso(const CdfEmObject& emObject, double cone, 
		                                   double zvertex){
    double phi = phoPhi(emObject);
    double eta = phoEta(emObject,zvertex);
    return   TrackingVariables::trkTrackIsolation(phi, eta, zvertex, cone);

  }

    // Track isolation in cone, needs to be reviewed
  double phoTrackIso(const CdfEmObject& emObject, 
		     const CdfTrackView& trackView, double cone, 
		     double zvertex){
    double phi = phoPhi(emObject);
    double eta = phoEta(emObject,zvertex);
    return TrackingVariables::trkTrackIsolation(phi, eta, zvertex, trackView, cone);

  }

    // Track isolation in a cone, same as above, but also return
    // number of tracks in the cone and Pt of max Pt track
  double phoTrackIsoVar(const CdfEmObject& emObject, double cone,
			double zvertex, int& nIsoTrk, double& maxIsoTrk) {
    double phi = phoPhi(emObject);
    double eta = phoEta(emObject,zvertex);
    nIsoTrk = -1;
    maxIsoTrk = -99.0;
    CdfTrackView_h trview;
    if(CdfTrackView::defTracks(trview,"PROD") != CdfTrackView::OK) {
      if(CdfTrackView::defTracks(trview,"USRP") != CdfTrackView::OK) {
	if(CdfTrackView::defTracks(trview,"UNKN") != CdfTrackView::OK) {
	  if(CdfTrackView::defTracks(trview,"L3") != CdfTrackView::OK) 
	    return -99.0;
	}
      }
    }
    return TrackingVariables::trkTrackIsolationInnerCore(phi, eta, zvertex, 
	       *trview, NULL, cone, 5.0, 0.0, 0, nIsoTrk, maxIsoTrk);

  }

  double phoTrackIsoVar(const CdfEmObject& emObject, 
			const CdfTrackView& trackView, double cone,
			double zvertex, int& nIsoTrk, double& maxIsoTrk) {
    double phi = phoPhi(emObject);
    double eta = phoEta(emObject,zvertex);
    nIsoTrk = -1;
    maxIsoTrk = -99.0;
    return TrackingVariables::trkTrackIsolationInnerCore(phi, eta, zvertex, 
	       trackView, NULL, cone, 5.0, 0.0, 0, nIsoTrk, maxIsoTrk);

  }


  // The number hits expected in the Silicon
  // the return value is (nSiHitsEx) | (nSiPhiHitsEx<<8) | (nSiHitsCone<<16)
  int phoSiExpected(const CdfEmObject& emObject) {
    return 0;
  }


  double phoSinTheta(const CdfEmObject& emObject, double zvertex) {
    
    HepPoint3D p;
    if(phoHasShowerMax(emObject)) {   // then there was a CES cluster
      p = phoShowerMaxPosition(emObject);
    } else {
      p = phoCentroidPosition(emObject);
    }
    p = p - HepPoint3D(0.,0.,zvertex);
    
    double c = p.cosTheta();   // no sinTheta()!
    double arg = 1.0-c*c;
    return ( arg >0 ? sqrt(arg) : 0.0 );
  }


  // the shower centroid according to the EM towers
  // use the ShowerMax position accessor for best value
  HepPoint3D phoCentroidPosition(const CdfEmObject& emObject) {
    const EmCluster* emc = emObject.getEmCluster();

    if(phoDetector(emObject)==PEM) {
      double p = emc->emPhi();
      double z = 0.5 * (TLZPESL0 + TLZPESL1);
      if(emc->side()==0) z = -z;
      double etaDet;
      if(emc->pem3x3FitTowers() != 0) {
	etaDet = emc->pem3x3FitDetEta();
      } else {
	etaDet = emc->emEtaDetector();
      }
      double cot = sinh(etaDet);
      double r = (cot==0? 0: z/cot);
      return HepPoint3D(r*cos(p),r*sin(p),z);
    } else {
      double z = emc->zCentroid();
      double p = emc->emPhi();
      double r = TLYCES;          // CES radius
      return HepPoint3D(r*cos(p),r*sin(p),z);
    }
  }

  // returns true if there are x and z CES clusters
  bool phoHasShowerMax(const CdfEmObject& emObject) {

    if(phoDetector(emObject)==PEM) {
      return (phoPes2dCluster(emObject)!=NULL);
    } else {
     return ( phoCesStripCluster(emObject)!=NULL && 
	      phoCesWireCluster(emObject)!=NULL );
    }
  }

  // shower position based on CES
  HepPoint3D phoShowerMaxPosition(const CdfEmObject& emObject) {

    if(!phoHasShowerMax(emObject)) return HepPoint3D(-999.0,-999.0,-999.0);
    if ( phoDetector(emObject)== PEM) {
      return phoPlugHepPoint3d(emObject);
    } else {
      double phi = phoCesWireCluster(emObject)->globalCoord();
      double localx = phoCesX(emObject);
      double r = sqrt(TLYCES*TLYCES + localx*localx); // TLYCES = CES radius
      double z = phoCesZ(emObject);
      return HepPoint3D(r*cos(phi), r*sin(phi), z);
    }

  }

    // CES local coordinates
  double phoCesX(const CdfEmObject& emObject) {
    const CesCluster* cesptr = 
                    phoCesWireCluster(emObject);
    if(cesptr==NULL) return -99.0;
    return cesptr->localCoord();
  }

    // CES local coordinates
  double phoCesZ(const CdfEmObject& emObject) {
    const CesCluster* cesptr = 
                phoCesStripCluster(emObject);
    if(cesptr==NULL) return -9999.0;

    /*Check and see which side of the detector the cluster is on, 0 -> eta<0 */
    /* 1 -> eta>0 and return the appropriate value for CesZ */
    if (cesptr->side()) {
    return cesptr->localCoord();
    } else {
      return -cesptr->localCoord();
    }
  }

    // CES energy (taken from wire view only)
  double phoCesE(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesWireCluster(emObject);
    if(cesptr==NULL) return -99.0;
    return cesptr->energy();
  }

    // CES strip energy
  double phoCesStripE(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesStripCluster(emObject);
    if(cesptr==NULL) return -99.0;
    double e = cesptr->energy();
    if(AbsEnv::instance()->monteFlag()) return e;
    double x = phoCesX(emObject);
    double z = phoCesZ(emObject);
    if(x<-98.0 || z<-9998.0) return e;
    z = fabs(z);
    const EmCluster* clus = emObject.getEmCluster();
    int wedge = clus->seedPhiAddr();
    int side =  clus->side();
    double corr = phoStripToWireEnergyCorr(wedge,side,x,z); 
    return corr*cesptr->energy();
  }


  //CES energy corrected, for now just a nominal factor
  double phoCesECorrected(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesWireCluster(emObject);
    if (cesptr==NULL) return -99.0;
    return (cesptr->energy())*1.3;
  }

    // CES chi**2
  double phoCesChiSq(const CdfEmObject& emObject) {
    const CesCluster* cesptrWire = phoCesWireCluster(emObject);
    if(cesptrWire==NULL) return -99.0;
    const CesCluster* cesptrStrip = phoCesStripCluster(emObject);
    if(cesptrStrip==NULL) return -99.0;
    //return (cesptrStrip->chi2_strip()+cesptrWire->chi2_strip())/2.0;
    double avgChi = (cesptrStrip->chi2()+cesptrWire->chi2())/2.0;
    // this correction comes from Electron/emobj_alg
    double corrChi = scaleCesChisq( avgChi, phoRawEnergy(emObject) );
    return corrChi;
  }


    // CES strip only chi**2
  double phoCesStrpChiSq(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesStripCluster(emObject);
    if(cesptr==NULL) return -99.0;
    double chi = cesptr->chi2();
    // this correction comes from Electron/emobj_alg
    double corrChi = scaleCesChisq( chi, phoRawEnergy(emObject) );
    return corrChi;
  }

    // CES wire only chi**2
  double phoCesWireChiSq(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesWireCluster(emObject);
    if(cesptr==NULL) return -99.0;
    double chi = cesptr->chi2();
    // this correction comes from Electron/emobj_alg
    double corrChi = scaleCesChisq( chi, phoRawEnergy(emObject) );
    return corrChi;
  }


    // the best ces wire cluster
  const CesCluster* phoCesWireCluster(const CdfEmObject& emObject){
    return phoCesCluster(emObject, CesCluster::WIRE, 1 );
  }

  const CesCluster* phoCesStripCluster(const CdfEmObject& emObject){
    return phoCesCluster(emObject, CesCluster::STRIP, 1 );
  }

    // find best wire or strip cluster, select by second argument
  const CesCluster* phoCesCluster(const CdfEmObject& emObject, 
		       CesCluster::CesView stripOrWire, int order){

    const CesCluster* cesptr = NULL;
    const CesCluster* cesptr2 = NULL;
    double e;
    double emax = -1;
    double emax2 = -1;

    // If no matching CES clusters, return a null link
    const CdfEmObject::Ces_CollType& assocCes = emObject.matchingCesClusters( "UNBIASED" );
    if (0 == assocCes.contents().size()) return cesptr;
    
    CdfEmObject::Ces_CollType::const_iterator ciIndxLinksCes(assocCes.contents().begin());
    CdfEmObject::Ces_CollType::const_iterator lastIndxLink(assocCes.contents().end());
    
    CdfEmObject::EmCluster_link emclst(emObject.matchingEmCluster());
    
    // We will require that the cluster selected as largest energy one also
    // be consistent with the CEM EmCluster seed tower.
    // The tolerance in Z is set by a global constant defined near top of this
    // file.  The tolerance in X is just that the CES cluster be in the same
    // half wedge as the seed tower since wires are split into two half wedges.
    for ( ; ciIndxLinksCes != lastIndxLink; ++ciIndxLinksCes) {
      const CdfEmObject::Ces_link pCesCluster = *ciIndxLinksCes;
      if (pCesCluster->view() == stripOrWire) {
 
       // Seed tower / CES cluster consistency check for wires
       if (stripOrWire == CesCluster::WIRE) {
	 if (0 == pCesCluster->half_wedge()) {
	   if ((0 == emclst->side()  &&
		emclst->seedEtaAddr() < MinEtaCEMW + 5)  ||
	       (1 == emclst->side()  &&
		emclst->seedEtaAddr() > MinEtaCEME + 4)) continue;
	 } else {
	   if ((0 == emclst->side()  &&
		emclst->seedEtaAddr() > MinEtaCEMW + 4)  ||
	       (1 == emclst->side()  &&
		emclst->seedEtaAddr() < MinEtaCEME + 5)) continue;
	 }
	 
         // Seed tower / CES cluster consistency check for strips
       } else {
	 TowerKey  seedTowerKey(emclst->towerIds().front());
	 const TowerGeometry*  pTowerGeom =
	   (CdfDetector::instance()->getCalorDetector()->getTowerGeometry());
	 double zMin = fabs(TLYCES * pTowerGeom->cotThetaMin( seedTowerKey ));
	 double zMax = fabs(TLYCES * pTowerGeom->cotThetaMax( seedTowerKey ));
	 // Have to swap zMin and zMax on West side
	 if (0 == emclst->side()) {
	   double zTemp = zMin;
	   zMin = zMax;
	   zMax = zTemp;
	 }
	 double zStrip = pCesCluster->fitted_position();
	 if (zStrip < zMin - 5.0  ||
	     zStrip > zMax + 5.0) continue;
       }
 
       e = pCesCluster->energy();

       if (e > emax) {
	 emax2 = emax;
	 cesptr2 = cesptr;
         cesptr = &(*pCesCluster);
         emax = e;
       } else if(e > emax2) {
         cesptr2 = &(*pCesCluster);
         emax2 = e;
       } // if curr>largest
      }  // requested strip or wire
      
    }

    return (order == 1 ? cesptr : cesptr2);
  }

  // variables for the second largest CES cluster
  // energy is reported separately since one may be missing
  double phoCesWireE2(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesCluster(emObject, CesCluster::WIRE, 2);
    if(cesptr==NULL) return -99.0;
    return cesptr->energy();
  }

  double phoCesStripE2(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesCluster(emObject, CesCluster::STRIP, 2);
    if(cesptr==NULL) return -99.0;
    double e = cesptr->energy();
    if(AbsEnv::instance()->monteFlag()) return e;
    double x = phoCesX(emObject);
    double z = phoCesZ2(emObject);
    if(x<-98.0 || z<-9998.0) return e;
    z = fabs(z);
    const EmCluster* clus = emObject.getEmCluster();
    int wedge = clus->seedPhiAddr();
    int side =  clus->side();
    double corr = phoStripToWireEnergyCorr(wedge,side,x,z); 
    return corr*cesptr->energy();
  }

  double phoCesX2(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesCluster(emObject, CesCluster::WIRE, 2);
    if(cesptr==NULL) return -99.0;
    return cesptr->localCoord();
  }

  double phoCesZ2(const CdfEmObject& emObject) {
    const CesCluster* cesptr = phoCesCluster(emObject, CesCluster::STRIP, 2);
    if(cesptr==NULL) return -9999.0;
    return cesptr->localCoord();
  }


// Run II CES strip-to-wire corrections
// implementation by  Anton Anastassov
double phoStripToWireEnergyCorr(int wedge, int side, double locX, double locZ)
{  
  const float parStrip[96][5]=
    {
      { 0.98927E+00,-0.72789E+00, 0.83092E+00, 0.15376E+01,-0.16665E+01},
      { 0.94156E+00, 0.11623E+01,-0.57364E+01, 0.99365E+01,-0.54566E+01},
      { 0.83466E+00, 0.13723E+01,-0.63057E+01, 0.96061E+01,-0.45566E+01},
      { 0.10428E+01,-0.58101E+00,-0.59273E+00, 0.34244E+01,-0.23008E+01},
      { 0.10752E+01,-0.18350E+01, 0.87342E+01,-0.13444E+02, 0.64237E+01},
      { 0.92379E+00, 0.15513E+01, 0.10889E+01,-0.55406E+01, 0.29358E+01},
      { 0.10323E+01, 0.24206E+00,-0.39890E+01, 0.91388E+01,-0.54926E+01},
      { 0.10167E+01, 0.83698E+00, 0.26662E+00,-0.16309E+01, 0.49134E+00},
      { 0.10412E+01, 0.66202E+00,-0.23451E+01, 0.43932E+01,-0.28348E+01},
      { 0.99916E+00,-0.10245E+01, 0.15430E+01, 0.36053E+00,-0.87852E+00},
      { 0.10817E+01,-0.81794E+00, 0.27638E+01,-0.23328E+01, 0.21600E+00},
      { 0.11385E+01,-0.18343E+01, 0.29822E+01,-0.12336E+01,-0.11167E+00},
      { 0.98685E+00,-0.47050E-02, 0.36977E-01, 0.77217E+00,-0.85445E+00},
      { 0.99318E+00,-0.52613E+00, 0.55471E+00, 0.14689E+01,-0.15045E+01},
      { 0.10200E+01,-0.76092E+00, 0.12957E+01, 0.10990E+00,-0.71795E+00},
      { 0.84253E+00,-0.24040E+00,-0.41351E+01, 0.92013E+01,-0.47557E+01},
      { 0.10234E+01,-0.19477E+01, 0.71961E+01,-0.97014E+01, 0.44851E+01},
      { 0.10477E+01,-0.18070E+01, 0.76925E+01,-0.10112E+02, 0.41983E+01},
      { 0.11158E+01,-0.33147E+00, 0.56600E+01,-0.11357E+02, 0.59802E+01},
      { 0.10545E+01,-0.13508E+01, 0.89643E+01,-0.14547E+02, 0.68950E+01},
      { 0.92812E+00, 0.34348E-01, 0.47719E+01,-0.10308E+02, 0.56518E+01},
      { 0.11070E+01,-0.28456E+01, 0.11661E+02,-0.16174E+02, 0.72460E+01},
      { 0.10405E+01,-0.19008E+00,-0.17872E+01, 0.40321E+01,-0.20663E+01},
      { 0.95256E+00, 0.22067E+01,-0.79974E+01, 0.89501E+01,-0.32218E+01},
      { 0.90762E+00, 0.42242E+01,-0.30835E+01,-0.66207E+01, 0.57449E+01},
      { 0.99587E+00, 0.16539E+01, 0.87925E+00,-0.57528E+01, 0.31765E+01},
      { 0.94610E+00,-0.42430E+00, 0.49536E+01,-0.79265E+01, 0.33986E+01},
      { 0.93755E+00, 0.32687E+01,-0.74024E+01, 0.52183E+01,-0.10626E+01},
      { 0.10378E+01, 0.40056E+00, 0.32893E+01,-0.64032E+01, 0.25354E+01},
      { 0.10155E+01, 0.50655E+00,-0.37420E+01, 0.63944E+01,-0.32463E+01},
      { 0.10935E+01,-0.12344E+01, 0.30010E+01,-0.37717E+01, 0.19661E+01},
      { 0.95837E+00, 0.27483E+01,-0.82024E+01, 0.86695E+01,-0.31957E+01},
      { 0.10377E+01, 0.73794E+00,-0.65297E+01, 0.13752E+02,-0.81819E+01},
      { 0.11332E+01,-0.34887E+01, 0.10744E+02,-0.12517E+02, 0.51309E+01},
      { 0.10031E+01,-0.12909E+01, 0.51972E+01,-0.75828E+01, 0.36672E+01},
      { 0.94969E+00,-0.57607E+00, 0.35291E+01,-0.52264E+01, 0.23135E+01},
      { 0.96809E+00,-0.18516E+00, 0.88589E+00,-0.26505E+00,-0.38087E+00},
      { 0.10353E+01, 0.15004E+01, 0.67913E+01,-0.18299E+02, 0.10082E+02},
      { 0.10512E+01,-0.28709E+01, 0.94787E+01,-0.10432E+02, 0.38206E+01},
      { 0.11321E+01,-0.18613E+00, 0.35865E+01,-0.65991E+01, 0.30551E+01},
      { 0.98576E+00, 0.35470E+00,-0.20281E+01, 0.51963E+01,-0.36502E+01},
      { 0.10382E+01,-0.97322E+00, 0.44857E+01,-0.59782E+01, 0.24103E+01},
      { 0.95063E+00,-0.66303E+00, 0.19444E+01,-0.10799E+01,-0.21045E+00},
      { 0.10673E+01,-0.42104E+00, 0.38712E+01,-0.65495E+01, 0.30592E+01},
      { 0.10659E+01, 0.12027E+01, 0.13632E+01,-0.57787E+01, 0.31707E+01},
      { 0.86085E+00, 0.13700E+01,-0.71213E+00,-0.21700E+01, 0.16738E+01},
      { 0.99506E+00, 0.20430E+00, 0.71564E+01,-0.14611E+02, 0.72151E+01},
      { 0.10149E+01,-0.14388E+01, 0.82624E+01,-0.12790E+02, 0.59559E+01},
      { 0.10552E+01,-0.13250E+01, 0.67305E+01,-0.47884E+01,-0.23068E+01},
      { 0.10144E+01,-0.15314E+01, 0.34527E+01, 0.16407E+01,-0.47991E+01},
      { 0.10020E+01,-0.12889E+00, 0.23416E+01,-0.27817E+01, 0.45152E+00},
      { 0.10412E+01,-0.85230E+00, 0.40170E+01,-0.73757E+01, 0.45427E+01},
      { 0.11416E+01,-0.64846E+00, 0.11791E+02,-0.26512E+02, 0.16991E+02},
      { 0.11044E+01,-0.94712E+00,-0.12905E+01, 0.86584E+01,-0.75167E+01},
      { 0.10913E+01,-0.72469E+00, 0.69306E-01, 0.48082E+01,-0.51368E+01},
      { 0.11550E+01,-0.42281E+01, 0.23547E+02,-0.45375E+02, 0.27691E+02},
      { 0.10790E+01,-0.10642E+00,-0.85369E+00, 0.73909E+01,-0.75295E+01},
      { 0.97801E+00, 0.27104E+00, 0.34896E+01,-0.10481E+02, 0.73478E+01},
      { 0.96482E+00, 0.21995E+01,-0.10250E+02, 0.20356E+02,-0.14114E+02},
      { 0.11462E+01,-0.21821E+01, 0.85516E+01,-0.11115E+02, 0.43069E+01},
      { 0.97531E+00,-0.18999E+00, 0.57105E+00,-0.11403E+01, 0.64286E+00},
      { 0.10287E+01, 0.63503E+00,-0.33713E+01, 0.10022E+02,-0.84370E+01},
      { 0.97172E+00, 0.18030E+01,-0.61123E+01, 0.13892E+02,-0.10930E+02},
      { 0.10929E+01,-0.72392E+00, 0.31122E+01,-0.32497E+01,-0.56932E-01},
      { 0.11048E+01,-0.16933E+01, 0.13012E+02,-0.25227E+02, 0.14388E+02},
      { 0.11459E+01,-0.76466E+00, 0.41916E+01,-0.45145E+01, 0.94634E-01},
      { 0.11132E+01,-0.57658E+00, 0.20561E+00, 0.19348E+01,-0.18563E+01},
      { 0.10681E+01,-0.90554E+00, 0.66167E+01,-0.12129E+02, 0.61939E+01},
      { 0.10233E+01, 0.37978E+00, 0.80626E-01, 0.10365E+01,-0.18212E+01},
      { 0.10592E+01,-0.21141E+01, 0.12367E+02,-0.20017E+02, 0.94538E+01},
      { 0.98210E+00, 0.87828E+00,-0.27029E+00,-0.32694E+01, 0.33547E+01},
      { 0.99625E+00, 0.15012E+01,-0.11209E+02, 0.26334E+02,-0.18372E+02},
      { 0.98738E+00, 0.18468E+00, 0.37655E+01,-0.10171E+02, 0.72753E+01},
      { 0.98944E+00,-0.37252E+00, 0.69208E-01, 0.40313E+01,-0.42588E+01},
      { 0.96699E+00,-0.20336E+00, 0.39630E+01,-0.88459E+01, 0.49565E+01},
      { 0.96384E+00, 0.22442E+00, 0.56968E+01,-0.16496E+02, 0.11431E+02},
      { 0.10549E+01,-0.12063E+01, 0.10323E+02,-0.15726E+02, 0.70324E+01},
      { 0.92561E+00, 0.20679E+00, 0.20544E+01,-0.79255E+01, 0.52314E+02},
      { 0.11443E+01,-0.84855E+00, 0.52197E+01,-0.10425E+02, 0.64171E+01},
      { 0.10687E+01,-0.14030E+01, 0.65939E+01,-0.10504E+02, 0.56451E+01},
      { 0.99060E+00, 0.67261E+00,-0.37280E+01, 0.64528E+01,-0.35660E+01},
      { 0.92246E+00, 0.23544E+01,-0.72280E+01, 0.73911E+01,-0.22879E+01},
      { 0.95429E+00, 0.14165E+01,-0.90442E+01, 0.23096E+02,-0.17504E+02},
      { 0.99243E+00,-0.52551E+00,-0.82068E+00, 0.16205E+01,-0.45873E-01},
      { 0.11231E+01, 0.33134E+00,-0.45179E+01, 0.11757E+02,-0.88601E+01},
      { 0.10490E+01,-0.14828E+01, 0.11809E+01, 0.50873E+01,-0.64088E+01},
      { 0.10309E+01,-0.28839E+00,-0.28975E+01, 0.76593E+01,-0.53057E+01},
      { 0.96774E+00, 0.37221E+00,-0.41721E+00,-0.13831E+01, 0.16910E+01},
      { 0.96776E+00, 0.22864E+00,-0.19311E+01, 0.41923E+01,-0.27102E+01},
      { 0.10984E+01,-0.16618E+01, 0.74613E+01,-0.10938E+02, 0.46438E+01},
      { 0.10580E+01,-0.47709E+00, 0.85925E+01,-0.19906E+02, 0.12149E+02},
      { 0.78031E+00, 0.23712E+01,-0.10094E+02, 0.18109E+02,-0.11093E+02},
      { 0.96844E+00, 0.36185E+00, 0.69727E+01,-0.16945E+02, 0.99676E+01},
      { 0.10329E+01,-0.13939E+01, 0.15143E+02,-0.31062E+02, 0.17660E+02},
      { 0.10048E+01,-0.73044E+00, 0.69556E+01,-0.15443E+02, 0.10001E+02},
      { 0.88110E+00, 0.19613E+01,-0.84979E+01, 0.19005E+02,-0.13941E+02}
    };
  const float parWire[96][5]=
    {
      { 0.49670E-01,-0.91118E+00, 0.46899E+01,-0.80364E+01, 0.42591E+01},
      { 0.35005E-01,-0.11172E+01, 0.50649E+01,-0.75352E+01, 0.35958E+01},
      { 0.37137E-02,-0.58745E+00, 0.28059E+01,-0.47251E+01, 0.26803E+01},
      {-0.34112E-01, 0.62521E-01, 0.52929E+00,-0.13322E+01, 0.79745E+00},
      { 0.67750E-01,-0.12134E+01, 0.61749E+01,-0.10298E+02, 0.52755E+01},
      {-0.25271E+00, 0.64957E-01, 0.28302E+01,-0.49669E+01, 0.25929E+01},
      { 0.44928E-01,-0.83844E+00, 0.32431E+01,-0.42590E+01, 0.17902E+01},
      { 0.26135E-01,-0.19073E+01, 0.83196E+01,-0.12774E+02, 0.67398E+01},
      { 0.56553E-01,-0.12240E+01, 0.51999E+01,-0.77802E+01, 0.38366E+01},
      {-0.11769E-01,-0.10194E+01, 0.43258E+01,-0.54243E+01, 0.21783E+01},
      {-0.37892E-01,-0.44243E+00, 0.89656E+00, 0.97671E+00,-0.14198E+01},
      { 0.18833E+00,-0.24430E+01, 0.72778E+01,-0.82853E+01, 0.33935E+01},
      { 0.18336E+00,-0.22808E+01, 0.79046E+01,-0.11016E+02, 0.53809E+01},
      {-0.86884E-01, 0.24352E+00,-0.61872E+00, 0.57860E+00, 0.13375E+00},
      { 0.23566E+00,-0.21555E+01, 0.62023E+01,-0.78864E+01, 0.37320E+01},
      { 0.11706E-01,-0.15301E+00,-0.15750E-01, 0.51718E+00,-0.29615E+00},
      {-0.98427E+00, 0.82864E+01,-0.23823E+02, 0.30009E+02,-0.13604E+02},
      {-0.10939E+00, 0.76631E+00,-0.35289E+01, 0.63686E+01,-0.34485E+01},
      {-0.22394E-01,-0.26912E+01, 0.12135E+02,-0.18431E+02, 0.96556E+01},
      { 0.70632E-01,-0.98542E+00, 0.28823E+01,-0.33571E+01, 0.15033E+01},
      {-0.68300E-01,-0.56080E+00, 0.52631E+01,-0.10861E+02, 0.65494E+01},
      {-0.12780E+00, 0.74661E+00,-0.61872E+00,-0.39715E+00, 0.29984E+00},
      { 0.14863E+00,-0.10075E+01, 0.19391E+01,-0.15218E+01, 0.44613E+00},
      {-0.79754E-01,-0.57837E+00, 0.52696E+01,-0.96799E+01, 0.51624E+01},
      { 0.53998E-01,-0.73894E+00, 0.36979E+01,-0.69622E+01, 0.41173E+01},
      {-0.56415E-01,-0.30018E+00, 0.25088E+01,-0.44932E+01, 0.24680E+01},
      { 0.33093E-01, 0.40520E-01,-0.11378E+01, 0.20191E+01,-0.89417E+00},
      { 0.13594E+00,-0.14767E+01, 0.14648E+01, 0.43981E+01,-0.49277E+01},
      { 0.26410E+00,-0.18941E+01, 0.57817E+01,-0.81457E+01, 0.39611E+01},
      {-0.22174E+00, 0.20635E+00, 0.47960E+00, 0.40228E+00,-0.70950E+00},
      { 0.79609E-01,-0.14764E+01, 0.58775E+01,-0.83890E+01, 0.39841E+01},
      { 0.16332E-01,-0.12509E+01, 0.62444E+01,-0.96726E+01, 0.47295E+01},
      {-0.96606E-01,-0.83100E+00, 0.60561E+01,-0.10151E+02, 0.51561E+01},
      { 0.56823E-01,-0.26450E+00, 0.33857E-01, 0.43756E+00,-0.22616E+00},
      { 0.17152E-02, 0.10256E-01,-0.14976E+01, 0.39518E+01,-0.24781E+01},
      {-0.60160E-01, 0.26912E+00,-0.27672E+00, 0.27495E+00,-0.25464E+00},
      {-0.11676E+00, 0.64737E+00,-0.29013E+01, 0.59247E+01,-0.36055E+01},
      {-0.28350E-01, 0.48997E+00,-0.23575E+01, 0.40377E+01,-0.22016E+01},
      { 0.19359E-01,-0.42799E+00, 0.10681E+01,-0.10967E+01, 0.56403E+00},
      { 0.33851E+00,-0.19537E+01, 0.55933E+01,-0.82174E+01, 0.41421E+01},
      {-0.70672E-01,-0.21541E-01, 0.35882E+01,-0.86409E+01, 0.52287E+01},
      { 0.10793E+00,-0.11017E+01, 0.33883E+01,-0.43105E+01, 0.19561E+01},
      { 0.76288E-01, 0.13369E+00,-0.19848E+01, 0.31249E+01,-0.13138E+01},
      { 0.32135E+00,-0.23189E+01, 0.68914E+01,-0.88888E+01, 0.38163E+01},
      {-0.65295E-01,-0.22657E+01, 0.11106E+02,-0.17636E+02, 0.95284E+01},
      { 0.37850E+00,-0.35580E+01, 0.10732E+02,-0.14789E+02, 0.76042E+01},
      { 0.14216E+00,-0.10566E+01, 0.18925E+01,-0.10598E+01, 0.10152E+00},
      { 0.20978E+00,-0.10150E+01, 0.14906E+01,-0.90718E+00, 0.13818E+00},
      { 0.31781E-01,-0.65668E+00, 0.30078E+01,-0.48525E+01, 0.25356E+01},
      { 0.56025E-01,-0.13198E+01, 0.61353E+01,-0.95233E+01, 0.46985E+01},
      {-0.48593E-01, 0.48500E+00,-0.20200E+01, 0.33211E+01,-0.17544E+01},
      { 0.72377E-01,-0.97329E+00, 0.44274E+01,-0.68816E+01, 0.32948E+01},
      { 0.48156E-01,-0.33353E+00,-0.33753E+01, 0.11695E+02,-0.84011E+01},
      {-0.33645E-01, 0.33773E+00,-0.86371E+00, 0.91303E+00,-0.37800E+00},
      { 0.62599E-01,-0.19353E+00, 0.64924E+00,-0.17067E+01, 0.12222E+01},
      {-0.18209E+00, 0.14527E+01,-0.36761E+01, 0.40082E+01,-0.16048E+01},
      { 0.14681E+00,-0.18888E+01, 0.52428E+01,-0.57565E+01, 0.24464E+01},
      { 0.25977E-01,-0.13965E+01, 0.62148E+01,-0.88135E+01, 0.40205E+01},
      { 0.33375E+00,-0.31867E+01, 0.97910E+01,-0.12086E+02, 0.50876E+01},
      { 0.29760E+00,-0.21336E+01, 0.57442E+01,-0.72264E+01, 0.33061E+01},
      { 0.13306E+00,-0.15375E+01, 0.60835E+01,-0.88990E+01, 0.41637E+01},
      {-0.31583E-01, 0.35776E+00,-0.10046E+01, 0.74813E+00, 0.27622E-02},
      {-0.27540E+00, 0.29228E+01,-0.73562E+01, 0.71700E+01,-0.26330E+01},
      {-0.25149E-01, 0.28187E+00,-0.11292E+00,-0.16144E+01, 0.16276E+01},
      {-0.18448E+00, 0.11969E+01,-0.25972E+01, 0.28842E+01,-0.13467E+01},
      { 0.17909E+00,-0.14645E+00,-0.16447E+01, 0.30201E+01,-0.15633E+01},
      { 0.74148E-01, 0.15719E+00,-0.21571E+01, 0.33400E+01,-0.13435E+01},
      { 0.38023E+00,-0.34372E+01, 0.91859E+01,-0.10375E+02, 0.43518E+01},
      { 0.70042E-01,-0.80325E+00, 0.51045E+01,-0.10898E+02, 0.67736E+01},
      {-0.31468E+00, 0.22350E+01,-0.46672E+01, 0.44691E+01,-0.18224E+01},
      {-0.54437E-01,-0.92869E+00, 0.44661E+01,-0.64237E+01, 0.31807E+01},
      { 0.34238E-01,-0.73981E+00, 0.48708E+01,-0.10146E+02, 0.62433E+01},
      {-0.59331E-01, 0.52402E-01, 0.11166E+01,-0.25705E+01, 0.15180E+01},
      { 0.32658E-01,-0.31759E+00, 0.15426E+01,-0.30981E+01, 0.19325E+01},
      {-0.91996E-01, 0.26008E+00,-0.21287E+01, 0.51086E+01,-0.30285E+01},
      {-0.92926E-01,-0.83337E+00, 0.50068E+01,-0.74741E+01, 0.35463E+01},
      {-0.60207E-01,-0.17143E+01, 0.93422E+01,-0.14487E+02, 0.71258E+01},
      { 0.18851E+00,-0.13169E+01, 0.38534E+01,-0.49623E+01, 0.21305E+01},
      { 0.81822E-01,-0.14274E+01, 0.57815E+01,-0.82225E+01, 0.38021E+01},
      {-0.16992E-01,-0.74765E+00, 0.46636E+01,-0.81165E+01, 0.43273E+01},
      {-0.11806E+00, 0.12006E+01,-0.45541E+01, 0.75283E+01,-0.42320E+01},
      { 0.24968E-01,-0.16607E+01, 0.63237E+01,-0.76593E+01, 0.30619E+01},
      { 0.13903E+00,-0.16721E+01, 0.56669E+01,-0.81946E+01, 0.42843E+01},
      { 0.10918E+00,-0.16271E+01, 0.65202E+01,-0.95666E+01, 0.46138E+01},
      {-0.14424E+00, 0.80510E+00,-0.33992E+01, 0.68048E+01,-0.41327E+01},
      {-0.25867E-01,-0.58949E+00, 0.25023E+01,-0.30289E+01, 0.12188E+01},
      { 0.17979E+00,-0.17392E+01, 0.51650E+01,-0.64943E+01, 0.29590E+01},
      { 0.31521E-01,-0.15899E+01, 0.76829E+01,-0.11760E+02, 0.57133E+01},
      {-0.90371E-01, 0.19170E+00,-0.11854E+00, 0.32643E+00,-0.23790E+00},
      { 0.61850E-01,-0.18795E+01, 0.85936E+01,-0.12948E+02, 0.62536E+01},
      {-0.60892E-01, 0.56664E+00,-0.25711E+01, 0.41328E+01,-0.19931E+01},
      { 0.38760E-02,-0.17341E-01,-0.11125E+01, 0.31157E+01,-0.20168E+01},
      { 0.45335E+00,-0.34381E+01, 0.99217E+01,-0.12772E+02, 0.57586E+01},
      { 0.55769E+00,-0.41962E+01, 0.11988E+02,-0.16556E+02, 0.84184E+01},
      { 0.75082E-01,-0.12301E+01, 0.46392E+01,-0.66870E+01, 0.33269E+01},
      { 0.16416E+00,-0.12534E+01, 0.28012E+01,-0.28846E+01, 0.12503E+01}
    };


  // Do not apply corrections for MC events

  const int   nFitPar = 5;

  const float xFid    =  21.5;
  const float zFidMin =   9.0;
  const float zFidMax = 220.0;

  float zMin = 6.16;
  float zMax = 121.183;

  int  halfBarrel = 0;   // half-barrel: near->0, far->1
  if (locZ > zMax) {
    halfBarrel = 1;
    zMin       = zMax;
    zMax       = zMin + 59.0*2.007;
  }

  int halfWedgeId = wedge + side*24 + halfBarrel*48;   // global CES module ID: numbered in the range 0-95

  float x = -locX;   // sign convention used in the fits is different from CES cluster local coordinates
  if (x < -xFid)   x = -xFid; 
  if (x >  xFid)   x =  xFid;

  float z =  locZ;
  if (z < zFidMin) z =  zFidMin;
  if (z > zFidMax) z =  zFidMax;
  if (halfWedgeId==77 && z>160.0) z = 160.0;      // special case: East side, wedge 5, "far" half-barrel

  float xLocRel = (x + 23.27744) / (2.0*23.27744);
  float zLocRel = (z - zMin) / (zMax - zMin);

  float fX  = parWire [halfWedgeId][0];
  float fZ  = parStrip[halfWedgeId][0];
  float argX = 1.0;
  float argZ = 1.0;

  for (int i=1; i<nFitPar; ++i) {
    argX *= xLocRel;
    argZ *= zLocRel;
    fX   += argX * parWire [halfWedgeId][i];
    fZ   += argZ * parStrip[halfWedgeId][i];
  }


  float swMulti = 1.0;
  if (fX + fZ > 0.0) swMulti = 1.0/(fX + fZ);

  return swMulti;

}



  // a number (1+eps), includes all final energy calibrations available
  // you must divide raw energy by this correction
  double phoEnergyCorrection(const CdfEmObject& emObject) {
        // but no corrections defined yet
    if(!phoHasShowerMax(emObject)) return 1.0;
    if ( phoDetector(emObject)== PEM) {
      float pemcorr = 1.0;
      PemCorrAlg::GetPemCorr(emObject,pemcorr);
      return pemcorr;
    } else {
      double cemcorr = 1.0;
      CemCorrAlg::GetCemCorr( emObject, "USEENERGY", cemcorr);
      return cemcorr;
    }
  }


} // namespace PhotonVariables