#include <float.h>
#include <iostream>
#include <math.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>

#include "TPhotonUtil.hh"
#include "TExtrap.hh"

#include "Stntuple/obj/TStnTrack.hh"
#include "Stntuple/obj/TStnVertexBlock.hh"
#include "Stntuple/obj/TStnHeaderBlock.hh"
#include "Stntuple/obj/TStnJetBlock.hh"
#include "Stntuple/obj/TStnElectronBlock.hh"
#include "Stntuple/obj/TStnPhotonBlock.hh"
#include "Stntuple/obj/TStnJetBlock.hh"
#include "Stntuple/obj/TStnJet.hh"
#include "Stntuple/obj/TStnElectronBlock.hh"
#include "Stntuple/obj/TPhoenixElectronBlock.hh"
#include "Stntuple/obj/TStnClusterBlock.hh"
#include "Stntuple/obj/TCalDataBlock.hh"
#include "Stntuple/obj/TCesDataBlock.hh"
#include "Stntuple/obj/TCprDataBlock.hh"
#include "Stntuple/obj/TStnMetBlock.hh"
#include "Stntuple/obj/TStnMuonBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "Stntuple/obj/TStnTriggerBlock.hh"
#include "Stntuple/obj/TStnEmTimingBlock.hh"
#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/obj/TStnDBManager.hh"
#include "Stntuple/alg/TStntuple.hh"

//#include "HighLevelObjects/PhotonBackgroundComputer.hh"

using std::cout;
using std::endl;

ClassImp(TPhotonUtil)

Int_t TPhotonUtil::L2ClusterMatch(Double_t eta, Double_t phi, 
				TStnTriggerBlock* trigBlock,
				int pass, int clusOrIso) 
{

  // convert seed trigger tower to eta
  float etatr[24] = 
    {-3.0,-2.3,-1.9,-1.6,-1.4,-1.2,-1.1,-0.9,-0.7,-0.5,-0.3,-0.1,
     0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.2, 1.4, 1.6, 1.9, 2.3, 3.0};

  int iclu = -1;
  int iiso = -1;

  if(trigBlock==NULL) return -1;
  TTl2d& tl2d = *(trigBlock->Tl2d());


  float dphi,deta,dr,drbest;
  drbest = 1000.0;
  if(clusOrIso==0) {
    int n = tl2d.NL2Clusters();
    for(int i=0; i<n; i++) {
      TTl2dCluster& clus =  *(tl2d.Cluster(i));
      if( clus.Pass()==pass) {
	dphi = fabs(phi - clus.Phi());
	if(dphi>M_PI) dphi = 2*M_PI - dphi;
	deta = fabs(eta - etatr[clus.IEta()]);
	dr = sqrt(dphi*dphi + deta*deta);
	if(dr<0.4 && dr<drbest) {
	  iclu = i;
	  drbest=dr;
	}
      }
    }
    return iclu;
  } else {
    int ni= tl2d.NL2IsoClusters();
    for(int i=0; i<ni; i++) {
      TTl2dClusterIso& isoc =  *(tl2d.ClusterIso(i));
      if( isoc.ClusteringPath()==pass) {
	dphi = fabs(phi - (isoc.Phi()+0.5)*2.0*M_PI/24.0);
	if(dphi>M_PI) dphi = 2.0*M_PI - dphi;
	deta = fabs(eta - etatr[isoc.Eta()]);
	dr = sqrt(dphi*dphi + deta*deta);
	if(dr<04. && dr<drbest) {
	  iiso = i;
	  drbest=dr;
	}
      }
    }
    return iiso;
  }

}

