#ifndef _PhotonVariables_HH_
#define _PhotonVariables_HH_ 1
////////////////////////////////////////////////////////////////////
//   File:  photonVariables.hh
//   Authors: Mike Kirby, Ray Culbertson
//   Purpose: High level analysis photon variables defined as functions
//   Created: 10/6/2000
//
//   Revision: 0.1
//   Revision: 1.0 Add CES and CPR bg weights, Ray Culbertson, 2/17/01
//   Revision: 2.0 Add PES quantities Mike Kirby, 3/19/01
//   Revision: 3.0 Add tower indices, Ray, 7/27/01
//   Revision: 4.0 Add hadron timing, Ray, 
//   Revision: 5.0 correction to CES chi2, Ray, 8/1/01
//   Revision: 6.0 Add CES variables, Ray, 8/11/01
//   Revision: 7.0 Add PES, PPR variables, Ray, 8/18/01
//   Revision: 8.0 Add deltaz cut to tracking iso, Ray, 10/??/01
//   Revision: 9.0 Fix LShr, add 3x3 chi2, Ray, 11/30/01
//   Revision:10.0 remove ces-seeded from stat word
//   Revision:11.0 fix trigger tower Et for z vertex
//   Revision:12.0 Fix ces and cpr weights, Ray, 2/20/01
//   Revision:15.0 tweak IsoLCor, add PES variables, Ray, 10/21/03
////////////////////////////////////////////////////////////////////

//#include "ElectronObjects/CdfEmObject.hh"
#include "CalorGeometry/CalConstants.hh"
#include "CalorObjects/Pes2dCluster.hh"
#include "CalorObjects/CesCluster.hh"
#include "CLHEP/Geometry/Point3D.h"
#include "CLHEP/Vector/LorentzVector.h"
// SLL Oct 06, 2003
#include "Electron/emobj_alg.hh"
class CesCluster;
class CprCluster;
class HepPoint3D;
class HepLorentzVector;
class CdfEmObject;
class CprWireCollectionMaker;

namespace PhotonVariables {

    // set lastVersion to the last version you know of to print changes
  int phoVariablesVersion(int lastVersion=999);

    // bits set for each cut passed
  int phoCutWord(const CdfEmObject& emObject, const double zvertex =0);

    // status word = 100, 200... + photon Et
  double phoStatus(const CdfEmObject& emObject, const double zvertex =0);


    // Detector region (CEM=0, PEM=1)  "CalorGeometry/CalConstants.hh"
  Detector phoDetector(const CdfEmObject& emObject);

    //the Et directly from z=0 and EmCluster centroid (for trigger)
  double phoTowerEt(const CdfEmObject& emObject);

    //the Et directly from the EmCluster method
  double phoRawEt(const CdfEmObject& emObject, const double zvertex =0 );

    //the Et corected for final offline calibrations
  double phoCorrectedEt(const CdfEmObject& emObject, 
			const double zvertex =0);

    //four vector using corrected energy
  HepLorentzVector phoFourVector(const CdfEmObject& emObject, 
			       const double zvertex = 0);
    // raw Energy
  double phoRawEnergy(const CdfEmObject& emObject);

    // corrected Energy
  double phoCorrectedEnergy(const CdfEmObject& emObject);

    // the phi index of the seed tower
  int phoSeedPhi(const CdfEmObject& emObject);
    // the eta index of the seed tower
  int phoSeedEta(const CdfEmObject& emObject);

    // the phi using only towers
  double phoTowerPhi(const CdfEmObject& emObject);
    // the eta using only towers
  double phoTowerEta(const CdfEmObject& emObject);

    // the tower usage pattern
    // CEM: bit0 = tow below seed, bit 1=seed, bit 2= above seed (count 0-51)  
  int phoTowerPattern(const CdfEmObject& emObject);

    // shower max phi, if avail, tower centroid otherwise
  double phoPhi(const CdfEmObject& emObject);

    // showermax event eta, if avail, tower centroid otherwise
  double phoEta(const CdfEmObject& emObject, double zvertex = 0);
    // showermax deta, if avail, tower centroid otherwise
  double phoDteta(const CdfEmObject& emObject);

    // had/em ("3-tower")
  double phoHadEm(const CdfEmObject& emObject);

    // had/em from trigger tower ("2-tower")
  double phoHadEmT(const CdfEmObject& emObject);

    // LShr
  double phoLshr(const CdfEmObject& emObject, double zvertex = 0, 
		 int nTowers = 3);

    // cone isolation energies
  double phoCo4(const CdfEmObject& emObject, const double zvertex);
  double phoCo7(const CdfEmObject& emObject, const double zvertex);

    // pt of highest pt track (order=1 is highest, 2 is next highest)
  double phoPt(const CdfEmObject& emObject, int order =1 );
    // number of tracks pointing at cluster
  int phoN3D(const CdfEmObject& emObject);

    // Track isolation in a cone
  double phoTrackIso(const CdfEmObject& emObject, double cone = 0.4,
		     double zvertex =0);

