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

#include "TPhotonUtil.hh"
#include "TEmTimeCorrection.hh"
#include "TExtrap.hh"
#include "phoconst.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/TGenpBlock.hh"
#include "Stntuple/obj/TObspBlock.hh"
#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/obj/TStnDBManager.hh"
#include "Stntuple/alg/TStntuple.hh"
#include <map>

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

ClassImp(TPhotonUtil)

TEmTimeCorrection* TPhotonUtil::fEmTimeCorr = NULL;

//_____________________________________________________________________________
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");
  if(!fPhotonBlock) {
    printf("TPhotonUtil::CorrectPhotonEnergy will not work without you registering PhotonBlock\n");
    return 1;
  }

  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::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);
    float emtime=-999.0;
    for(int ilink=0; ilink<nlink; ilink++) {
      int l = links->Index(i,ilink);
      int iphi = l & 0x3F;
      int ieta = (l>>8) & 0x3F;
      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%7d%7d%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");
  if(!fZVertexBlock) {
    printf("TPhotonUtil::SelectVertex will not work without you registering ZVertexBlock\n");
    return false;
  }
  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");
  if(!fJetBlock) {
    printf("TPhotonUtil::MatchJet will not work without you registering MatchJet\n");
    return NULL;
  }

  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");
  if(!fElectronBlock) {
    printf("TPhotonUtil::MatchElectron will not work without you registering ElectronBlock\n");
    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->TowPhi(),pho->TowEta(),
				 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) {
    printf("TPhotonUtil::MatchPhoenixElectrons will not work without you registering Phoenix_Electrons block\n");
    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->TowPhi(),pho->TowEta(),
				 ele->EmClusPhi(),ele->DetEta());
    if(dr<dmin) {
      mele = ele;
      dmin = dr;
    }
  }

  // report only perfect matches
  if(dmin>0.01) return NULL;

  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");
  if(!fPhElectronBlock) {
    printf("TPhotonUtil::MatchPhoenixElectrons will not work without you registering Phoenix_Electrons block\n");
    return NULL;
  }
  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;
    }
  }

  // report only perfect matches
  if(dmin>0.01) return NULL;

  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");
  if(!fPhotonBlock) {
    printf("TPhotonUtil::MatchPhoton will not work without you registering PhotonBlock\n");
    return NULL;
  }
  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->TowPhi(),pho->TowEta(),
				 ele->EmClusPhi(),ele->DetEta());
    if(dr<dmin) {
      mpho = pho;
      dmin = dr;
    }
  }
  return mpho;

}


