////////////////////////////////////////////////////////////////////
//
//     File: PhotonVariables2.cc
//     Authors: Mike Kirby, Ray Culbertson
//     Description: High level analysis photon variables defined as functions
//     Created: 1/31/2001
//
//     Revision: 0.0
//
////////////////////////////////////////////////////////////////////

//#include <cstdio>
//#include <math.h>
#include "Rtypes.h"
#include "HighLevelObjects/PhotonVariables.hh"
#include "CalorGeometry/CalParameters.hh"
#include "ElectronObjects/CdfEmObject.hh"
#include "CalorObjects/CalData.hh"
#include "CalorObjects/CesCluster.hh"
#include "CalorObjects/Pes2dCluster.hh"
#include "CLHEP/Geometry/Point3D.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Tracks/track_sort.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Tracks/SimpleReconstructedTrack.hh"
#include "Trajectory/Helix.hh"

namespace PhotonVariables { 


  const Pes2dCluster* phoPes2dCluster(const CdfEmObject& emObject, int order)
  {

    //get the unbiased, highest energy PES cluster

    const Pes2dCluster* pesptr  = NULL;    
    const Pes2dCluster* pesptr2 = NULL;

    std::vector<CdfEmObject::Pes_link>::const_iterator ci,first,last;
    first = emObject.matchingPes2dClusters().contents().begin();
    last = emObject.matchingPes2dClusters().contents().end();
    if (first==last) {  //there are no clusters
      return NULL;
    }  

    double emax = -1.0;
    double emax2 = -1.0;
    double e;
    for (ci=first; ci != last; ci++){
      //e = (*ci)->energy(); // doesn't work
      e = ((*ci)->pesClusterU().energy() + (*ci)->pesClusterV().energy() )/2.0;
      if (e>emax){
	emax2 = emax;
	emax = e;
	pesptr2 = pesptr;
 	pesptr = (*ci).operator->();
      } else if (e>emax2){
	
	emax2 = e;
	pesptr2 = (*ci).operator->();
      }
    }
    return (order ==1 ? pesptr : pesptr2);
  }

  const PesCluster* phoPesClusterU(const CdfEmObject& emObject)
  {
    const Pes2dCluster* p2 = phoPes2dCluster(emObject);
    if(p2==NULL) return NULL;
    return &(p2->pesClusterU());
  }

  const PesCluster* phoPesClusterV(const CdfEmObject& emObject)
  {
    const Pes2dCluster* p2 = phoPes2dCluster(emObject);
    if(p2==NULL) return NULL;
    return &(p2->pesClusterV());
  }

