////////////////////////////////////////////////////////////////////
//
//     File: PhotonVariables3.cc
//     Authors: Mike Kirby, Ray Culbertson
//     Description: High level analysis photon variables defined as functions
//     Created: 2/1/2001
//
//     Revision: 0.0
//     Revision: 1.0 Add CES and CPR bg weights, Ray Culbertson, 2/17/01
//
////////////////////////////////////////////////////////////////////

#include "HighLevelObjects/PhotonVariables.hh"
#include "HighLevelObjects/PhotonBackgroundComputer.hh"
#include "HighLevelObjects/TrackingVariables.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "ElectronObjects/CdfEmObject.hh"
#include "AbsEnv/AbsEnv.hh"
#include "RawDataBanks/CESD_StorableBank.hh"
#include "ShowerMax/CesDefs.hh"
#include "Calor/CprWirePar.hh"
#include "Calor/CprClustPar.hh"
#include "Calor/CprClusterMaker.hh"
#include "Calor/CprWireCollectionMaker.hh"
#include "CalorObjects/CalData.hh"
#include "Electron/emobj_alg.hh"
#include "MuonObjects/CdfMuonView.hh"
#include "MuonObjects/MuonStubColl.hh"
#include "JetObjects/CdfJetColl.hh"
#include "MetObjects/CdfMet.hh"
#include "Electron/IsoCorrAlg.hh"

IsoCorrAlg::Error plugisocor(const CdfEmObject& emobj, double& leakf);

namespace PhotonVariables {

  // this is Run I sliding cut on the energy of the second 
  // CES wire cluster in the wedge
  double phoCesSlide(const CdfEmObject& emObject){
    if(phoDetector(emObject)!=CEM) return -1;
    double e = phoCorrectedEnergy(emObject);
    if(e<18.0) {
      return 0.14*e;
    } else {
      return  2.4 + 0.01*e;
    }
  }

  // this returns the Run I correction to isolation energy
  // due to phi-crack leakage, from Peter Wilson
  // used:  CorIsoEtCone0.4 = IsoEtCone0.4 - LCor*PhotonCorrectedEt
  double phoLCor(const CdfEmObject& emObject) {
    if(!phoHasShowerMax(emObject)) return 0.0;
    if(phoDetector(emObject)==CEM) {
      //From fit to z->ee data (Run 1B) 10/13/97
      const double p10=0.000280, p11=0.407;
      double xces = phoCesX(emObject);
      double corr = 1.0/(1.0+p10*exp(p11*(48.4-fabs(xces))));
      // samll additional correction which is not a fn of x
      double par3 = 0.002 + 0.003*(phoCorrectedEt(emObject)/400.0);
      return corr + par3;
    } else {
      double leakf;
      IsoCorrAlg::Error status = plugisocor(emObject, leakf);
      return (status == IsoCorrAlg::OK ? leakf : 0.0);
    }
  }

  // this returns the Run I correction to isolation energy
  // due to multiple underlying events, from Peter Wilson
  // user substracts this energy: CorIsoEtCone0.4 = IsoEtCone0.4 - LCor;
  double phoVCor(const CdfEmObject& emObject) {
    //these constants are from Fall 2004 
    // revision of the ZVertexModule, standard for gen5
    int nVert=TrackingVariables::trkNGoodZVertex();

    if(nVert<=1) return 0.0;
    return  0.3563 * (nVert-1);
  }

  // this returns the Run I cone 0.4 isolation energy
  // with VCor and Lcor corrections (see above)
  double phoCo4PJW(const CdfEmObject& emObject, const double zvertex) {
    double isoBase = phoCo4(emObject, zvertex);
    double isoCorr = isoBase - phoVCor(emObject) - 
      phoLCor(emObject)*phoCorrectedEt(emObject, zvertex);
    return (isoCorr>0.0 ? isoCorr : 0.0);
  }

