////////////////////////////////////////////////////////////////////
//
//     File: MuonVariables2.cc
//     Authors: Abe DeBenedetti and Henry Frisch
//     Description: High level analysis muon variables defined as functions
//     Created: 7/26/2001
//
//     Revision History: none so far
//
////////////////////////////////////////////////////////////////////

//#include <cstdio>
#include <math.h>
#include <cmath>
#include <algorithm>
//#include "Rtypes.h"
#include "HighLevelObjects/MuonVariables.hh"
#include "MuonObjects/CdfMuon.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "CLHEP/Vector/LorentzVector.h"
#include "Extrapolator/MExtrapolator.hh"
#include "VertexObjects/VertexColl.hh"
#include "AbsEnv/AbsEnv.hh"
#include "VertexObjects/Vertex.hh"
#include "HighLevelObjects/TrackingVariables.hh"
#include "CalorObjects/CalData.hh"
#include "MuonUser/MuonAnalysisTools.hh"

using namespace std;
namespace MuonVariables {

  // Pt taken from the best track
  double muoPt(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().perp();
  } 
  
  // Px taken from the best track
  double muoPx(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().px();
  } 
  
  // Py taken from the best track
  double muoPy(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().py();
  } 

  // Pz taken from the best track
  double muoPz(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().pz();
  } 

  // E taken from the best track
  double muoE(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().e();
  } 

  // Charge (+=1; -=-1): inherited from SimpleReconstructedTrack 
  double muoCharge(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->charge();
  } 
  
  // This is event eta, because it comes from the track
  double muoEventEta(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().eta();
  } 

  // This is phi at r=0, because it comes from the track
  double muoPhi0(const CdfMuon& muonObject) {
    return muonObject.fourMomentum().phi();
  } 

  //Detector Eta
    // Needs input argument giving muon Type:
    // Type of muon taken from CdfMuon
    //typedef enum { CMU = 0,
    //               CMP = 1,
    //               CMX = 2,
    //               BMU = 3,
    //               CSP = 4,
    //               CSX = 5,
    //               BSU = 6,
    //               TSU = 7 } System;
  //
  //spacePoint() should be in Global Coordinates taken from the Stub 
  double muoDetectorEta(const CdfMuon& muonObject,int system) {
      return muonObject.stub(system)->spacePoint().eta();
  } 

  //Detector Phi
  //spacePoint() should be in Global Coordinates taken from the Stub 
  double muoDetectorPhi(const CdfMuon& muonObject,int system) {
    return muonObject.stub(system)->spacePoint().phi();
  } 

  //Stub X
  double muoDetectorX(const CdfMuon& muonObject,int system) {
      return muonObject.stub(system)->spacePoint().x();
  } 

  //Stub Y
  double muoDetectorY(const CdfMuon& muonObject,int system) {
      return muonObject.stub(system)->spacePoint().y();
  } 

  //Stub Z
  double muoDetectorZ(const CdfMuon& muonObject,int system) {
      return muonObject.stub(system)->spacePoint().z();
  } 


  //Type of muon taken from CdfMuon
  //typedef enum { CMU = 0,
  //               CMP = 1,
  //               CMX = 2,
  //               BMU = 3,
  //               CSP = 4,
  //               CSX = 5,
  //               BSU = 6,
  //               TSU = 7 } System;
  //bit 0 corresponds to CMU et cetra. So a CMU-CMP-CSP is 10011 or 19
  unsigned int muoType(const CdfMuon& muonObject) {
    return muonObject.muonType();
  } 

  //Does it has have a system (eg. CMU)?
  bool muoHasSystem(const CdfMuon& muonObject,int system) {
    return ((muoType(muonObject) >> system) & 1);
  }


  //Wedge or region of first hit (0-31)
  unsigned int muoWedge(const CdfMuon& muonObject,int system){
    const MuonStub * stub = muonObject.stub(system);
    if (!stub) return 0x1f;
    MuonStub::AssHitCIter FirstHit= stub->beginAssHits();
    if (!FirstHit->getHit()) return (unsigned int)-1;
    return  FirstHit->getHit()->getCode().getWedge();//identical to getRegion
  }