//_____________________________________________________________________________
Int_t  TPhotonUtil::MatchGenp(TStnEvent* event, TStnPhoton* pho, int dump)
{
  int rc = -1;
  int debug=dump;

  if(pho==NULL) return rc;

  TStnHeaderBlock* fHeaderBlock = (TStnHeaderBlock*)
    event->GetDataBlock("HeaderBlock");
  if(fHeaderBlock->McFlag()==0) return rc;

  int ientry = event->GetCurrentTreeEntry();
  TStnPhotonBlock* fPhotonBlock = (TStnPhotonBlock*)
    event->GetDataBlock("PhotonBlock");
  if(!fPhotonBlock) {
    printf("TPhotonUtil::MatchGenp will not work without you registering PhotonBlock\n");
    return rc;
  }
  fPhotonBlock->GetEntry(ientry);

  TGenpBlock* fGenpBlock = (TGenpBlock*)
    event->GetDataBlock("GenpBlock");
  if(!fGenpBlock) {
    printf("TGenpUtil::MatchGenp will not work without you registering GenpBlock\n");
    return rc;
  }
  fGenpBlock->GetEntry(ientry);

  TObspBlock* fObspBlock = (TObspBlock*)
    event->GetDataBlock("ObspBlock");
  if(!fObspBlock) {
    printf("TGenpUtil::MatchGenp will not work without you registering ObspBlock\n");
    return rc;
  }
  fObspBlock->GetEntry(ientry);

  float pphi = pho->Phi();
  float peta = pho->DetEta();
  if(dump>0) {
    printf("photon %7.3f %7.3f \n",peta,pphi);
  }

  int nobs = fObspBlock->NParticles();
  int nvobs = fObspBlock->NVertices();
  int ngen = fGenpBlock->NParticles();

  // find the highest-pt thing pointing at the photon
  int ipt = -1;
  float etmax = -1.0;
  for(int j=0; j<nobs; j++) {
    TObspParticle* po = fObspBlock->Particle(j);
    float eta = po->Momentum()->PseudoRapidity();
    float phi = po->Momentum()->Phi();
    if(phi<0.0) phi+=TMath::TwoPi();
    float et = po->Momentum()->Pt();

    int iv = po->VertexNumber()-1;
    
    if(dump>5) {
      printf("%3d %5d %8.3f %8.3f %8.3f %3d ",j,po->PdgCode(),et,eta,phi,iv);
    }

    if(iv>=0 && iv<nvobs) {
      TObsvVertex* vtx = fObspBlock->Vertex(iv);
      float zorig = vtx->Z();
      float deta = TStntuple::DetEtaFromEta(eta,zorig);
      if(phi<0.0) phi+=TMath::TwoPi();
      float dr = TStntuple::DeltaR(pphi,peta,phi,deta);
      if(dump>5) printf(" %7.3f %7.3f ",zorig,deta);
      if(dr<0.2 && et>etmax) {
	ipt = j;
	etmax = et;
	if(dump>5) printf(" max\n");
      } else {
	if(dump>5) printf("\n");
      }
    } else {
      if(dump>5) printf("\n");
    }
  }

  if(ipt<0) { // did not find anything pointing to the cluster
    if(dump>0) printf("returning with -1\n");
    return -1;
  }

  // trace the leading particle back to the most useful genp index

  TObspParticle* po = fObspBlock->Particle(ipt);
  int ig0 = po->GenpNumber();
  if(ig0<0 || ig0>ngen) {
    if(dump>0) printf("returning with -2\n");
    return -2;
  }


  bool done = false;
  int imo,ig,ipdg,ipdg0;
  ig = ig0; // needed for first pass

  while(!done) {

    ig0 = ig; // move up the decay chain (no effect in first pass)
    TGenParticle* p0 = fGenpBlock->Particle(ig0);
    if(!p0) {
      if(dump>0) printf("returning with -3\n");
      return -3;
    }
    int ipdg0 = p0->GetPdgCode();

    int imo = p0->GetFirstMother();
    if(imo<0 || imo>ngen) { //no more to check
      if(dump>0) printf("returning in loop with ig0=%d\n",ig0);
      return ig0;
    }
    TGenParticle* p = fGenpBlock->Particle(imo);
    ipdg = p->GetPdgCode();
    ig = imo;

    if(dump>0) printf("%3d %3d %3d %3d \n",ig,ipdg,ig0,ipdg0);

    if(ipdg0==22 && ipdg==22) { // photon with mother photon
      if(dump>0) printf("returning photon mother with ig0=%d\n",ig0);
      return ig0;
    }
    if(abs(ipdg0)==11 && abs(ipdg)==11) { // ele with mother ele
      if(dump>0) printf("returning ele with mother ele with ig0=%d\n",ig0);
      return ig0;
    }

    done = true;
    if(ipdg0==22 && abs(ipdg)==11) { // photon from electron
      if(dump>5) printf("found photon from electron\n");
      done = false;
    }
    if(abs(ipdg0)==11 && ipdg==22) { // electron from photon
      if(dump>5) printf("found electron from photon\n");
      done = false;
    }
    if(abs(ipdg0)==22 && abs(ipdg)==111) { // photon from pi0
      if(dump>5) printf("found photon from pi0\n");
      done = false;
    }
    if(abs(ipdg0)==111 && abs(ipdg)==221) { // pi0 from eta
      if(dump>0) printf("returning with pi0 from eta ig0=%d\n",ig);
      return ig;
    }
    if(abs(ipdg0)==22 && abs(ipdg)==221) { // photon from eta
      if(dump>0) printf("returning with photon from eta ig0=%d\n",ig);
      return ig;
    }
    if(abs(ipdg0)==111 && abs(ipdg)==310) { // pi0 from Ks
      if(dump>0) printf("returning with pi0 from Ks ig0=%d\n",ig);
      return ig;
    }

  }

  if(dump>0) printf("returning end loop with ig0=%d\n",ig0);
  return ig0;

}



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

  if(pele==NULL) return NULL;

  int ientry = event->GetCurrentTreeEntry();
  TPhoenixElectronBlock* fPhoenixElectronBlock = (TPhoenixElectronBlock*)
    event->GetDataBlock("Phoenix_Electrons");
  if(!fPhoenixElectronBlock) {
    printf("TPhotonUtil::PhoenixTrack will not work without you registering Phoenix_Electrons block\n");
    return NULL;
  }
  fPhoenixElectronBlock->GetEntry(ientry);

  TStnTrackBlock* fPhoenixSiTrackBlock = (TStnTrackBlock*)
    event->GetDataBlock("PROD@PhoenixSI_Tracking");
  if(!fPhoenixSiTrackBlock) {
    printf("TPhotonUtil::MatchPhoenixElectrons will not work without you registering PROD@PhoenixSI_Tracking block\n");
    return NULL;
  }
  fPhoenixSiTrackBlock->GetEntry(ientry);

  TStntuple::SetTrackPointersPh(fPhoenixElectronBlock,fPhoenixSiTrackBlock);
  return pele->Track();


}