//_____________________________________________________________________________
Int_t  TPhotonUtil::SetPhotonZ(TStnPhoton* pho, 
				Double_t z, 
				TClonesArray* trackList) 
{
// Correct photon quantites for a new z vertex.  Track list is needed for
// Track iso, can be retrieved from TrackBlock or set to NULL (iso=999.0)
// note that BoxIso, TowEt, TowEta are not corrected by definition

  double hemi = (pho->fEveta>0.0? 1.0 : -1.0);
  double sint =  pho->fSth;
  double cott = hemi*sqrt(1.0-sint*sint)/sint;

  double corr,newsint,newcott;
  if(pho->fDetector==0) {  // central
    double dz = (183.9*cott+pho->fZv) - z;
    newcott = dz/183.9;
    newsint = 183.9/sqrt(183.9*183.9 + dz*dz);
    corr = newsint/sint;
  } else {  // plug
    double r = (hemi*185.4-pho->fZv)/cott;
    double dz = hemi*185.4-z;
    newcott = dz/r;
    newsint = r/sqrt(r*r+dz*dz);
    corr = newsint/sint;
  }

  pho->fSth    = newsint;
  pho->fEt     *= corr;
  pho->fEtc    *= corr;
  pho->fCo4    *= corr;
  pho->fCo4PJW *= corr;
  pho->fCo7    *= corr;
  pho->fEveta  = log(1.0/newsint + newcott);

  double e = pho->fMomentum.E();
  double cost  = hemi*sqrt(1.0-newsint*newsint);
  double phi = pho->fPhi;
  pho->fMomentum = TLorentzVector(e*cos(phi)*newsint,
				  e*sin(phi)*newsint, e*cost, e);

  pho->fSumpt4 = TrkIsoSum0(pho->fPhi, pho->fEveta, z, 
			    trackList, 0.4, 5.0, NULL);

  pho->fZv = z;

  //this is difficult to correct, code-management-wise
  // set to 0 for now
  pho->fCprwht = 0.0;
  //if(pho->fCprwht>-98.0 && abs(pho->fCprwht)>1.e-6) {
  //  // if cpr wht was valid, then need to recompute based on the new sinTheta
  //  pho->fCprwht =  PhotonBackgroundComputer::cprWeight(
  //				   pho->fCescprq,pho->fEt,pho->fSth);
  //}

  return 0;
}

//_____________________________________________________________________________
Int_t  TPhotonUtil::SetAllPhotonZ(TStnEvent* event, Double_t z) 
{

  if(z<-998.0) z=0.0;

  int ientry = event->GetCurrentTreeEntry();

  TStnPhotonBlock* fPhotonBlock = 
    (TStnPhotonBlock*)event->GetDataBlock("PhotonBlock");
  fPhotonBlock->GetEntry(ientry);

  bool corr=false;
  for(int i=0; i<fPhotonBlock->NPhotons(); i++) {
    TStnPhoton* photon = fPhotonBlock->Photon(i);
    if( fabs(photon->VertexZ()-z)>0.001 ) corr=true;
  }

  if(corr) {
    TStnTrackBlock* fTrackBlock = 
      (TStnTrackBlock*)event->GetDataBlock("TrackBlock");
    fTrackBlock->GetEntry(ientry);
    for(int i=0; i<fPhotonBlock->NPhotons(); i++) {
      TStnPhoton* photon = fPhotonBlock->Photon(i);
      TStntuple::SetPhotonZ(photon, z, fTrackBlock->TrackList() );
    }
  }
  return 0;
}


//_____________________________________________________________________________
Int_t  TPhotonUtil::CorrectPhotonEnergy(TStnEvent* event) 
{

  int ientry = event->GetCurrentTreeEntry();

  TStnPhotonBlock* fPhotonBlock = 
    (TStnPhotonBlock*)event->GetDataBlock("PhotonBlock");
  if(!fPhotonBlock) {
    printf("TPhotonUtil::CorrectPhotonEnergy will not work without Regitering PhotonBlock\n");
    return 1;
  }
  fPhotonBlock->GetEntry(ientry);

  TStnHeaderBlock* fHeaderBlock = 
    (TStnHeaderBlock*)event->GetDataBlock("HeaderBlock");
  bool qMc = fHeaderBlock->McFlag(); 
  int run = fHeaderBlock->RunNumber(); 

  for(int i=0; i<fPhotonBlock->NPhotons(); i++) {
    TStnPhoton* pho = fPhotonBlock->Photon(i);
    float corr = 1.0;
    if(pho->Detector() == 0) { // central
      if(qMc) {
	  corr *= 0.996;
      } else {
      }
    } else { // plug
      // "face" corrections should already be in Etc and Momentum
      // add PPR energy
      corr = (pho->E()+pho->PprE())/pho->E();
      // leakage correction
      corr *= 1.0 + pho->fLcor;
      // final scale factors for 5.3.3
      if(qMc) {
	if(fabs(pho->DetEta())>1.78) {
	  corr *= 0.998;
	} else {
	  corr *= 0.996;
	}
      } else {
	if(pho->DetEta()<-1.78) {
	  corr *= 1.007;
	} else if(pho->DetEta()<0.0) {
	  corr *= 1.015;
	} else if(pho->DetEta()<1.78) {
	  corr *= 1.020;
	} else {
	  corr *= 1.010;
	}
      }
    }

    pho->fMomentum *= corr;
    pho->fEtc *= corr;
    // note that Et() and E() remain uncorrected
  }

  return 0;
}


