//-----------------------------------------------------------------------------
//  algorithms working on different STNTUPLE objects
//  Nov 26 2000 P.Murat
//-----------------------------------------------------------------------------

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

#include "Stntuple/alg/TStntuple.hh"
#include "Stntuple/alg/TStnMuonID.hh"

#include "Stntuple/data/TStnTypes.hh"

#include "Stntuple/loop/TStnModule.hh"

#include "Stntuple/obj/TStnTau.hh"
#include "Stntuple/obj/TStnTrack.hh"
#include "Stntuple/obj/TStnVertex.hh"
#include "Stntuple/obj/TStnPi0.hh"
#include "Stntuple/obj/TStnTriggerTable.hh"

#include <Stntuple/obj/TStnHeaderBlock.hh>
#include <Stntuple/obj/TStnJetBlock.hh>
#include <Stntuple/obj/TStnElectronBlock.hh>
#include <Stntuple/obj/TPhoenixElectronBlock.hh>
#include <Stntuple/obj/TStnPhotonBlock.hh>
#include <Stntuple/obj/TStnMetBlock.hh>
#include <Stntuple/obj/TStnMuonBlock.hh>
#include <Stntuple/obj/TStnPi0Block.hh>
#include <Stntuple/obj/TStnTauBlock.hh>
#include <Stntuple/obj/TStnTrackBlock.hh>
#include <Stntuple/obj/TStnTriggerBlock.hh>
#include <Stntuple/obj/TCalDataBlock.hh>

#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/obj/TStnErrorLogger.hh"

#include "Stntuple/alg/TStntuple.hh"

#include "Stntuple/obj/TStnDBManager.hh"
#include "Stntuple/data/TXftTrack.hh"

#include "HighLevelObjects/PhotonBackgroundComputer.hh"
#include "CalorGeometry/CalParameters.hh"
#include "Stntuple/obj/TStnTofMatch.hh"

class CdfTrack;

ClassImp(TStntuple)


TStntuple*        TStntuple::fgInstance       = 0;

Int_t             TStntuple::fgRunNumber = 0;
Int_t             TStntuple::fgTriggerBucket;
Float_t           TStntuple::fgEventVertex     = -999.0;

TStnBeamPos*      TStntuple::fgCotBeamPos      = 0;
TStnBeamPos*      TStntuple::fgSvxBeamPos      = 0;
TStnBeamPos*      TStntuple::fgBeamPos         = 0;
TStnBeamPosBlock* TStntuple::fgCotBeamPosBlock = 0;
TStnBeamPosBlock* TStntuple::fgSvxBeamPosBlock = 0;
TStnBeamPosBlock* TStntuple::fgBeamPosBlock    = 0;

//_____________________________________________________________________________
TStntuple::TStntuple() {
  fL1TrigEt = 0;
  fL2TrigEt = 0;
  fL3TrigEt = 0;
}

//_____________________________________________________________________________
TStntuple::~TStntuple() {
}

//_____________________________________________________________________________
TStntuple*  TStntuple::Instance() { 
  static Cleaner cleaner;
  return (fgInstance) ? fgInstance : (fgInstance = new TStntuple());
}

//------------------------------------------------------------------------------
TStntuple::Cleaner::Cleaner() {
}

//------------------------------------------------------------------------------
TStntuple::Cleaner::~Cleaner() {
  if (TStntuple::fgInstance) {
    delete TStntuple::fgInstance;
    TStntuple::fgInstance = 0;
  }
}

//_____________________________________________________________________________
Int_t TStntuple::Init(Int_t RunNumber) {

  int rc = 0;
  //----------------------------------------------------------------------------
  // Did we do this alread?
  //----------------------------------------------------------------------------
  if (fgRunNumber == RunNumber)
    return rc;

  //----------------------------------------------------------------------------
  // Pointer to our data base manager
  //----------------------------------------------------------------------------
  TStnDBManager* dbm = TStnDBManager::Instance();
  // rc = dbm->GetTriggerBucket(RunNumber,fgTriggerBucket);
  // if (rc < 0)
  //   return rc;

//----------------------------------------------------------------------------
// Get access to the beam positions (average is first in the block)
//----------------------------------------------------------------------------
  fgCotBeamPosBlock = (TStnBeamPosBlock*) dbm->GetTable("CotBeamPosBlock");
  fgSvxBeamPosBlock = (TStnBeamPosBlock*) dbm->GetTable("SvxBeamPosBlock");
//---------------------------------------------------------------------------
// Assign average beam position and take care of backward compatibility
//---------------------------------------------------------------------------
  if (fgCotBeamPosBlock->NBeamPos() > 0)
    fgCotBeamPos = fgCotBeamPosBlock->BeamPos(0);
  else
    fgCotBeamPos = (TStnBeamPos*) dbm->GetTable("CotBeamPos");

  if (fgSvxBeamPosBlock->NBeamPos() > 0)
    fgSvxBeamPos = fgSvxBeamPosBlock->BeamPos(0);
  else
    fgSvxBeamPos = (TStnBeamPos*) dbm->GetTable("SvxBeamPos");
//----------------------------------------------------------------------------
// Assign the best beam position version (svx best, cot next best)
//
//  checking X0() is a kludge, however status codes are not always correct - 
//  see run 142227 in root://fcdfdata030.fnal.gov
//     /cdf/scratch/ewk/bhel08-hpte/bhel08-hpte.0068.PR6094.0.s.01
//----------------------------------------------------------------------------
  if ((fgSvxBeamPos->Status() >= 0) && (TMath::Abs(fgSvxBeamPos->X0()) < 98.)){
    fgBeamPosBlock = fgSvxBeamPosBlock;
    fgBeamPos      = fgSvxBeamPos;
  }
  else {                        // SVX beam position not defined -> COT
    fgBeamPosBlock = fgCotBeamPosBlock;
    fgBeamPos      = fgCotBeamPos;
  }

  //----------------------------------------------------------------------------
  // Assign the run number to make sure we do not repeat this unnecessarily
  //----------------------------------------------------------------------------
  fgRunNumber = RunNumber;

  return rc;
}

//_____________________________________________________________________________
Float_t TStntuple::CalIso(TCalDataBlock* calData, 
			  TStnTrack* trk, Float_t cone) {
  if(trk==NULL || calData==NULL) return -999.0;
  float detEta = TStntuple::DetEtaFromCot(trk->Lam0(),trk->Z0());
  return TStntuple::CalIso(calData,detEta,trk->Phi0(),trk->Z0(),cone);
}

//_____________________________________________________________________________
Float_t TStntuple::CalIso(TCalDataBlock* calData, Float_t detEta, 
			  Float_t phi, Float_t z, Float_t cone) {

  int nt = calData->NTowers();
  float sum = 0.0;
  for(int i=0; i<nt; i++) {
    TCalTower* t = calData->Tower(i);
    float phit;
    float locCOG = 0.5;
    if( (t->IEta()>=14 && t->IEta()<=37) ) locCOG = t->EmCOGPhi();

    if ((t->IEta()>7 && t->IEta()<14) ||(t->IEta()>37 && t->IEta()<44)) {
      phit=( 7.5 * (t->IPhi() +locCOG) )/360.0*2*M_PI;
    } else {
      phit=(15.0 * (t->IPhi() +locCOG) )/360.0*2*M_PI;
    }
    float theta = (t->Theta())*M_PI/180.0;
    float etad = -log(tan(theta/2));
    float deleta = etad-detEta;

    float delphi = fabs(phit-phi);
    if(delphi>M_PI) delphi = 2*M_PI-delphi;

    float dr=sqrt(deleta*deleta+delphi*delphi);
    if (dr<cone) {
      float sint = TStntuple::SinFromDetEta(etad, z);
      float etnew = t->Energy()*sint;
      if(etnew>0.1) sum += etnew;
    }
  }
  return sum;
}


//_____________________________________________________________________________
Float_t TStntuple::DetEtaFromCot(Float_t cott, Float_t z) {
  if(cott>1000.0) return 1000.0;
  if(cott<-1000.0) return -1000.0;

  float rces = TLYCES;
  float zces = cott*rces + z;
  float cott0 = zces/rces;
  float detEta = log(sqrt(1.0+cott0*cott0)+cott0);
  if(fabs(detEta)<1.1) return detEta;
  // plug
  float zpes = (cott>0.0? 1.0: -1.0)*(TLZPESL0+TLZPESL1)/2;
  float r = (zpes-z)/cott;
  cott0 = zpes/r;
  detEta = log(sqrt(1.0+cott0*cott0)+cott0);
  return detEta;
}


//_____________________________________________________________________________
Float_t TStntuple::CotFromDetEta(Float_t detEta, Float_t z) {
  if(detEta>1000.0) return 1000.0;
  if(detEta<-1000.0) return -1000.0;
  float cott0 = (exp(detEta)-exp(-detEta))/2.0;
  if(fabs(z)<1.0e-4) return cott0;
  if(fabs(detEta)<1.1) {
    float rces = TLYCES;
    float dz = cott0*rces - z;
    return dz/rces;
  } else {
    float zpes = (detEta>0.0? 1.0: -1.0)*(TLZPESL0+TLZPESL1)/2;
    float r = zpes/cott0;
    float dz = zpes-z;
    return dz/r;
  }
}

//_____________________________________________________________________________
Float_t TStntuple::SinFromDetEta(Float_t detEta, Float_t z) {
  float cott = CotFromDetEta(detEta, z);
  return 1.0/sqrt(1.0+cott*cott);
}


//_____________________________________________________________________________
Float_t TStntuple::CesZPesRFromDetEta(Float_t detEta) {
  if(detEta>1000.0) return 1000.0;
  if(detEta<-1000.0) return -1000.0;
  float cott = (exp(detEta)-exp(-detEta))/2.0;

  if(fabs(detEta)<1.1) {
    return cott*TLYCES;
  } else {
    float zpes = (detEta>0.0? 1.0: -1.0)*(TLZPESL0+TLZPESL1)/2;
    return zpes/cott;
  }
}

//_____________________________________________________________________________
Int_t TStntuple::CenTowerFromZ(Float_t cesz) {
  // boundaries of towers at CES radius
  const float bound[11]={4.22, 24.16, 48.32, 72.48, 96.64,120.8,
			 144.97,169.12,193.28,217.45,245.96};
  float az = fabs(cesz);
  int w = 1;
  while(az>bound[w] && w<10) w++;
  w--;
  if(cesz<0) w=-w-1;
  w += 26;

  return w;

}

//_____________________________________________________________________________
// ieta is 0 to 51
// boundary is 0=tower center (defaut), 1 high boundary, -1 low boundary
Float_t TStntuple::TowerTheta(Int_t ieta, Int_t boundary) {

  if(ieta<0 || ieta>51) return -999.0;

  float thetaMin,thetaMax,theta;
  int ind;
  if(ieta<(TOWER_NETA/2)) {
    ind = ieta;
    thetaMin = (180-TOWER_THETA[ind+1])*DEGRAD;
    thetaMax = (180-TOWER_THETA[ind])*DEGRAD;
  } else {
    ind = TOWER_NETA - 1 - ieta;
    thetaMin = TOWER_THETA[ind]*DEGRAD;
    thetaMax = TOWER_THETA[ind+1]*DEGRAD;
  }
   
  if(boundary==-1) return thetaMin;
  if(boundary== 1) return thetaMax;
  theta = (thetaMax + thetaMin)/2;

  return theta;
}


//_____________________________________________________________________________
// ieta is 0 to 51
Float_t TStntuple::EmTowerRadius(Int_t ieta) {

  if(ieta<0 || ieta>51) return -999.0;
  float theta = TowerTheta(ieta);

  int ind = (ieta<(TOWER_NETA/2) ? ieta : TOWER_NETA - 1 - ieta);

  float rem;
  if (TETTYP[ind]<=2) {           // CEM
    rem = TLYCEM + TLDXDRCEM*(TLRTOT - TLRCOIL/sin(theta));
  } else if (TETTYP[ind]<=7) {      // PEM
    float zem = EmTowerZ(ieta);
    rem = zem*tan(theta);
  }

  return rem;
}


//_____________________________________________________________________________
// ieta is 0 to 51
Float_t TStntuple::EmTowerZ(Int_t ieta) {

  if(ieta<0 || ieta>51) return -999.0;

  int ind = (ieta<(TOWER_NETA/2) ? ieta : TOWER_NETA - 1 - ieta);
  float theta = TowerTheta(ieta);

  float zem;
  if (TETTYP[ind]<=2) {           // CEM
    float rem = TLYCEM + TLDXDRCEM*(TLRTOT - TLRCOIL/sin(theta));
    zem = rem/tan(theta);
  } else if (TETTYP[ind]<=7) {      // PEM
    float fcos = -fabs(cos(theta));
    zem = TLZPEM + TLDXDRPEM*(TLRTOT-TLREND/fcos)*fcos;
    if (TETTYP[ind]==3) {    // one short PEM tower
      if (zem>=TLZTOW10) zem = TLZTOW10;
    }
    if(ieta<(TOWER_NETA/2)) zem = -zem;
  }
  return zem;
}

//_____________________________________________________________________________
Float_t  TStntuple::DeltaR(Float_t phi0, Float_t eta0, 
			   Float_t phi1, Float_t eta1) {
  float dphi = fabs(phi0-phi1);
  if(dphi>M_PI) dphi = 2*M_PI- dphi;
  float deta = eta0-eta1;
  return sqrt(dphi*dphi + deta*deta);
}

//_____________________________________________________________________________
Double_t TStntuple::TrackPt(Double_t C0, Double_t BField) {
  if(C0>-1.e-10 && C0<0.0) C0=-1.e-10;
  if(C0< 1.e-10 && C0>0.0) C0= 1.e-10;
  Double_t pt =  0.5/C0*0.0029979*BField;
  return pt;
}


