////////////////////////////////////////////////////////////////////
//
//     File: MuonVariables.cc
//     Authors: Abe DeBenedetti, routines writen by Vladimir Rekovic
//     Description: High level analysis muon variables defined as functions
//     Created: 7/26/2001
//
//     Revision: none so far
//
////////////////////////////////////////////////////////////////////

//#include <cstdio>
//#include <math.h>
//#include "Rtypes.h"

#include "HighLevelObjects/MuonVariables.hh"
#include "MuonObjects/CdfMuon.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "CLHEP/Vector/LorentzVector.h"

//#include "CalorGeometry/CalParameters.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"

namespace MuonVariables {

  // function returnig reference to structure CdfMuon::cmpMatch
  // holding delta and chi-square data.
  const MuonMatchData& cmpMatch(const CdfMuon& muonObject) {
    return muonObject.cmp();
  } 

  // function accessing dphi, memeber of structure cmpMatch
  float cmpMatchDphi(const CdfMuon& muonObject) {
    return muonObject.cmp().dphi();
  }

  // function accessing drphi, memeber of structure cmpMatch
  float cmpMatchDrphi(const CdfMuon& muonObject) {
    return muonObject.cmp().drphi();
  }

  // function accessing chsqX, memeber of structure cmpMatch
  float cmpMatchChsqX(const CdfMuon& muonObject) {
    return muonObject.cmp().chsqX();
  }
  
  // function accessing chsqXPosition, memeber of structure cmpMatch
  float cmpMatchChsqXPosition(const CdfMuon& muonObject) {
    return muonObject.cmp().chsqXPosition();
  }


  // function returning reference to structure CdfMuon::cmuMatch
  // holding delta and chi-square data.
  const MuonMatchData& cmuMatch(const CdfMuon& muonObject) {
    return muonObject.cmu();
  } 

  // function accessing dphi, memeber of structure cmuMatch
  float cmuMatchDphi(const CdfMuon& muonObject) {
    return muonObject.cmu().dphi();
  }

  // function accessing drphi, memeber of structure cmuMatch
  float cmuMatchDrphi(const CdfMuon& muonObject) {
    return muonObject.cmu().drphi();
  }

  // function accessing dtheta, memeber of structure cmuMatch
  float cmuMatchDtheta(const CdfMuon& muonObject) {
    return muonObject.cmu().dtheta();
  }

  // function accessing dz, memeber of structure cmuMatch
  float cmuMatchDz(const CdfMuon& muonObject) {
    return muonObject.cmu().dz();
  }

  // function accessing chsqZ, memeber of structure cmuMatch
  float cmuMatchChsqZ(const CdfMuon& muonObject) {
    return muonObject.cmu().chsqZ();
  }

  // function accessing chsqX, memeber of structure cmuMatch
  float cmuMatchChsqX(const CdfMuon& muonObject) {
    return muonObject.cmu().chsqX();
  }

  // function accessing chsqZPosition, memeber of structure cmuMatch
  float cmuMatchChsqZPosition(const CdfMuon& muonObject) {
    return muonObject.cmu().chsqZPosition();
  }

  // function accessing chsqXPosition, memeber of structure cmuMatch
  float cmuMatchChsqXPosition(const CdfMuon& muonObject) {
    return muonObject.cmu().chsqXPosition();
  }

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


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


  //Tracks


  // best track
   CdfTrack_clnk getBestTrack(const CdfMuon& muonObject) {
    return muonObject.bestTrack();
  } 

  // access to best track parameters
  double bestTrackCurvature(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->curvature();
  }

  //const CdfMuon& testMuon;
  //double myTrackLink = getBestTrack(testMuon)->curvature();