//_____________________________________________________________________________
Int_t  TPhotonUtil::CorrectPhotonIsolation(TStnEvent* event) 
{

  int ientry = event->GetCurrentTreeEntry();

  TStnPhotonBlock* fPhotonBlock = 
    (TStnPhotonBlock*)event->GetDataBlock("PhotonBlock");
  if(!fPhotonBlock) {
    printf("TPhotonUtil::CorrectPhotonEnergy will not work without Regitering PhotonBlock\n");
    return 1;
  }
  fPhotonBlock->GetEntry(ientry);

  float z;
  int nv12;
  TPhotonUtil::SelectVertex(event, z, nv12, 12);

  for(int i=0; i<fPhotonBlock->NPhotons(); i++) {
    TStnPhoton* pho = fPhotonBlock->Photon(i);
    pho->fVcor = (nv12>1? 0.3563*(nv12-1) : 0.0 );
    pho->fCo4PJW = pho->fCo4 - pho->fLcor*pho->fEtc - pho->fVcor;
  }

  return 0;
}



//_____________________________________________________________________________
Double_t  TPhotonUtil::TrkIsoSum1(TStnTrack*    track,
			        TClonesArray* trackList,
			        Double_t      cone,
			        Double_t      dzCut,
			        TStnTrack*    skipTrack)
{
  return TPhotonUtil::TrkIsoSum0( track->Phi0(), track->Eta(), track->Z0(),
			        trackList, cone, dzCut, skipTrack);
}

//_____________________________________________________________________________
Double_t  TPhotonUtil::TrkIsoSum0(Double_t phi, Double_t eta, Double_t z, 
			        TClonesArray* trackList,
			        Double_t      cone,
			        Double_t      dzCut,
			        TStnTrack*    skipTrack)
{
  // perform sum of Pt of tracks in cone in phi-eta, centered at z.
  // all tracks have to pass delta-z and cone cuts
  // method skips the indicated track in the trackList

  double  phiTk, etaTk, zTk, ptTk;
  double  dphi, deta, dr, dz, sum_pt;
  int     ntrk;


  if(trackList==NULL) return 999.0;

  sum_pt = 0;
  ntrk   = trackList->GetEntriesFast();

  for (int i=0; i<ntrk; i++) {
    TStnTrack* ti = (TStnTrack*) trackList->UncheckedAt(i);
    if (ti != skipTrack) {
      TLorentzVector* mi = ti->Momentum();
      phiTk = ti->Phi0();
      etaTk = ti->Eta();
      ptTk  = ti->Pt();
      zTk   = ti->Z0();

      dphi = TVector2::Phi_mpi_pi(phi-phiTk);
      deta = eta - etaTk;
      dr   = sqrt(dphi*dphi+deta*deta);
      dz   = fabs(z - zTk);

      if (dr < cone && dz < dzCut ) {
	sum_pt += ptTk;
      }

    }
  }

  return sum_pt;

}