//-----------------------------------------------------------------------------
Double_t TStntuple::LarrysPt(TStnTrack* Track, Int_t Mode, Int_t OffVer) {

  // run number is required
  // should be filled at BeginRun with TStntuple::Init(runNumber)
  if( fgRunNumber<138425 ) {
    printf("TStntuple::LarrysPt: Error: called without setting run number - failing\n");
    return 0.0;
  }
  double pt, pinv, phi, cot, corrPt;

  pt   = Track->Charge()*Track->Pt();
  pinv   = 1.0/pt;

  if (Mode == 1) {
    //----------------------------------------------------------------------
    // correction for COT-only beam-constrained tracks , do not do anything for the
    // tracks which do not have COT-only fits
    //----------------------------------------------------------------------
    if (Track->Chi2Cot() < 1e5) {
      phi = Track->BcPhi0();
      cot = Track->BcLam0();
      pt   = TrackPt(Track->BcC0());
    } else {
      return pt;
    }
  }  else if (Mode == 2) {  
    //----------------------------------------------------------------------
    // correction for COT+SVX, use unconstrained parameters
    //----------------------------------------------------------------------
    phi = Track->Phi0();
    cot = Track->Lam0();
    pt   = TrackPt(Track->C0());
  }


  if ( OffVer == 5 ) {

    if ( Mode == 1 ) {
      pinv = pinv + 0.00020*sin(phi+3.4);
      pinv = pinv + 0.00022*sin(3*phi+0.9);
      pinv = pinv - (0.000026 + 0.000072*cot - 0.00024*cot*cot);
      pinv = pinv - 0.0002*cot*sin(phi-0.9) - 0.0002*cot*cot*sin(phi-4.1);
    } else {
      pinv = pinv + 0.000195*sin(phi+3.73);
      pinv = pinv + 0.00022*sin(3*phi+0.9);
      pinv = pinv - (0.000039+0.000082*cot - 0.00023*cot*cot);
      pinv = pinv - 0.00023*cot*sin(phi-0.9) - 0.0002*cot*cot*sin(phi-4.1);
    }
    
  }  else if (OffVer == 6) {
    
    if ( fgRunNumber<=213000 ) {

      if ( Mode == 1 ) {
	pinv = pinv - (0.000048 + 0.000085*cot - 0.000185*cot*cot - 0.000185*sin(phi+3.4));
	pinv = pinv + 0.000235*sin(3*phi+0.84); 
	pinv = pinv - 0.00020*cot*sin(phi-0.9) - 0.00024*cot*cot*sin(phi-4.1); 
      } else {
	pinv = pinv - (0.000071 + 0.000095*cot - 0.000195*cot*cot + 0.000260*sin(phi+0.15));
	pinv = pinv + 0.000215*sin(3*phi+0.84); 
	pinv = pinv - 0.00020*cot*sin(phi-0.9) - 0.00024*cot*cot*sin(phi-4.1);
      }    
      
    } else if ( fgRunNumber>213000 && fgRunNumber<226086 ) { 
      
      if ( Mode == 1 ) {
	pinv = pinv - (0.000054 + 0.000085*cot - 0.000185*cot*cot + 0.00029*sin(phi-0.385));
	pinv = pinv + 0.000235*sin(3*phi+0.84); 
	pinv = pinv - 0.00035*cot*sin(phi-1.1) - 0.00040*cot*cot*sin(phi-4.3); 
      } else {
	pinv = pinv - (0.000080 + 0.000100*cot - 0.000185*cot*cot + 0.00033*sin(phi-0.58));
	pinv = pinv + 0.000235*sin(3*phi+0.84); 
	pinv = pinv - 0.00035*cot*sin(phi-1.1) - 0.00040*cot*cot*sin(phi-4.3);
      }    
      
    } else if ( fgRunNumber>=226086 ) { 
      
      if ( Mode == 1 ) {
	pinv = pinv - (0.000040 + 0.000104*cot - 0.000213*cot*cot + 0.00030*sin(phi-0.27));
	pinv = pinv + 0.000235*sin(3*phi+0.84); 
	pinv = pinv - 0.00050*cot*sin(phi-0.9) - 0.00024*cot*cot*sin(phi-4.1); 
      } else {
	pinv = pinv - (0.000085 + 0.000120*cot - 0.000213*cot*cot + 0.00030*sin(phi-0.27));
	pinv = pinv + 0.000235*sin(3*phi+0.84); 
	pinv = pinv - 0.00050*cot*sin(phi-0.9) - 0.00024*cot*cot*sin(phi-4.1);
      }
    }
  } else {
    printf("TStntuple::LarrysPt: Error: don't know offline version %d\n",OffVer);
  }
  corrPt = (fabs(pinv)>1.e-10 ? 1/pinv : pt );
  return corrPt;

}


//_____________________________________________________________________________
Int_t TStntuple::GoodHptElectron(const TStnElectron* Ele) {
  // no isolation requirements, just make sure electron has a track

  if ( (Ele->IsCentral()           ) &&
       (Ele->TrackNumber()   >  0  ) &&
       (Ele->Et()            > 20. ) &&
       (Ele->TrackPt()       > 10. ) && 
       (Ele->D0()            <  1. ) && 
       (Ele->HadEm()         <  0.1)    ) {
    return 1;
  }
  else return 0;
}


//_____________________________________________________________________________
Int_t TStntuple::GoodHptMuon(const TStnMuon* Mu) {
  // 

  int mu_is_central = Mu->Detector() & 0x3;

  if ( (mu_is_central              ) &&
       (Mu->Momentum()->Pt()  > 20.) &&
       (Mu->EmEnergy()        <  2.) &&
       (Mu->HadEnergy()       <  6.) &&
       (Mu->D0()              <  1.)    ) {
    return 1;
  }
  else return 0;
}


//_____________________________________________________________________________
Int_t TStntuple::GoodHptTau(const TStnTau* Tau) {
  // neet to check D0 of the seed track, not caching it so far

  if ( (Tau->IsCentral()            ) &&
       (Tau->Momentum()->Pt() > 20. ) &&
       (Tau->SeedTrackPt()    > 4.5 ) &&
       (Tau->NTracks1030()   == 0   ) &&
       (Tau->Mass()           < 5.0 ) && 
       (Tau->TrackMass()      < 1.8 )    ) {
    return 1;
  }
  else return 0;
}


//_____________________________________________________________________________
Double_t TStntuple::LeptonZ0(TStnLepton* Lepton, TStnTrackBlock* TrackBlock) 
{
  // returns Z0 of a [supposedly - good and high-Pt!] lepton with the 
  // primary goal to use it for recalculating the MET

  int              it;
  TStnTrack*       trk;
  TStnErrorLogger* el;

  if      (Lepton->Type() == TStnLepton::kElectron) {
//-----------------------------------------------------------------------------
// electron: find its track and take Z0
//-----------------------------------------------------------------------------
    TStnElectron* ele = (TStnElectron*) Lepton;
    it  = ele->TrackNumber();
    if (it >= 0) {
      trk = TrackBlock->Track(it);
      return trk->Z0();
    }
    else {
      el = TrackBlock->GetEvent()->GetErrorLogger();
      el->Report(-109,Form("TStntuple::LeptonZ0: electron with it=-1"));
      return 0;
    }
  }
  else if (Lepton->Type() == TStnLepton::kMuon    ) {
//-----------------------------------------------------------------------------
// muon: find its track and take Z0
//-----------------------------------------------------------------------------
    TStnMuon* muo = (TStnMuon*) Lepton;
    it  = muo->TrackNumber();
    if (it >= 0) {
      trk = TrackBlock->Track(it);
      return trk->Z0();
    }
    else {
      el = TrackBlock->GetEvent()->GetErrorLogger();
      el->Report(-110,Form("TStntuple::LeptonZ0: muon with it=-1"));
      return 0;
    }
  }
  else if (Lepton->Type() == TStnLepton::kTau     ) {
//-----------------------------------------------------------------------------
// tau: its Z0 is Z0 of the seed track, use it
//-----------------------------------------------------------------------------
    TStnTau* tau = (TStnTau*) Lepton;
    return tau->Zv();
  }

  return TStnDataBlock::kUndefined;
}

//_____________________________________________________________________________
Int_t TStntuple::CorrectMetForMuons(TStnEvent*      Event, 
				    TVector2&       CorrMet, 
				    const char*     MuonBranch,
				    Int_t           Mode) 
{
  // corrects RAW missing Et for muons:
  // Mode = 0: only for high-Pt 'CMUO'
  // Mode = 1: plus for high-Pt 'CMIO'
  // Mode = 2: plus for low-Pt  'CMIO'
  // only mode = 0 is implemented for the moment, 
  // don't do anything for Ht either
  // input: MET branch and the used muon branch are supposed to be read in

  float       phi;
  TStnMuonID  muon_id;

  if (Mode != 0) {
    printf (" <TStntuple::CorrectMetForMuons>: only Mode=0 is implemented\n");
    return -1;
  }
//-----------------------------------------------------------------------------
// start from raw missing Et (fMet[0])
//-----------------------------------------------------------------------------
  TStnMuonBlock* muon_block = (TStnMuonBlock*) Event->GetDataBlock(MuonBranch);
  TStnMetBlock*  met_block  = (TStnMetBlock* )  Event->GetDataBlock("MetBlock");
  
  phi  =  met_block->fMetphi[0];
  TVector2 met(met_block->Met(0)*cos(phi),met_block->Met(0)*sin(phi));

  CorrectMetForMuons(muon_block,&muon_id,&met,&CorrMet,Mode);

  return 0;
}

//_____________________________________________________________________________
Int_t TStntuple::CorrectMetForMuons(TStnMuonBlock*  MuonBlock, 
				    TStnMuonID*     MuonID,
				    TVector2*       Met,
				    TVector2*       CorrMet, 
				    Int_t           Mode) 
{
  // corrects RAW missing Et for muons:
  // Mode = 0: only for high-Pt 'CMUO'
  // Mode = 1: plus for high-Pt 'CMIO'
  // Mode = 2: plus for low-Pt  'CMIO'
  // only mode = 0 is implemented for the moment, 
  // don't do anything for Ht either
  // correct only for muons passing ID cuts defined by MuonID

  double  metx, mety, dmetx, dmety, pt, ecal, p;
  int     nmuons, id_word;

  metx   = Met->X();
  mety   = Met->Y();
				// loop over the muons, note that track
				// pointer should be set at run-time by the user
  dmetx  = 0.;
  dmety  = 0.;

  nmuons = MuonBlock->NMuons();
  for (int i=0; i<nmuons; i++) {
    TStnMuon* mu   = MuonBlock->Muon(i);
    TStnTrack* trk = mu->Track();
    if      (Mode == 0) id_word = MuonID->IDWord(mu);
    else if (Mode == 1) id_word = MuonID->LooseIDWord(mu);
    else {
      id_word = -1;
      printf (" <TStntuple::CorrectMetForMuons>: only Mode=0 is implemented\n");
    }

    if (id_word == 0) {
				// muon passes today's "good muon" cuts

      //      phi    = mu->Momentum()->Phi();
      pt     = TMath::Abs(fgInstance->LarrysPt(trk,1));
      ecal   = mu->CalEnergy();
      p      = mu->Momentum()->P()*pt/mu->Momentum()->Pt();
      dmetx += mu->Momentum()->Px()*(1-ecal/p);
      dmety += mu->Momentum()->Py()*(1-ecal/p);
    }
  }

  CorrMet->Set(metx-dmetx,mety-dmety);
  return 0;
}



//_____________________________________________________________________________
Int_t TStntuple::CorrectMet(TStnEvent*      event, 
			    TVector2*       corrMet, 
			    Double_t*       corrHt) 
{
  // calculates corrected missing Et ant total energy for given `event'. 
  // Calculation starts from MET corrected for high-Pt
  // CMUO, high-Pt CMIO amd minimum-ionizing CMUO with Pt>10 GeV
  // calculates MET accounting for corrected jet energies and also for
  // CEM/PEM corrections
  // this routine returns different answer from stntuple/get_metc.F because
  // when taking into account jet corrections it uses only the correction
  // for non-linear response:
  // - when calculating MET one shouldn't correct jets for out-of-cone energy 
  //   which is already included into MET - maybe with wrong scale factor, 
  //   though.. 
  // - also one shouldn't subtract the energy of underlying event from the jet
  //   one assumes that the underlying event in average has zero MET...

				// don't correct MET for jets with Et < 8 
				// in cone 0.4
  float const kEtMin = 8.;
  float phi, metx, mety, dmetx, dmety, dx, dy;

				// start from raw MET 
				// and high Pt CMIO's (is this correct ?

  TStnMetBlock* met = (TStnMetBlock*) event->GetDataBlock("MetBlock");

  phi  =  met->fMetphi[0];
  metx = -met->fMet[0]*cos(phi);
  mety = -met->fMet[0]*sin(phi);

  corrMet->Set(0.,0.);
  *corrHt = 0.;

  dmetx  = 0.;
  dmety  = 0.;

  TStnElectronBlock* edata;
  TStnPhotonBlock*   pdata;
  TStnJetBlock*      jdata;

  edata = (TStnElectronBlock*) event->GetDataBlock("ElectronBlock");
  pdata = (TStnPhotonBlock*)   event->GetDataBlock("PhotonBlock");
  jdata = (TStnJetBlock*)      event->GetDataBlock("JetBlock");

  TStnJet*        jet;
  TStnElectron*   e;
  TStnPhoton*     photon;

  for (int i=0; i<jdata->NJets(); i++) {
    jet = jdata->Jet(i);

    if (jet->Et() > kEtMin) {
      if (jet->ObjectType() == TStnTypes::kJet) {
 	dx = jet->Momentum()->Px()*(jet->Corfd()-1.);
 	dy = jet->Momentum()->Py()*(jet->Corfd()-1.);
 	if (fabs(jet->DetEta()) < 2.) {
 	  *corrHt += jet->Et()*jet->Corfm();
	}
      }
      else if (jet->ObjectType() == TStnTypes::kElectron) {
//-----------------------------------------------------------------------------
//  electrons
//-----------------------------------------------------------------------------
 	e        = edata->Electron(jet->ObjectNumber());
 	dx       = e->Momentum()->Px()*(1-1./e->Etcor());
 	dy       = e->Momentum()->Py()*(1-1./e->Etcor());
 	*corrHt += e->Et()*e->Etcor();
       }
       else if (jet->ObjectType() == TStnTypes::kPhoton) {
 	photon   = pdata->Photon(jet->ObjectNumber());
 	dx       = photon->Momentum()->Px()*(1-photon->Et()/photon->Etc());
 	dy       = photon->Momentum()->Py()*(1-photon->Et()/photon->Etc());
 	*corrHt += photon->Etc();
       }
      else {
					// this shouldn't happen

	printf(" << Stntuple::CorrectMissingEt >> : shouldn't happen!\n");
	dx = 0; 
	dy = 0;
      }
      dmetx += dx;
      dmety += dy;
    }
  }
					// corrected missing Et is here

//   std::cout << " dmetx, dmety = " << dmetx << " " << dmety << std::endl;

  corrMet->Set(-metx-dmetx,-mety-dmety);
  *corrHt += corrMet->Mod();

  return 0;
}