  // Ces energy for all wires in the wedge, call with
  // double cesWireArr[64]
  // the first 32 are the wires closest to |z|=0, 32-63 are the ones 
  // further from zero.  Order is decreasing phi on west, increasing on east.
  // double cesStripArr[128]   order is from low |z| to high
  // return value is 0 for all OK
  int phoCesStripAndWireE(const CdfEmObject& emObject, 
			  double* cesWireArr, double* cesStripArr) {
 
    for(int i=0; i<64; i++) cesWireArr[i] = 0;

    int emBarrel = emObject.getEmCluster()->side();
    int emWedge  = emObject.getEmCluster()->seedPhiAddr();

    AbsEvent* anEvent = AbsEnv::instance()->theEvent();

    EventRecord::ConstIterator iter(anEvent, "CESD_StorableBank");

    if (! iter.is_valid()) return 1;

    ConstHandle<CESD_StorableBank> cesd(iter);

    int wedge, barrel, stripOrWire, stripNo, energy;

    for (CESD_StorableBank::ConstBankIter blk(cesd); blk.is_valid(); ++blk) {
      for (CESD_StorableBank::ConstBlockIter ch(blk); ch.is_valid(); ++ch) {

        wedge       = cesd->get_module(ch);
        barrel      = cesd->get_we   (blk);
        stripOrWire = cesd->get_strip  (ch);
        stripNo     = cesd->get_stripNo(ch);
	energy      = cesd->get_data(ch);

	if (wedge==emWedge && barrel==emBarrel) {
	  if (stripOrWire == 0) {
	    cesWireArr[stripNo]  = float(energy-PEDESTAL_CES)/GEV_TO_ADC_WIRE;
	  } else {
	    cesStripArr[stripNo] = float(energy-PEDESTAL_CES)/GEV_TO_ADC_STRIP;
	  }
	}
      }
    }
    return 0;
  }


  // the CPR local x and charge of the highest energy cluster in this wedge
  // return value is 0 if all OK
  // cpr local x is in +phi sense
  int phoCprUnbiased(const CdfEmObject& emObject, 
			double& CprX, double& CprQ) {
    CprX = -99;
    CprQ = -99;
    const CprCluster* clus = phoCprUnbiasedCluster(emObject);
    if(clus==NULL) return 1;
    CprQ =  clus->energy();
    CprX =  clus->localCoord();

    return 0;
  }

  // find the CPR unbiased cluster
  const CprCluster* phoCprUnbiasedCluster(const CdfEmObject& emObject) {
    double CprQ = -99;
    const CprCluster* p = NULL;
    int emBarrel = emObject.getEmCluster()->side();
    int emWedge  = emObject.getEmCluster()->seedPhiAddr();

    AbsEvent* anEvent = AbsEnv::instance()->theEvent();

    EventRecord::ConstIterator iCprClusterColl( anEvent, "CprClusterColl");
    for ( ; iCprClusterColl.is_valid(); ++iCprClusterColl) {
      CprClusterColl_ch   cprClusterColl_ch( iCprClusterColl );
      if(cprClusterColl_ch->description() == "SeedBased_CprClusterCollection"){
	CprClusterColl::const_iterator ri;
	for (ri  = cprClusterColl_ch->contents().begin();
	     ri != cprClusterColl_ch->contents().end();  ++ri) {
	  if( ri->side()==emBarrel && ri->module()==emWedge) {
	    if(ri->energy() > CprQ) {
	      CprQ =  ri->energy();
	      p = &(*ri);
	    }
	  }
	}  // end loop over clusters

      }  // end if on collection description
    } // end loop over CPR collections

    return p;
  }