//_____________________________________________________________________________
Int_t TPhotonUtil::DumpPhotons(TStnEvent* event, float photonEt) {
  // remove electrons, photons etc  from the list of event cone 0.4 jets
  // assuming we know the branches with the jets, electrons and photons

	 			// loop over all the jets in cone 0.4
  TStnJet       *jet;
  TStnPhoton    *photon;
  TCalTower     *tower;

  TStnHeaderBlock*   hdata;
  TStnPhotonBlock*   pdata;
  TStnJetBlock*      jdata;
  TStnMetBlock*      mdata;
  TStnMuonBlock*     udata;
  TStnElectronBlock*  edata;
  //TStnClusterBlock*  cdata;
  TCalDataBlock*     caldata;
  //TCesDataBlock*     cesdata;
  //TCprDataBlock*     cprdata;
  TStnVertexBlock*     vdata;
  TStnEmTimingBlock*     emtdata;

  hdata = (TStnHeaderBlock*)   event->GetDataBlock("HeaderBlock");
  pdata = (TStnPhotonBlock*)   event->GetDataBlock("PhotonBlock");
  jdata = (TStnJetBlock*)      event->GetDataBlock("JetBlock");
  mdata = (TStnMetBlock*)      event->GetDataBlock("MetBlock");
  udata = (TStnMuonBlock*)     event->GetDataBlock("MuonBlock");
  edata = (TStnElectronBlock*)     event->GetDataBlock("ElectronBlock");
  //cdata = (TStnClusterBlock*)  event->GetDataBlock("ClusterBlock");
  caldata = (TCalDataBlock*)   event->GetDataBlock("CalDataBlock");
  vdata = (TStnVertexBlock*)   event->GetDataBlock("ZVertexBlock");
  emtdata = (TStnEmTimingBlock*)   event->GetDataBlock("EmTimingBlock");

  int ientry =  event->GetCurrentTreeEntry();

  if(hdata) hdata->Print();

  float met=999,metphi=99.;
  if(mdata) {
    met = mdata->Met(0);
    metphi = mdata->MetPhi(0);
    cout <<"Met, MetPhi " << met << "  " << metphi << endl;
  } else {
    cout <<"Could not find Met block" << endl;
  }
  if(vdata) {
    printf("Number of Z vertices: %i\n",vdata->NVertices());
    for(int i=0; i<vdata->NVertices(); i++) {
      TStnVertex* vtx = vdata->Vertex(i);
      printf("z= %5.2f, class=%3i, Sumpt=%5.1f Ntrk=%3i \n",
	     vtx->Z(),vtx->VClass(),vtx->SumPt(),vtx->NTracks());
    }
  } else {
    cout <<"Could not find Z Vertex block" << endl;
  }

    
  caldata->GetEntry(ientry);
  if(emtdata!=NULL && emtdata->IsValid()) emtdata->GetEntry(ientry);

  for (int i=0; i<pdata->fNPhotons; i++) {
    photon   = pdata->Photon(i);
    if(photon->Etc()>photonEt) {
    cout << "*************** Photon " << i << " ************* " << endl;
    cout << "Detector, Status, Zv \t\t" 
	 << (photon->fDetector==0? "CEM" : "PEM")
	 << "  " << photon->fStat << "  " << photon->fZv << endl;
    cout << "Et, Etcorr, E  \t\t\t" << photon->fEt 
	 << "  " << photon->fEtc << "  " << photon->fE << endl;
    cout << "Momentum  \t\t\t" << photon->fMomentum.X()
	 << "  " << photon->fMomentum.Y()
	 << "  " << photon->fMomentum.Z()
	 << "  " << photon->fMomentum.T() << endl;
    cout << "Iso 0.4, 0.7, Lcor, L+Vcor \t" << photon->fCo4 
	 << "  " << photon->fCo7 << "  " << photon->EIso4(1) 
	 << "  " << photon->EIso4(2)<< endl;
    cout << "DetEta, EvEta, SinTheta, Phi \t" << photon->fDteta
	 << "  " << photon->fEveta 
	 << "  " << photon->fSth
	 << "  " << photon->fPhi << endl;
    cout << "Had/Em, SlidingCut, LShr \t" << photon->fHadem 
	 << "  " << 0.055+0.00045*photon->fE
	 << "  " << photon->fLshr << endl;
    cout << "N3D, Pt1 Pt2, Sumpt4 \t\t" << photon->fN3d
	 << "  " << photon->fPt
	 << "  " << photon->fPt2
	 << "  " << photon->fSumpt4 << endl;
    if(photon->fDetector==0) {
      cout << "CesX CesZ WireE \t\t" << photon->fCesx
	   << "  " << photon->fCesz 
	   << "  " << photon->fCese << endl;
      cout << "ChiAvg, Wire, Strip \t\t" << photon->fChi
	   << "  " << photon->fChiwir
	   << "  " << photon->fChistr << endl;
      cout << "2nd Ces X, Z, WirE, StrE slide " << photon->fWwir2
	   << "  " << photon->fStr2
	   << "  " << photon->fWire2
	   << "  " << photon->fStre2 
	   << "  " << photon->fCesslide << endl;
      cout << "Unbiased CprX CprQ \t\t" << photon->fCpr5ps
	   << "  " << photon->fCpr5ph << endl;
      cout << "Ces seeded CprX CprQ \t\t" << photon->fCescprx
	   << "  " << photon->fCescprq << endl;
      cout << "Ces, Cpr Weight \t\t" 
	   << photon->fCeswht
	   << "  " << photon->fCprwht << endl;

    } else {
      cout << "3x3 chi2 \t\t\t" << photon->fChi3x3 << endl;
      cout << "PES E, X, Y \t\t\t" 
	   << "  " << photon->fPesE
	   << "  " << photon->fPesX
	   << "  " << photon->fPesY << endl;
	cout << "2nd PES E, PPR E\t\t" 
	   << "  " << photon->fPesE2
	   << "  " << photon->fPprE << endl;
    }

    cout << "L3 Et,Eta,Phi,Iso \t\t" 
	 << photon->fTowEt
	 << "  " << photon->fTowEta
	 << "  " << photon->fTowPhi 
	 << "  " << photon->fTowIso << endl;
    cout << "2-tower Had/Em \t\t\t" << photon->fHademT << endl;
    cout << "HTDC time \t\t\t" << photon->fTime << endl;
    cout << "Likelihood \t\t\t" << photon->fLike << endl;

    cout << "Towers:  Seed phi, eta \t" 
	 << "  " << photon->PhiSeedIndex()
	 << "  " << photon->EtaSeedIndex() << endl;

    cout << "Cosmic Stubs:  all, costub, 30deg: \t" 
	 << "  " << photon->NStub()
	 << "  " << photon->NCosStub()
	 << "  " << photon->NCosStubPho() << endl;

    bool ccut = true;
    if(photon->Sidewedges()<=2) {
      if( (photon->Seedwedge() - 1.5*photon->Sidewedges())>=6 ) ccut = false;
    } else if(photon->Sidewedges()<=12) {
      if( (photon->Seedwedge() - 0.5*photon->Sidewedges())>=8 ) ccut = false;
    } else {
      if(photon->Seedwedge()>=14) ccut = false;
    }
    cout << "Cosmic:  seed, sidewedges, cut: \t" 
	 << "  " << photon->Seedwedge()
	 << "  " << photon->Sidewedges()
	 << "  " << ccut << endl;

    cout << "Cosmic towers:  contig, East, West: \t" 
	 << "  " << photon->Ncontig()
	 << "  " << photon->HaloEast()
	 << "  " << photon->HaloWest() << endl;

    cout << "IPhi  IEta  EmEt  HadEt  Time0  Time1 EmTime  Assym" << endl;

    const TStnLinkBlock*  links = pdata->TowerLinkList();
    int nlink = links->NLinks(i);
    for(int ilink=0; ilink<nlink; ilink++) {
      int l = links->Index(i,ilink);
      int iphi = l & 0x3F;
      int ieta = (l>>8) & 0x3F;
      float emtime=-99.0;
      if(emtdata!=NULL && emtdata->IsValid()) {
	for(int iemt=0; iemt<emtdata->NHits(); iemt++) {
	  if(emtdata->IPhi(iemt) == iphi &&
	     emtdata->IEta(iemt) == ieta) {
	    emtime = emtdata->Time(iemt);
	  }
	}
      }

      //int istt = (l>>16) & 0xFFFF;
      //cout << "  " << iphi << "    " << ieta;
      printf("  %2d   %2d",iphi,ieta);
      tower=NULL;
      if(caldata) tower=caldata->Tower(ieta,iphi);
      if(tower) {
	
	float assym = float(tower->EmPmt(0)-tower->EmPmt(1))/float(tower->EmPmt(0)+tower->EmPmt(1));
	printf("%7.2f%7.2f%7.1f%7.1f%7.1f%7.3f\n",
	       tower->EmEnergy()*photon->fSth,tower->HadEnergy()*photon->fSth,
	       tower->HadTdc(0),tower->HadTdc(1),emtime,assym);
      } else {
	cout << "  can't find tower in CalData" << endl;
      }

    } // loop over links


    float dphi,deta,dr,jphi,dphim;
    if(jdata) {
      cout << endl;
      cout << " jet    Et      phi   dteta      dphi      deta      dr    dphi_met" 
	   << endl;
      for (int j=0; j<jdata->NJets(); j++) {
	jet = jdata->Jet(j);
	jphi = jet->Momentum()->Phi();
	if(jphi<0.0) jphi+=2.0*M_PI;
	dphi = fabs(jphi - photon->fPhi);
	if(dphi>M_PI) dphi = 2.0*M_PI-dphi;
	dphim = fabs(jphi - metphi);
	if(dphim>M_PI) dphim = 2.0*M_PI-dphim;
	deta = jet->DetEta() - photon->fDteta;
	dr = sqrt(dphi*dphi + deta*deta);
	cout << j 
	     << "  " << jet->Et()
	     << "  " << jphi 
	     << "  " << jet->DetEta()
	     << "  " << dphi 
	     << "  " << deta 
	     << " " << dr 
	     << " " << dphim 
	     << endl;
      }
    } else {
      cout << "No Jet block found" << endl;
    }


    }  // if Et passes cut
  } // loop over photons

  if(udata) {
    udata->Print();
  } else {
    cout << "No Muon block found" << endl;
  }
  if(edata) {
    edata->Print();
  } else {
    cout << "No Electron block found" << endl;
  }

  return 0;
}


