////////////////////////////////////////////////////////////////////////////////
//
//     File: ElectronVariables.cc
//     Authors: David Goldstein, Martin Von Der Mey
//     Description: High level analysis electron variables defined as functions
//     Created: 10/24/2000
//
//     Revision: 0.1
//
////////////////////////////////////////////////////////////////////////////////

#include "HighLevelObjects/ElectronVariables.hh"
#include "CalorGeometry/CalConstants.hh"
#include "CLHEP/Vector/LorentzVector.h"

//---------------
// C++ Headers --
//---------------

#include <iostream>
#include <vector>
#include <string>
#include <cstdio>
#include <math.h>
#include <assert.h>
#include <iomanip>
#include "cern_i/packlib.h"

//-------------------------------
// Collaborating Class Headers --
//-------------------------------

#include "BaBar/Cdf.hh"
#include "AbsEnv/AbsEnv.hh"
#include "Edm/EventRecord.hh"
#include "Edm/ConstHandle.hh"
#include "Edm/Handle.hh"
#include "RawDataBanks/PESD_StorableBank.hh"
#include "RawDataBanks/CESD_StorableBank.hh"
#include "Trybos/TRY_Bank_Number.hh"
#include "AbsEnv/AbsEnv.hh"
#include "CLHEP/Units/PhysicalConstants.h"
// #include "ElectronObjects/CdfEmObject.hh"
// #include "ElectronObjects/CdfEmObjectColl.hh"
#include "CalorObjects/Pes2dCluster.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "MetObjects/CdfMet.hh"
#include "TrackingObjects/Tracks/track_cut.hh"
#include "Electron/CemCorrAlg.hh"
#include "Electron/PemCorrAlg.hh"
#include "Electron/emobj_alg.hh"

#ifdef __GNUG__
#pragma implementation
#endif

#include "Calor/Cal2.hh"


/////////////////////////////////////////////////////////////////////////////
// Below are electron accessor functions.  They can be called by including
// "HighLevelObjects/ElectronVariables.hh" in your code and then typing 
// ElectronVariables::eleOne_of_the_functions( input_parameter, etc... )
/////////////////////////////////////////////////////////////////////////////

namespace ElectronVariables {


// Returns detector region (CEM=0, PEM=1 currently) 

Detector eleDetector(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  Detector region = cluster->region();
  return region;
}

// Returns detector eta...
  // current algorithm: in CEM, if emObject has track, use CES clusters 
  // closest to track, else use highest energy clusters
  // In PEM, use closest PES cluster to PEM centroid
  // Failing everything else, use cluster centroid
  double eleEtaDetector(const CdfEmObject& em_object) {
    const EmCluster* cluster = em_object.getEmCluster();
    if (cluster->region() == CEM) {
      CdfEmObject::Ces_link cesstr = em_object.bestMatchingCesCluster(CesCluster::STRIP,
								      "DEFAULT");
      CdfEmObject::Ces_link ceswir = em_object.bestMatchingCesCluster(CesCluster::WIRE,
								      "DEFAULT");
      if (cesstr.is_nonnull() && ceswir.is_nonnull()) {
	double localx = ceswir->localCoord();
	double y = sqrt(TLYCES*TLYCES + localx*localx);
	double z = cesstr->side() ? cesstr->localCoord() : -cesstr->localCoord();
	return HepPoint3D(0,y,z).eta();
      }
    } else { // PEM
      CdfEmObject::Pes_link pes = em_object.bestMatchingPes2dCluster();
      if (pes.is_nonnull()) {
	return pes->plug2dClusterLocation().eta();
      }
    }
  return cluster->emEtaDetector();
  }
  

// Returns event eta...

double eleEtaEvent(const CdfEmObject& em_object) {
  HepLorentzVector mom = eleFourMomentum(em_object);
  if (mom.e() > 0) {
    return mom.pseudoRapidity();
  }
  const EmCluster* cluster = em_object.getEmCluster();
  double event_eta = cluster->emEtaEvent();
  return event_eta;
}

// Returns phi...

double elePhi(const CdfEmObject& em_object) {
  bool phiFound = false;
  double phi;
  HepLorentzVector mom = eleFourMomentum(em_object);
  if (mom.e() > 0) {
    phi = mom.phi();
    phiFound = true;
  }
  if (! phiFound) {
    phi = eleStripPhi(em_object);
    if (phi != -99999.) {
      phiFound = true;
    }
  }
  if (! phiFound) {
    phi = em_object.matchingEmCluster()->emPhi();
  }
  phi = (phi < 0) ? phi + TWOPI : phi;
  return phi;
}