  // This returns the CPR cluster seeded from the photon's CES cluster
  // two variables here in order to not have to do the clustering twice
  // nstrips is the number of the Cpr strips to consider (not implemented)
  // cpr local x is in +phi sense
  int phoCprCesSeeded(const CdfEmObject& emObject, double& CprX, 
		      double& CprQ, const int nstrips, const double zvertex,
		      CprWireCollectionMaker* pMaker) {

    CprWireCollectionMaker* pMakerLoc = NULL;
    CprQ = -99.0;
    CprX = -99.0;
    if(!phoHasShowerMax(emObject) || phoDetector(emObject)!=CEM) return 1;

    // values such as the numbers of wires to cluster is in here
    // but not settable
    CprWirePar  wirePar(NULL);
    CprClustPar clustPar(NULL);

    if(pMaker!=NULL) {
      pMakerLoc = pMaker;
    } else {
      pMakerLoc = new CprWireCollectionMaker();
      pMakerLoc->initJob(false);
      pMakerLoc->initRun(false);
      AbsEvent* anEvent = AbsEnv::instance()->theEvent();
      pMakerLoc->initCollection(anEvent);
    }

    const CesCluster* pCluWire = phoCesWireCluster(emObject);
    const CesCluster* pCluStrp = phoCesStripCluster(emObject);
    // method uses vector iterators
    vector<CesCluster> wireclus;  wireclus.push_back(*pCluWire);
    vector<CesCluster> strpclus;  strpclus.push_back(*pCluStrp);
    vector<CesCluster>::const_iterator is=strpclus.begin();
    vector<CesCluster>::const_iterator iw=wireclus.begin();

    HepPoint3D primvtx(0,0,zvertex);

    const vector<CprWire> cluWires = 
               pMakerLoc->get_ces_based_wires(&wirePar, is, iw, primvtx);

    int nWires = cluWires.size();
    if (nWires<1) {
      if(pMaker==NULL) delete pMakerLoc;
      return 2;
    }

    CprClusterMaker doCluster;
    doCluster.initCprClusterMaker();

    doCluster.setSeedWire(pMakerLoc->get_SeedWireNo());

    CdfTrack_clnk emptyTrackLink;
    CprCluster myCluster = 
      doCluster.do_cluster_wires( &cluWires, &clustPar, emptyTrackLink);
    if (!doCluster.is_valid()) {
      if(pMaker==NULL) delete pMakerLoc;
      return 3;
    }
    CprX = myCluster.localCoord();
    CprQ = myCluster.energy();

    if(pMaker==NULL) delete pMakerLoc;

    return 0;
  }

  // these return the CES and CPR background weights
  // the CES weight is a positive number for events with avg chi^2 <4,
  // and negative for events with 4<chi^2<20
  // the CPR weight is positive for event without a cpr clusters
  // and neagitve otherwise.  There are details to using these such as the
  // ID and fiducial cuts on the photon.  returns -99 for plug photons
  double phoCesBgWeight(const CdfEmObject& emObject, const double zvertex) {
    if(!phoHasShowerMax(emObject) || phoDetector(emObject)!=CEM) return -99.0;
    double et    = phoCorrectedEt(emObject,zvertex);
    double chisq = phoCesChiSq(emObject);
    return PhotonBackgroundComputer::cesWeight(et,chisq);
  }

  double phoCprBgWeight(const CdfEmObject& emObject, const double zvertex) {

    if(phoDetector(emObject)!=CEM) return -99.0;
    double cprx,cprq;
    int stat = phoCprCesSeeded(emObject, cprx, cprq, 5, zvertex);
    if(stat==1) {    // no CES cluster so can't find CPR cluster
      return -99.0;
    } else if(stat==2 || stat==3) {      // cluster not found b/c no hits
      cprq = 0.0;
    } else if(stat==0) {
      // sucessful
    } else {
      std::cerr <<"phoCprBgWeight found an unknown CPR error code" << std::endl;
    }
    double et  = phoCorrectedEt(emObject,zvertex);
    double sth = phoSinTheta(emObject,zvertex);
    return PhotonBackgroundComputer::cprWeight(cprq,et,sth);
  }