//_____________________________________________________________________________
Double_t TPhotonUtil::ExpectedHadEm(Double_t energy) {

  // following the PDG Review of Particle 
  // Properties 2000 Sec. 23 pg 21 at http://pdg.lbl.gov/
  double e_critical = 0.100; // 100 MeV -- this is rough, but close
  double y = energy/e_critical;
  // This value is chosen to set tmax and the expected 
  // had fraction in roughly the right place
  double b = 0.60; 
  // location of shower max, in units of radiation lengths
  double tmax = log(y); 
  double a = b*tmax+1;

  // for CDF.  This number oscillates with theta 
  // -- see my radiation length plot created with cdfSim
  double emCalorimeterRadiationLengths = 18.3; 

// the fraction of energy expected to fall in the hadronic calorimeter
// double ans = 1.-incompleteGammaFunction(a,b*emCalorimeterRadiationLengths); 
  double ans = 1.-TMath::Gamma(a,b*emCalorimeterRadiationLengths); 

  return(ans);

}

//_____________________________________________________________________________
Double_t TPhotonUtil::CprTrack(TStnPhoton* pho, TStnTrackBlock* tblk) {

  if(pho->Detector()!=0) return 999.0;
  if(fabs(pho->ZCes())>250.0) return 999.0;
  if(fabs(pho->XCes())> 30.0) return 999.0;

  float zv=0;
  float cesz = pho->ZCes();
  if(fabs(pho->VertexZ())<150.0) zv = pho->VertexZ();
  float cprz = 0.91*(cesz-zv) + zv;
  int module;
  if(cprz<0.0) {
    module = (cprz<-122.0? 0 : 1 );
  } else {
    module = (cprz< 122.0? 2 : 3 );
  }
  //cprx is in +phi sense, opposite to Ces x
  float cprx = -0.91*pho->XCes();

  float dxmin = 99.0;
  TExtrap extrap;
  for(int i=0; i<tblk->NTracks(); i++){
    extrap.SetParameters(tblk->Track(i));
    if(extrap.GetCprWedge()==pho->PhiSeedIndex()) {
      int modulet;
      float crpzt = extrap.GetCprLocalZ();
      if(crpzt<0.0) {
	modulet = (crpzt<-122.0? 0 : 1 );
      } else {
	modulet = (crpzt< 122.0? 2 : 3 );
      }
      if(modulet==module) {
	float cprxt = extrap.GetCprLocalX();
	if(fabs(cprxt-cprx)< fabs(dxmin)) dxmin = cprxt - cprx;
      }
    }
  }

  return dxmin;

}