  double eleStripPhi(const CdfEmObject& em_object) {
    double phi;
    bool phiFound = false;
    if (eleDetector(em_object) == CEM) {
      CdfEmObject::Ces_link cesstr = em_object.bestMatchingCesCluster(CesCluster::STRIP,
								      "DEFAULT");
      CdfEmObject::Ces_link ceswir = em_object.bestMatchingCesCluster(CesCluster::WIRE,
								      "DEFAULT");
      if (cesstr.is_nonnull() && ceswir.is_nonnull()) {
	phi = ceswir->globalCoord();
	phiFound = true;
      }
    } else {
      CdfEmObject::Pes_link pes = em_object.bestMatchingPes2dCluster();
      if (pes.is_nonnull()) {
	phi = pes->plug2dClusterLocation().phi();
	phiFound = true;
      }
    }
    if (phiFound) {
      phi = (phi < 0) ? phi + TWOPI : phi;
      return phi;
    } else {
      return -99999.;
    }
  }

// Returns z-vertex...

double eleZVertex(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster(); 
  double z_vertex = cluster->vertexZ();
  return z_vertex;
}

// Returns cluster centroid Z-coordinate...

double eleZCentroid(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double z_centroid = cluster->zCentroid();
  return z_centroid;
}

// Returns raw energy of cluster in EM calorimeter...

double eleRawEmEnergy(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double raw_energy = cluster->emEnergy();
  return raw_energy;
}

// Returns EM calorimeter Et...

double eleEmEt(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double transverse_energy = cluster->emEt();
  return transverse_energy;
}

// Returns corrected Et...

double eleCorrectedEt(const CdfEmObject& em_object) {
  HepLorentzVector mom = eleFourMomentum(em_object);
  if (mom.e() > 0) {
    return mom.et();
  }
  const EmCluster* cluster = em_object.getEmCluster();
  double corrected_et = cluster->emEt();

  if(em_object.getEmCluster()->region()==CEM) {
    HepDouble cemcorr = 1.0;
    if (CemCorrAlg::GetCemCorr( em_object, "USETRACK", cemcorr) 
	== CemCorrAlg::OK) {
      corrected_et = corrected_et/cemcorr;
    }
  } else if(em_object.getEmCluster()->region()==PEM) {
    float pemcorr = 1.0;
    if (PemCorrAlg::GetPemCorr( em_object, pemcorr) == PemCorrAlg::OK) 
      corrected_et /= pemcorr;
  }

  return corrected_et;
}

// Returns raw energy of cluster in Had calorimeter...

double eleRawHadEnergy(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double raw_energy = cluster->hadEnergy();
  return raw_energy;
}

// Returns transverse Had energy...

double eleHadEt(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double transverse_energy = cluster->hadEt();
  return transverse_energy;
}

// Returns total (raw) energy (EM cal. + HAD cal.)...

double eleTotalRawEnergy(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double raw_energy = cluster->totalEnergy();
  return raw_energy;
}

// Returns total (transverse) energy (EM Et + HAD Et)...

double eleTotalEt(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double transverse_energy = cluster->totalEt();
  return transverse_energy;
}

// Returns (hadronic energy / electromagnetic energy)...
// ...method called uses total em energy & total 
// had energy as of 11/01/00...

double eleHadEm(const CdfEmObject& em_object) {
  const EmCluster* cluster = em_object.getEmCluster();
  double had_over_em = cluster->hadEm(); 
  return had_over_em;
}

// Returns total isolation energy (had + em) in cone of .4

double eleTotalIso4(const CdfEmObject& em_object) {  
  const EmCluster* cluster = em_object.getEmCluster();
  double iso_energy = cluster->totalIsolationEt4();
  return iso_energy;
}

// Returns total isolation energy (had + em) in cone of .7

double eleTotalIso7( const CdfEmObject& em_object ) {  
  const EmCluster* cluster = em_object.getEmCluster();
  double iso_energy = cluster->totalIsolationEt7();
  return iso_energy;
}

// Returns Em isolation energy in cone of .4

double eleEmIso4(const CdfEmObject& em_object) {  
  const EmCluster* cluster = em_object.getEmCluster();
  double iso_energy = cluster->emIsolationEt4();
  return iso_energy;
}

// Returns Em isolation energy in cone of .7

double eleEmIso7(const CdfEmObject& em_object) {  
  const EmCluster* cluster = em_object.getEmCluster();
  double iso_energy = cluster->emIsolationEt7();
  return iso_energy;
}

// Returns Had isolation energy in cone of .4

double eleHadIso4(const CdfEmObject& em_object) {  
  const EmCluster* cluster = em_object.getEmCluster();
  double iso_energy = cluster->hadIsolationEt4();
  return iso_energy;
}

// Returns Had isolation energy in cone of .7

double eleHadIso7(const CdfEmObject& em_object) {  
  const EmCluster* cluster = em_object.getEmCluster();
  double iso_energy = cluster->hadIsolationEt7();
  return iso_energy;
}

// Returns Lshr...

double eleLshr(const CdfEmObject& em_object) {
  double lshr = em_object.lshr();
  return lshr;
}

// Returns Lshr2 (two-tower Lshr?)...

double eleLshr2(const CdfEmObject& em_object) {
  double lshr2 = em_object.lshr(2);
  return lshr2;
}

// Returns pt of the track associated with the electron.  This is
// currently (11/01/00) the max pt track of those matching the em
// cluster.  -check CdfEmObject.hh/.cc to make sure...
// If there is no track for the CdfEmObject in question, -9999.0
// will be returned.

double eleTrackPt(const CdfEmObject& em_object) {
  double no_track = -9999.0;
  CdfEmObject::Track_link  seed_track = em_object.maxPtTrack();
  if (seed_track) {
  double track_pt = seed_track->pt();
  return track_pt;
  }
  else {
  return no_track;
  }
}