  double phoCesProbGam(const CdfEmObject& emObject, const double zvertex) {
    if(!phoHasShowerMax(emObject) || phoDetector(emObject)!=CEM) return -99.0;
    double et    = phoCorrectedEt(emObject,zvertex);
    double chisq = phoCesChiSq(emObject);
    return PhotonBackgroundComputer::phoEff(et);
  }
  double phoCesProbBgd(const CdfEmObject& emObject, const double zvertex) {
    if(!phoHasShowerMax(emObject) || phoDetector(emObject)!=CEM) return -99.0;
    double et    = phoCorrectedEt(emObject,zvertex);
    double chisq = phoCesChiSq(emObject);
    return PhotonBackgroundComputer::backEf(et);
  }
  double phoCprProbGam(const CdfEmObject& emObject, const double zvertex) {
    // currently this can't be done as a function
    // for Stntuple, done in StntupleMakerModule
    return -99.0;
  }
  double phoCprProbBgd(const CdfEmObject& emObject, const double zvertex) {
    // currently this can't be done as a function
    // for Stntuple, done in StntupleMakerModule
    return -99.0;
  }
  // quantity needed for Cpr method
  double phoSumetNoJets(){

    AbsEvent* anEvent = AbsEnv::instance()->theEvent();

    CdfMet_ch met;
    CdfMet::Error err = CdfMet::find(met);
    if (err == CdfMet::ERROR) return -999.0;
    double sumet = met->etSum();

    vector<std::string> processNames(3);
    processNames[0] = anEvent->process_name();
    processNames[1] = "PROD";
    processNames[2] = "L3";


    //get list of cone 0.7 jets
    //first the newer collection name with the "-cone0.7"

    // first current process, then PROD, then L3, then any
    EventRecord::ConstIterator jetCIter;
    for (int i = 0; i < 3; i++) {
      StorableObject::SelectByProcessName procsel( processNames[i] );
      EventRecord::ConstIterator jetCIter1(anEvent, 
	       StorableObject::SelectByClassName("CdfJetColl") &&
	       StorableObject::SelectByDescription("JetCluModule-cone0.4") &&
					      procsel);
      if (jetCIter1.is_valid()) {
	jetCIter = jetCIter1;
	break;
      }
    }
    // now look for one with any process name
    if (jetCIter.is_invalid()) {
      EventRecord::ConstIterator jetCIter1(anEvent, 
	      StorableObject::SelectByClassName("CdfJetColl") &&
	      StorableObject::SelectByDescription("JetCluModule-cone0.4"));
      jetCIter = jetCIter1;
    }

    CdfJetColl::const_iterator jetIter;
    ConstHandle<CdfJetColl> jetColl(jetCIter);
    std::string usedproc("PROD");
    if(jetColl.is_nonnull()) usedproc = jetColl->process_name();


    if (jetColl.is_nonnull()) {
      for(jetIter=jetColl->contents().begin();
	  jetIter!=jetColl->contents().end(); jetIter++) {
	if(jetIter->et()>10.0) sumet -= jetIter->et();
      }
    }
    
    return sumet;

  }


  // distance between expected Cpr hit and nearest track
  double phoCprTrack(const CdfEmObject& emObject, const double zvertex) {

    double cesx = phoCesX(emObject);
    double cesz = phoCesZ(emObject);
    if(fabs(cesz)>250.0) return 0.0; // failures will fail the cut
    if(fabs(cesx)> 30.0) return 0.0;
    int emWedge  = emObject.getEmCluster()->seedPhiAddr();
    
    double cprz = 0.91*(cesz-zvertex) + zvertex;
    int module;
    if(cprz<0.0) {
      module = (cprz<-122.0? 0 : 1 );
    } else {
      module = (cprz< 122.0? 2 : 3 );
    }
    
    //cprx is in +phi sense, opposite to Ces x
    double cprx = -0.91*cesx;
    
    CdfTrackView_h trview;
    if(CdfTrackView::defTracks(trview) != CdfTrackView::OK) return 0.0;
    
    CprWireCollectionMaker mycol;
    
    double dxmin = 99.0;
    for(CdfTrackView::const_iterator itrack=trview->contents().begin();
        itrack!=trview->contents().end(); ++itrack) {
      const CdfTrack_clnk & theTrack = *itrack;
      int stat = mycol.extrapolate_track(theTrack);
      if(stat==1) {
	double x =-mycol.get_Xlocal(); // extrap returns x in -phi sense
	double z = mycol.get_Zlocal(); // extrap returns local z
	if(mycol.get_Side()==0) z = -z;
	if(mycol.get_Wedge()==emWedge) {
	  int modulet;
	  if(z<0.0) {
	    modulet = (z<-122.0? 0 : 1 );
	  } else {
	    modulet = (z< 122.0? 2 : 3 );
	  }
	  if(modulet==module) {
	    if(fabs(x-cprx)< fabs(dxmin)) dxmin = x - cprx;
	  }
	}
      }
    } // end track loop
    return dxmin;
    
  }



  // returns emObject time, averaged over all firing TDCs
  // -9999 if no TDC fired for the cluster
  double phoTime(const CdfEmObject& emObject) {
    int ncha, nwha, npha;
    double chatime = phoDetTime(emObject, CHA, ncha);
    double whatime = phoDetTime(emObject, WHA, nwha);
    double phatime = phoDetTime(emObject, PHA, npha);
    
    if (ncha + nwha + npha == 0) 
      return -9999;
    else {
      if (ncha > 0) return chatime;
      if (npha > 0) return phatime;
      // if we get here we have only whatime
      return whatime;
    }
  }