  //Bit map of CMU/CMP/CMX Layer hits
  unsigned int muoLayerBitmap(const CdfMuon& muonObject,int system){
   const MuonStub * stub = muonObject.stub(system);
   if (!stub) return 0;
   unsigned int map=0;
   for (MuonStub::AssHitCIter ihit= stub->beginAssHits();
	 ihit !=stub->endAssHits(); ihit++){
     int layer = ihit->getHit()->getCode().getLayer();


     if ((layer>=0)&&(layer<=8))    // for safety, expect <=4
	map = map |(1 << layer);
   }
   return map;
  }
  
  //Number of hits in the system for the stub
  int muoSystemHits(const CdfMuon& muonObject,int system){
    return muonObject.stub(system)->numAssHits();
  }

  //array of DigiCodes for up to 4 hits
  std::vector<MuonDigiCode> muoHitDigiCodes(const CdfMuon& muonObject,int system){
    const MuonStub * stub = muonObject.stub(system);
    std::vector<MuonDigiCode> codes(4);
    if (!stub) return codes;
    int i=0;
    for (MuonStub::AssHitCIter ihit= stub->beginAssHits(); 
         ihit !=stub->endAssHits(); ihit++){
      codes[i++]=ihit->getHit()->getCode();
    }
    return codes;
  }

  //array of Widths for up to 4 hits
  std::vector<unsigned int> muoHitWidths(const CdfMuon& muonObject,int system){
    const MuonStub * stub = muonObject.stub(system);
    std::vector<unsigned int> widths (4);
    if (!stub) return widths;
    int i=0;
    for (MuonStub::AssHitCIter ihit= stub->beginAssHits();
         ihit !=stub->endAssHits(); ihit++){
      widths[i++]=ihit->getHit()->getTDCWidth();
    }
    return widths;
  }

  //Tracking Algorithm for the best track
  int muoTrackAlgo(const CdfMuon& muonObject){
    return muonObject.bestTrack()->algorithm().value();
  }


  //X at the point of closest approach to beam 
  double muoTrackX0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->d0()*sin(muonObject.bestTrack()->phi0());
  } 

  //Y at the point of closest approach to beam
  double muoTrackY0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->d0()*cos(muonObject.bestTrack()->phi0());
  } 

  //Z at the point of closest approach to beam
  double muoTrackZ0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->z0();
  } 

  //According to Pasha, ) is positively curved and, i.e. dphi/dr>0, 
  //if y is positive then impact parameter is positive
  double muoTrackD0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->d0();
  } 

  //Phi at the point of closest approach to beam
  double muoTrackPhi0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->phi0();
  } 

  //Chi^2 of track
  double muoTrackChi2(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->chi2();
  } 

  //Degrees of Freedom from trackfit
  unsigned int muoTrackdof(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->getHelixFit().getNumDegreesOfFreedom();
  } 

  //Track Covariance Matrix
  const HepSymMatrix & muoTrackErrorMatrix(const CdfMuon& muonObject){
    return muonObject.bestTrack()->getHelixFit().getErrorMatrix();
  }

  //Extrapolation of the Track extrapolated to the radius of the Stub
  MExtrapolator muoTrackExtrapolation(const CdfMuon& muonObject,int system){

    const MuonStub * stub= muonObject.stub(system); //Muon Stub for given system

    MExtrapolator extr(CdfDetector::instance());    //Used to extrapolate the track
    HepPoint3D stubPoint = stub->spacePoint();      //Location of the Stub
                                                    //give the extrapolator to the momentum of the track 
    extr.setInitialMomentum(muonObject.bestTrack()->momentum());
                                                    //Get the Helix of the track
    Helix helix = muonObject.bestTrack()->getTrajectory();
                                                    //Give the extrapolator the position of helix at its 
                                                    //closest approch to the beam
    extr.setInitialPosition(helix.getPosition(helix.getPathLengthAtRhoEquals(muonObject.bestTrack()->d0())));
                                                    //Give the extrapolator the charge
    extr.setParticleCharge(int(muonObject.bestTrack()->charge()));
    if (muonObject.bestTrack()->charge() == -1) {   //Tell the Extrapolator what kind of particle it is
      extr.setParticleType(AbsExtrapolator::MUM);
    }
    else {
      extr.setParticleType(AbsExtrapolator::MUP);
    }
    
    extr.extrapolateToRadius(stubPoint.perp());     //Extrapolate to position of the system
    return  extr;
  }

  //3 point of the Track extrapolated to the radius of the Stub
  HepPoint3D muoExtrTrackPoint(const CdfMuon& muonObject,int system){
    return muoTrackExtrapolation(muonObject,system).getCurrentPosition();
  }


  // X position of the Track extrapolated to the radius of the Stub
  double muoExtrTrackX(const CdfMuon& muonObject,int system){
    return muoExtrTrackPoint(muonObject,system).x() ;
  }

  // Y position of the Track extrapolated to the radius of the Stub
  double muoExtrTrackY(const CdfMuon& muonObject,int system){
    return muoExtrTrackPoint(muonObject,system).y() ;
  }

  // Z position of the Track extrapolated to the radius of the Stub
  double muoExtrTrackZ(const CdfMuon& muonObject,int system){
    return muoExtrTrackPoint(muonObject,system).z() ;
  }

  //3 Momentum of the Track extrapolated to the radius of the Stub
  HepVector3D muoExtrTrackMomentum(const CdfMuon& muonObject,int system){
    return muoTrackExtrapolation(muonObject,system).getCurrentMomentum();
  }

 