//_____________________________________________________________________________
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;
}

// Below are CP2 related methods


//____________________________________________________________________________
int   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(zCPR)-z_middle[ipad];
     float delxDist = (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;
}

//===========================================================================
//
//      New Algorithm  :cp2Sum[]={4pad, 1+1pad, 1+2+4pad};
//
//===========================================================================

// more generic
// photon
//_____________________________________________________________________________
void TPhotonUtil::GetCP2(TStnPhoton* pho, int* ibest, int* padE, 
			 int* cp2Sum) {

  // fetch the cpr charge from 4 pads and indices of 4 pads stored in the bank
  const int NPAD=4;
  for(int i=0; i < NPAD; i++)
    {
      ibest[i]= pho->Cpr2Ind(i); 
      padE[i] = pho->Cpr2E(i);
    }
  const int NALG=3;
  for(int i=0; i < NALG; i++)cp2Sum[i]=0;

  // this algorithm returns the same best bad as pho->Cp2Ind(i)
  float cp2z = TPhotonUtil::GetCP2Z(pho);
  float delx = -999;
  float delz = -999;
  int temp= TPhotonUtil::cp2_close_pad(cp2z, pho->Phi(), delx, delz);

  int hitPad = ibest[0];

  int ph4=0, ph12=0, phfinal=0;


  // ----------------------------------------
  // calculate cp2 hits using 4 pad algorithm
  //-----------------------------------------

  for(int i=0; i< NPAD; i++)
    {
      if ( padE[i] > 0)
	ph4 += padE[i];
    }

  //-------------------------------------------
  // calculate cp2 hits using 1+2 pad algorithm
  //-------------------------------------------

  // within the center
  if ( fabs(delx) <=4.25 && fabs(delz) <=4.25) { 
    if ( padE[0] > 0) 
      ph12 += padE[0];
  } 
  // close to one or two sides
  else
    {
      for(int i=0; i< 2; i++)
	{
	  if ( padE[i] > 0) 
	    ph12 += padE[i];
	}
    }


  // --------------------------------------
  // calculate using 1+2+4pad algorithm
  //---------------------------------------

  // within the center, every pad is the same and 
  // special treatment for the corner and boundary pads

  if (  
      (fabs(delx) <=4.25 && fabs(delz) <=4.25) ||
      (hitPad ==  0 && delx <  4.25 && delz <  4.25)   || 
      (hitPad ==  2 && delx > -4.25 && delz <  4.25)   || 
      (hitPad == 51 && delx <  4.25 && delz > -4.25)   || 
      (hitPad == 53 && delx > -4.25 && delz > -4.25)   ||
      (hitPad % 3 ==0 && delx < 0 && fabs(delz) <= 4.25) || 
      (hitPad % 3 ==2 && delx > 0 && fabs(delz) <= 4.25) || 
      (hitPad == 1 && fabs(delx) <= 4.25 && delz < 0)  ||
      (hitPad == 52 && fabs(delx) <= 4.25 && delz > 0)  
     )
    { 

      if ( padE[0] > 0) phfinal += padE[0];

    }

  // close to one side only
  else if( 
  	  (fabs(delx) >  4.25 && fabs(delz)<= 4.25) ||
  	  (fabs(delx) <= 4.25 && fabs(delz)>  4.25)
  	  )
    {
       for(int pad=0; pad < 2; pad ++)
	 if ( padE[pad] > 0)phfinal += padE[pad];
    }

// close to two sides but are boundary pads
  else if( (hitPad % 3 == 0 &&  delx < -4.25 && fabs(delz) > 4.25) || 
           (hitPad % 3 == 2 &&  delx >  4.25 && fabs(delz) > 4.25) || 
           (hitPad >= 0 && hitPad <=2  && fabs(delx) > 4.25 && delz < -4.25) ||
	   (hitPad >=51 && hitPad <=53 && fabs(delx) > 4.25 && delz >  4.25) 

	   )

     {
       for(int pad=0; pad < 2; pad ++)
	 if ( padE[pad] > 0)phfinal += padE[pad];

     }
  else
    { 
      for(int pad=0; pad<4; pad++)
	{
	  if(padE[pad]>0) phfinal += padE[pad];
	}
    }

  cp2Sum[0]=ph4;
  cp2Sum[1]=ph12;
  cp2Sum[2]=phfinal;

  return;
}