  // access to best track parameters
  double bestTrackD0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->d0();
  }

  // access to best track parameters
  double bestTrackZ0(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->z0();
  }

  // access to best track parameters
  double bestTrackLambda(const CdfMuon& muonObject) {
    return muonObject.bestTrack()->lambda();
  }


  // second best track
   CdfTrack_clnk getSecondBestTrack(const CdfMuon& muonObject) {
    return muonObject.secondBestTrack();
  } 

  // access to best track parameters
  double secondBestTrackCurvature(const CdfMuon& muonObject) {
    return muonObject.secondBestTrack()->curvature();
  }

  //const CdfMuon& testMuon;
  //double myTrackLink = getBestTrack(testMuon)->curvature();


  // access to best track parameters
  double secondBestTrackD0(const CdfMuon& muonObject) {
    return muonObject.secondBestTrack()->d0();
  }

  // access to best track parameters
  double secondBestTrackZ0(const CdfMuon& muonObject) {
    return muonObject.secondBestTrack()->z0();
  }

  // access to best track parameters
  double secondBestTrackLambda(const CdfMuon& muonObject) {
    return muonObject.secondBestTrack()->lambda();
  }



  // four Momentum
  HepLorentzVector getFourMomentum(const CdfMuon& muonObject) {
    return muonObject.fourMomentum();
  } 

  
  // test function
  bool isItCmu(const CdfMuon& muonObject) {
    return muonObject.hasCmu();
  } 

  // test function
  bool isItCmp(const CdfMuon& muonObject) {
    return muonObject.hasCmp();
  } 
  // test function
  bool isItCmx(const CdfMuon& muonObject) {
    return muonObject.hasCmx();
  } 
  // test function
  bool isItImu(const CdfMuon& muonObject) {
    return muonObject.hasBmu();
  } 



  /*
    // set lastVersion to the last version you know of to print changes
  int phoVariablesVersion(int lastVersion) {
    int currentVersion = 2;
    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;
    }
    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)>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(etiso>4.0) word |= (1<<12);
    if(phoTrackIso(emObject,0.4,zvertex)>2.0) word |= (1<<13);
    if(phoTrackIso(emObject,0.4,zvertex)>5.0) word |= (1<<14);
    double cesslide = phoCesSlide(emObject);
    bool qces = false;
    if( phoCesStripE2(emObject)<cesslide && 
	(phoCesWireE2(emObject)<cesslide || 
	 fabs( phoCesX(emObject) - phoCesX2(emObject))<7.0) )
	qces = true;
    if(!qces) word |= (1<<15);
    double az = fabs(phoCesZ(emObject));
    if(az<14.0 || az>217.0) word |= (1<<16);
    if(phoPt(emObject)>1.0) word |= (1<<17);
    if(fabs(phoCesX(emObject))>21.5) word |= (1<<18);
    double CprX,CprQ;
    int stat = phoCprCesSeeded(emObject, CprX, CprQ, 5, zvertex);
    if(stat==0 && CprQ>500.0) word |= (1<<19);
    if(!phoHasShowerMax(emObject)) word |= (1<<20);

    //skip a few for central

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

    return word;
  }

    // status word = 100, 200... + 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();
  }


  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
  double 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 sumPt = 0;

    Hep3Vector p = phoFourVector(emObject,zvertex).vect();

    CdfTrackView_h trview;
    // this should be the best tracklist available, i.e. use SVX if it exists
    CdfTrackView::Error trackStatus = CdfTrackView::defTracks(trview);
    if(trackStatus == CdfTrackView::OK) {
      const CdfTrackView::CollType& tracks = trview->contents();
      //select all tracks within cone
      CdfTrackView* pTrkInCone = new CdfTrackView(tracks.begin(),tracks.end(),
					  track_cut::InConeAround(p,cone));

      //select within cone and within 5cm to vertex and 2 ax layers with 4 hits
      //and two stereo layers with 2 hits
      //CdfTrackView* pTrkInCone = 
      //new CdfTrackView(tracks.begin(),tracks.end(),
      //HasCOThits(0,0,2,2,4,2)&&InConeAround(p,cone)
      //&&Z0WrtZInRange(zvertex,5.0));
      const CdfTrackView::CollType& conetracks = pTrkInCone->contents();
      for(CdfTrackView::const_iterator itrack=conetracks.begin();
	  itrack!=conetracks.end(); ++itrack) {
	sumPt += (*itrack)->pt();
      }
      pTrkInCone->destroy();
    }              // found a track list
    return sumPt;

  }
  */

} // namespace PhotonVariables