//_____________________________________________________________________________
Int_t TStntuple::InitJets(TStnEvent* event) {
  // remove electrons, photons etc  from the list of event cone 0.4 jets
  // assuming we know the branches with the jets, electrons and photons
  // electrons/photons and such are supposed to have ID word defined !
  // (make sure this is the case!)

				// loop over all the jets in cone 0.4
  float const kConeSize = 0.4;

  Float_t drmin;

  float   dr;

  int nmatches;

  TStnJet       *jet, *best_jet;
  TStnElectron  *e;
  TStnPhoton    *photon;

  int   kCemLoose    = 3981; // 0x0f33
  float kLeptonPtCut = 25.;

  TStnElectronBlock* edata;
  TStnPhotonBlock*   pdata;
  TStnJetBlock*      jdata;

  edata = (TStnElectronBlock*) event->GetDataBlock("ElectronBlock");
  pdata = (TStnPhotonBlock*)   event->GetDataBlock("PhotonBlock");
  jdata = (TStnJetBlock*)      event->GetDataBlock("JetBlock");

				// do not rely on the previous info

  for (int j=0; j<jdata->NJets(); j++) {
    jdata->Jet(j)->SetObjectType(TStnTypes::kJet);
    jdata->Jet(j)->SetObjectNumber(j);
  }
				// loop over the electrons and remove the
				// corresponding jets from the jet list
  int  status;

  for (int i=0; i<edata->fNElectrons; i++) {
    e = edata->Electron(i);
				// check that it really is an electron, not
				// just EM cluster
				// unlike Toback, don't check conversion status:
				// if we have high-Pt electron from conversion,
				// correct for it
    status = 0;
    if      ( (e->IDWord()          == 0           ) &&
	      (e->Etcor()            > kLeptonPtCut)    ) {
      status = 1;
    }
    else if ( (e->DetectorCode()     == 0           ) && 
	      (e->Etcor()             > kLeptonPtCut) && 
	      (e->IDWord()&kCemLoose == 0           )    ) {
      status = 2;
    }

    if (status > 0) {
      //      double ele_eta = e->fEveta;
      double ele_phi = e->Momentum()->Phi();
      drmin   = kConeSize;
				// loop over the jets in cone 0.4
      best_jet = NULL;
      for (int j=0; j<jdata->NJets(); j++) {
	jet = jdata->Jet(j);
				// here I assume (and this is not something given) 
				// that jet and electron have the same common
				// vertex, so subtracting eta's makes sense

 	dr      = jet->Momentum()->DrEtaPhi(*e->Momentum());
 	if (dr < drmin) {
 	  best_jet  = jet;
 	  drmin     = dr;
	  nmatches   += 1;
	  if (nmatches > 1) {
	  }
 	}
      }
				// see if there is a jet-electron match
      if      (nmatches == 0) {
				// something is basically wrong
      }
      else if (nmatches == 1) {
				// this jet is an electron
	best_jet->SetObjectNumber(i);
	best_jet->SetObjectType(TStnTypes::kElectron);
      }
    }
  }

//------------------------------------------------------------------------------
// photons: whatever EM cluster we have in the calorimeter, we don't want to 
// correct it for HAD response, so
// loop over ALL the photons and for each photon find the closest jet 
// (among those which do not have an electron match)
//------------------------------------------------------------------------------
  for (int i=0; i<pdata->fNPhotons; i++) {

    photon   = pdata->Photon(i);
    drmin    = kConeSize;
    nmatches = 0;
    best_jet = NULL;
				// check if it is a good photon, status is only
				// checked for being non-zero
    status   = 0;
    if ( (photon->StatusCode() > 3025.) && (photon->StatusCode() < 4000.) ) {

				// Dave's good plug photon
      status = 1;
    }
    else if ((photon->StatusCode() > 5025.) && (photon->StatusCode() < 6000.)){

				// Jeff`s good photon with cone isolation
      status = 2;
    }
    else if ((photon->StatusCode() > 6025.) && (photon->StatusCode() < 7000.)){

				// Jeff`s good photon with box and cone 
				// isolation
      status = 3;
    }

    if (status > 0) {
				// loop over the jets in cone 0.4

      for (int j=0; j<jdata->NJets(); j++) {
	jet = jdata->Jet(j);
	if (jet->Type() == TStnTypes::kJet) {

				// here I asume (and this is not something given) 
				// that jet and photon have the same common
				// vertex, so eta's could be subtracted

	  dr      = photon->Momentum()->DrEtaPhi(*jet->Momentum());
	  if (dr < drmin) {
	    best_jet  = jet;
	    drmin     = dr;
	    nmatches   += 1;
	    if (nmatches > 1) {
	      // std::cout << "e-moe! nmatches = " << nmatches << std::endl;
	    }
	  }
	}
      }
    }

    if (nmatches == 1) {
				// best_jet is a photon

      best_jet->SetObjectNumber(i);
      best_jet->SetObjectType(TStnTypes::kPhoton);
    }
				// see if there is a jet-photon match
    else if (nmatches > 1) {
				// something is basically wrong - more than
				// 1 jet matches to 25 GeV photon

      printf(" init_jets: something is wrong: # photon-jet matches > 1\n");
    }
  }
  return 0;
}


//_____________________________________________________________________________
Double_t  TStntuple::TrkIso(TStnTrack*    track, 
			    TClonesArray* trackList,
			    Double_t      cone,
			    Double_t      dz,
			    Double_t      ptMin) 
{
  // calculate `track' isolation: (sumPt of other tracks)/trackPt
  // from `trackList' within `cone' in eta-phi from the `track'. 
  // Tracks with Pt > `ptMin' are required to have |Ztrk-Z0|<`dz'
  // where Z0 is z-coordinate of the `track' in the point of closest
  // approach to the origin, thus cutting against multiple interactions


  double  phi0, eta0, dphi, deta, sum_pt, z0;
  double  pt, dr, trk_iso;
  int     ntrk, ok;

  sum_pt = 0;

  phi0   = track->Momentum()->Phi();
  eta0   = track->Momentum()->Eta();
  z0     = track->Z0();
  ntrk   = trackList->GetEntriesFast();

  for (int i=0; i<ntrk; i++) {
    TStnTrack* ti = (TStnTrack*) trackList->UncheckedAt(i);
    if (ti != track) {
      ok = 1;
      TLorentzVector* mi = ti->Momentum();
      pt  = mi->Pt();
      if (pt > ptMin) {
				// high-Pt track, check Z-coordinate

	if (fabs(ti->Z0()-z0) >= dz) {
	  ok = 0;
	}
      }
      if (ok) {
	dphi = TVector2::Phi_mpi_pi(phi0-ti->Momentum()->Phi());
	deta = eta0-ti->Momentum()->Eta();
	dr   = sqrt(dphi*dphi+deta*deta);
	if (dr < cone) {
	  sum_pt += pt;
	}
      }
    }
  }
				// finally calculate isolation

  trk_iso = sum_pt/track->Momentum()->Pt();
  return trk_iso;
}


Int_t TStntuple::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.NClusterWords();
    for(int i=0; i<n; i++) {
      TTl2dCluster& clus =  *(tl2d.Cluster(i));
      if( clus.Pass()==pass) {
	dphi = fabs(phi - (clus.Phi()+0.5)*2.0*M_PI/24.0);
	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.NIsoClusterWords();
    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  TStntuple::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;
}

//_____________________________________________________________________________
TStnTrack*  TStntuple::SeedTrack(const TStnTau*  Tau, 
				 TStnTauBlock*   TauBlock,
				 TStnTrackBlock* TrackBlock) 
{
  // returns seed track for a given tau

  int          iseed;
  TStnTrack*   track;

  int i  = Tau->Number();
  iseed  = TauBlock->TrackNumber(i,0);
  if (iseed < 0) {
    printf("<< TStntuple::SeedTrack ERROR: a tau w/o a seed track\n");
    track = 0;
  }
  else {
    track = TrackBlock->Track(iseed);
  }
  return track;
}

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

//_____________________________________________________________________________
Double_t  TStntuple::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;

}



//------------------------------------------------------------------------------
//  this is a primitive algorithm, reconstructing primary vertices in Z
//  using only COT tracks (hmmm - am I using this assumption? )
//------------------------------------------------------------------------------
int TStntuple::VprimReco(TClonesArray* trackList, 
			 TClonesArray* vertexList,
			 double        maxChi2,
			 int           minTracks,
			 int           searchMode) {
  TStnTrack*     track;
  TStnVertex*    vertex;
  TStnVertex*    closest_vertex;
  double         chi2;
  double         chi2_min;
  int            vertex_number;
  int            ntracks, nv;

  ntracks = trackList->GetEntriesFast();

  vertexList->Clear();
  nv   = 0;
  for (int i=0; i<ntracks; i++) {
    track = (TStnTrack*) trackList->UncheckedAt(i);
    if (track->Is3D()) {
                                        // use only 3D tracks
                                        // loop over the existing vertices and
      chi2_min = DBL_MAX;
      for (int iv=0; iv<nv; iv++) {
	vertex = (TStnVertex*) vertexList->UncheckedAt(iv);

                                        // calculate track chi**2 in the vertex
	if      (searchMode == 0) {
                                        // during the 1st 3D reconstruction pass
                                        // look only at the distance between
                                        // track and vertex

	  chi2 = fabs(track->Z0()-vertex->Z());
	}
	else if (searchMode == 1) {
                                        // at the 2nd 3D reconstruction pass 
                                        // take errors into account

	  chi2 = TStntuple::TrackChi2RZ(track,vertex);
	}
	if (chi2 < chi2_min) {
                                        // new closest vertex for the track
	  closest_vertex  = vertex;
	  chi2_min        = chi2;
	}
      }

      if (track->Pt() > 0.2) {
//------------------------------------------------------------------------------
//  to define vertices use only tracks with Pt > 0.2
//------------------------------------------------------------------------------
	if (chi2_min < maxChi2) {

                                        // chi**2 of the track in the closest
                                        // vertex is small enough, add track to
                                        // this vertex
	  closest_vertex->Add(track);
	}
	else {
                                        // track is too far from the vertex,
                                        // closest to it, create new vertex

	  closest_vertex = new ((*vertexList)[nv]) TStnVertex();
	  closest_vertex->SetNumber(nv);
	  closest_vertex->Init(track);
	  nv++;
	}
      }
                                        // end of loop over the tracks
                                        // last step: for each track define 
                                        // number of the corresponding vertex in
                                        // the vertex list
      if (closest_vertex) {
	vertex_number = closest_vertex->Number();
      }
      else {
	vertex_number = -1;
      }
      track->SetVertexNumber(vertex_number);
    }
  }

//------------------------------------------------------------------------------
//  vertices are found, leave only those, which have N(tracks) > N(min) = 3...
//  order found vertices in descending ntracks
//------------------------------------------------------------------------------

  if (vertexList->GetEntriesFast()) {
    vertexList->Sort();
    int n_removed = 0;
    for (int iv=0; iv<nv; iv++) {
      TStnVertex* v = (TStnVertex*) vertexList->UncheckedAt(iv);
      if (v->NTracks() < minTracks) {
	vertexList->RemoveAt(iv);
	n_removed++;
      }
    }
    if (n_removed > 0) vertexList->Compress();
  }
  
//------------------------------------------------------------------------------
//  second step : calculate vertex chi**2
//------------------------------------------------------------------------------
//   double  track_chi2, dz;
//   TStnTrack* t;

//   for (int iv=0; iv<nv; iv++) {
//     vertex = (TStnVertex*) vertexList->UncheckedAt(iv);
//     chi2   = 0.;
//     ntr    = vertex->NTracks();
//     for (it=0; it<ntr; it++) {
//       itt  = vertex->TrackNumber(it);
//       t    = (TStnTrack*) trackList->UncheckedAt(itt);

//                                         // all the tracks in the list are 3D
//                                         // and all have Pt > 0.2

//       dz         = t->ZAtGivenS(0.)-vertex->Z();
//       track_chi2 = dz*dz/t->sigmaZ2AtGivenS(0.);
//       t->SetVertexChi2RZ(track_chi2);
//       chi2      += track_chi2;
//     }
//     chi2         = chi2/(vertex->NTracks()-0.9999999);
//     vertex->SetChi2(chi2);
//   } 
                                        // end of procedure
  return 0;
}