  // returns emObject time, averaged over all firing TDCs in detector det
  // -9998 if det != CHA, WHA or PHA
  // -9999 if no timing in det
  double phoDetTime(const CdfEmObject& emObject, Detector det, int& ntow) {
    double sum_times = 0;
    double itmp;

    ntow = 0;

    if ((det != CHA) && (det != WHA) && (det != PHA)) return -9998;
    
    const EmCluster& emCluster = *(emObject.getEmCluster());
    const EmCluster::TowerIdList& towids = emCluster.towerIds();
    CalData_ch caldata_hndl;

    CalData::Error result = CalData::find(caldata_hndl);

    if (result != CalData::OK) return -9999; 

    
    for (int i = 0; i < towids.size(); i++) {
      const CalTower& tower_ref = *caldata_hndl->tower(towids[i]);
      if (tower_ref.includesDetector(det)) {
	CellKey ck(tower_ref.iEta(), tower_ref.iPhi(), det);
	itmp = tower_ref.time(ck,0);
	if (itmp >-250 && itmp < 300 && 
            tower_ref.nTDCHits(ck) > 0 && tower_ref.nTDCHits(ck) <= 2) {
	  sum_times += itmp;
	  ntow++;
	}
      }
    }
    
    if (ntow == 0)
      return -9999;
    else
      return sum_times / ntow;
  }

  // returns a likelihood that this is consistent with a single EM cluster
  // Stephen Levy: October 6, 2003
  // Provide argument for likelihood version.
  // See Electron/Electron/emobj_alg.hh for details.
  double phoLikelihood(const CdfEmObject& emObject, LikelihoodVersion ver) {
    if(!phoHasShowerMax(emObject)) return -9.0;
    if(emObject.getEmCluster()->region()==CEM) {
      return centralEmIdLikelihood(emObject,ver);
    } else {
      return plugEmIdLikelihood(&emObject);
    }

  }