//_____________________________________________________________________________
float TPhotonUtil::GetTOFZ(TStnPhoton* pho)
{
  return TPhotonUtil::GetLocalZ(pho, PhoConst::TOF_RADIUS);
}


//_____________________________________________________________________________
float TPhotonUtil::GetCOTZ(TStnPhoton* pho)
{
  return TPhotonUtil::GetLocalZ(pho, PhoConst::COTOUTER_RADIUS);
}





//_____________________________________________________________________________
float TPhotonUtil::GetCP2Z(TStnPhoton* pho)
{
  return TPhotonUtil::GetLocalZ(pho, PhoConst::CP2_RADIUS);
}


//_____________________________________________________________________________
float TPhotonUtil::GetLocalZ(TStnPhoton* pho, float Rdetector)
{
  float localz = (pho->ZCes() - pho->VertexZ())*(Rdetector)/(PhoConst::CES_RADIUS) + pho->VertexZ();
  
  return localz;

}



//_____________________________________________________________________________

float TPhotonUtil::GetCP2X(TStnPhoton* pho)
{
  return (pho->XCes()* PhoConst::CP2_RADIUS/PhoConst::CES_RADIUS);
}




//_____________________________________________________________________________

bool TPhotonUtil::isCP2BadChan(int side, int wedge, int pad)
{

  int wedgeGlobal = side*24 + wedge;
  
  if(wedgeGlobal == 0 && pad == 2)return true;

  if(wedgeGlobal ==31 && pad == 5)return true;

  if(wedgeGlobal ==39 && pad == 3)return true;

  if(wedgeGlobal ==44 && pad == 2)return true;
  

  return false;


}

//_____________________________________________________________________________
bool TPhotonUtil::isCP2Fid(TStnPhoton* pho)
{
  float cp2x = TPhotonUtil::GetCP2X(pho);
  float cp2z = TPhotonUtil::GetCP2Z(pho);

  if(fabs(cp2z) <  2.25) return false;
  if(fabs(cp2z) >  227.25) return false;  
  if(fabs(cp2x) >  18.5)return false;

  float delx=-999, delz=-999;

  int hitPad = TPhotonUtil::cp2_close_pad(cp2z, pho->Phi(), delx, delz);

  int side = cp2z<0? 0 : 1;
  int  wedge = pho->PhiSeedIndex();  

  if(isCP2BadChan(side, wedge, hitPad))return false;

  if(hitPad < 0)return false;
  
  return true;

}


//_____________________________________________________________________________
Int_t TPhotonUtil::LoadEmTimeCorr(char* fn)
{
  if(fEmTimeCorr) delete fEmTimeCorr;
  fEmTimeCorr = 
    new TEmTimeCorrection(TEmTimeCorrection::READM, fn);
  bool rc = fEmTimeCorr->readInFile();

  if(!rc) return 1;
  return 0;
}