//_____________________________________________________________________________
double TStntuple::TrackChi2RZ(TStnTrack* track, TStnVertex* v) {
  printf(" Dummy TStntuple::TrackChi2RZ called ! \n");
  return -1;
}


//_____________________________________________________________________________
double TStntuple::CorrectedD0(TStnTrack* Track) {
  // returns track impact parameter corrected for the beam offset
  // assume that the constants for the corresponding run have already been 
  // prefetched from the DB
  // if beam position is not defined, then do not correctD0
  if(Track==NULL) return -999.0;
  Double_t corrected_d0 = Track->D0();
  if (fgBeamPos->Status() != -1) {
    double beamX, beamY;
    if (Track->NSvxHits() <= 0 && fgCotBeamPos->Status() >= 0) {
      beamX = fgCotBeamPos->X0() + Track->Z0()*fgCotBeamPos->DxDz();
      beamY = fgCotBeamPos->Y0() + Track->Z0()*fgCotBeamPos->DyDz();
    } else {
      beamX = fgBeamPos->X0() + Track->Z0()*fgBeamPos->DxDz();
      beamY = fgBeamPos->Y0() + Track->Z0()*fgBeamPos->DyDz();
    }
    corrected_d0 -= -beamX*sin(Track->Phi0()) + beamY*cos(Track->Phi0());
  }

  return corrected_d0;
}


//_____________________________________________________________________________
double TStntuple::CorrectedD0(TStnMuon* Muon) {
  // returns track impact parameter corrected for the beam offset
  // assume that the constants for the corresponding run have already been 
  // prefetched from the DB
   TStnTrack* Track = Muon->Track();
   if (Track != NULL) {
     return CorrectedD0(Track);
   } else {
     Double_t corrected_d0 = Muon->D0();
     const TLorentzVector* mom = Muon->Momentum();
     double beamX= fgBeamPos->X0() + Muon->Z0()*fgBeamPos->DxDz();
     double beamY = fgBeamPos->Y0() + Muon->Z0()*fgBeamPos->DyDz();
     corrected_d0 -= -beamX*sin(mom->Phi()) + beamY*cos(mom->Phi());
     return corrected_d0;
   }
}


//_____________________________________________________________________________
double TStntuple::CorrectedD0(TStnElectron* Electron) {
  // returns impact parameter of electron track corrected for the beam offset
   TStnTrack* Track = Electron->Track();
   if (Track != NULL) {
     return CorrectedD0(Track);
   } else {
     Double_t corrected_d0 = Electron->D0();
     const TLorentzVector* mom = Electron->Momentum();
     double beamX = fgBeamPos->X0() + Electron->Z0()*fgBeamPos->DxDz();
     double beamY = fgBeamPos->Y0() + Electron->Z0()*fgBeamPos->DyDz();
     corrected_d0 -= -beamX*sin(mom->Phi()) + beamY*cos(mom->Phi());
     return corrected_d0;
   }

}

//_____________________________________________________________________________
Double_t TStntuple::XftPt(Int_t PtBin, Int_t IsolBit, Int_t ShortBit) {
  // 
  return TXftTrack::TrackPt(PtBin,IsolBit,ShortBit);
}

//_____________________________________________________________________________
void TStntuple::DefinePi0Momentum(TStnPi0* Pi0, double Z0) {
  // 

  TLorentzVector* mom = (TLorentzVector*) Pi0->Momentum();

  double p      = mom->P();
  double pt     = mom->Pt();
  double cosphi = mom->Px()/pt;
  double sinphi = mom->Py()/pt;

  Pi0->SetZ0(Z0);
//-----------------------------------------------------------------------------
// check Z-coordinate - correct for a bug of filling routine... 
// (which is already fixed)
//-----------------------------------------------------------------------------
  if (Pi0->IEta() <= 25) {
    if (Pi0->ZCes() > 0) Pi0->SetCesCoordinates(Pi0->XCes(),-Pi0->ZCes());
  }
  double dz     = Pi0->ZCes()-Pi0->Z0();
  double sqr    = sqrt(dz*dz+TLYCES*TLYCES);
  double sinth  = TLYCES/sqr;
  double costh  = dz/sqr;
				// 4th component doesn't change
  mom->SetX(p*cosphi*sinth);
  mom->SetY(p*sinphi*sinth);
  mom->SetZ(p*costh);
}

//_____________________________________________________________________________
Float_t TStntuple::CmwInt(Int_t id, Float_t E) {
  /*
   The energy-dependent number of interaction lengths in a central wedge, 
   for different hadron species.  The probability of the meson passing 
   through the wedge without interacting is P = exp(-cmwint/sin(theta)),
   id = 0, 1, 2, 3, 4 ==>> pi, K+, K-, p, p-bar respectively.
   From CMWINT.FOR, David A. Smith (University of Illinois), 14 June 1988
   rlc: for E<2.3 now return value for 2.3
  */

// On average, minimum ionizing particles will loose 1.3 GeV
// as they traverse a central wedge.
   float dE = 1.3;

   float par[5][6] = {
8.492184,-4.713234,2.967638,-0.9973104,0.1679028,-1.0990458E-02,//pi
1.717775,2.161341,-0.5330036,-9.2531711E-02,5.5526044E-02,-5.6797615E-03,//K+
6.676070,-1.831558,1.326447,-0.6193235,0.1336822,-1.0318661E-02,//K-
3.813816,1.441914,0.2904211,-0.4477635,0.1200408,-9.9918647E-03,//p
12.43994,-3.136523,0.7221342,-0.1178043,1.2072344E-02,-4.0744245E-04};//pbar

   if(id<0 || id>4) return -1.0;
   if(E<=1.3 || E>150.0) return -1.0;
   if(E<2.3) E=2.3;

// The fit is a 5th order (6 parameters) polynomial in log(E).
// Integrate the fit from E to E-dE and divide by dE to get
// the effective number of interaction lengths.

   float x = log(E);
   float xf = log(E-dE);
   float sum1 = 0.;
   float sum2 = 0.;
   
   for(int i=0; i<6; i++) {
     int ii = i+1;
     sum1 += par[id][i]*pow(x,ii)/float(ii);
     sum2 += par[id][i]*pow(xf,ii)/float(ii);
   }
   return (sum2 - sum1)/(xf-x);
}


//_____________________________________________________________________________
TStnTrack* TStntuple::TauSeedTrack(TStnTau*        Tau, 
				   TStnTauBlock*   TauBlock,
				   TStnTrackBlock* TrackBlock) {
  //
  
  int          iseed;
  TStnTrack*   track(NULL);

  int i  = Tau->Number();
  iseed  = TauBlock->TrackNumber(i,0);
  if (iseed >= 0) {
    track = TrackBlock->Track(iseed);
  }
  return track;
}



//_____________________________________________________________________________
TStnTrack* TStntuple::ElectronTrack(TStnElectron*      Ele, 
				    TStnTrackBlock*    TrackBlock) {
  //
  TStnTrack*   track(NULL);

  int it  = Ele->TrackNumber();
  if (it >= 0) {
    track = TrackBlock->Track(it);
  }
  return track;
}



//_____________________________________________________________________________
Int_t TStntuple::CheckL3TriggerPath(TStnTriggerBlock* TriggerBlock, 
				    TObjArray*        List) {
  //
  int passed = 1;
  int nt     = List->GetEntriesFast();

  if (nt > 0) {
					// check requested bits, require
					// at least one of them to pass
    passed = 0;
    TTl3d* tl3d = TriggerBlock->Tl3d();

    for (int i=0; i<nt; i++) {  
      TStnTrigger* trig = (TStnTrigger*) List->UncheckedAt(i);
      if (tl3d->PathPassed(trig->Bit()) != 0) passed = 1;
      if (passed) goto PASSED;
    }
  }
 PASSED:;
  return passed;
}


//_____________________________________________________________________________
Int_t TStntuple::PrintListOfPassedTriggers(TStnTriggerBlock* TriggerBlock, 
					   Int_t  Level, Int_t printAll) 
{
  //
  TStnDBManager*   dbm = TStnDBManager::Instance();
  TStnTriggerTable* tt = (TStnTriggerTable*) dbm->GetTable("TriggerTable");

  int banner_printed = 0;

  if ((Level == 1) || (Level == 0)) {
    TTl2dL1Decision* l1_decision = TriggerBlock->Tl2d()->L1Decision();
    for (int i=0; i<64; i++) {
      int passed = l1_decision->L1DecisionBit(i);
      if (passed || printAll>0) {
	for (int i1=0; i1<tt->NL1Triggers(); i1++) {
	  const TStnTrigger* trig = tt->GetTrigger(1,i1);
	  if (trig && (trig->Bit() == i)) {
	    if (! banner_printed) {
	      trig->Print("banner");
	      banner_printed = 1;
	    }
	    if(printAll>0) {
	      if(passed) {
		printf(" on ");
	      } else {
		printf("off ");
	      }
	    }
	    trig->Print("data");
	    break;
	  }
	}
      }
    }
  }

  if ((Level == 2) || (Level == 0)) {
//-----------------------------------------------------------------------------
// L2 : look at unprescaled bits
//-----------------------------------------------------------------------------
    TTl2d*           tl2d = TriggerBlock->Tl2d();
    TTl2dL2Decision* l2_decision = tl2d->L2Decision();
    int              nl2bits     = tl2d->NL2TriggerBitWords()*32;
    for (int i=0; i<nl2bits; i++) {
      int passed = l2_decision->PrescaledBit(i);
      if (passed || printAll>0) {
	for (int i2=0; i2<tt->NL2Triggers(); i2++) {
	  const TStnTrigger* trig = tt->GetTrigger(2,i2);
	  if (trig && (trig->Bit() == i)) {
	    if (! banner_printed) {
	      trig->Print("banner");
	      banner_printed = 1;
	    }
	    if(printAll>0) {
	      if(passed) {
		printf(" on ");
	      } else {
		printf("off ");
	      }
	    }
	    trig->Print("data");
	    break;
	  }
	}
      }
    }
  }

  if ((Level == 3) || (Level == 0)) {
    TTl3d* tl3d = TriggerBlock->Tl3d();
    Int_t npaths = tl3d->NPaths();
    for (int i = 0; i<npaths; i++) {
      int passed = tl3d->PathPassed(i);
      if (passed || printAll>0) {
	const TStnTrigger* trig = tt->GetTrigger(3,i);
	if (! banner_printed) {
	  trig->Print("banner");
	  banner_printed = 1;
	}
	if(printAll>0) {
	  if(passed) {
	    printf(" on ");
	  } else {
	    printf("off ");
	  }
	}
	trig->Print("data");
      }
    }
  }
  return 0;
}


//_____________________________________________________________________________
Int_t TStntuple::InitListOfL3Triggers(TStnModule* Module) {
  // given list of L3 triger names fill list of L3 triggers corresponding to
  // these bits for a given run

  TStnDBManager*   dbm = TStnDBManager::Instance();
  TStnTriggerTable* tt = (TStnTriggerTable*) dbm->GetTable("TriggerTable");
//-----------------------------------------------------------------------------
//  level 1,2,3 are handled the same way
//-----------------------------------------------------------------------------
  TString     name;
  TString     pattern;

  TObjArray* list_of_l3_names = Module->GetListOfL3TrigNames();
  TObjArray* list_of_triggers = Module->GetListOfL3Triggers ();

  list_of_triggers->Clear();
  int nt = list_of_l3_names->GetEntriesFast();

  for (int i=0; i<nt; i++) {
    TObjString* pat = (TObjString*) list_of_l3_names->UncheckedAt(i);
    for (int j=0; j<tt->NTriggers(3); j++) {

					// keep "human" numbering : 1,2,3

      const TStnTrigger* trig = tt->GetTrigger(3,j);
	
      name = trig->Name();
      name.ToUpper();
      pattern = pat->GetString();
      pattern.ToUpper();
      char* q = strstr(name.Data(),pattern.Data());
      if (q) {
				// match

	list_of_triggers->Add((TObject*) trig);
      }
    }
  }
  return 0;
}


//_____________________________________________________________________________
Int_t TStntuple::CheckTrigger(TStnTriggerBlock* TriggerBlock, 
					   Int_t level, char* trigger) 
{
  //
  TStnDBManager*   dbm = TStnDBManager::Instance();
  TStnTriggerTable* tt = (TStnTriggerTable*) dbm->GetTable("TriggerTable");
  TString tname(trigger);
  int star = 0;
  if(tname[tname.Length()-1]=='*') {
    tname = tname(0,tname.Length()-1);
    star = 1;
  }

  if (level == 1) {
    TTl2dL1Decision* l1_decision = TriggerBlock->Tl2d()->L1Decision();
    for (int i1=0; i1<tt->NL1Triggers(); i1++) {
      const TStnTrigger* trig = tt->GetTrigger(1,i1);
      if(trig) {
	const TString& tn = trig->Name();
	if (tn == tname || (star==1 && tn.BeginsWith(tname))) {
	  return l1_decision->L1DecisionBit(i1);
	}
      }
    }
    return -1;
  }

  if ((level == 2) || (level == 12)) {
    TTl2d*           tl2d = TriggerBlock->Tl2d();
    TTl2dL2Decision* l2_decision = tl2d->L2Decision();
    TTl2dL2DecisionUnPrescaled* l2_decisionu = tl2d->L2DecisionUnPrescaled();

    for (int i2=0; i2<tt->NL2Triggers(); i2++) {
      const TStnTrigger* trig = tt->GetTrigger(2,i2);
      const TString& tn = trig->Name();
      if (tn == tname || (star==1 && tn.BeginsWith(tname))) {
	int passed;
	if(level==2) {
	  passed = l2_decision->PrescaledBit(i2);
	} else{
	  passed = l2_decisionu->UnprescaledBit(i2);
	} 
	return passed;
      }
    }
    return -1;
  }

  if (level == 3) {
    TTl3d* tl3d = TriggerBlock->Tl3d();
    Int_t npaths = tl3d->NPaths();
    for (int i = 0; i<npaths; i++) {
      const TStnTrigger* trig = tt->GetTrigger(3,i);
      const TString& tn = trig->Name();
      if (tn == tname || (star==1 && tn.BeginsWith(tname))) {
	return tl3d->PathPassed(i);
      }
    }
    return -1;
  }

  return -1;
}