  // returns number of muon stubs not associated with jets or tracks
  // nPhoMatchStub is stubs within 30deg of photon and not with track or jet
  // nStub is all stubs
  // input argument muonSystem is  bit flag for the allowed muon systems:
  // bit 0=CMU, 1=CMP, 2=CMX, 3=BMU
  int phoCosStub(const CdfEmObject& emObject, int& nCosStubPho, 
		 int& nStub, int muonSystem) {

  float dphiCut = 0.5236; // 30deg
    int nCosStub = 0;
    nCosStubPho = 0;
    nStub = 0;

    AbsEvent* anEvent = AbsEnv::instance()->theEvent();

    // get full muons
    CdfMuonView_h muonView_handle;
    CdfMuonView::defMuons(muonView_handle);
    CdfMuonView::const_iterator muIter;

    vector<std::string> processNames(3);
    processNames[0] = anEvent->process_name();
    processNames[1] = "PROD";
    processNames[2] = "L3";


    //get list of cone 0.7 jets
    //first the newer collection name with the "-cone0.7"

    // first current process, then PROD, then L3, then any
    EventRecord::ConstIterator jetCIter;
    for (int i = 0; i < 3; i++) {
      StorableObject::SelectByProcessName procsel( processNames[i] );
      EventRecord::ConstIterator jetCIter1(anEvent, 
	       StorableObject::SelectByClassName("CdfJetColl") &&
	       StorableObject::SelectByDescription("JetCluModule-cone0.7") &&
					      procsel);
      if (jetCIter1.is_valid()) {
	jetCIter = jetCIter1;
	break;
      }
    }
    // now look for one with any process name
    if (jetCIter.is_invalid()) {
      EventRecord::ConstIterator jetCIter1(anEvent, 
	      StorableObject::SelectByClassName("CdfJetColl") &&
	      StorableObject::SelectByDescription("JetCluModule-cone0.7"));
      jetCIter = jetCIter1;
    }


    if (jetCIter.is_invalid()) {
      // now try to find the older collection desc w/o "-cone0.7"
      // first current process, then PROD, then L3, then any
      for (int i = 0; i < 3; i++) {
	StorableObject::SelectByProcessName procsel( processNames[i] );
	EventRecord::ConstIterator jetCIter1(anEvent, 
	      StorableObject::SelectByClassName("CdfJetColl") &&
	      StorableObject::SelectByDescription("JetCluModule") &&
					     procsel);
	if (jetCIter1.is_valid()) {
	  jetCIter = jetCIter1;
	  break;
	}
      }
      // now look for one with any process name
      if (jetCIter.is_invalid()) {
	EventRecord::ConstIterator jetCIter1(anEvent, 
	       StorableObject::SelectByClassName("CdfJetColl") &&
	       StorableObject::SelectByDescription("JetCluModule"));
	jetCIter = jetCIter1;
      }
      
    } // end if for trying older name


    CdfJetColl::const_iterator jetIter;
    ConstHandle<CdfJetColl> jetColl(jetCIter);
    std::string usedproc("PROD");
    if(jetColl.is_nonnull()) usedproc = jetColl->process_name();

    // phi of emObject
    float phip = emObject.getEmCluster()->emPhi();

    // get list of stubs
    EventRecord::ConstIterator imustubs;
    StorableObject::SelectByProcessName procsel(usedproc);
    imustubs = EventRecord::ConstIterator(anEvent, 
		    StorableObject::SelectByClassName("MuonStubColl") && 
					    procsel);

    while(imustubs.is_valid()) {

      ConstHandle<MuonStubColl> hSColl(imustubs);
      for(MuonStubColl::container_type::const_iterator
	        istub=hSColl->contents().begin();
	  istub != hSColl->contents().end(); istub++) {
	
	if( (istub->system()==MuonDigiCode::CMU && (muonSystem&1)>0 ) || 
	    (istub->system()==MuonDigiCode::CMP && (muonSystem&2)>0 ) ||
	    (istub->system()==MuonDigiCode::CMX && (muonSystem&4)>0 ) ||
	    (istub->system()==MuonDigiCode::BMU && (muonSystem&8)>0 )
	    ) {
	  if(istub->numAssHits()>=2) {
	    
	    nStub++;

	    float stubPhi = istub->spacePoint().phi();
	    if(stubPhi<0.0) stubPhi += 2*M_PI;
	    
	    bool match=false;
	    
	    if (muonView_handle.is_nonnull()) {
	      // loop over muons, reject skip stubs matched to tracks
	      for (muIter = muonView_handle->contents().begin(); muIter != 
		     muonView_handle->contents().end(); ++muIter ) {
		const CdfMuon& muon = **muIter;
		if( &(*muon.stub(istub->system()))
		    ==(&*istub) ) match=true;
	      }
	    }
	    if (jetColl.is_nonnull()) {
	      for(jetIter=jetColl->contents().begin();
		  jetIter!=jetColl->contents().end(); jetIter++) {
		float dphip = fabs(jetIter->phi() - phip);
		if(dphip>M_PI) dphip = 2*M_PI - dphip;
		if(dphip>dphiCut) { // if this jet is not the photon
		  float dphij = fabs(jetIter->phi() - stubPhi);
		  if(dphij>M_PI) dphij = 2*M_PI - dphij;
		  if(dphij<dphiCut) match=true;  // stub matched to a jet
		}
	      }
	    }
	    
	    if(!match) {
	      nCosStub++;
	      float dphip = fabs(phip - stubPhi);
	      if(dphip>M_PI) dphip = 2*M_PI - dphip;
	      if(dphip<dphiCut) nCosStubPho++;
	    }
	    
	  } // endif 2 hits
	} // endif system
      } // stub loop
      
      imustubs++;
    } // container loop

    return nCosStub;
  }