//____________________________________________________________________________
Int_t   TPhotonUtil::PhotonTowerTransients(TStnEvent* event) {

  int ientry = event->GetCurrentTreeEntry();

  TStnPhotonBlock* fPhotonBlock = NULL;
  fPhotonBlock = (TStnPhotonBlock*)event->GetDataBlock("PhotonBlock");
  if(!fPhotonBlock) {
    printf("TStntuple::PhotonTowerTransients will not work without registering PhotonBlock\n");
    return 1;
  }
  fPhotonBlock->GetEntry(ientry);
  int np = fPhotonBlock->NPhotons();
  if(np<=0) return 0;

  TCalDataBlock* fCalDataBlock = NULL;
  fCalDataBlock = (TCalDataBlock*)event->GetDataBlock("CalDataBlock");
  if(!fCalDataBlock) {
    printf("TStntuple::PhotonTowerTransients will not work without registering CalDataBlock\n");
    return 1;
  }
  fCalDataBlock->GetEntry(ientry);


  for(int ip = 0; ip<np; ip++) {
    TStnPhoton* pho = fPhotonBlock->Photon(ip);
    TStnLinkBlock* links = fPhotonBlock->TowerLinkList();
    for(int i=0; i<links->NLinks(ip); i++) {
      int key = links->Index(ip,i);
      int phit = key & 0x3F;
      int etat = (key>>8) & 0x3F;
      TCalTower* tow = fCalDataBlock->Tower(etat,phit);
      if(!tow) {
	printf("TStntuple::PhotonTowerTransients: Error - could not find Cal tower p=%2d eta=%2d for tau=%d\n",phit,etat,ip);
	return 1;
      }
      pho->SetCalTower(i,tow);
    } // end loop over towers

  } // end photon loop
  
}


//____________________________________________________________________________
Int_t   TPhotonUtil::PhotonTimeTransients(TStnEvent* event, Int_t useCorr) {

  int ientry = event->GetCurrentTreeEntry();

  TStnPhotonBlock* fPhotonBlock = NULL;
  fPhotonBlock = (TStnPhotonBlock*)event->GetDataBlock("PhotonBlock");
  if(!fPhotonBlock) {
    printf("TStntuple::PhotonTimeTransients will not work without registering PhotonBlock\n");
    return 1;
  }
  fPhotonBlock->GetEntry(ientry);
  int np = fPhotonBlock->NPhotons();
  if(np<=0) return 0;

  TStnEmTimingBlock* fEmTimingBlock = NULL;
  fEmTimingBlock = (TStnEmTimingBlock*)event->GetDataBlock("EmTimingBlock");
  if(!fEmTimingBlock) {
    printf("TStntuple::PhotonTimeTransients will not work without registering EmTimingBlock\n");
    return 1;
  }
  fEmTimingBlock->GetEntry(ientry);

  int run=0;
  if(useCorr!=0) {
    if(fEmTimeCorr==NULL) {
      printf("TStntuple::PhotonTimeTransients run-dependent correction requested,\n but no file loaded\n");
      return 1;
    }
    TStnHeaderBlock* fH = (TStnHeaderBlock*)event->GetDataBlock("HeaderBlock");
    run=fH->RunNumber();
  }


  for(int ip = 0; ip<np; ip++) {
    TStnPhoton* pho = fPhotonBlock->Photon(ip);
    TStnLinkBlock* links = fPhotonBlock->TowerLinkList();
    pho->SetEmTime(-999.0);
    for(int i=0; i<links->NLinks(ip); i++) {
      int key = links->Index(ip,i);
      int phit = key & 0x3F;
      int etat = (key>>8) & 0x3F;
      int ehit = fEmTimingBlock->FindHit(etat,phit);
      if(ehit<0) { // no EM hit
	pho->SetEmTimeInd(i,-1);
	pho->SetEmTimeTow(i,-999.0);
      } else {
	pho->SetEmTimeInd(i,ehit);
	float t = fEmTimingBlock->Time(ehit);
	if(useCorr!=0) {
	  float corr = fEmTimeCorr->getRunCorrection(run);
	  if(corr>-998.0) t = t -corr;
	}
	pho->SetEmTimeTow(i,t);
	if(etat == pho->EtaSeedIndex() && phit == pho->PhiSeedIndex()) {
	  pho->SetEmTime(t);
	}
      }

    } // end loop over towers

  } // end photon loop
  
}