//_____________________________________________________________________________
Int_t TStntuple::SetJetTriggerThresholds(const char* TriggerName) {
 
  if (strcmp(TriggerName,"JET_20") == 0) {
    fL1TrigEt =   5.;
    fL2TrigEt =  15.;
    fL3TrigEt =  20.;
  }
  else if (strcmp(TriggerName,"JET_50") == 0) {
    fL1TrigEt =   5.;
    fL2TrigEt =  40.;
    fL3TrigEt =  50.;
  }
  else if (strcmp(TriggerName,"JET_70") == 0) {
    fL1TrigEt =  10.;
    fL2TrigEt =  60.;
    fL3TrigEt =  70.;
  }
  else if (strcmp(TriggerName,"JET_100") == 0) {
    fL1TrigEt =  10.;
    fL2TrigEt =  90.;
    fL3TrigEt = 100.;
  }
  else {
    printf(" TStntuple::SetJetTriggerThresholds: unknown trigger %s\n",
	   TriggerName);
  }
  return 0;
}



//____________________________________________________________________________
void TStntuple::SetTrackPointers(TStnMuonBlock* fMuonBlock, 
				 TStnTrackBlock* fTrackBlock) {
  int ntrk = fTrackBlock->NTracks();
  for(int i=0; i<fMuonBlock->NMuons(); i++) {
    TStnMuon* muo = fMuonBlock->Muon(i);
    int it = muo->TrackNumber();
    if(it>=0 && it<ntrk) {
      muo->SetTrack(fTrackBlock->Track(it));
    } else {
      muo->SetTrack(NULL);
    }
  }
}

//____________________________________________________________________________
void TStntuple::SetTrackPointers(TStnElectronBlock* fElectronBlock, 
				 TStnTrackBlock* fTrackBlock) {
  int ntrk = fTrackBlock->NTracks();
  for(int i=0; i<fElectronBlock->NElectrons(); i++) {
    TStnElectron* ele = fElectronBlock->Electron(i);
    int it = ele->TrackNumber();
    if(it>=0 && it<ntrk) {
      ele->SetTrack(fTrackBlock->Track(it));
    } else {
      ele->SetTrack(NULL);
    }
  }
}

//____________________________________________________________________________


void TStntuple::SetTrackPointersPh(
		   TPhoenixElectronBlock* fPhoenixElectronBlock, 
		   TStnTrackBlock* fPhoenixSiTrackBlock) {
  int nTrkBlk = fPhoenixSiTrackBlock->NTracks();
  if( nTrkBlk<=0) return;
  for(int ip=0; ip<fPhoenixElectronBlock->NElectrons(); ip++) {
    TStnElectron* pele = fPhoenixElectronBlock->Electron(ip);
    
    TStnTrack* best=NULL;
    // 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<nTrkBlk) {
	  TStnTrack* trk = fPhoenixSiTrackBlock->Track(itrk);
	  float chi2 = trk->Chi2()/trk->NSvxHits();
	  if(best==NULL) {
	    best = trk;
	    bestChi2 = chi2;
	  } else {
	    bool better = (trk->Algorithm()==4 && best->Algorithm()==11)
	      || ( (trk->Algorithm()==best->Algorithm()) && 
		   (chi2 < bestChi2)    );
	    if(better) {
	      best = trk;
	      bestChi2 = chi2;
	    }
	  }
	}
      }
    }
    
    pele->SetTrack(best);
    
  }
  return;
}


//____________________________________________________________________________
void TStntuple::SetTrackPointers(TStnTauBlock* fTauBlock, 
				 TStnTrackBlock* fTrackBlock) {
  int ntrk = fTrackBlock->NTracks();
  for(int i=0; i<fTauBlock->NTaus(); i++) {
    TStnTau* tau = fTauBlock->Tau(i);
    int it = fTauBlock->TrackNumber(i,0);
    if(it>=0 && it<ntrk) {
      tau->SetSeedTrack(fTrackBlock->Track(it));
    } else {
      tau->SetSeedTrack(NULL);
    }
  }
}

  
//____________________________________________________________________________
Int_t TStntuple::CorrectElectronEnergy(TStnEvent* event, const double zVtx) {

  int ientry = event->GetCurrentTreeEntry();
  
  TStnElectronBlock* fElectronBlock = 
    (TStnElectronBlock*)event->GetDataBlock("ElectronBlock");
  if(!fElectronBlock) {
    printf("TStntuple::CorrectElectronEnergy will not work without registering ElectronBlock\n");
    return 1;
  }
  fElectronBlock->GetEntry(ientry);
  
  TStnHeaderBlock* fHeaderBlock = 
    (TStnHeaderBlock*)event->GetDataBlock("HeaderBlock");
  bool qMc = fHeaderBlock->McFlag(); 
  int run = fHeaderBlock->RunNumber();
  
  for(int i=0; i<fElectronBlock->NElectrons(); i++) {
    TStnElectron* ele = fElectronBlock->Electron(i);

    float faceCorr = ele->EtFaceCor();
    float corr;
    if(ele->DetectorCode()==0) {
      corr = TStntuple::CorrectCemEnergy( qMc, run, faceCorr);
    } else {
      float pem3x3Eta = ele->Pem3x3Eta();
      float rawE = ele->EmE();
      float leakcorr = ele->LeakCorr()*ele->EmE();
      float pprE = ele->PprEnergy();
      corr = TStntuple::CorrectPemEnergy( qMc, run, faceCorr, pem3x3Eta, 
					  rawE, leakcorr, pprE, zVtx);
      //printf("new cor=%f\n",corr);
    }
    ele->SetTweak(corr/faceCorr);
  }
  
  return 0;
}  
  