//   //this is purely track based and has nothing to do with the stub
//   //We think this is useless but aren't sure
   
//   // the difference in theta between the position and Momentum
//   // of the Track extrapolated to the radius of the Stub
//   double muoExtrTrackDeltaTheta(const CdfMuon& muonObject,int system){
//     double theta = muoExtrTrackPoint(muonObject,system).theta();
//     double ptheta = muoExtrTrackMomentum(muonObject,system).theta();
//     return theta - ptheta;
//   }
  
//   // the difference in Phi between the position and Momentum
//   // of the Track extrapolated to the radius of the Stub
//   double muoExtrTrackDeltaPhi(const CdfMuon& muonObject,int system){
//     double phi = muoExtrTrackPoint(muonObject,system).phi();
//     double pphi = muoExtrTrackMomentum(muonObject,system).phi();
//     return phi - pphi;
//   }

  

// These are not implemented in MExtrapolator
//   //sigma x of the Track extrapolated to the radius of the Stub
//   HepVector3D muoExtrTrackSigmaX(const CdfMuon& muonObject,int system){
//     return muoTrackExtrapolation(muonObject,system).getSigmaXIntercept();
//   }

//   //sigma x slope of the Track extrapolated to the radius of the Stub
//   HepVector3D muoExtrTrackSigmaXSlope(const CdfMuon& muonObject,int system){
//     return muoTrackExtrapolation(muonObject,system).getSigmaXSlope();
//   }

//   //sigma z of the Track extrapolated to the radius of the Stub
//   HepVector3D muoExtrTrackSigmaZ(const CdfMuon& muonObject,int system){
//     return muoTrackExtrapolation(muonObject,system).getSigmaZIntercept();
//   }