  HepPoint3D phoPlugHepPoint3d(const CdfEmObject& emObject) {

    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr !=NULL){
      return pesptr->plug2dClusterLocation();
    } else {
      return HepPoint3D(-999.,-999.,-999.);
    }
  }


  double phoPesZ(const CdfEmObject& emObject) {
   
    double globalZ = -999.0;
    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr != NULL) {
      globalZ = pesptr->plug2dClusterZ();
    }

    return globalZ;
  }


  double phoPes2dX( const CdfEmObject& emObject) {

    double localX = -999.0;
    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr != NULL) {
      localX = pesptr->plug2dClusterX();
    }
    return localX;
  }

  double phoPes2dY( const CdfEmObject& emObject) {

    double localY = -999.0;
    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr != NULL) {
      localY = pesptr->plug2dClusterY();
    }
    return localY;
  }
    

  double phoPesE(const CdfEmObject& emObject) {

    double pesE = -99.0;
    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr != NULL) {
      //pesE = pesptr->energy();  //doesn't currently work
      pesE= ( pesptr->pesClusterU().energy() + 
	      pesptr->pesClusterV().energy() )/2.0;
    }

    return pesE;
  }
    
  double phoPem3x3ChiSq(const CdfEmObject& emObject) {
    if(phoDetector(emObject)!=PEM) return 99.0;
    return emObject.getEmCluster()->pem3x3Chisq();
  }

  // 5/9 chi2 of individual clusters
  double phoPesChi2ClusterU(const CdfEmObject& emObject){
    const PesCluster* p = phoPesClusterU(emObject);
    if(p==NULL) return -99.0;
    return p->pesProfileRatio5by9();
  }
  double phoPesChi2ClusterV(const CdfEmObject& emObject){
    const PesCluster* p = phoPesClusterV(emObject);
    if(p==NULL) return -99.0;
    return p->pesProfileRatio5by9();
  }

  // Delta-R PES-PEM
  double phoPesDeltaR(const CdfEmObject& emObject) {
    if(phoDetector(emObject)!=PEM) return 99.0;
    if(!phoHasShowerMax(emObject)) return 99.0;
    double pesx = phoPes2dX(emObject);
    double pesy = phoPes2dY(emObject);
    double pemeta = emObject.getEmCluster()->pem3x3FitDetEta();
    double pemphi = emObject.getEmCluster()->pem3x3FitPhi();
    double cott = (exp(pemeta) - exp(-pemeta))/2.0;
    double z = 0.5 * (TLZPESL0 + TLZPESL1);
    if(emObject.getEmCluster()->side()==0) z = -z;
    double pemr = z/cott;
    double pemx = pemr*cos(pemphi);
    double pemy = pemr*sin(pemphi);
    double rsq = pow(pesx-pemx,2)+pow(pesy-pemy,2);
    if(rsq<1.e-10) {
      return 99.0;
    } else {
      return sqrt(rsq);
    }
  }

  double phoPesErrorX(const CdfEmObject& emObject) {

    double errorX = -999.0;

    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr != NULL) {
      errorX = pesptr->plug2dClusterErrX();
    }

    return errorX;
  }

  double phoPesErrorY(const CdfEmObject& emObject) {

    double errorY = -999.0;

    const Pes2dCluster* pesptr = phoPes2dCluster(emObject,1);

    if (pesptr != NULL) {
      errorY = pesptr->plug2dClusterErrY();
    }

    return errorY;
  }

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

    eUMax2 = -99.0;
    eVMax2 = -99.0;

    const PesCluster* u1p = phoPesClusterU(emObject);
    const PesCluster* v1p = phoPesClusterV(emObject);
    if(u1p==NULL || v1p==NULL) return -99.0;

    std::vector<CdfEmObject::Pes_link>::const_iterator ci,first,last;
    first = emObject.matchingPes2dClusters().contents().begin();
    last = emObject.matchingPes2dClusters().contents().end();
    if (first==last) return -99.0;

    // loop over 2-d clusters, exclude the photon's 1d clusters look
    // for next highest E
    for (ci=first; ci != last; ci++){
      const PesCluster* uTestp = &((*ci)->pesClusterU());
      const PesCluster* vTestp = &((*ci)->pesClusterV());
      if( uTestp != u1p && uTestp->energy()>eUMax2 ) 
	eUMax2 = uTestp->energy();
      if( vTestp != v1p && vTestp->energy()>eVMax2 ) 
	eVMax2 = vTestp->energy();
    }
    return (eUMax2> eVMax2 ? eUMax2 : 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) {
    ePprTowMax = -99.0;

    if(phoDetector(emObject) != PEM) return -99.0;
    const EmCluster& clust = *(emObject.getEmCluster());
    const EmCluster::TowerIdList&  tlist = clust.towerIds();

    // get access to the tower energies
    CalData_ch data; 
    CalData::Error result = CalData::find(data); 
    if(result != CalData::OK) return -99.0; 

    // loop over towers, examine PPR energies
    double emax = -99.0;
    double etot = 0.0;
    for(EmCluster::TowerIdList::const_iterator it=tlist.begin(); 
	                                      it<tlist.end(); ++it) { 
      TowerKey tk = clust.towerKey(*it);
      const CalTower* tower = data->tower(tk.iEta(),tk.iPhi()); 
      if(tower->includesDetector(PPR)) {
	float Etow = tower->pprEnergy();
	etot += Etow;
	if(Etow>emax) emax = Etow;

      }  // end if has PPR

    } // end loop over towers in cluster

    ePprTowMax = emax;
    return etot;
  }



} // namespace PhotonVariables