//_____________________________________________________________________________
Bool_t  TPhotonUtil::SelectVertex(TStnEvent* event, Float_t& z, Int_t& nv12, 
				  int minClass, bool useVxprimToo)
{

  int ientry = event->GetCurrentTreeEntry();

  TStnVertexBlock* fZVertexBlock = 
    (TStnVertexBlock*)event->GetDataBlock("ZVertexBlock");
  fZVertexBlock->GetEntry(ientry);

  nv12 = 0;
  z = 0.0;
  int vclassMax = -1;

  float ptmax = -1;
  int nv = fZVertexBlock->NVertices();
  for(int iv=0; iv<nv; iv++) {
    float pt = fZVertexBlock->Vertex(iv)->SumPt();
    int vclass = fZVertexBlock->Vertex(iv)->VClass();
    if(pt>ptmax) {
      ptmax = pt;
      z = fZVertexBlock->Vertex(iv)->Z();
      vclassMax = vclass;
    }
    if(vclass>=12) nv12++;
  }

  if(ptmax<0.0 && useVxprimToo) {
    TStnVertexBlock* fVertexBlock = 
      (TStnVertexBlock*)event->GetDataBlock("VertexBlock");
    fVertexBlock->GetEntry(ientry);
    int nv = fVertexBlock->NVertices();
    for(int iv=0; iv<nv; iv++) {
      int nt = fVertexBlock->Vertex(iv)->NTracks();
      if(nt>ptmax) {
	ptmax = nt;
	z = fVertexBlock->Vertex(iv)->Z();
	vclassMax = fZVertexBlock->Vertex(iv)->VClass();
      }
    }
  }

  if(ptmax<0.0 || vclassMax<minClass || fabs(z)>60.0) {
    z = 0.0;
    return false;
  }

  return true;
}


//_____________________________________________________________________________
TStnJet*  TPhotonUtil::MatchJet(TStnEvent* event, TStnPhoton* pho)
{

  TStnJet* mjet=NULL;
  if(pho==NULL) return mjet;

  int ientry = event->GetCurrentTreeEntry();
  TStnJetBlock* fJetBlock = (TStnJetBlock*)
    event->GetDataBlock("JetBlock");

  float dmin=999.0;
  for(int i=0; i<fJetBlock->NJets(); i++) {
    TStnJet* jet = fJetBlock->Jet(i);
    float dr = TStntuple::DeltaR(pho->Phi(),pho->DetEta(),
				 jet->Phi(),jet->DetEta());
    if(dr<dmin) {
      mjet = jet;
      dmin = dr;
    }
  }
  return mjet;

}