  // for a central emObject, returns:
  // seedwedge: number of central towers above threshold in same wedge as
  //            emObject seed tower
  // sidewedges: number of central towers above threshold in wedges on each
  //             side of seed tower wedge
  // ncontig: number of contiguous towers above threshold in same wedge as
  //          seed tower (there is allowed to be a gap of one tower)
  // threshold: E cutoff to consider a tower occupied
  void phoBeamHaloOccupancy(const CdfEmObject& emObject, int& seedwedge,
			    int& sidewedges, int& ncontig, float threshold) {
    const CalTower* tow;
    seedwedge = 0;
    sidewedges = 0;
    ncontig = 0;
    if (TOWER_TYPE[emObject.getEmCluster()->seedEtaAddr()] > 2) {
      return;
    }
    int seedphi = emObject.getEmCluster()->seedPhiAddr();
    int phiplusone = (seedphi == 23) ? 0 : seedphi + 1;
    int phiminusone = (seedphi == 0) ? 23 : seedphi - 1;

    CalData_ch data; 
    CalData::Error result = CalData::find(data); 
    if(result != CalData::OK) return; 

    int ieta;

    for (ieta = 16; ieta < 36; ieta++) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->energy() > threshold) {
	seedwedge++;
      }
      tow = data->tower(ieta, phiplusone);
      if (tow && tow->energy() > threshold) {
	sidewedges++;
      }
      tow = data->tower(ieta, phiminusone);
      if (tow && tow->energy() > threshold) {
	sidewedges++;
      }
    }

    bool iscontig = true;
    ieta = emObject.getEmCluster()->seedEtaAddr();
    while (iscontig && ieta < 36) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->energy() > threshold) {
	iscontig = true;
	ncontig++;
	ieta++;
      } else {
	iscontig = false;
	ieta++;
      }
    }

    iscontig = true;
    int contup = 0;
    while (iscontig && ieta < 36) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->energy() > threshold) {
	iscontig = true;
	contup++;
	ieta++;
      } else {
	iscontig = false;
      }
    }

    iscontig = true;
    ieta = emObject.getEmCluster()->seedEtaAddr() - 1;
    while (iscontig && ieta > 15) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->energy() > threshold) {
	iscontig = true;
	ncontig++;
	ieta--;
      } else {
	iscontig = false;
	ieta--;
      }
    }

    iscontig = true;
    int contdown = 0;
    while (iscontig && ieta > 15) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->energy() > threshold) {
	iscontig = true;
	contdown++;
	ieta--;
      } else {
	iscontig = false;
      }
    }
    
    // allow a one tower gap
    ncontig += contup > contdown ? contup : contdown;

  }


  //hadron tower occupancy in regions where a beam halo muon would cross
  // east and west are the number of hit towers on each side
  void phoBHHadOccupancy(const CdfEmObject& emObject, int& east, int& west) {
    west = 0; east = 0;
    if (TOWER_TYPE[emObject.getEmCluster()->seedEtaAddr()] > 2) return;
    
    int seedphi = emObject.matchingEmCluster()->seedPhiAddr();
    const CalTower* tow;
    
    CalData_ch data; 
    CalData::Error result = CalData::find(data); 
    if(result != CalData::OK) return; 
    
    for (int ieta = 15; ieta > 13; ieta--) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->pemEnergy() == 0 && tow->whaEnergy() > 0) {
	west++;
      }
    }
    for (int ieta = 13; ieta > 11; ieta--) {
      tow = data->tower(ieta, seedphi*2);
      if (tow && tow->pemEnergy() == 0 && tow->phaEnergy() > 0) {
	west++;
      }
      tow = data->tower(ieta, seedphi*2 + 1);
      if (tow && tow->pemEnergy() == 0 && tow->phaEnergy() > 0) {
	west++;
      }
    }
    for (int ieta = 35; ieta < 38; ieta++) {
      tow = data->tower(ieta, seedphi);
      if (tow && tow->cemEnergy() == 0 && tow->pemEnergy() == 0 && tow->whaEnergy() > 0) {
	east++;
      }
    }
  }

  int phoNTracksInitPoint(const CdfEmObject& emObject, std::string processName) {
    CdfTrackView_h trview;
    if(CdfTrackView::defTracks(trview,processName) != CdfTrackView::OK) {
      return -99;
    } else {
      return phoNTracksInitPoint(emObject, trview);
    }
  }

  int phoNTracksInitPoint(const CdfEmObject& emObject, CdfTrackView_h trview) {
    if (trview.is_null()) {
      return -99;
    }
    // algorithm: loop over tracks, for each check if straight line
    // with phi=phi0, lambda=lambda, z0=z0 would hit seed tower

    int retval = 0;
    TowerGeometry geo;
    TowerKey seedkey(emObject.matchingEmCluster()->seedEtaAddr(),
                     emObject.matchingEmCluster()->seedPhiAddr());

    for(CdfTrackView::const_iterator itrack=trview->contents().begin();
        itrack!=trview->contents().end(); ++itrack) {
      float phi0 = TrackingVariables::trkPhi0(&**itrack);
      float detlambda = TrackingVariables::trkLambda(&**itrack) +
        TrackingVariables::trkZ0(&**itrack)/TLYCES;
      if (phi0 > geo.phiMin(seedkey) && phi0 < geo.phiMax(seedkey)
          && detlambda > geo.cotThetaMin(seedkey)
          && detlambda < geo.cotThetaMax(seedkey)) {
        retval++;
      }
    }
    return retval;
  }




} // namespace PhotonVariables