#include "TROOT.h"
#include "TFolder.h"

#include <BaBar/Experiment.hh>
#include <Edm/EventRecord.hh>
#include <Edm/ConstHandle.hh>


#include <evt/Event.hh>

#include "Stntuple/mod/InitStntupleDataBlocks.hh"
#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/obj/TStnPhotonBlock.hh"
#include "Stntuple/obj/TStnClusterBlock.hh"
#include "Stntuple/obj/TStnLinkBlock.hh"
#include "ElectronObjects/CdfEmObjectColl.hh"
#include "HighLevelObjects/PhotonVariables.hh"
#include "HighLevelObjects/TrackingVariables.hh"
#include "Stntuple/obj/TStnDBManager.hh"
#include "Stntuple/data/TStnBeamPos.hh"
#include "Stntuple/mod/StntupleSiExpected.hh"
#include "TrackingSI/Utils/SiExpected.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Utils/HelixFitting.hh"

//_____________________________________________________________________________
Int_t StntupleInitPhotonBlock(TStnDataBlock* block, AbsEvent* event, int mode) 
{
  // initialize photon data block starting from event (from default collection
  // of EM objects. Provide interface for non-default initialization - `opt'
  // may have various meanings, but don't use it right away
  // using the functions defined within the HighLevelObjects package to fill
  // Photon variables.  See HighLevelObjects/photonVariables.hh

  using namespace PhotonVariables;
  using namespace TrackingVariables;


  static TFolder*         fol    = NULL;
  static StntupleSiExpected siexpected("PhotonBlockSiExpected");
  

  //---------------------------------------------------------------------------
  //  initialize static variables
  //--------------------------------------------------------------------------
  if (! fol) {
    fol = (TFolder*) gROOT->GetRootFolder()->FindObject("StntupleFolder");
    fol->Add((TObject*) &siexpected);
  }
  siexpected.BeginRun();
  siexpected.BeginEvent(event);

  double CprX, CprQ;
  double eUMax2,eVMax2,ePprTowMax;

  TStnPhotonBlock* data = (TStnPhotonBlock*) block;
  data->Clear();

  // pick up default primary vertex
  // ZVertexColl then vxprim
  double zVertex = trkZVertexZ();
  if(zVertex< -200.0) zVertex = trkPrimaryVertexZ();
  if(zVertex< -200.0) zVertex = 0.0;

  CdfEmObjectColl::const_handle h;
  
  if (CdfEmObjectColl::find(h)) {
    const CdfEmObjectColl::EmObjectList& list = h->contents();


    CdfEmObjectColl::const_iterator  it;
    for (it = list.begin(); it != list.end(); ++it) {
      TStnPhoton* photon = data->NewPhoton();

      const CdfEmObject& pho = *it;

      photon->fVersion  = phoVariablesVersion();

      HepLorentzVector momentum = phoFourVector(pho,zVertex);
      HepLorentzVector* mom = &momentum;
  
      photon->fMomentum.SetX(mom->px());
      photon->fMomentum.SetY(mom->py());
      photon->fMomentum.SetZ(mom->pz());
      photon->fMomentum.SetT(mom->t());
      photon->fWrd      = phoCutWord(pho, zVertex);
      double status = phoStatus(pho, zVertex);
      photon->fStat     = status;
      photon->fZv       = zVertex;
      int det = phoDetector(pho);
      photon->fDetector = det;
      photon->fEt       = phoRawEt(pho, zVertex);
      photon->fEtc      = phoCorrectedEt(pho, zVertex);
      photon->fE        = phoRawEnergy(pho);
      photon->fCesslide = phoCesSlide(pho);
      photon->fCo4      = phoCo4(pho, zVertex);
      photon->fCo4PJW   = phoCo4PJW(pho, zVertex);
      photon->fCo7      = phoCo7(pho, zVertex);
      photon->fSumpt4  	= phoTrackIso(pho, 0.4, zVertex);  // cone, vertex
      photon->fCesx     = phoCesX(pho);
      photon->fCesz     = phoCesZ(pho);
      photon->fChi      = phoCesChiSq(pho);
      photon->fChi3x3   = phoPem3x3ChiSq(pho);

      photon->fChieta   = phoPesChi2ClusterU(pho);
      photon->fChiphi   = phoPesChi2ClusterV(pho);

      photon->fChistr   = phoCesStrpChiSq(pho);
      photon->fChiwir   = phoCesWireChiSq(pho);
      photon->fPhi      = phoPhi(pho);
      photon->fDteta    = phoDteta(pho);
      photon->fEveta    = phoEta(pho, zVertex);
      photon->fSth      = phoSinTheta(pho, zVertex);
      photon->fCese     = phoCesE(pho);
      photon->fCeses    = phoCesStripE(pho);
      photon->fHadem    = phoHadEm(pho);
      photon->fLshr     = phoLshr(pho, zVertex);
      photon->fStr2     = phoCesZ2(pho);
      photon->fStre2    = phoCesStripE2(pho);
      photon->fWwir2    = phoCesX2(pho);
      photon->fWire2    = phoCesWireE2(pho);
      photon->fPt       = phoPt(pho,1);
      photon->fPt2      = phoPt(pho,2);
      photon->fN3d      = phoN3D(pho);
      phoCprUnbiased(pho, CprX, CprQ);
      photon->fCpr5ps = CprX; 
      photon->fCpr5ph = CprQ;

      //phoCprCesSeeded(pho, CprX, CprQ, 5, zVertex);
      //photon->fCescprx = CprX;
      //photon->fCescprq = CprQ;
      //photon->fCprwht = phoCprBgWeight(pho,zVertex);
      // filling these is now done in StntupleMakerModule
      // in order to deal with the Cpr database access
      // which shouldn't be done on an event-by-event basis
      photon->fCescprx = -99.0;
      photon->fCescprq = -99.0;
      photon->fCprwht = -999.0;
      photon->fTkcprx = phoCprTrack(pho, zVertex);

      photon->fCeswht = phoCesBgWeight(pho, zVertex);
      photon->fCespg = phoCesProbGam(pho, zVertex);
      photon->fCespb = phoCesProbBgd(pho, zVertex);
      photon->fCprpg = 0.0; // filled in StntupleMakerModule
      photon->fCprpb = 0.0;
      photon->fTime  = phoTime(pho);
      photon->fVcor  = phoVCor(pho);
      photon->fLcor  = phoLCor(pho);

      photon->fPesE = phoPesE(pho);
      photon->fPesX = phoPes2dX(pho);
      photon->fPesY = phoPes2dY(pho);
      photon->fPesE2= phoPesE2(pho,eUMax2,eVMax2);
      photon->fPprE = phoPprEnergy(pho,ePprTowMax);

      photon->fTowInd = (phoSeedPhi(pho)<<24) | 
	                (phoSeedEta(pho)<<16) | phoTowerPattern(pho);
      photon->fTowEt = phoTowerEt(pho);
      photon->fTowPhi = phoTowerPhi(pho);
      photon->fTowEta = phoTowerEta(pho);
      photon->fTowIso = phoCo4(pho,0.0);
      photon->fHademT = phoHadEmT(pho);

      photon->fLike   = phoLikelihood(pho);

      int nCosStub=0,nCosStubPho=0,nStub=0;
      nCosStub        = phoCosStub(pho,nCosStubPho,nStub,3);
      photon->fCosStub = (nCosStub<<16) | (nCosStubPho<<8) | nStub;

      int nIsoTrk;
      double maxIsoTrk;
      phoTrackIsoVar(pho,0.4,zVertex,nIsoTrk,maxIsoTrk);
      photon->fNIsoTrk = nIsoTrk;
      photon->fMaxIsoTrk = maxIsoTrk;

      int seedwedge,sidewedges,ncontig,east,west;
      phoBeamHaloOccupancy(pho, seedwedge, sidewedges, ncontig);
      phoBHHadOccupancy(pho, east, west);
      photon->fHaloVar = (seedwedge<<24) | (sidewedges<<18) | (ncontig<<12) |
	(east<<6) | west;

      photon->fPemPesDr = phoPesDeltaR(pho);
      photon->fSpare1 = 0.0;

      CdfEmObject::Cpr2EList     elist = pho.cpr2PadEnergies();
      CdfEmObject::Cpr2PadNoList ilist = pho.cpr2PadNumbers();
      int ncpr2 = elist.size();
      if(ilist[0]==0 && ilist[1]==0 && ilist[2]==0 && ilist[3]==0) ncpr2=0;
      photon->fCpr2I = (ncpr2<<24);
      photon->fCpr2E0 = photon->fCpr2E1 = 0;
      if(ncpr2>0) { photon->fCpr2E0 |= (elist[0]&0xFFFF); 
                           photon->fCpr2I  |= (ilist[0]); }
      if(ncpr2>1) { photon->fCpr2E0 |= (elist[1]<<16); 
                           photon->fCpr2I  |= (ilist[1]<<6); }
      if(ncpr2>2) { photon->fCpr2E1 |= (elist[2]&0xFFFF); 
                           photon->fCpr2I  |= (ilist[2]<<12); }
      if(ncpr2>3) { photon->fCpr2E1 |= (elist[3]<<16); 
                           photon->fCpr2I  |= (ilist[3]<<18); }

      // save pointers for links
      photon->fWirep = phoCesWireCluster(pho);
      photon->fStrpp = phoCesStripCluster(pho);
      photon->fPesUp = phoPesClusterU(pho);
      photon->fPesVp = phoPesClusterV(pho);
      photon->fCprcp = phoCprUnbiasedCluster(pho);
      photon->fCdfEmObject = &pho;

      // fill in the counters
      if(det==0) (data->fNcencl)++;
      if(status>2000.0) (data->fNphoidl)++;
      if(status>3000.0) (data->fNphoisol)++;
      if(status>4000.0) (data->fNphoid)++;
      if(status>5000.0) (data->fNphoiso)++;

    }
  }

  //now sort them according to status word
  data->PhotonList()->Sort();

  return 0;
}