    // Track isolation in a cone, specified track view
  double phoTrackIso(const CdfEmObject& emObject, 
		     const CdfTrackView& trackView, double cone = 0.4,
		     double zvertex =0);

    // 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);

    // Track isolation in a cone, same as above, but also return
    // number of tracks in the cone and Pt of max Pt track
    // specified track view
  double phoTrackIsoVar(const CdfEmObject& emObject, 
			const CdfTrackView& trackView, double cone,
			double zvertex, int& nIsoTrk, double& maxIsoTrk);

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

    // Shower position based on towers centroids
    // for best value, use CES position below
  HepPoint3D phoCentroidPosition(const CdfEmObject& emObject);

  // returns true if there are x and z CES or u and v PES clusters
  bool phoHasShowerMax(const CdfEmObject& emObject);

  // CES or PES position
  HepPoint3D phoShowerMaxPosition(const CdfEmObject& emObject);
  double phoSinTheta(const CdfEmObject& emObject, double zvertex =0);

  // CES local coordinates
  // note that we have flipped Z as necessary so that it is CDF z coords
  double phoCesX(const CdfEmObject& emObject);
  double phoCesZ(const CdfEmObject& emObject);

  // CES energy (taken from wire cluster only)
  double phoCesE(const CdfEmObject& emObject);
  double phoCesECorrected(const CdfEmObject& emObject);
  double phoCesStripE(const CdfEmObject& emObject);

  // CES chi**2 (returns -99 if strip or wire is missing)
  // this is the average of strip and wire chi**2
  double phoCesChiSq(const CdfEmObject& emObject);
  double phoCesStrpChiSq(const CdfEmObject& emObject);
  double phoCesWireChiSq(const CdfEmObject& emObject);

  const CesCluster* phoCesWireCluster(const CdfEmObject& emObject);
  const CesCluster* phoCesStripCluster(const CdfEmObject& emObject);

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

  // variables for the second largest CES cluster
  // energy is reported separately since one may be missing
  double phoCesWireE2(const CdfEmObject& emObject);
  double phoCesStripE2(const CdfEmObject& emObject);
  double phoCesX2(const CdfEmObject& emObject);
  double phoCesZ2(const CdfEmObject& emObject);

  // for CES, given a position, return a factor
  // which should multiply the strip e to make it 
  // calibrated to the wire e
  // x is ces cluster x (-phi), and z is fabs(CDF z)
  double phoStripToWireEnergyCorr(int wedge, int side, 
				  double locX, double locZ);
  

    // a number (1+eps), includes all final energy calibrations available
  double phoEnergyCorrection(const CdfEmObject& emObject);


  // These are Run I quantities, need to be tuned for Run II

  // this is Run I sliding cut on the energy of the second 
  // CES wire cluster in the wedge
  double phoCesSlide(const CdfEmObject& emObject);

  // this returns the Run I correction to isolation energy
  // due to phi-crack leakage, from Peter Wilson
  // used:  CorIsoEtCone0.4 = IsoEtCone0.4 - LCor*PhotonCorrectedEt
  double phoLCor(const CdfEmObject& emObject);

  // this returns the Run I correction to isolation energy
  // due to multiple underlying events, from Peter Wilson
  // user substracts this energy: CorIsoEtCone0.4 = IsoEtCone0.4 - LCor;
  double phoVCor(const CdfEmObject& emObject);

  // this returns the Run I cone 0.4 isolation energy
  // with VCor and Lcor corrections (see above)
  double phoCo4PJW(const CdfEmObject& emObject, const double zvertex);


  // Ces energy for all wires in the wedge, call with
  // double cesWireArr[64]
  // the first 32 are the wires closest to |z|=0, 32-63 are the ones 
  // further from zero.  Order is decreasing phi on west, increasing on east.
  // double cesStripArr[128]   order is from low |z| to high
  // return value is 0 for all OK
  int phoCesStripAndWireE(const CdfEmObject& emObject, 
			  double* cesWireArr, double* cesStripArr);

  // the CPR local x and charge of the highest energy cluster in this wedge
  // return value is 0 if all OK
  // cpr local x is in +phi sense
  int phoCprUnbiased(const CdfEmObject& emObject, 
		     double& CprX, double& CprQ);

  // return the cluster itself
  const CprCluster* phoCprUnbiasedCluster(const CdfEmObject& emObject);

  // This returns the CPR cluster seeded from the photon's CES cluster
  // two variables here in order to not have to do the clustering twice
  // nstrips is the size of the Cpr cluster to consider
  // return value is 0 if all OK
  // if the CprWireCollectionMaker is left null , then it will be made 
  // for each event, which requires db connection.  For faster access,
  // once:
  // CprWireCollectionMaker maker;
  // begin job:
  // maker.initJob(false)
  // begin run:
  // maker.initRun(false)
  // per event:
  // maker.initCollection(AbsEvent* anEvent);
  // this should really be done in production
  // cpr local x is in +phi sense
  int phoCprCesSeeded(const CdfEmObject& emObject, double& CprX, 
		       double& CprQ, const int nstrips=5, 
		       const double zvertex = 0,
		       CprWireCollectionMaker* pMaker = NULL);