  // Returns 4-momentum (HepLorentzVector) with direction taken from track
  // and energy taken from em calorimeter.  CEM corrections are used if
  // applicable.
  // If there is no track for the CdfEmObject in question, a
  // HepLorentzVector with all entries = -9999.0 will be returned.
  
HepLorentzVector eleFourMomentum(const CdfEmObject& em_object) {
  const HepLorentzVector no_track = HepLorentzVector(-9999.0, -9999.0, -9999.0, -9999.0);
  HepLorentzVector twotow = emobj_alg::CEMFourMomentum::TwoTower ( em_object, false );
  CdfEmObject::Track_link  seed_track = em_object.maxPtTrack();
  if (seed_track) {
    HepDouble cemcorr = 1.0;
    bool corvalid = false;

    if (em_object.getEmCluster()->region() == CEM) {
      CemCorrAlg::Error err = CemCorrAlg::GetCemCorr(em_object, "USETRACK", cemcorr);
      if (err == CemCorrAlg::OK) {
        corvalid = true;
      } else {
        errlog.setSubroutine("ElectronVariables::eleFourMomentum");
        ERRLOG(ELwarning, "USETRACK correction failed")
          << "USETRACK correction failed, using uncorrected energy for electron four-vector"
          << endmsg;
        corvalid = false;
      }
    }
    return corvalid ? twotow/cemcorr : twotow;
  } else {
    return no_track;
  }
}

// Returns the charge...
// If there is no track for the CdfEmObject in question, -9999.0
// will be returned.

double eleCharge(const CdfEmObject& em_object) {
  double no_track = -9999.0;
  CdfEmObject::Track_link  seed_track = em_object.maxPtTrack();
  if (seed_track) {
  double charge = seed_track->charge();
  return charge;
  }
  else {
  return no_track;
  }
}

// Returns (Energy of cluster in EM cal.)/(magnitude of track mom.)...
// If there is no track for the CdfEmObject in question, -9999.0
// will be returned.

double eleEOverP(const CdfEmObject& em_object) {
  double no_track = -9999.0;
  CdfEmObject::Track_link  seed_track = em_object.maxPtTrack();
  if (seed_track) {
  const EmCluster* cluster = em_object.getEmCluster();
  double em_energy = cluster->emEnergy();
  const Hep3Vector p = seed_track->momentum();
  const HepDouble mag_of_p = p.Hep3Vector::mag();
  double e_over_p = em_energy/mag_of_p ;
  return e_over_p ;
  }
  else {
  return no_track;
  }
}

// Returns xces (wire info.)...
// If wire cluster is null, returns -8888.0

double eleXces(const CdfEmObject& em_object) {
  double null_wire_cluster = -8888.0;
  ValueVectorLink<CesClusterColl> wire_cluster;
  wire_cluster = em_object.bestMatchingCesCluster(CesCluster::WIRE);
  if (wire_cluster.is_nonnull()) {
    double xces = wire_cluster->fitted_position();
    return xces;
  }
  else {
    return null_wire_cluster;
  }
}

// Returns chi-squared (with emobj_alg correction) for wire cluster...
// If wire cluster is null, returns -8888.0

double eleChiW(const CdfEmObject& em_object) {
  double null_wire_cluster = -8888.0;
  ValueVectorLink<CesClusterColl> wire_cluster;
  wire_cluster = em_object.bestMatchingCesCluster(CesCluster::WIRE);
  if (wire_cluster.is_nonnull()) {
    return scaleCesChisq(wire_cluster->chi2(), eleRawEmEnergy(em_object));
  }
  else {
  return null_wire_cluster;
  }
}

// Returns delta-x...
// If wire cluster is null, returns -8888.0
// If there is no track for the CdfEmObject in question,
// -9999.0 will be returned.

double eleDeltaX(const CdfEmObject& em_object) {
  double null_wire_cluster = -8888.0;
  double no_track = -9999.0;
  ValueVectorLink<CesClusterColl> wire_cluster;
  wire_cluster = em_object.bestMatchingCesCluster(CesCluster::WIRE);
  CdfEmObject::Track_link  seed_track = em_object.maxPtTrack();
  if (wire_cluster.is_null()) {
    return null_wire_cluster;
  }
  else {
    if (seed_track) {
      return emobj_alg::TrackLocalCoord::xdiff(em_object);
      /*      double track_x = em_object.maxPtTrackLocalCoord().x();
      double wire_coord = wire_cluster->fitted_position();
      double delta_x = track_x - wire_coord;
      return delta_x; */
    }
    else {
    return no_track;
    }
  }
}

// Returns zces (strip info.)...
// If strip cluster is null, returns -8888.0

double eleZces(const CdfEmObject& em_object) {
  double null_strip_cluster = -8888.0;
  ValueVectorLink<CesClusterColl> strip_cluster;
  strip_cluster = em_object.bestMatchingCesCluster(CesCluster::STRIP);
  if (strip_cluster.is_nonnull()) {
    double zces = strip_cluster->fitted_position();
    return zces;
  }
  else {
    return null_strip_cluster;
  }
}

// Returns chi-squared (with emobj_alg correction) for strip cluster...
// If strip cluster is null, returns -8888.0

double eleChiS(const CdfEmObject& em_object) {
  double null_strip_cluster = -8888.0;
  ValueVectorLink<CesClusterColl> strip_cluster;
  strip_cluster = em_object.bestMatchingCesCluster(CesCluster::STRIP);
  if (strip_cluster.is_nonnull()) {
    return scaleCesChisq(strip_cluster->chi2(), eleRawEmEnergy(em_object));
  }
  else {
  return null_strip_cluster;
  }
}

// Returns delta-z...
// If strip cluster is null, returns -8888.0
// If there is no track for the CdfEmObject in question,
// -9999.0 will be returned.

double eleDeltaZ(const CdfEmObject& em_object) {
  double null_strip_cluster = -8888.0;
  double no_track = -9999.0;
  ValueVectorLink<CesClusterColl> strip_cluster;
  strip_cluster = em_object.bestMatchingCesCluster(CesCluster::STRIP);
  CdfEmObject::Track_link  seed_track = em_object.maxPtTrack();
  if (strip_cluster.is_null()) {
    return null_strip_cluster;
  }
  else {
    if (seed_track) {
      return emobj_alg::TrackLocalCoord::zdiff(em_object);
      /*      double track_z = em_object.maxPtTrackLocalCoord().z();
      double wire_coord = strip_cluster->fitted_position();
      if(strip_cluster->side()==0)  wire_coord = -wire_coord;
      double delta_z = track_z - wire_coord;
      return delta_z; */
    }
    else {
    return no_track;
    }
  }
}

//########################################################################
// Fiducial routine for electrons.
// Returns an integer according to
// the following code:             -1 = error (null strip/wire cluster)
//                                  0 = not fiducial
//                                  1 = fiducial
//                                  2 = fiducial (tower 9)

// int eleFid(const CdfEmObject& em_object) {

//    int error = -1;
//    int status;
//    float xces;
//    float zces;
//    float plug_eta;

//    const EmCluster* cluster = em_object.getEmCluster();
//    float phi = cluster->emPhi();

//    // first check if EM object is central or plug
//    int region = cluster->region();

//    if (region == 0) { // if electron is central.

//      ValueVectorLink<CesClusterColl> strip_cluster;
//      strip_cluster = em_object.bestMatchingCesCluster(CesCluster::STRIP);
//      ValueVectorLink<CesClusterColl> wire_cluster;
//      wire_cluster = em_object.bestMatchingCesCluster(CesCluster::WIRE);

//      // check for null Ces info.
//      if (strip_cluster.is_nonnull()) {
//        if (strip_cluster->side()) {   // get zces info.
//          zces = strip_cluster->fitted_position();
//        }
//        else {
//          zces = -strip_cluster->fitted_position();
//        }
//      }
//      else { 
//        return error; }

//      if (wire_cluster.is_nonnull()) {
//        xces = wire_cluster->fitted_position();
//      }
//      else { 
//        return error; }

//      // now check and assign fiducial status
//      if ( (abs(xces)>21) || (abs(zces)<9) || (abs(zces)>230) || 
//        // chimney condition: 0.77 < eta < 1.0 and 75 < phi < 90.
//        // taking radius of CES strips at 184.15 cm.
//        ( (phi>1.31) && (phi<1.57) && (zces>156.23) && (zces<216.41) ) ) {
//        status = 0;  // set to return non-fiducial
//      }
//      else {
//        if (abs(zces)<217) {
//          status = 1;  // set to return fiducial
//        } 
//        else {
//          status = 2;  // set to return fiducial (tower9)
//        }
//      }
//   }

//   if (region == 1) { // if electron is plug.

//     CdfEmObject::Pes_link pes_cluster; //ValueVectorLink<PesClusterColl>
//     pes_cluster = em_object.bestMatchingPes2dCluster();

//     if (pes_cluster.is_nonnull()) {
//       plug_eta = pes_cluster->plug2dClusterEta();
//     }
//     else { 
//       return error; }

//     if ((plug_eta>1.2) && (plug_eta<2.6)) {
//       status = 1;  // set to return fiducial
//     }
//     else {
//       status =2;  // set to return non-fiducial
//     }
//   }
//  return status;
//}
//########################################################################


// likelihood that object is really a single EM cluster
// Stephen Levy: October 6, 2003
// Provide argument for likelihood version.
// See Electron/Electron/emobj_alg.hh for details.

double eleLikelihood(const CdfEmObject& em_object, LikelihoodVersion ver) {
  if(em_object.getEmCluster()->region()==CEM) {
    return centralEmIdLikelihood(em_object,ver);
  } else {
    return plugEmIdLikelihood(&em_object);
  }
}


//*********************************************************************
// Tracking isolation.  This takes five input parameters:  Limits are
// currently set as follows:
//  
//               (1) must be a POINTER to a CdfTrack. -Taken as seed track.
//               (2) must be a REFERENCE to a CdfTrackView.   
//               (3)  0.0 <= cone_size <= 5.0  [eta-phi space]
//               (4)  0.0 <= dz_max <= 20.0    [cm]
//               (5)  0.0 <= pt_min <= 100.0   [GeV]
//
// Returns -9.0 if an input parameter is unacceptable.  Otherwise returns
// sum of pt of tracks within specified cone (excluding pt of seed track).
//
// -If using your own code and are unclear about selecting a CdfTrackView, 
// here is an example which will work, provided you have already defined
// a pointer to a specific CdfEmObject -- type:
//
// const CdfTrackView& yourTracks = yourpointertoemobject->matchingTracks();
// 
// -you could then pass `yourtracks' as the 2nd input parameter...

double eleTiso(const CdfTrack* seed_track,     
               const CdfTrackView& track_list, 
	       double cone_size,                
               double dz_max,                   
               double pt_min ) {

  double error = -9.0;

  // test input parameters and return error code if any are unacceptable...

   if (cone_size <= 0.0 || cone_size >= 5.0 || dz_max <= 0.0 ||
       dz_max >= 20.0 || pt_min <= 0.0 || pt_min >= 100.0) {

     std::cout << "One or more input parameters are out of bounds." << std::endl;
     std::cout << "See comments in ElectronVariables(.cc or .hh) for" << std::endl;
     std::cout << "limits on input parameters." << std::endl;

     return error;
   }    

  else {
    double tisosumpt = 0.0;
    double seed_zvtx = seed_track->z0();
    double seed_eta  = log(tan((atan(seed_track->lambda())+M_PI/2.0)/2.0));
    double seed_phi  = seed_track->phi0();
    // print out cone_size
    //    std::cout << "using cone size = " << cone_size << std::endl;
    // print out dz_max
    //std::cout << "setting max. delta z = " << dz_max << std::endl;
    // print out pt_min
    //std::cout << "setting min. pt value = " << pt_min << std::endl;
  //     print out seed_zvtx
    //std::cout << "seed track z-vertex = " << seed_zvtx << std::endl;
    //     print out seed_eta
    //std::cout << "seed track eta = " << seed_eta << std::endl;
   //    print out seed_phi
    // std::cout << "seed track phi = " << seed_phi << std::endl;
  
    for ( CdfTrackView::const_iterator i = track_list.contents().begin();
	i != track_list.contents().end() ; ++i ) {
    // print out z0 and pt for ith track
      //std::cout << "track z0 = " << (*i)->z0() << std::endl;
      //std::cout << "track pt = " << (*i)->pt() << std::endl;
  
      if ( (abs(((*i)->z0())-seed_zvtx)) <= dz_max &&
	 (abs((*i)->pt()) >= pt_min) ) {
      double eta_i = log(tan((atan((*i)->lambda())+M_PI/2.0)/2.0));
      double phi_i = (*i)->phi0();
    // print out eta_i and phi_i
      //std::cout << "eta_i = " << eta_i << std::endl;
      //std::cout << "phi_i = " << phi_i << std::endl;
  
        if ( (sqrt( (eta_i-seed_eta)*(eta_i-seed_eta) +
			 (phi_i-seed_phi)*(phi_i-seed_phi) ) ) <= cone_size &&
	   (sqrt( (eta_i-seed_eta)*(eta_i-seed_eta) +             
			 (phi_i-seed_phi)*(phi_i-seed_phi) ) ) > 0.0 ) {
	  tisosumpt += (*i)->pt();
        }

      }

    }
   return tisosumpt;
  }
}
//******************************************************************

//*********************************************************************
// ###THIS TISO IS OVERLOADED, TAKES SLIGHTLY DIFFERENT INPUT FORMS###
// Tracking isolation.  This takes five input parameters:  Limits are
// currently set as follows:
//  
//               (1) must be a CdfEmObject::Track_link
//               (2) must be a CdfTrackView_h.   
//               (3)  0.0 <= cone_size <= 5.0  [eta-phi space]
//               (4)  0.0 <= dz_max <= 20.0    [cm]
//               (5)  0.0 <= pt_min <= 100.0   [GeV]
//
// Returns -9.0 if an input parameter is unacceptable.  Otherwise returns
// sum of pt of tracks within specified cone (excluding pt of seed track).
//
// -If using your own code and are unclear about selecting a CdfTrackView, 
// here is an example which will work, provided you have already defined
// a pointer to a specific CdfEmObject -- type:
//
// const CdfTrackView& yourTracks = yourpointertoemobject->matchingTracks();
// 
// -you could then pass `yourtracks' as the 2nd input parameter...

double eleTiso(CdfEmObject::Track_link seed_track,     
               CdfTrackView_h track_list, 
	       double cone_size,                
               double dz_max,                   
               double pt_min ) {

  double error = -9.0;

  // test input parameters and return error code if any are unacceptable...

   if (cone_size <= 0.0 || cone_size >= 5.0 || dz_max <= 0.0 ||
       dz_max >= 20.0 || pt_min <= 0.0 || pt_min >= 100.0) {

     std::cout << "One or more input parameters are out of bounds." << std::endl;
     std::cout << "See comments in ElectronVariables(.cc or .hh) for" << std::endl;
     std::cout << "limits on input parameters." << std::endl;

     return error;
   }    

  else {
    double tisosumpt = 0.0;
    double seed_zvtx = seed_track->z0();
    double seed_eta  = log(tan((atan(seed_track->lambda())+M_PI/2.0)/2.0));
    double seed_phi  = seed_track->phi0();
    // print out cone_size
    //std::cout << "using cone size = " << cone_size << std::endl;
    // print out dz_max
    //std::cout << "setting max. delta z = " << dz_max << std::endl;
    // print out pt_min
    //std::cout << "setting min. pt value = " << pt_min << std::endl;
  //     print out seed_zvtx
    //std::cout << "seed track z-vertex = " << seed_zvtx << std::endl;
    //     print out seed_eta
    //std::cout << "seed track eta = " << seed_eta << std::endl;
   //    print out seed_phi
    // std::cout << "seed track phi = " << seed_phi << std::endl;
  
    for ( CdfTrackView::const_iterator i = track_list->contents().begin();
	i != track_list->contents().end() ; ++i ) {
    // print out z0 and pt for ith track
      //std::cout << "track z0 = " << (*i)->z0() << std::endl;
      //std::cout << "track pt = " << (*i)->pt() << std::endl;
  
      if ( (abs(((*i)->z0())-seed_zvtx)) <= dz_max &&
	 (abs((*i)->pt()) >= pt_min) ) {
      double eta_i = log(tan((atan((*i)->lambda())+M_PI/2.0)/2.0));
      double phi_i = (*i)->phi0();
    // print out eta_i and phi_i
      //std::cout << "eta_i = " << eta_i << std::endl;
      //std::cout << "phi_i = " << phi_i << std::endl;
  
        if ( (sqrt( (eta_i-seed_eta)*(eta_i-seed_eta) +
			 (phi_i-seed_phi)*(phi_i-seed_phi) ) ) <= cone_size &&
	   (sqrt( (eta_i-seed_eta)*(eta_i-seed_eta) +             
			 (phi_i-seed_phi)*(phi_i-seed_phi) ) ) > 0.0 ) {
	  tisosumpt += (*i)->pt();
        }

      }

    }
   return tisosumpt;
  }
}
//******************************************************************

////////////////////////////////////////////////////////////////
// this has been...namespace ElectronVariables...brought to you
// by the letter P and the number 9
////////////////////////////////////////////////////////////////
}