//_____________________________________________________________________________
Int_t StntuplePhotonBlockLinks(TStnDataBlock* block,AbsEvent* event,int mode) 
{
  // set links between the photons and the objects in the other blocks
  // return 0 if everything is OK

  //SiExpected makes db calls so have only one permanent copy 
  static TFolder*         fol       = NULL;
  static SiExpected*      siexpected= NULL;

  if (! fol) {
    fol     = (TFolder*) gROOT->GetRootFolder()->FindObject("StntupleFolder");
    StntupleSiExpected*  stnsiexp= (StntupleSiExpected*)  
                        fol->FindObject("PhotonBlockSiExpected");
    if(stnsiexp) {
      siexpected = stnsiexp->GetSiExpected();
    } else {
      std::cout << "Could not find PhotonBlockSiExpected" << std::endl;
    }
  }


  char message[200];

  Int_t rc = 0;

  TStnPhotonBlock* pho_block = (TStnPhotonBlock*) block;
  Int_t npho = pho_block->NPhotons();
  if(npho<=0) return rc;
  
  TStnPhoton*        pho;
  const CdfEmObject*   emObjectp;
  
  TStnClusterBlock* clu_block = (TStnClusterBlock* ) 
    block->GetEvent()->GetDataBlock("ClusterBlock");


  // find the vertex position
  double vx, vy;
  bool vrtOk = true;
  double vz = pho_block->Photon(0)->VertexZ();
  if(vz<-998.0) vrtOk = false;

  TStnDBManager* dbm = TStnDBManager::Instance();
  TStnBeamPos* svx_beam_pos = (TStnBeamPos*) dbm->GetTable("SvxBeamPos");
  TStnBeamPos* cot_beam_pos = (TStnBeamPos*) dbm->GetTable("CotBeamPos");
  if(svx_beam_pos && svx_beam_pos->Status()==0 ) {
    vx = svx_beam_pos->X0() + vz*svx_beam_pos->DxDz();
    vy = svx_beam_pos->Y0() + vz*svx_beam_pos->DyDz();
  } else if (cot_beam_pos && cot_beam_pos->Status()==0) {
    vx = cot_beam_pos->X0() + vz*cot_beam_pos->DxDz();
    vy = cot_beam_pos->Y0() + vz*cot_beam_pos->DyDz();
  } else {
    vrtOk = false;
  }
  
  for (int ipho=0; ipho<npho; ipho++) {
    pho = pho_block->Photon(ipho);
    
    if (clu_block) {
      int ic;
      bool found;

      if(pho->Detector()==0) {
	if (pho->GetCesWireCluster()) {
	  found = false;
	  for (ic=0; ic<clu_block->NCesUnbClusters(); ic++) {
	    TCesCluster* clu = clu_block->CesUnbCluster(ic);
	    if (pho->GetCesWireCluster() == clu->GetCesCluster()) {
	      pho->SetCesWireCluster(ic);
	      found = true;
	    }
	  }
	  if(!found) {
	    sprintf(message,"no CES wire cluster matches photon %d",ipho); 
	    ERRLOG(ELwarning,message) << endmsg;
	  }
	} else {
	  pho->SetCesWireCluster(0xFFFF);
	}
	
	if (pho->GetCesStrpCluster()) {
	  found = false;
	  for (ic=0; ic<clu_block->NCesUnbClusters(); ic++) {
	    TCesCluster* clu = clu_block->CesUnbCluster(ic);
	    if (pho->GetCesStrpCluster() == clu->GetCesCluster()) {
	      pho->SetCesStrpCluster(ic);
	      found = true;
	    }
	  }
	  if(!found) {
	    sprintf(message,"no CES strip cluster matches photon %d",ipho); 
	    ERRLOG(ELwarning,message) << endmsg;
	  }
	} else {
	  pho->SetCesStrpCluster(0xFFFF);
	}

	if (pho->GetCprCluster()) {
	  found = false;
	  for (ic=0; ic<clu_block->NCprUnbClusters(); ic++) {
	    TCprCluster* clu = clu_block->CprUnbCluster(ic);
	    if (pho->GetCprCluster() == clu->GetCprCluster()) {
	      pho->SetCprCluster(ic);
	      found = true;
	    }
	  }
	  if(!found) {
	    sprintf(message,"no PES U cluster matches photon %d",ipho); 
	    ERRLOG(ELwarning,message) << endmsg;
	  }
	} else {
	  pho->SetCprCluster(0xFFFF);
	}

      } else { // central or else plug
      
	if (pho->GetPesUCluster()) {
	  found = false;
	  for (ic=0; ic<clu_block->NPesClusters(); ic++) {
	    TPesCluster* clu = clu_block->PesCluster(ic);
	    if (pho->GetPesUCluster() == clu->GetPesCluster()) {
	      pho->SetPesUCluster(ic);
	      found = true;
	    }
	  }
	  if(!found) {
	    sprintf(message,"no PES U cluster matches photon %d",ipho); 
	    ERRLOG(ELwarning,message) << endmsg;
	  }
	} else {
	  pho->SetPesUCluster(0xFFFF);
	}
	
	if (pho->GetPesVCluster()) {
	  found = false;
	  for (ic=0; ic<clu_block->NPesClusters(); ic++) {
	    TPesCluster* clu = clu_block->PesCluster(ic);
	    if (pho->GetPesVCluster() == clu->GetPesCluster()) {
	      pho->SetPesVCluster(ic);
	      found = true;
	    }
	  }
	  if(!found) {
	    sprintf(message,"no PES U cluster matches photon %d",ipho); 
	    ERRLOG(ELwarning,message) << endmsg;
	  }
	} else {
	  pho->SetPesVCluster(0xFFFF);
	}

      } // if central or else plug
     
    } // if cluster block exists
    

   
//----------------------------------------------------------------------
// fill tower link block
//----------------------------------------------------------------------

    emObjectp=  pho->GetCdfEmObject();

    const std::vector<TowerKey>& towlist = 
                emObjectp->getEmCluster()->towerIds();
   
    int ntowers = towlist.size();
    for (int itower=0; itower<ntowers; itower++) {
      const TowerKey pierre_key =  towlist[itower];

      int tower_key = (((pierre_key.iPhi() & 0x3f  )      ) |
		       ((pierre_key.iEta() & 0x3f  ) <<  8) |
		       ((pierre_key.isValid()) << 16)    );

      pho_block->TowerLinkList()->Add(ipho,tower_key);
    }


//----------------------------------------------------------------------
// fill SiExpected
//----------------------------------------------------------------------
    pho->fSiExpected = 0;

    if(vrtOk) {
      siexpected->setCountHitsCone( true );
      siexpected->setCone( 0.4 );
      std::vector<SiExpected::Info> exphits;
      exphits.clear();
      int nSiHitsCone = 0;
      int nSiHitsEx = 0;
      int nSiPhiHitsEx = 0;

      HepVector3D momentum(pho->Momentum()->Px(),
			   pho->Momentum()->Py(),pho->Momentum()->Pz());
      // give it a huge momentum so tracking doesn't bend it
      momentum *= 1.e4;
      HepPoint3D vertex(vx,vy,vz);
      Helix helix(momentum,vertex,1.0,14.0);
      HepSymMatrix cov(5);
      cov[0][0] = pow(0.01,2);   // lambda
      cov[1][1] = pow(1.0e-8,2); // c, make it small
      cov[2][2] = pow(0.05,2);   // z
      cov[3][3] = pow(0.005,2);  // d
      cov[4][4] = pow(0.005,2);  // phi
      
      CdfTrack* trk = new CdfTrack();
      HelixFit fit(0.,1,helix.getParameters(),cov);
      trk->setHelixFit(fit);

      if ( siexpected->findExpected ( *trk, exphits )) {
	for ( std::vector<SiExpected::Info>::iterator
		expHitIter=exphits.begin(), expHitIterend = exphits.end () ;
	      expHitIter != expHitIterend ; ++expHitIter )
	  {
	    if ( (*expHitIter).integrated ) {
	      nSiHitsEx++;
	      if ((*expHitIter).isPhi) nSiPhiHitsEx++;
	      nSiHitsCone+=(*expHitIter).nHitsCone;
	    }
	  }
      }

      if(nSiHitsCone>4096) nSiHitsCone=4096;
      pho->fSiExpected = 
	( (nSiHitsEx) | (nSiPhiHitsEx<<8) | (nSiHitsCone<<16) );

      trk->destroy();

    } // endif vrtOk

  } // end loop over photons

  return rc;
}