//_____________________________________________________________________________
TStnElectron*  TPhotonUtil::MatchElectron(TStnEvent* event, TStnPhoton* pho)
{

  TStnElectron* mele=NULL;
  if(pho==NULL) return mele;

  int ientry = event->GetCurrentTreeEntry();
  TStnElectronBlock* fElectronBlock = (TStnElectronBlock*)
    event->GetDataBlock("ElectronBlock");
  fElectronBlock->GetEntry(ientry);

  float dmin=999.0;
  for(int i=0; i<fElectronBlock->NElectrons(); i++) {
    TStnElectron* ele = fElectronBlock->Electron(i);
    float dr = TStntuple::DeltaR(pho->Phi(),pho->DetEta(),
				 ele->EmClusPhi(),ele->DetEta());
    if(dr<dmin) {
      mele = ele;
      dmin = dr;
    }
  }
  return mele;

}


//_____________________________________________________________________________
TStnElectron*  TPhotonUtil::MatchPhoenixElectron(TStnEvent* event, 
						 TStnPhoton* pho)
{

  TStnElectron* mele=NULL;
  if(pho==NULL) return mele;

  int ientry = event->GetCurrentTreeEntry();
  TPhoenixElectronBlock* fElectronBlock = (TPhoenixElectronBlock*)
    event->GetDataBlock("Phoenix_Electrons");
  if(!fElectronBlock) return NULL;

  fElectronBlock->GetEntry(ientry);

  float dmin=999.0;
  for(int i=0; i<fElectronBlock->NElectrons(); i++) {
    TStnElectron* ele = fElectronBlock->Electron(i);
    float dr = TStntuple::DeltaR(pho->Phi(),pho->DetEta(),
				 ele->EmClusPhi(),ele->DetEta());
    if(dr<dmin) {
      mele = ele;
      dmin = dr;
    }
  }
  return mele;

}


//_____________________________________________________________________________
TStnElectron*  TPhotonUtil::MatchPhoenixElectron(TStnEvent* event, 
						 TStnElectron* ele)
{

  TStnElectron* mele=NULL;
  if(ele==NULL) return mele;

  int ientry = event->GetCurrentTreeEntry();
  TPhoenixElectronBlock* fPhElectronBlock = (TPhoenixElectronBlock*)
    event->GetDataBlock("Phoenix_Electrons");
  fPhElectronBlock->GetEntry(ientry);

  float dmin=999.0;
  for(int i=0; i<fPhElectronBlock->NElectrons(); i++) {
    TStnElectron* phele = fPhElectronBlock->Electron(i);
    float dr = TStntuple::DeltaR(ele->EmClusPhi(),  ele->DetEta(),
				 phele->EmClusPhi(),phele->DetEta());
    if(dr<dmin) {
      mele = phele;
      dmin = dr;
    }
  }
  return mele;

}


//_____________________________________________________________________________
TStnPhoton*  TPhotonUtil::MatchPhoton(TStnEvent* event, TStnElectron* ele)
{

  TStnPhoton* mpho=NULL;
  if(ele==NULL) return mpho;

  int ientry = event->GetCurrentTreeEntry();
  TStnPhotonBlock* fPhotonBlock = (TStnPhotonBlock*)
    event->GetDataBlock("PhotonBlock");
  fPhotonBlock->GetEntry(ientry);

  float dmin = 999.0;
  for(int i=0; i<fPhotonBlock->NPhotons(); i++) {
    TStnPhoton* pho = fPhotonBlock->Photon(i);
    float dr = TStntuple::DeltaR(pho->Phi(),pho->DetEta(),
				 ele->EmClusPhi(),ele->DetEta());
    if(dr<dmin) {
      mpho = pho;
      dmin = dr;
    }
  }
  return mpho;

}



//_____________________________________________________________________________
Int_t  TPhotonUtil::PhoenixHits(TStnEvent* event, TStnElectron* pele,
				Int_t& nHit, Int_t& nHitPhi)
{

  nHit = 0;
  nHitPhi = 0;
  if(pele==NULL) {
    return -10;
  }

  TStnTrack*  trk = TPhotonUtil::PhoenixTrack(event, pele);

  if(trk==NULL) {
    return -9;
  }

  nHit = trk->NSvxHits();
  nHitPhi = trk->NSvxRPhiHits();

  return 0;

}