//____________________________________________________________________________
Int_t TStntuple::CorrectPhotonEnergy(TStnEvent* event, const double zVtx) {

  int ientry = event->GetCurrentTreeEntry();
  
  TStnPhotonBlock* fPhotonBlock = 
    (TStnPhotonBlock*)event->GetDataBlock("PhotonBlock");
  if(!fPhotonBlock) {
    printf("TStntuple::CorrectPhotonEnergy will not work without registering 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);

    // only correction is the face corr, important to not use accessors
    float faceCorr = (fabs(pho->fEt)>1e-3 && fabs(pho->fEt)>1e-3 ? 
		      pho->fEtc/pho->fEt : 1.0 );
    float corr;
    if(pho->Detector()==0) {
      corr = TStntuple::CorrectCemEnergy( qMc, run, faceCorr);
    } else {
      float pem3x3Eta = pho->DetEta();
      float rawE = pho->EmE();
      float leakcorr = pho->fLcor*pho->Etc(); // leakage in energy
      float pprE = pho->PprE();
      corr = TStntuple::CorrectPemEnergy( qMc, run, faceCorr, pem3x3Eta, 
					  rawE, leakcorr, pprE, zVtx);
    }
    // tweak is everything after face corr
    pho->fTweak = corr/faceCorr;
  }
  
  return 0;
}  
  

//
// methods called by CorrectPhotonEnergy and CorrectElectronEnergy
// faceCorr is inverse of the offline sense = electronBlock Ecorr variable
// returns a factor to multiply raw energy
//____________________________________________________________________________
Float_t TStntuple::CorrectCemEnergy(bool qMc, int run, float faceCorr) {

  if(faceCorr<=0.0) faceCorr=1.;
  //-1 value occurs if there is a mapping failure, 
  //in this case dont apply correction
  
  //the energy scale factor to make the
  //a gausian fit between 86-98 peak at 91 GeV
  //see note 8274 for details
  float energyScaleFactor;

  float corrFac = 1.0; //the complete energy correction factor

  // all set to 1.0, 
  // CDF Note 8614 and Daryl Hare, 4/2007
  if(qMc) {
    // gen6 MC
    energyScaleFactor = 1.0;
  } else {
    energyScaleFactor = 1.0;
  }
  corrFac = energyScaleFactor*faceCorr;
    
  return corrFac;
}  
  

// methods called by CorrectPhotonEnergy and CorrectElectronEnergy
// faceCor is the in the offline sense = inverse of electronBlock Ecorr variable
// leakcorr is an energy
// returns a factor to multiply raw energy
//____________________________________________________________________________
Float_t TStntuple::CorrectPemEnergy(bool qMc, int run, float faceCorr, float pem3x3Eta, 
				  float rawE, float leakCorr, float pprE, float zVtx) {

  //note that you would *think* BcE() and EmE()*EtCor() 
  //would be the same thing but currently (06/05/06) 
  //  if its a plug electron and it doesnt have a 
  //pes cluster, BcE() is set to zero which occurs in
  //  StdPemElectron in ElectronUser 

  if(faceCorr<=0) faceCorr=1.;
  //-1 value occurs if there is a mapping failure, 
  //in this case dont apply correction

  float energyScaleFactor;
  float corrFac = 1.0; //the complete energy correction factor

  // CDF Note 8614 and Daryl Hare, 4/2007
  // added P11 and 12 8/2007
  if(qMc) {
    // gen6 MC
    energyScaleFactor = (pem3x3Eta<0? 0.9971 : 0.9965 );
  } else {
    if(run<=212133) {
      energyScaleFactor = (pem3x3Eta<0? 1.0140 : 1.0135);
    } else if(run<=222426) {
      energyScaleFactor = (pem3x3Eta<0? 1.0130 : 1.0230);
    } else if(run<=228596) {
      energyScaleFactor = (pem3x3Eta<0? 1.0257 : 1.0268);
    } else if(run<=233111) {
      energyScaleFactor = (pem3x3Eta<0? 1.0163 : 1.0147);
    } else if(run<=237795) {
      energyScaleFactor = (pem3x3Eta<0? 1.0128 : 1.0182);
    } else if(run<=241664) {
      energyScaleFactor = (pem3x3Eta<0? 1.0098 : 1.0177); // P12
    } else if(run<=246231) {
      energyScaleFactor = (pem3x3Eta<0? 1.0190 : 1.0163); //P13
    } else  {
      energyScaleFactor = (pem3x3Eta<0? 1.0190 : 1.0163);  //place holder
    }
  }

  //get the lateral leakage corrections for eta index 4 and 15 towers 
  //(eta ranges 3.00-3.50 and 1.1-1.2) which are discused in note 6181
  //note that we use a Stntuplised copy of the orginal 
  // AC++ namespace PemCorrAlg to do this

  float lateralLeakage=1.0,lowEtaHighZLeakage=1.0;
  if(!qMc) {
    int corr_status =  TStntuple::PEMfid3by3(pem3x3Eta,zVtx,lateralLeakage);
    if( corr_status != 0 ) lateralLeakage = 1.0;
    TStntuple::PEMLowEtaHighZLeakage(pem3x3Eta,zVtx,lowEtaHighZLeakage);
  }

  // the complete uncorrected energy:
  // 2x2 + leakage energy outside 2x2 cluster (note 6167) + pprEnergy
  if( fabs(pprE)>1000.0 ) pprE = 0.0; // protect against 1e32 flag value 
  float corrE = rawE + leakCorr + pprE;

  //the correction factors applied
  corrE *= energyScaleFactor/(lateralLeakage*lowEtaHighZLeakage*lowEtaHighZLeakage)*faceCorr;
  //printf("new %f %f %f %f %f %f\n",rawE,leakCorr,pprE,energyScaleFactor,lowEtaHighZLeakage,faceCorr);  
  corrFac = corrE/rawE; //now obtaining the correction factor
    
  return corrFac;
}  
  

//____________________________________________________________________________
void TStntuple::GetPemRadius(const float detEta, float& rPes, int& side)
{
//----------------------------------------------------------------------------
// Calcuates the radius and detector side (w/e = 0/1) at the PES z midpoint
// from detector eta
//----------------------------------------------------------------------------

  const float zpem = 0.5*(184.8 + 186.0);		// z at PES midpont
  const float tanby2 = exp(-fabs(detEta));		// tan(theta/2)

  rPes = zpem * 2.0*tanby2/(1.0 - tanby2*tanby2);
  side = (detEta < 0.0) ? 0 : 1;

  return;
}

//____________________________________________________________________________
//returns the low eta high Z leakage when the electron doesnt fully pass through
//the entire tower. Apply by dividing the energy by the square of this number
//(dont ask me why, thats Willis's convention)
void TStntuple::PEMLowEtaHighZLeakage(const float detEta,
				      const float eventZ0,float &leakLowEta)
{
  leakLowEta = 1.0;
  int side;
  float rPes  =  0.0;				// Radius in cm at PES z
  TStntuple::GetPemRadius( detEta, rPes, side );
  
  float delR = 138.8 - rPes;
  if(delR<40.0) leakLowEta = (90.61 + 0.05275*(delR-40.0))/90.61;
 
}


//____________________________________________________________________________

Int_t TStntuple::PEMfid3by3(const float detEta, const float zVtx, float& leakage)
{

  //----------------------------------------------------------------------------
  // Calcuates a leakage correction term for eta towers 4 and 15 and fiducial 
  // recomendations for these towers.
  // The leakage corrections for tower  4 are ESTIMATES and are for Zvtx = 0.
  // The leakage corrections for tower 15 are fits from CP ee data.
  //----------------------------------------------------------------------------

  leakage = 1.0;				// Default
  
  int side;
  float rPes  =  0.0;				// Radius in cm at PES z
  TStntuple::GetPemRadius( detEta, rPes, side );

  const float eta = fabs(detEta);
  if ( eta > 1.40 && eta < 3.00 ) return 0;
  
  //
  // Estimate the energy leakage for eta tower 4. Lateral leakage profile
  // in eta tower 5, from TB data, run 1009
  //

  if ( rPes < 18.5054 ) {
    float dr = rPes - 11.66;
    if ( dr < 0.0 ) return 1;
    leakage -= (0.308195*exp(-3.2447*dr)+0.191805*exp(-0.57457*dr));

  }  else {

    //
    // Leakage for eta tower 15 is from CP ee data. zExit is the z at the exit
    // radius, and z0 is the z at the front of the PEM. The 36.82 degree outer
    // edge of the PEM is defined by tan1 = tan(36.82deg) and the z at which it
    // intersects the outer radius, z1. zPes is the z of the PES midpoint.
    //
    const float rPem = 138.80;                  // Outer radius, |z| > |z1|
    if (rPes > rPem) return 1;

    const float z0   = 177.09;                  // |z| front edge in cm
    const float z1   = 186.53;

    const float zPes = ( detEta > 0.0 ) ?  
                          0.5*(184.8 + 186.0) :
                         -0.5*(184.8 + 186.0);

    float tanV  = rPes / (zPes - zVtx);    // tan(theta) of particle
    float zExit = zVtx + rPem/tanV;
    if (fabs(zExit) < z1) {
      float tan1 = ( detEta > 0.0 ) ? 0.74864 : -0.74864;
      zExit = (fabs(tanV-tan1) < 1e-5) ? z0 : zVtx/(1.0-tan1/tanV);
    }
    float dL = (fabs(zExit) - z0) * sqrt(tanV*tanV + 1.0);
    if (dL < 40.0) {
      float ddL = dL - 40.0;
      float corr = 1.0 - (5.964196e-4 + 5.723980e-04*ddL)*ddL;
      leakage = corr * corr;
    }
  }

  //
  // Rough fiducial recommendations. Eta tower 4: 3.8cm from edge (97% lateral
  // containment. Eta tower 15: 8cm from edge (85% total containment). Any
  // closer, the position measurements get unreliable. The inner beam hole
  // edge is r = 11.66cm, and the outer edge is r = 138.8cm.
  //

  if ( rPes <  15.6 ) return 1;
  if ( rPes > 130.8 ) return 1;

  return 0;
    
}



// fill in some missing Tau variables
// if (mode&0x1) then also fill TowerEt
// if (mode&0x2) then use TauBlock instead of PROD@CdfTauCollection-LeptonTrack
//____________________________________________________________________________
Int_t   TStntuple::TauTransients(TStnEvent* event, int mode) {

  int ientry = event->GetCurrentTreeEntry();

  TStnTauBlock* fTauBlock = NULL;
  if( (mode&0x2) > 0 ) {
    fTauBlock = (TStnTauBlock*)event->GetDataBlock("TauBlock");
  } else { // default
    fTauBlock = (TStnTauBlock*)event->
      GetDataBlock("PROD@CdfTauCollection-LeptonTrack");
  }
  if(!fTauBlock) {
    printf("TStntuple::TauTransients will not work without registering TauBlock\n");
    return 1;
  }
  fTauBlock->GetEntry(ientry);
  int nt = fTauBlock->NTaus();
  if(nt<=0) return 0;

  TStnTrackBlock* fTrackBlock = 
    (TStnTrackBlock*)event->GetDataBlock("TrackBlock");
  if(!fTrackBlock) {
    printf("TStntuple::TauTransients will not work without registering TrackBlock\n");
    return 1;
  }
  fTrackBlock->GetEntry(ientry);

  SetTrackPointers(fTauBlock, fTrackBlock);

  TStnJetBlock* fJetBlock = 
    (TStnJetBlock*)event->GetDataBlock("JetBlock");
  if(!fJetBlock) {
    printf("TStntuple::TauTransients will not work without registering JetBlock\n");
    return 1;
  }
  fJetBlock->GetEntry(ientry);

  // if requested, also fill highest towers
  TCalDataBlock* fCalDataBlock = NULL;
  if( mode & 0x1 ) {
    fCalDataBlock = (TCalDataBlock*)event->GetDataBlock("CalDataBlock");
    if(!fCalDataBlock) {
      printf("TStntuple::TauTransients will not work without registering CalDataBlock\n");
      return 1;
    }
    fCalDataBlock->GetEntry(ientry);
  }

  for(int it = 0; it<nt; it++) {
    TStnTau* tau = fTauBlock->Tau(it);
    TStnTrack* trk = tau->SeedTrack();
    float zces=-999.0,d0=-999.0,zeta=-999.0;
    if(!trk) {
      printf("TauTransients found TStnTau with no track\n");
    } else {
      float theta = trk->Momentum()->Theta();
      if(fabs(theta-TMath::PiOver2())<0.0001) {
	zces = trk->Z0();
      } else {
	zces = trk->Z0() + 184.0/tan(theta);
      }
      d0 = CorrectedD0(trk);
    }
    if(tau->SumPt10()<=0.0 || tau->CalEt()<=0.0) {
      printf("TauTransients found TStnTau with SumPt10==%f and CalEt=%f\n",
	     tau->SumPt10(),tau->CalEt());
    } else {
      zeta = tau->CalEt()/tau->SumPt10()*
	(1.0 - tau->EtEm()/tau->CalEt());
    }
    tau->SetSeedTrackZCes(zces);
    tau->SetSeedTrackD0(d0);
    tau->SetZeta(zeta);
    tau->SetJetNumber(-1);
    float minDr = 999.0;
    float tPhi =  TVector2::Phi_0_2pi(tau->VisMomentum()->Phi());
    for(int j=0; j<fJetBlock->NJets(); j++) {
      TStnJet* jet = fJetBlock->Jet(j);
      float dr =  DeltaR(tPhi, tau->DetEta(),
			 jet->Phi(), jet->DetEta());
      if(dr<minDr) {
	tau->SetJetNumber(j);
	minDr = dr;
      }
    }

    if(fabs(fgEventVertex)<200.0 && fabs(tau->SeedTrackZ())<200.0) {
      tau->SetDZVrt(tau->SeedTrackZ()-fgEventVertex);
    } else {
      tau->SetDZVrt(-999.0);
    }

    if( mode & 0x1 ) {  // fill towerEt
      TStnLinkBlock* links = fTauBlock->TowerLinkList();
      tau->fTowerEt[0] = 0.0;
      tau->fTowerEt[1] = 0.0;
      for(int i=0; i<links->NLinks(it); i++) {
	int key = links->Index(it,i);
	int phit = key & 0x3F;
	int etat = (key>>8) & 0x3F;
	TCalTower* tow = fCalDataBlock->Tower(etat,phit);
	if(!tow) {
	  printf("TStntuple::TauTransients: Error - could not find Cal tower p=%2d eta=%2d for tau=%d\n",phit,etat,it);
	  return 1;
	}
	float e = tow->Energy();
	if(e > tau->fTowerEt[0]) {
	  tau->fTowerEt[1] = tau->fTowerEt[0];
	  tau->fTowerEt[0] = e;
	} else if(e > tau->fTowerEt[1]) {
	  tau->fTowerEt[1] = e;
	}
      }
    } // end if fill cal

  } // end tau loop
  
}

//____________________________________________________________________________
// find the z vertex position indicated by the tracks
// associated with a jet
Float_t TStntuple::JetZpos(TStnJet* jet, TStnJetBlock* jbl, 
			   TStnTrackBlock* trbl, Int_t& atracks) {
  Int_t jn = jet->Number();
  Int_t ntracks = jbl->NTracks(jn);
  float jz = -999.;
  atracks = 0;
  if(ntracks < 1) return jz;

  float window = 5;
  float jz_back = jz;
  for(int iter = 0; iter < 3; iter++){

    atracks = 0;
    jz = -999.;
    for(int i = 0; i < ntracks; i++){
      Int_t itrk = jbl->TrackNumber(jn,i);
      TStnTrack *trk = trbl->Track(itrk);
      float zpos = trk->Z0();

      if(atracks == 0 ){ jz = zpos; atracks++;}
      else if(iter > 0 && fabs(jz_back-zpos) > window) continue;
      else{
	jz = jz + zpos;
	atracks++;
      }
    }
    if (atracks>0) jz = jz/atracks;
    jz_back = jz;
  }
  return jz;

}


//____________________________________________________________________________
// ADD the return value to the raw tower timing
// if distance to vertex is smaller than distance to z=0, correction is positive
Float_t TStntuple::EmTimingTOFCorrection(Float_t z, Int_t ieta) {

  double deltaT = -999.;
  int side = (ieta <  TENETA/2 ? -1 : 1);

  double rem = EmTowerRadius(ieta) ;
  double zem = EmTowerZ(ieta);

  double trad = sqrt(rem*rem + zem*zem);
 
  double dztr = zem - z;
  double rtrack =  sqrt(dztr*dztr + rem*rem);
  deltaT = 1./29.9792*(trad - rtrack);
  return deltaT;

}


Int_t TStntuple::LastLayerCOT(TStnTrack* trk) 
{
  float slr[8] = {46.774,58.534,70.295,82.055,93.815,
		  105.575,117.335,129.096};
  float wireSpacing = 0.7112;
  float c35 = cos(35.0*TMath::Pi()/180.);
  float zMax = 155.0;

  double trkd0 = trk->D0();
  double cu = trk->C0();
  double z0 = trk->Z0();
  double cott = trk->Lam0();

  double  outerCircleR = ( fabs(cu) >1.0e-30 ? 
			   fabs(trkd0 + 0.99 / cu) : 1.0e10 );

  int ilr=0;
  int exitLayer = 96;
  for (int isl=0; isl<8; isl++) {
    for (int iwr=0; iwr<12; iwr++, ilr++) {

      double r = slr[isl] + (iwr-5.5)*wireSpacing*c35;
      double chord = (r * r - trkd0 * trkd0) / (1. + 2. * cu * trkd0);
      double l2d = 0.;
      if (chord <= 0.) {
	return 0;
      }
      chord = sqrt(chord);
      double absc = fabs(cu);
      double cr = absc * chord;
      if (cr > 1.) {
	l2d = TMath::Pi() / (2. * absc);
      } else if (cu == 0.) {
	l2d = r;
      } else {
	l2d = asin(cr) / absc;
      }

      if ((z0 + cott * l2d) > zMax ||
	  (z0 + cott * l2d) < -zMax ||
	  outerCircleR < r) {
	exitLayer = ilr;
	break;
      }
    }
  }
  return exitLayer-1;
}


//____________________________________________________________________________
// Run II CES strip-to-wire corrections
// implementation by  Anton Anastassov
double TStntuple::swMultiplier(int wedge, int side, 
			       float locX, float locZ, int mcflag)
{  
  const float parStrip[96][5]=
    {
      { 0.98927E+00,-0.72789E+00, 0.83092E+00, 0.15376E+01,-0.16665E+01},
      { 0.94156E+00, 0.11623E+01,-0.57364E+01, 0.99365E+01,-0.54566E+01},
      { 0.83466E+00, 0.13723E+01,-0.63057E+01, 0.96061E+01,-0.45566E+01},
      { 0.10428E+01,-0.58101E+00,-0.59273E+00, 0.34244E+01,-0.23008E+01},
      { 0.10752E+01,-0.18350E+01, 0.87342E+01,-0.13444E+02, 0.64237E+01},
      { 0.92379E+00, 0.15513E+01, 0.10889E+01,-0.55406E+01, 0.29358E+01},
      { 0.10323E+01, 0.24206E+00,-0.39890E+01, 0.91388E+01,-0.54926E+01},
      { 0.10167E+01, 0.83698E+00, 0.26662E+00,-0.16309E+01, 0.49134E+00},
      { 0.10412E+01, 0.66202E+00,-0.23451E+01, 0.43932E+01,-0.28348E+01},
      { 0.99916E+00,-0.10245E+01, 0.15430E+01, 0.36053E+00,-0.87852E+00},
      { 0.10817E+01,-0.81794E+00, 0.27638E+01,-0.23328E+01, 0.21600E+00},
      { 0.11385E+01,-0.18343E+01, 0.29822E+01,-0.12336E+01,-0.11167E+00},
      { 0.98685E+00,-0.47050E-02, 0.36977E-01, 0.77217E+00,-0.85445E+00},
      { 0.99318E+00,-0.52613E+00, 0.55471E+00, 0.14689E+01,-0.15045E+01},
      { 0.10200E+01,-0.76092E+00, 0.12957E+01, 0.10990E+00,-0.71795E+00},
      { 0.84253E+00,-0.24040E+00,-0.41351E+01, 0.92013E+01,-0.47557E+01},
      { 0.10234E+01,-0.19477E+01, 0.71961E+01,-0.97014E+01, 0.44851E+01},
      { 0.10477E+01,-0.18070E+01, 0.76925E+01,-0.10112E+02, 0.41983E+01},
      { 0.11158E+01,-0.33147E+00, 0.56600E+01,-0.11357E+02, 0.59802E+01},
      { 0.10545E+01,-0.13508E+01, 0.89643E+01,-0.14547E+02, 0.68950E+01},
      { 0.92812E+00, 0.34348E-01, 0.47719E+01,-0.10308E+02, 0.56518E+01},
      { 0.11070E+01,-0.28456E+01, 0.11661E+02,-0.16174E+02, 0.72460E+01},
      { 0.10405E+01,-0.19008E+00,-0.17872E+01, 0.40321E+01,-0.20663E+01},
      { 0.95256E+00, 0.22067E+01,-0.79974E+01, 0.89501E+01,-0.32218E+01},
      { 0.90762E+00, 0.42242E+01,-0.30835E+01,-0.66207E+01, 0.57449E+01},
      { 0.99587E+00, 0.16539E+01, 0.87925E+00,-0.57528E+01, 0.31765E+01},
      { 0.94610E+00,-0.42430E+00, 0.49536E+01,-0.79265E+01, 0.33986E+01},
      { 0.93755E+00, 0.32687E+01,-0.74024E+01, 0.52183E+01,-0.10626E+01},
      { 0.10378E+01, 0.40056E+00, 0.32893E+01,-0.64032E+01, 0.25354E+01},
      { 0.10155E+01, 0.50655E+00,-0.37420E+01, 0.63944E+01,-0.32463E+01},
      { 0.10935E+01,-0.12344E+01, 0.30010E+01,-0.37717E+01, 0.19661E+01},
      { 0.95837E+00, 0.27483E+01,-0.82024E+01, 0.86695E+01,-0.31957E+01},
      { 0.10377E+01, 0.73794E+00,-0.65297E+01, 0.13752E+02,-0.81819E+01},
      { 0.11332E+01,-0.34887E+01, 0.10744E+02,-0.12517E+02, 0.51309E+01},
      { 0.10031E+01,-0.12909E+01, 0.51972E+01,-0.75828E+01, 0.36672E+01},
      { 0.94969E+00,-0.57607E+00, 0.35291E+01,-0.52264E+01, 0.23135E+01},
      { 0.96809E+00,-0.18516E+00, 0.88589E+00,-0.26505E+00,-0.38087E+00},
      { 0.10353E+01, 0.15004E+01, 0.67913E+01,-0.18299E+02, 0.10082E+02},
      { 0.10512E+01,-0.28709E+01, 0.94787E+01,-0.10432E+02, 0.38206E+01},
      { 0.11321E+01,-0.18613E+00, 0.35865E+01,-0.65991E+01, 0.30551E+01},
      { 0.98576E+00, 0.35470E+00,-0.20281E+01, 0.51963E+01,-0.36502E+01},
      { 0.10382E+01,-0.97322E+00, 0.44857E+01,-0.59782E+01, 0.24103E+01},
      { 0.95063E+00,-0.66303E+00, 0.19444E+01,-0.10799E+01,-0.21045E+00},
      { 0.10673E+01,-0.42104E+00, 0.38712E+01,-0.65495E+01, 0.30592E+01},
      { 0.10659E+01, 0.12027E+01, 0.13632E+01,-0.57787E+01, 0.31707E+01},
      { 0.86085E+00, 0.13700E+01,-0.71213E+00,-0.21700E+01, 0.16738E+01},
      { 0.99506E+00, 0.20430E+00, 0.71564E+01,-0.14611E+02, 0.72151E+01},
      { 0.10149E+01,-0.14388E+01, 0.82624E+01,-0.12790E+02, 0.59559E+01},
      { 0.10552E+01,-0.13250E+01, 0.67305E+01,-0.47884E+01,-0.23068E+01},
      { 0.10144E+01,-0.15314E+01, 0.34527E+01, 0.16407E+01,-0.47991E+01},
      { 0.10020E+01,-0.12889E+00, 0.23416E+01,-0.27817E+01, 0.45152E+00},
      { 0.10412E+01,-0.85230E+00, 0.40170E+01,-0.73757E+01, 0.45427E+01},
      { 0.11416E+01,-0.64846E+00, 0.11791E+02,-0.26512E+02, 0.16991E+02},
      { 0.11044E+01,-0.94712E+00,-0.12905E+01, 0.86584E+01,-0.75167E+01},
      { 0.10913E+01,-0.72469E+00, 0.69306E-01, 0.48082E+01,-0.51368E+01},
      { 0.11550E+01,-0.42281E+01, 0.23547E+02,-0.45375E+02, 0.27691E+02},
      { 0.10790E+01,-0.10642E+00,-0.85369E+00, 0.73909E+01,-0.75295E+01},
      { 0.97801E+00, 0.27104E+00, 0.34896E+01,-0.10481E+02, 0.73478E+01},
      { 0.96482E+00, 0.21995E+01,-0.10250E+02, 0.20356E+02,-0.14114E+02},
      { 0.11462E+01,-0.21821E+01, 0.85516E+01,-0.11115E+02, 0.43069E+01},
      { 0.97531E+00,-0.18999E+00, 0.57105E+00,-0.11403E+01, 0.64286E+00},
      { 0.10287E+01, 0.63503E+00,-0.33713E+01, 0.10022E+02,-0.84370E+01},
      { 0.97172E+00, 0.18030E+01,-0.61123E+01, 0.13892E+02,-0.10930E+02},
      { 0.10929E+01,-0.72392E+00, 0.31122E+01,-0.32497E+01,-0.56932E-01},
      { 0.11048E+01,-0.16933E+01, 0.13012E+02,-0.25227E+02, 0.14388E+02},
      { 0.11459E+01,-0.76466E+00, 0.41916E+01,-0.45145E+01, 0.94634E-01},
      { 0.11132E+01,-0.57658E+00, 0.20561E+00, 0.19348E+01,-0.18563E+01},
      { 0.10681E+01,-0.90554E+00, 0.66167E+01,-0.12129E+02, 0.61939E+01},
      { 0.10233E+01, 0.37978E+00, 0.80626E-01, 0.10365E+01,-0.18212E+01},
      { 0.10592E+01,-0.21141E+01, 0.12367E+02,-0.20017E+02, 0.94538E+01},
      { 0.98210E+00, 0.87828E+00,-0.27029E+00,-0.32694E+01, 0.33547E+01},
      { 0.99625E+00, 0.15012E+01,-0.11209E+02, 0.26334E+02,-0.18372E+02},
      { 0.98738E+00, 0.18468E+00, 0.37655E+01,-0.10171E+02, 0.72753E+01},
      { 0.98944E+00,-0.37252E+00, 0.69208E-01, 0.40313E+01,-0.42588E+01},
      { 0.96699E+00,-0.20336E+00, 0.39630E+01,-0.88459E+01, 0.49565E+01},
      { 0.96384E+00, 0.22442E+00, 0.56968E+01,-0.16496E+02, 0.11431E+02},
      { 0.10549E+01,-0.12063E+01, 0.10323E+02,-0.15726E+02, 0.70324E+01},
      { 0.92561E+00, 0.20679E+00, 0.20544E+01,-0.79255E+01, 0.52314E+02},
      { 0.11443E+01,-0.84855E+00, 0.52197E+01,-0.10425E+02, 0.64171E+01},
      { 0.10687E+01,-0.14030E+01, 0.65939E+01,-0.10504E+02, 0.56451E+01},
      { 0.99060E+00, 0.67261E+00,-0.37280E+01, 0.64528E+01,-0.35660E+01},
      { 0.92246E+00, 0.23544E+01,-0.72280E+01, 0.73911E+01,-0.22879E+01},
      { 0.95429E+00, 0.14165E+01,-0.90442E+01, 0.23096E+02,-0.17504E+02},
      { 0.99243E+00,-0.52551E+00,-0.82068E+00, 0.16205E+01,-0.45873E-01},
      { 0.11231E+01, 0.33134E+00,-0.45179E+01, 0.11757E+02,-0.88601E+01},
      { 0.10490E+01,-0.14828E+01, 0.11809E+01, 0.50873E+01,-0.64088E+01},
      { 0.10309E+01,-0.28839E+00,-0.28975E+01, 0.76593E+01,-0.53057E+01},
      { 0.96774E+00, 0.37221E+00,-0.41721E+00,-0.13831E+01, 0.16910E+01},
      { 0.96776E+00, 0.22864E+00,-0.19311E+01, 0.41923E+01,-0.27102E+01},
      { 0.10984E+01,-0.16618E+01, 0.74613E+01,-0.10938E+02, 0.46438E+01},
      { 0.10580E+01,-0.47709E+00, 0.85925E+01,-0.19906E+02, 0.12149E+02},
      { 0.78031E+00, 0.23712E+01,-0.10094E+02, 0.18109E+02,-0.11093E+02},
      { 0.96844E+00, 0.36185E+00, 0.69727E+01,-0.16945E+02, 0.99676E+01},
      { 0.10329E+01,-0.13939E+01, 0.15143E+02,-0.31062E+02, 0.17660E+02},
      { 0.10048E+01,-0.73044E+00, 0.69556E+01,-0.15443E+02, 0.10001E+02},
      { 0.88110E+00, 0.19613E+01,-0.84979E+01, 0.19005E+02,-0.13941E+02}
    };
  const float parWire[96][5]=
    {
      { 0.49670E-01,-0.91118E+00, 0.46899E+01,-0.80364E+01, 0.42591E+01},
      { 0.35005E-01,-0.11172E+01, 0.50649E+01,-0.75352E+01, 0.35958E+01},
      { 0.37137E-02,-0.58745E+00, 0.28059E+01,-0.47251E+01, 0.26803E+01},
      {-0.34112E-01, 0.62521E-01, 0.52929E+00,-0.13322E+01, 0.79745E+00},
      { 0.67750E-01,-0.12134E+01, 0.61749E+01,-0.10298E+02, 0.52755E+01},
      {-0.25271E+00, 0.64957E-01, 0.28302E+01,-0.49669E+01, 0.25929E+01},
      { 0.44928E-01,-0.83844E+00, 0.32431E+01,-0.42590E+01, 0.17902E+01},
      { 0.26135E-01,-0.19073E+01, 0.83196E+01,-0.12774E+02, 0.67398E+01},
      { 0.56553E-01,-0.12240E+01, 0.51999E+01,-0.77802E+01, 0.38366E+01},
      {-0.11769E-01,-0.10194E+01, 0.43258E+01,-0.54243E+01, 0.21783E+01},
      {-0.37892E-01,-0.44243E+00, 0.89656E+00, 0.97671E+00,-0.14198E+01},
      { 0.18833E+00,-0.24430E+01, 0.72778E+01,-0.82853E+01, 0.33935E+01},
      { 0.18336E+00,-0.22808E+01, 0.79046E+01,-0.11016E+02, 0.53809E+01},
      {-0.86884E-01, 0.24352E+00,-0.61872E+00, 0.57860E+00, 0.13375E+00},
      { 0.23566E+00,-0.21555E+01, 0.62023E+01,-0.78864E+01, 0.37320E+01},
      { 0.11706E-01,-0.15301E+00,-0.15750E-01, 0.51718E+00,-0.29615E+00},
      {-0.98427E+00, 0.82864E+01,-0.23823E+02, 0.30009E+02,-0.13604E+02},
      {-0.10939E+00, 0.76631E+00,-0.35289E+01, 0.63686E+01,-0.34485E+01},
      {-0.22394E-01,-0.26912E+01, 0.12135E+02,-0.18431E+02, 0.96556E+01},
      { 0.70632E-01,-0.98542E+00, 0.28823E+01,-0.33571E+01, 0.15033E+01},
      {-0.68300E-01,-0.56080E+00, 0.52631E+01,-0.10861E+02, 0.65494E+01},
      {-0.12780E+00, 0.74661E+00,-0.61872E+00,-0.39715E+00, 0.29984E+00},
      { 0.14863E+00,-0.10075E+01, 0.19391E+01,-0.15218E+01, 0.44613E+00},
      {-0.79754E-01,-0.57837E+00, 0.52696E+01,-0.96799E+01, 0.51624E+01},
      { 0.53998E-01,-0.73894E+00, 0.36979E+01,-0.69622E+01, 0.41173E+01},
      {-0.56415E-01,-0.30018E+00, 0.25088E+01,-0.44932E+01, 0.24680E+01},
      { 0.33093E-01, 0.40520E-01,-0.11378E+01, 0.20191E+01,-0.89417E+00},
      { 0.13594E+00,-0.14767E+01, 0.14648E+01, 0.43981E+01,-0.49277E+01},
      { 0.26410E+00,-0.18941E+01, 0.57817E+01,-0.81457E+01, 0.39611E+01},
      {-0.22174E+00, 0.20635E+00, 0.47960E+00, 0.40228E+00,-0.70950E+00},
      { 0.79609E-01,-0.14764E+01, 0.58775E+01,-0.83890E+01, 0.39841E+01},
      { 0.16332E-01,-0.12509E+01, 0.62444E+01,-0.96726E+01, 0.47295E+01},
      {-0.96606E-01,-0.83100E+00, 0.60561E+01,-0.10151E+02, 0.51561E+01},
      { 0.56823E-01,-0.26450E+00, 0.33857E-01, 0.43756E+00,-0.22616E+00},
      { 0.17152E-02, 0.10256E-01,-0.14976E+01, 0.39518E+01,-0.24781E+01},
      {-0.60160E-01, 0.26912E+00,-0.27672E+00, 0.27495E+00,-0.25464E+00},
      {-0.11676E+00, 0.64737E+00,-0.29013E+01, 0.59247E+01,-0.36055E+01},
      {-0.28350E-01, 0.48997E+00,-0.23575E+01, 0.40377E+01,-0.22016E+01},
      { 0.19359E-01,-0.42799E+00, 0.10681E+01,-0.10967E+01, 0.56403E+00},
      { 0.33851E+00,-0.19537E+01, 0.55933E+01,-0.82174E+01, 0.41421E+01},
      {-0.70672E-01,-0.21541E-01, 0.35882E+01,-0.86409E+01, 0.52287E+01},
      { 0.10793E+00,-0.11017E+01, 0.33883E+01,-0.43105E+01, 0.19561E+01},
      { 0.76288E-01, 0.13369E+00,-0.19848E+01, 0.31249E+01,-0.13138E+01},
      { 0.32135E+00,-0.23189E+01, 0.68914E+01,-0.88888E+01, 0.38163E+01},
      {-0.65295E-01,-0.22657E+01, 0.11106E+02,-0.17636E+02, 0.95284E+01},
      { 0.37850E+00,-0.35580E+01, 0.10732E+02,-0.14789E+02, 0.76042E+01},
      { 0.14216E+00,-0.10566E+01, 0.18925E+01,-0.10598E+01, 0.10152E+00},
      { 0.20978E+00,-0.10150E+01, 0.14906E+01,-0.90718E+00, 0.13818E+00},
      { 0.31781E-01,-0.65668E+00, 0.30078E+01,-0.48525E+01, 0.25356E+01},
      { 0.56025E-01,-0.13198E+01, 0.61353E+01,-0.95233E+01, 0.46985E+01},
      {-0.48593E-01, 0.48500E+00,-0.20200E+01, 0.33211E+01,-0.17544E+01},
      { 0.72377E-01,-0.97329E+00, 0.44274E+01,-0.68816E+01, 0.32948E+01},
      { 0.48156E-01,-0.33353E+00,-0.33753E+01, 0.11695E+02,-0.84011E+01},
      {-0.33645E-01, 0.33773E+00,-0.86371E+00, 0.91303E+00,-0.37800E+00},
      { 0.62599E-01,-0.19353E+00, 0.64924E+00,-0.17067E+01, 0.12222E+01},
      {-0.18209E+00, 0.14527E+01,-0.36761E+01, 0.40082E+01,-0.16048E+01},
      { 0.14681E+00,-0.18888E+01, 0.52428E+01,-0.57565E+01, 0.24464E+01},
      { 0.25977E-01,-0.13965E+01, 0.62148E+01,-0.88135E+01, 0.40205E+01},
      { 0.33375E+00,-0.31867E+01, 0.97910E+01,-0.12086E+02, 0.50876E+01},
      { 0.29760E+00,-0.21336E+01, 0.57442E+01,-0.72264E+01, 0.33061E+01},
      { 0.13306E+00,-0.15375E+01, 0.60835E+01,-0.88990E+01, 0.41637E+01},
      {-0.31583E-01, 0.35776E+00,-0.10046E+01, 0.74813E+00, 0.27622E-02},
      {-0.27540E+00, 0.29228E+01,-0.73562E+01, 0.71700E+01,-0.26330E+01},
      {-0.25149E-01, 0.28187E+00,-0.11292E+00,-0.16144E+01, 0.16276E+01},
      {-0.18448E+00, 0.11969E+01,-0.25972E+01, 0.28842E+01,-0.13467E+01},
      { 0.17909E+00,-0.14645E+00,-0.16447E+01, 0.30201E+01,-0.15633E+01},
      { 0.74148E-01, 0.15719E+00,-0.21571E+01, 0.33400E+01,-0.13435E+01},
      { 0.38023E+00,-0.34372E+01, 0.91859E+01,-0.10375E+02, 0.43518E+01},
      { 0.70042E-01,-0.80325E+00, 0.51045E+01,-0.10898E+02, 0.67736E+01},
      {-0.31468E+00, 0.22350E+01,-0.46672E+01, 0.44691E+01,-0.18224E+01},
      {-0.54437E-01,-0.92869E+00, 0.44661E+01,-0.64237E+01, 0.31807E+01},
      { 0.34238E-01,-0.73981E+00, 0.48708E+01,-0.10146E+02, 0.62433E+01},
      {-0.59331E-01, 0.52402E-01, 0.11166E+01,-0.25705E+01, 0.15180E+01},
      { 0.32658E-01,-0.31759E+00, 0.15426E+01,-0.30981E+01, 0.19325E+01},
      {-0.91996E-01, 0.26008E+00,-0.21287E+01, 0.51086E+01,-0.30285E+01},
      {-0.92926E-01,-0.83337E+00, 0.50068E+01,-0.74741E+01, 0.35463E+01},
      {-0.60207E-01,-0.17143E+01, 0.93422E+01,-0.14487E+02, 0.71258E+01},
      { 0.18851E+00,-0.13169E+01, 0.38534E+01,-0.49623E+01, 0.21305E+01},
      { 0.81822E-01,-0.14274E+01, 0.57815E+01,-0.82225E+01, 0.38021E+01},
      {-0.16992E-01,-0.74765E+00, 0.46636E+01,-0.81165E+01, 0.43273E+01},
      {-0.11806E+00, 0.12006E+01,-0.45541E+01, 0.75283E+01,-0.42320E+01},
      { 0.24968E-01,-0.16607E+01, 0.63237E+01,-0.76593E+01, 0.30619E+01},
      { 0.13903E+00,-0.16721E+01, 0.56669E+01,-0.81946E+01, 0.42843E+01},
      { 0.10918E+00,-0.16271E+01, 0.65202E+01,-0.95666E+01, 0.46138E+01},
      {-0.14424E+00, 0.80510E+00,-0.33992E+01, 0.68048E+01,-0.41327E+01},
      {-0.25867E-01,-0.58949E+00, 0.25023E+01,-0.30289E+01, 0.12188E+01},
      { 0.17979E+00,-0.17392E+01, 0.51650E+01,-0.64943E+01, 0.29590E+01},
      { 0.31521E-01,-0.15899E+01, 0.76829E+01,-0.11760E+02, 0.57133E+01},
      {-0.90371E-01, 0.19170E+00,-0.11854E+00, 0.32643E+00,-0.23790E+00},
      { 0.61850E-01,-0.18795E+01, 0.85936E+01,-0.12948E+02, 0.62536E+01},
      {-0.60892E-01, 0.56664E+00,-0.25711E+01, 0.41328E+01,-0.19931E+01},
      { 0.38760E-02,-0.17341E-01,-0.11125E+01, 0.31157E+01,-0.20168E+01},
      { 0.45335E+00,-0.34381E+01, 0.99217E+01,-0.12772E+02, 0.57586E+01},
      { 0.55769E+00,-0.41962E+01, 0.11988E+02,-0.16556E+02, 0.84184E+01},
      { 0.75082E-01,-0.12301E+01, 0.46392E+01,-0.66870E+01, 0.33269E+01},
      { 0.16416E+00,-0.12534E+01, 0.28012E+01,-0.28846E+01, 0.12503E+01}
    };


  // Do not apply corrections for MC events

  if (mcflag) return 1.0;


  const int   nFitPar = 5;

  const float xFid    =  21.5;
  const float zFidMin =   9.0;
  const float zFidMax = 220.0;

  float zMin = 6.16;
  float zMax = 121.183;

  int  halfBarrel = 0;   // half-barrel: near->0, far->1
  if (locZ > zMax) {
    halfBarrel = 1;
    zMin       = zMax;
    zMax       = zMin + 59.0*2.007;
  }

  int halfWedgeId = wedge + side*24 + halfBarrel*48;   // global CES module ID: numbered in the range 0-95

  float x = -locX;   // sign convention used in the fits is different from CES cluster local coordinates
  if (x < -xFid)   x = -xFid; 
  if (x >  xFid)   x =  xFid;

  float z =  locZ;
  if (z < zFidMin) z =  zFidMin;
  if (z > zFidMax) z =  zFidMax;
  if (halfWedgeId==77 && z>160.0) z = 160.0;      // special case: East side, wedge 5, "far" half-barrel

  float xLocRel = (x + 23.27744) / (2.0*23.27744);
  float zLocRel = (z - zMin) / (zMax - zMin);

  float fX  = parWire [halfWedgeId][0];
  float fZ  = parStrip[halfWedgeId][0];
  float argX = 1.0;
  float argZ = 1.0;

  for (int i=1; i<nFitPar; ++i) {
    argX *= xLocRel;
    argZ *= zLocRel;
    fX   += argX * parWire [halfWedgeId][i];
    fZ   += argZ * parStrip[halfWedgeId][i];
  }


  float swMulti = 1.0;
  if (fX + fZ > 0.0) swMulti = 1.0/(fX + fZ);

  return swMulti;

}

//------------------------------------------------------------------------------
Float_t TStntuple::PredictedTof(Float_t arclen, Float_t p, Float_t mass) 
{
  static Float_t const c(299792458. * 100.0/1.0e9);  // cm/nsec
  Float_t const inverseBeta(sqrt(1.0+(mass*mass)/(p*p)));
  return arclen * inverseBeta / c;
}

//------------------------------------------------------------------------------
Float_t TStntuple::PredictedTofMuon(TStnTofMatch* tm, int run) 
{   

  const Float_t parMPi115[6] = {0.3192,1.099,-0.006002,
				-0.00000001568,0.0000000005732,-3.706e-14};
  
  const Float_t parMPi1619[6] = {1.0,1.0,-0.006002,
				 -0.00000001568,0.0000000005732,-3.706e-14};

  Float_t zc;
  const Float_t mmu = 0.10568;

  if (run<=168766)
    zc = ZCorrMPi(parMPi115,tm->fTofZBar);
  else
    zc = ZCorrMPi(parMPi1619,tm->fTofZBar);
  
  return PredictedTof(tm->ArcLength(),tm->Momentum(),mmu) - 0.005 + zc; 
}

//------------------------------------------------------------------------------
Float_t TStntuple::PredictedTofPion(TStnTofMatch* tm, int run) 
{

  const Float_t parMPi115[6] = {0.3192,1.099,-0.006002,
				-0.00000001568,0.0000000005732,-3.706e-14};
  
  const Float_t parMPi1619[6] = {1.0,1.0,-0.006002,
				 -0.00000001568,0.0000000005732,-3.706e-14};

  const Float_t mpi = 0.13957;
  Float_t zc;

  if (run<=168766)
    zc = ZCorrMPi(parMPi115,tm->TofZBar());
  else
    zc = ZCorrMPi(parMPi1619,tm->TofZBar());
  
  return PredictedTof(tm->ArcLength(),tm->Momentum(),mpi) - 0.005 + zc; 
}

//------------------------------------------------------------------------------
Float_t TStntuple::PredictedTofKaon(TStnTofMatch* tm, int run) 
{

  const Float_t parK115[13] = {-0.277558,537.904,-0.0499686,0.0199611,
			       -0.195865,-0.404445,0.63299,0.6366,
			       13.5379,-4.691,-0.00875,0.000006526,
			       -0.0000000004516};
  
  const Float_t parK1619[13] = { 0.684301,2.88137,-0.0260869,
				 0.0145981,0.0345695,-6.13759,
				 0.775240,0.3646,4.78312,
				 -2.83314,-0.01374,0.000008436,
				 -0.0000000005318};

  const Float_t mk = 0.49368;

  Float_t ptc, zc;
  Float_t exptof = PredictedTof(tm->ArcLength(),tm->Momentum(),mk);
  Float_t pt = tm->Momentum()*(tm->ArcLength()/tm->PathLength());

  if (run<=168766){
    ptc = PtCorrKP(parK115,pt);
    zc  = ZCorrKP(parK115,pt,tm->TofZBar());
  }
  
  if (run>168766){
    ptc = PtCorrKP(parK1619,pt);
    if (tm->TofNPulses()>2)
      ptc = 0.5*ptc;
    zc  = ZCorrKP(parK1619,pt,tm->TofZBar());
  }
  return exptof + ptc + zc;
}

//------------------------------------------------------------------------------
Float_t TStntuple::PredictedTofProton(TStnTofMatch* tm, int run) 
{

  const Float_t parP115[13] ={ 1.74245,0.77923,-0.128467,-0.550532,
			       76.6338,-3.00108,0.7251,0.3745,
			       3.17085,-1.83036,-0.021526,0.0000102575,
			       -0.000000000526471};
  
  const Float_t parP1619[13] = { 0.617078,1.14603,-0.0200317,0.0406387,
				 1.53921,-1.53244,0.7903,0.2564,
				 2.621,-1.255,-0.02142,0.000009875,
				 -0.0000000005694};
  const Float_t mp = 0.93827;

  Float_t ptc, zc;
  Float_t exptof = PredictedTof(tm->ArcLength(),tm->Momentum(),mp);
  Float_t pt = tm->Momentum()*(tm->ArcLength()/tm->PathLength());

  if (run<=168766){
    ptc = PtCorrKP(parP115,pt);
    zc  = ZCorrKP(parP115,pt,tm->TofZBar());
  }
  
  if (run>168766){
    ptc = PtCorrKP(parP1619,pt);
    if (tm->TofNPulses()>2)
      ptc = 0.5*ptc;
    zc  = ZCorrKP(parP1619,pt,tm->TofZBar());
  }
  
  return exptof+ptc+zc;
}

//------------------------------------------------------------------------------

Float_t TStntuple::PtCorrKP(const Float_t par[13], Float_t pt) 
{
  Float_t tmp = pt-par[5];
  return par[0]*exp(-par[1]*pt)+par[2]/(pt*pt) - par[3]/pt - 
    par[4]/(tmp*tmp*tmp); 
}

//------------------------------------------------------------------------------

Float_t TStntuple::ZCorrKP(const Float_t par[13], Float_t pt, Float_t z) 
{
  Float_t tmp = z*(par[6]+pt*par[7]);
  return par[8]*exp(par[9]*pt)*(par[10]+par[11]*tmp*tmp+
				par[12]*tmp*tmp*tmp*tmp);
}

//------------------------------------------------------------------------------

Float_t TStntuple::ZCorrMPi(const Float_t par[6], Float_t z) 
{
  Float_t tmp = z*z*par[1]*par[1];
  return par[0]*(par[2]+par[3]*tmp+par[4]*tmp*tmp+par[5]*tmp*tmp*tmp);
}

//------------------------------------------------------------------------------