//   //sigma z slope of the Track extrapolated to the radius of the Stub
//   HepVector3D muoExtrTrackSigmaZSlope(const CdfMuon& muonObject,int system){
//     return muoTrackExtrapolation(muonObject,system).getSigmaZSlope();
//   }

  //Bit map of the SVX hits
  unsigned int muoSiLayerBitMap(const CdfMuon& muonObject){
    return muonObject.bestTrack()->usedSI();
  }

  //Bit map of the COT hits
  CT_LayerMask muoCTLayerBitMap(const CdfMuon& muonObject){
   return muonObject.bestTrack()->usedCT();
  }

  //Number of Axial Superlayer Hits
  int muoNAxialHits(const CdfMuon& muonObject){
    return muonObject.bestTrack()->numCTHitsAx();
  }

  //Number of Stereo Superlayer Hits
  int muoNStereoHits(const CdfMuon& muonObject){
    return muonObject.bestTrack()->numCTHitsSt();
  }


  //Closest Vertex
  const Vertex *  muoClosestVertex(const CdfMuon& muonObject){
    AbsEvent*  pEvent = AbsEnv::theEvent();
    EventRecord::ConstIterator 
      iter(pEvent,"VertexColl","DEFAULT_VXPRIM_VERTICES");

    if(iter.is_valid()) {

      ConstHandle<VertexColl> vertexset;
      vertexset=iter;
      const VertexColl::CollType contents = vertexset->contents();
      
      double z = muoTrackZ0(muonObject);
      const Vertex * gvertex = NULL;
      
      double dzmin = 999.0,dz;
      for(VertexColl::const_iterator vertex = contents.begin();
	         vertex != contents.end();++vertex)
	{
	  dz = fabs((*vertex)->z() - z);
	  if ( dz < dzmin ){
	    dzmin = dz;
	    gvertex= &(**vertex);
	  }
	}

      return gvertex;

    } else {
      return NULL;
    }
  }
  
  //Z of Closest Vertex
  double muoClosestVertexZ(const CdfMuon& muonObject){
    const Vertex* v = muoClosestVertex(muonObject);
    return (v==NULL? -999.0 : v->z());
  }

  //Ht of Closest Vertex
  double muoClosestVertexHt(const CdfMuon& muonObject){
    const Vertex* v = muoClosestVertex(muonObject);
    return (v==NULL? -999.0 : v->ht());
  }

  //Number of tracks associated with the ClosestVertex
  int muoClosestVertexNTracks(const CdfMuon& muonObject){
    const Vertex* v = muoClosestVertex(muonObject);
    return (v==NULL? -999 : int(v->nTracks()));
  }
  
  //Delta Z between the Closest vertex and the Z0 of the track
  double muoDelZClosestVertex(const CdfMuon& muonObject){
    double zv = muoClosestVertexZ(muonObject);
    return (zv<-998.0? -999.0 : muoTrackZ0(muonObject) - zv );
  }


  //Delta Z between the primary vertex and the Z0 of the track
  double muoDelZPrimaryVertex(const CdfMuon& muonObject){
    const Vertex* pv = TrackingVariables::trkPrimaryVertex();
    return ( pv==NULL? -999.0 : muoTrackZ0(muonObject) - pv->z());
  }

  //Best Track
  CdfTrack_clnk muoBestTrack(const CdfMuon& muonObject) {
    return muonObject.bestTrack();
  } 


  //Second Best Track
  CdfTrack_clnk muoSecondBestTrack(const CdfMuon& muonObject) {
    return muonObject.secondBestTrack();
  } 

  const MuonStub * muoStub(const CdfMuon& muonObject,int system){
    return muonObject.stub(system);
  }

  //Return a pointer to the closest stub if there is one in the same wedge 
  const MuonStub * muoClosestStub(const CdfMuon& muonObject,int system){
    vector<const MuonStub*> stubs = muonObject.muonStubs();
    const MuonStub * ourstub= muonObject.stub(system);
    
    const MuonStub * closestub;
    
    for (vector<const MuonStub*>::const_iterator is = 
	   stubs.begin(); is != stubs.end(); ++is ) {
      if ((*is)->system() == ourstub->system()){
	int wedge=(*is)->beginAssHits()->getHit()->getCode().getWedge();
	                                              //In current releases there is a
	                                              //MuonStub::wedge() routine
	if (wedge == muoWedge(muonObject,system)){
	  if ((*is) != ourstub){
	    if (!closestub){
	      closestub=(*is);
	    }else if ( fabs( (*is)->spacePoint().phi() - ourstub->spacePoint().phi())<
		       fabs( closestub->spacePoint().phi() - ourstub->spacePoint().phi())){
	      closestub=(*is);
	    }
	  }
	}
      }
    }
	
    return closestub;
  }

  HepVector3D muoStubDirection(const CdfMuon& muonObject,int system){
    return muonObject.stub(system)->direction();
  }

  double muoStubDirectionPhi(const CdfMuon& muonObject,int system){
    return muonObject.stub(system)->direction().phi();
  }

  double muoStubDirectionCotanTheta(const CdfMuon& muonObject,int system){
    return 1/tan(muonObject.stub(system)->direction().theta());
  }

  // delta X between the Track extrapolated to the radius of the Stub and the stub
  double muoDeltaX(const CdfMuon& muonObject,int system){
    return muonObject.matchData(system).drphi(); //in later versions there is a function dx()
  }

  // delta Z between the Track extrapolated to the radius of the Stub and the stub
  double muoDeltaZ(const CdfMuon& muonObject,int system){
    return muonObject.matchData(system).dz();
  }
  
  // delta phi between the Track extrapolated to the radius of the Stub and the stub
  double muoDeltaPhi(const CdfMuon& muonObject,int system){
    return muonObject.matchData(system).dphi();
  }

  // delta Theta between the Track extrapolated to the radius of the Stub and the stub
  double muoDeltaTheta(const CdfMuon& muonObject,int system){
    return muonObject.matchData(system).dtheta();
  }

  // Chi Squared X of the Track extrapolated to the radius of the Stub and the stub
  double muoChiSqX(const CdfMuon& muonObject,int system){
    return muonObject.matchData(system).chsqX();
  }
  
  // Chi Squared Z of the Track extrapolated to the radius of the Stub and the stub
  double muoChiSqZ(const CdfMuon& muonObject,int system){
    return muonObject.matchData(system).chsqZ();
  }
  
  //X error matrix of the Track extrapolated to the radius of the Stub
  //This is used as the error matrix for the chi squared calculation
  HepSymMatrix muoErrorMatrixX(const CdfMuon& muonObject,int system){
    return  muoTrackExtrapolation(muonObject,system).getXErrorMatrix();
  }

  //Z error matrix of the Track extrapolated to the radius of the Stub
  //This is used as the error matrix for the chi squared calculation
  HepSymMatrix muoErrorMatrixZ(const CdfMuon& muonObject,int system){
    return  muoTrackExtrapolation(muonObject,system).getZErrorMatrix();
  }

  
  // Et in cone of R = 0.4 around muon tower
  double muoCone4Et(const CdfMuon& muonObject){
    return muonObject.coneR4Et();
  }

  // Et in cone of R = 0.7 around muon tower
  double muoCone7Et(const CdfMuon& muonObject){
    return muonObject.coneR7Et();
  }

  // Em fraction in cone of R = 0.4 around muon tower
  double muoCone4EmFraction(const CdfMuon& muonObject){
    return muonObject.coneR4EmFraction();
  }

  // Em fraction in cone of R = 0.7 around muon tower
  double muoCone7EmFraction(const CdfMuon& muonObject){
    return muonObject.coneR7EmFraction();
  }


  // Standard track Iso
  double muoTrackPtCone4(const CdfMuon& muonObject){
    CdfTrack_clnk  seed_track = muonObject.bestTrack(); 
    return TrackingVariables::trkTrackIsolation(*seed_track);
  }

  // This is the Sum of the Pt of all track within a cone of .7, pt>.2 and 
  // < 10 cm from the Z0 of the track
  double muoTrackPtCone7(const CdfMuon& muonObject){
    CdfTrack_clnk  seed_track = muonObject.bestTrack(); 
    return TrackingVariables::trkTrackIsolation(*seed_track,0.7);
  }

  // Standard track Iso, for given track view
  double muoTrackPtCone4(const CdfMuon& muonObject, const CdfTrackView& trview){
    CdfTrack_clnk  seed_track = muonObject.bestTrack();
    return TrackingVariables::trkTrackIsolation(*seed_track,trview);
  }

  // This is the Sum of the Pt of all track within a cone of .7, pt>.2 and 
  // < 10 cm from the Z0 of the track, for given track view
  double muoTrackPtCone7(const CdfMuon& muonObject, const CdfTrackView& trview){
    CdfTrack_clnk  seed_track = muonObject.bestTrack(); 
    return TrackingVariables::trkTrackIsolation(*seed_track,trview,0.7);
  }

  
  // em energy
  float muoEmEnergy(const CdfMuon& muonObject){
    return muonObject.emEnergy();
  } 


  // hadron Energy
  float muoHadronEnergy(const CdfMuon& muonObject){
    return muonObject.hadronEnergy();
  } 

  // em Et
  float muoEmEt(const CdfMuon& muonObject){
    return muonObject.emEnergy()*sin(muonObject.bestTrack()->theta());
  } 


  // hadron Et
  float muoHadronEt(const CdfMuon& muonObject){
    return muonObject.hadronEnergy()*sin(muonObject.bestTrack()->theta());
  } 

  //Et of Neighboring towers
  float muoNeighborEt(const CdfMuon& muonObject){
    return muonObject.neighborEt();
  }


  float muoTime(const CdfMuon& muonObject) {
    CalData_ch cal_hndl;
    CalData::Error calStatus = CalData::find(cal_hndl);
    if (calStatus == CalData::OK){
      float maxHad = 0.5; //Set threshold HadE for central
      float tdcTime = 999999;
      float tmptdcTime = 999999;
      for (CdfMuon::TowerKeySet::const_iterator itow = muonObject.hadTowers().begin();
	   itow != muonObject.hadTowers().end(); ++itow) {
	if (itow->isValid()){
	  const CalTower* myTwr = cal_hndl->tower(itow->iEta(),itow->iPhi());
	  if (myTwr && myTwr->hadEnergy() > maxHad) {
	    if(myTwr->includesDetector(CHA) && myTwr->includesDetector(WHA)){
	      //take TDC value closest to zero
	      CellKey keyCHA(myTwr->iEta(),myTwr->iPhi(),CHA);
              CellKey keyWHA(myTwr->iEta(),myTwr->iPhi(),WHA);
              //dead or hot TDC
	      if(myTwr->nTDCHits(keyCHA)==0 || myTwr->nTDCHits(keyCHA)>2) {
		if(myTwr->nTDCHits(keyWHA)==1)
		  tdcTime = myTwr->time(keyWHA,0);
		else if(myTwr->nTDCHits(keyWHA)==2)
		  tdcTime = (abs(myTwr->time(keyWHA,0)) <
			     abs(myTwr->time(keyWHA,1))) ?
		    myTwr->time(keyWHA,0) : tmptdcTime = myTwr->time(keyWHA,1);             
	      }
	      //dead or hot TDC
	      else if(myTwr->nTDCHits(keyWHA)==0 || myTwr->nTDCHits(keyWHA)>2) {
		if(myTwr->nTDCHits(keyCHA)==1)
		  tdcTime = myTwr->time(keyCHA,0);
		else if(myTwr->nTDCHits(keyCHA)==2)
		  tdcTime = (abs(myTwr->time(keyCHA,0)) <
			     abs(myTwr->time(keyCHA,1))) ?
		    myTwr->time(keyCHA,0) : tmptdcTime = myTwr->time(keyCHA,1);
	      }
	      //take TDC of Tower with larger E
	      else if(myTwr->chaEnergy()>myTwr->whaEnergy()) {
		if(myTwr->nTDCHits(keyCHA)==1)
		  tdcTime = myTwr->time(keyCHA,0);
		else if(myTwr->nTDCHits(keyCHA)==2)
		  tdcTime = (abs(myTwr->time(keyCHA,0)) <
			     abs(myTwr->time(keyCHA,1))) ?
		    myTwr->time(keyCHA,0) : tmptdcTime = myTwr->time(keyCHA,1);
	      }
	      else {
		if(myTwr->nTDCHits(keyWHA)==1)
		  tdcTime = myTwr->time(keyWHA,0);
		else if(myTwr->nTDCHits(keyWHA)==2)
		  tdcTime = (abs(myTwr->time(keyWHA,0)) <
			     abs(myTwr->time(keyWHA,1))) ?
		    myTwr->time(keyWHA,0) : tmptdcTime = myTwr->time(keyWHA,1);
	      }
	    }
	    else if(myTwr->includesDetector(CHA)){
	      CellKey keyCHA(myTwr->iEta(),myTwr->iPhi(),CHA);
	      if(myTwr->nTDCHits(keyCHA)==1)
		tdcTime = myTwr->time(keyCHA,0);
	      else if(myTwr->nTDCHits(keyCHA)==2)
		tdcTime = (abs(myTwr->time(keyCHA,0)) <
			   abs(myTwr->time(keyCHA,1))) ?
		  myTwr->time(keyCHA,0) : tmptdcTime = myTwr->time(keyCHA,1);
	    }
	    else if(myTwr->includesDetector(WHA)){
	      CellKey keyWHA(myTwr->iEta(),myTwr->iPhi(),WHA);
	      if(myTwr->nTDCHits(keyWHA)==1)
		tdcTime = myTwr->time(keyWHA,0);
	      else if(myTwr->nTDCHits(keyWHA)==2)
		tdcTime = (abs(myTwr->time(keyWHA,0)) <
			   abs(myTwr->time(keyWHA,1))) ?
		  myTwr->time(keyWHA,0) : tmptdcTime = myTwr->time(keyWHA,1);
	    }
	    else if(myTwr->includesDetector(PHA) && myTwr->hadEnergy() > 3){
	      CellKey keyPHA(myTwr->iEta(),myTwr->iPhi(),PHA);
	      if(myTwr->nTDCHits(keyPHA)==1)
		tdcTime = myTwr->time(keyPHA,0);
	      else if(myTwr->nTDCHits(keyPHA)==2)
		tdcTime = (abs(myTwr->time(keyPHA,0)) <
			   abs(myTwr->time(keyPHA,1))) ?
		  myTwr->time(keyPHA,0) : tmptdcTime = myTwr->time(keyPHA,1);
	    }
	    if(tmptdcTime != 999999){ //make sure TDC info is valid
	      tdcTime = tmptdcTime;
	      maxHad = myTwr->hadEnergy();
	    }
	  }
	}
      }
      return tdcTime;
    } else {
      return 999999;
    }
  }

  //Returns a bitmask of which detector the MuonFiducialTool claims 
  // a muon should hit. See description of muoType for explanation of
  // bitmask.
  int muoFidType(const CdfMuon& muonObject) {
    MuonFiducialTool mutool;
    int fidtype = 0;
    double dist, xdist, zdist;
    if (mutool.isFiducial(muonObject, MuonDigiCode::CMU, dist, xdist, zdist)) {
      fidtype |= 0x1;
    }
    if (mutool.isFiducial(muonObject, MuonDigiCode::CMP, dist, xdist, zdist)) {
      if ((xdist < 0) && (zdist < -3)) {
	fidtype |= 0x2;
      }
    }
    if (mutool.isFiducial(muonObject, MuonDigiCode::CMX, dist, xdist, zdist)) {
      if ((xdist < 0) && (zdist < -3)) {
	fidtype |= 0x4;
      }
    }
    if (mutool.isFiducial(muonObject, MuonDigiCode::BMU, dist, xdist, zdist)) {
      if ((xdist < 0) && (zdist < -3)) {
	fidtype |= 0x8;
      }
    }
    /*
    if (mutool.isFiducial(muonObject, MuonDigiCode::CMP)) fidtype |= 0x2;
    if (mutool.isFiducial(muonObject, MuonDigiCode::CMX)) fidtype |= 0x4;
    if (mutool.isFiducial(muonObject, MuonDigiCode::BMU)) fidtype |= 0x8;
    //    if (mutool.isFiducial(muonObject, MuonDigiCode::CSP)) fidtype |= 0x10;
    //    if (mutool.isFiducial(muonObject, MuonDigiCode::CSX)) fidtype |= 0x20;
    //    if (mutool.isFiducial(muonObject, MuonDigiCode::BSU)) fidtype |= 0x40;
    //    if (mutool.isFiducial(muonObject, MuonDigiCode::TSU)) fidtype |= 0x80;
    */
    return fidtype;
  }

  // Likelihood code stolen from Bruce K. The end result is a muoLikelihood
  // function which tells you the likelihood that a CdfMuon is a real muon.


  vector<double> createVector(int n, double v[]) 
  {
    vector<double> ans;
    for(int i=0; i<n; i++)
      ans.push_back(v[i]);
    return(ans);
  }

  int getBin(vector<double> p, vector<vector<double> > edges)
  {
    int whichBin = 0, base = 1;
    assert(p.size()==edges.size());
    for(int i=0; i<p.size(); i++)
      {
	whichBin += base*(upper_bound(edges[i].begin(), edges[i].end(), p[i]) - edges[i].begin());
	base *= edges[i].size()+1;
      }
    return(whichBin);
  }
  
  double muonLikelihood_iso4(double iso4)
  {
    vector<double> p;
    p.push_back(iso4);
    
    double edges0[20] = { 0.0187173, 0.119364, 0.154063, 0.201518, 0.258561, 0.326921, 0.407762, 0.49741, 0.613965, 0.755347, 0.938163, 1.18288, 1.50894, 2.00428, 2.79375, 4.37551, 6.43208, 8.98571, 11.8521, 16.0836 };
    vector<vector<double> > edges = vector<vector<double> >(1);
    edges[0] = createVector(20,edges0);
    
    double l[21] = {   32.3831, 24.0487, 24.0768, 17.5741, 13.8094, 18.8514, 13.8094, 13.1035, 7.82005, 5.5229, 4.95121, 3.6952, 2.22512, 1.31549, 0.768936, 0.196835, 0.0527002, 0.0244818, 0.0117734, 0.0125392, 0.00623436 };
    
    return(log10(l[getBin(p, edges)]));
  }
  
  double muonLikelihood_tIso(double tIso)
  {
    vector<double> p;
    p.push_back(tIso);
    
    double edges0[13] = { 0.0985694, 0.376415, 0.511858, 0.678915, 0.908669, 1.24814, 1.80379, 2.79936, 4.30274, 6.53151, 9.41064, 14.0733, 23.725 };
    vector<vector<double> > edges = vector<vector<double> >(1);
    edges[0] = createVector(13,edges0);
    
    double l[14] = {   8.4709, 4.64937, 4.08554, 3.2654, 2.37204, 1.40995, 0.910167, 0.345624, 0.112021, 0.0611021, 0.0436616, 0.0395559, 0.0294865, 0.0596036 };
    
    return(log10(l[getBin(p, edges)]));
  }
  
  double muonLikelihood_hadE(double hadE)
  {
    vector<double> p;
    p.push_back(hadE);
    
    double edges0[27] = { 0.561405, 1.01, 1.18181, 1.3025, 1.38699, 1.46689, 1.53378, 1.59818, 1.65944, 1.71969, 1.78423, 1.84659, 1.91124, 1.98193, 2.05918, 2.14718, 2.24443, 2.3619, 2.50378, 2.70579, 3.01179, 3.63395, 5.50667, 9.86534, 13.5223, 17.1989, 22.6365 };
    vector<vector<double> > edges = vector<vector<double> >(1);
    edges[0] = createVector(27,edges0);
    
    double l[28] = {   1.48117, 3.13325, 5.65726, 10.8774, 6.57508, 7.36178, 6.94717, 7.14871, 7.36178, 7.14871, 5.52756, 7.82665, 5.06067, 6.40298, 4.75675, 7.37941, 4.31477, 4.57256, 3.58218, 3.27092, 2.34947, 1.39055, 0.397903, 0.0333279, 0.00541998, 0.0114016, 0.014482, 0.0333378 };
    
    return(log10(l[getBin(p, edges)]));
  }
  
  double muonLikelihood_emE(double emE)
  {
    vector<double> p;
    p.push_back(emE);
    
    double edges0[27] = { 0.169595, 0.216775, 0.235961, 0.250576, 0.263205, 0.273381, 0.284246, 0.29406, 0.304078, 0.314175, 0.324811, 0.336067, 0.348997, 0.362548, 0.378541, 0.39628, 0.418619, 0.447135, 0.489225, 0.547881, 0.633144, 0.778566, 1.08025, 1.79586, 4.23103, 8.56469, 17.841 };
    vector<vector<double> > edges = vector<vector<double> >(1);
    edges[0] = createVector(27,edges0);
    
    double l[28] = {   3.5266, 5.96344, 3.78241, 4.2364, 4.86603, 4.47989, 3.75904, 4.10256, 4.03639, 3.76371, 3.6394, 4.1658, 3.42435, 3.24006, 3.69833, 3.09359, 2.71479, 2.96912, 2.49112, 1.81963, 1.54289, 0.736024, 0.601465, 0.336765, 0.0437354, 0.0117834, 0.0106406, 0.0204915 };
    
    return(log10(l[getBin(p, edges)]));
  }
  
  double muonLikelihood(double iso4, double tIso, double hadE, double emE)
  {
    return( muonLikelihood_iso4(iso4) + muonLikelihood_tIso(tIso) + muonLikelihood_hadE(hadE) + muonLikelihood_emE(emE) );
  }
  
  double muonLikelihood(vector<double> vars)
  {
    return(muonLikelihood(vars[0], vars[1], vars[2], vars[3]));
  }
  

  int muoRegion(const CdfMuon& muonObject) {
    int ret = 0;
    if (MuonAnalysisTools::MuonIsInMiniskirt(muonObject)) ret |= 0x1;
    if (MuonAnalysisTools::MuonIsInBluebeam(muonObject)) ret |= 0x2;
    if (MuonAnalysisTools::MuonIsInKeystone(muonObject)) ret |= 0x4;
    return ret;
    
  }

}   // end of name space bracket