  // these return the CES and CPR background weights
  // the CES weight is a positive number for events with avg chi^2 <4,
  // and negative for events with 4<chi^2<20
  // the CPR weight is positive for event without a cpr clusters
  // and neagitve otherwise.  There are details to using these such as the
  // ID and fiducial cuts on the photon.  returns -99 for plug photons
  // or if no CES cluster or other failure
  double phoCesBgWeight(const CdfEmObject& emObject, const double zvertex = 0);
  double phoCprBgWeight(const CdfEmObject& emObject, const double zvertex = 0);
  double phoCesProbGam(const CdfEmObject& emObject, const double zvertex = 0);
  double phoCesProbBgd(const CdfEmObject& emObject, const double zvertex = 0);
  double phoCprProbGam(const CdfEmObject& emObject, const double zvertex = 0);
  double phoCprProbBgd(const CdfEmObject& emObject, const double zvertex = 0);
  // quantity needed for Cpr method
  double phoSumetNoJets();

  // distance between expected Cpr hit and nearest track
  double phoCprTrack(const CdfEmObject& emObject, const double zvertex = 0);

  // returns emObject hadron TDC time, averaged over all firing TDCs
  // -9999 if no TDC fired for the cluster
  double phoTime(const CdfEmObject& emObject);

  // returns emObject time, averaged over all firing TDCs in detector det
  // -9998 if det != CHA, WHA or PHA
  // -9999 if no timing in det
  double phoDetTime(const CdfEmObject& emObject, Detector det, int& ntow);

  // returns a likelihood that this is consistent with a single EM cluster
  double phoLikelihood(const CdfEmObject& emObject, LikelihoodVersion ver=UPDATED);

  // return value is number of muon stubs not associated with jets or tracks
  // nPhoMatchStub is stubs within 30deg of photon and not with track or jet
  // nStub is all stubs
  // input argument muonSystem is  bit flag for the allowed muon systems:
  // bit 0=CMU, 1=CMP, 2=CMX, 3=BMU
  int phoCosStub(const CdfEmObject& emObject, int& nCosStubPho, 
		 int& nStub, int muonSystem=7);

  // for a central emObject, returns:
  // seedwedge: number of central towers above threshold in same wedge as
  //            emObject seed tower
  // sidewedges: number of central towers above threshold in wedges on each
  //             side of seed tower wedge
  // ncontig: number of contiguous towers above threshold in same wedge as
  //          seed tower (there is allowed to be a gap of one tower)
  // threshold: E cutoff to consider a tower occupied
  void phoBeamHaloOccupancy(const CdfEmObject& emObject, int& seedwedge,
			    int& sidewedges, int& ncontig,
			    float threshold=0.1);

  //hadron tower occupancy in regions where a beam halo muon would cross
  // east and west are the number of hit towers on each side
  void phoBHHadOccupancy(const CdfEmObject& emObject, int& east, int& west);

  // 2-d cluster selected for photons, NULL if no cluster
  // order = 1 for best match, 2 for second best match
  const Pes2dCluster* phoPes2dCluster(const CdfEmObject& emObject, 
				      int order =1);

  const PesCluster* phoPesClusterU(const CdfEmObject& emObject);
  const PesCluster* phoPesClusterV(const CdfEmObject& emObject);

  //plug PES energy
  double phoPesE(const CdfEmObject& emObject);

  // chi2 of fit to shower shape in PEM towers
  double phoPem3x3ChiSq(const CdfEmObject& emObject);

  // 5/9 chi2 of individual clusters
  double phoPesChi2ClusterU(const CdfEmObject& emObject);
  double phoPesChi2ClusterV(const CdfEmObject& emObject);

  // Delta-R PES-PEM
  double phoPesDeltaR(const CdfEmObject& emObject);

  // space point of PES 2-d cluster, with variations
  HepPoint3D phoPlugHepPoint3d(const CdfEmObject& emObject);
  double phoPesZ(const CdfEmObject& emObject);

  double phoPes2dX(const CdfEmObject& emObject);
  double phoPes2dY(const CdfEmObject& emObject);
  double phoPesErrorX(const CdfEmObject& emObject);
  double phoPesErrorY(const CdfEmObject& emObject);

  //after excluding the U and V 1d clusters from the best matching 2d
  //cluster, loop over U and V 1d clusters associated to the photon
  //and return the energy of highest energy cluster found.
  //As arguments, also return the same selection for U and V separately
  double phoPesE2(const CdfEmObject& emObject, double& eUMax2, double& eVMax2);

  //among the towers associated to the cluster, return the total 
  //PPR energy.  As arguments, also return the tower with the highest energy
  double phoPprEnergy(const CdfEmObject& emObject, double& ePprTowMax);

  // Number of tracks whose initial parameters would make them hit the
  // photon seed tower if they did not curve
  int phoNTracksInitPoint(const CdfEmObject& emObject, std::string processName = "PROD");
  int phoNTracksInitPoint(const CdfEmObject& emObject, CdfTrackView_h trview);
}

#endif