//_____________________________________________________________________________
TStnTrack*  TPhotonUtil::PhoenixTrack(TStnEvent* event, TStnElectron* pele)
{

  TStnTrack* best=NULL;

  if(pele==NULL) return best;

  int ientry = event->GetCurrentTreeEntry();
  TPhoenixElectronBlock* fPhoenixElectronBlock = (TPhoenixElectronBlock*)
    event->GetDataBlock("Phoenix_Electrons");
  //TStnTrackBlock* fPhoenixTrackBlock = (TStnTrackBlock*)
  //  event->GetDataBlock("Phoenix_Tracking");
  TStnTrackBlock* fPhoenixSiTrackBlock = (TStnTrackBlock*)
    event->GetDataBlock("PROD@PhoenixSI_Tracking");

  //fPhoenixElectronBlock->GetEntry(ientry);
  fPhoenixSiTrackBlock->GetEntry(ientry);

  if( fPhoenixSiTrackBlock->NTracks()<=0) return best;

  int ip=-1;
  for(int i=0; i<fPhoenixElectronBlock->NElectrons(); i++){
    if(pele==fPhoenixElectronBlock->Electron(i)) ip=i;
  }
  if(ip<0) return best;

  // loop over the two charges
  float bestChi2 = 1e10;
  for(int iseedInd=0; iseedInd<2; iseedInd++) {
    // index into PhoenixTrack
    int iseed = pele->PhoenixSeedID(iseedInd);
    int nTrk = fPhoenixElectronBlock->NTracks(ip, iseedInd);
    for(int ind=0; ind<nTrk; ind++) {
      int itrk = fPhoenixElectronBlock->TrackNumber(ip, iseedInd, ind);
      if(itrk>=0 && itrk<fPhoenixSiTrackBlock->NTracks()) {
	TStnTrack* trk = fPhoenixSiTrackBlock->Track(itrk);
	float chi2 = trk->Chi2()/trk->NSvxHits();
	if(best==NULL) {
	  best = trk;
	  bestChi2 = chi2;
	} else {
	  bool better = (trk->Algorithm()==11 && best->Algorithm()==4)
	    || ( (trk->Algorithm()==best->Algorithm()) && 
		 (chi2 < bestChi2 )    );
	  if(better) {
	    best = trk;
	    bestChi2 = chi2;
	  }
	}
      }
    }
  }
  return best;

}

//_____________________________________________________________________________
void TPhotonUtil::PhoBgWeights(Float_t jetEt,
			       Float_t* bgEt, Float_t* bgWht, int mode)
{
  float gaWht[15] = {
    0.003174,0.008981,0.021651,0.044481,0.077871,0.116170,
    0.147681,0.159981,0.147681,0.116170,0.077871,0.044481,
    0.021651,0.008981,0.003174};
  

  float a=-1.52611e-01;
  float b=2.39690e+00;
  float c=4.03686e-01;

  float fakeRate = 1.0e-3*(exp(a*jetEt + b) + c);
  float scale = 1.0;
  if(mode==-1) {
    scale = exp(-5.95958e-04*jetEt + 3.47161e+00) - 3.12553e+01;
    if(scale<0.0) scale = 0.0;
  } else if (mode==1) {
    scale = exp( 2.43173e-02*jetEt - 8.33021e-01) + 6.91892e-01;
  }
  for(int i=0; i<15; i++) {
    float dx = 6.0/15;
    float x = -3.0+i*dx+dx/2.0;
    bgEt[i] = jetEt*(0.937 + x*0.0481);
    bgWht[i] = scale*gaWht[i]*fakeRate;
  }
  return;
}

//____________________________________________________________________________
Int_t TPhotonUtil::cp2_close_pad( Float_t zCPR, Float_t phiGlobal, 
				  Float_t &delx, Float_t &delz){

  const float z_middle[54] = {8.5,8.5,8.5,21.,21.,21.,33.5,33.5,33.5,
                              46.,46.,46.,58.5,58.5,58.5,71.,71.,71.,
                              83.5,83.5,83.5,96.,96.,96.,108.5,108.5,108.5,
                              121.,121.,121.,133.5,133.5,133.5,
                              146.,146.,146.,158.5,158.5,158.5,
                              171.,171.,171.,183.5,183.5,183.5,
                              196.,196.,196.,208.5,208.5,208.5,
                              221.,221.,221.};
 
  const float phi_middle[54] = {3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69,
                                3.31,7.50,11.69,3.31,7.50,11.69}; 
  int ipad;
  float delr, bestdelr=999.;  
 
  // Need to convert phiGlobal from radians to degrees:
  phiGlobal = 180.0 * phiGlobal / M_PI;
  int tmpphi = (int) phiGlobal;
  float local_phi = (tmpphi % 15) + (phiGlobal-tmpphi);
   
  int ibest = -1;
  for(ipad=0; ipad<54; ++ipad) {
     float delzDist = fabs(fabs(zCPR)-z_middle[ipad]);
     float delxDist = fabs((local_phi-phi_middle[ipad])*(12.5/4.13));
     delr = sqrt(delzDist*delzDist + delxDist*delxDist);
     if(delr<bestdelr) {
       delz = delzDist;
       delx = delxDist;
       bestdelr = delr;
       ibest = ipad;
     }
  }
  return ibest;
}