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

#include "evt/Event.hh"
#include <Bbos/bank_interface.hh>

#include <iostream>

#include <JetObjects/CdfJet.hh>
#include <JetObjects/CdfJetColl.hh>

#include <HighLevelObjects/JetVariables.hh>

#include <Stntuple/data/TStnTypes.hh>
#include <Stntuple/mod/InitStntupleDataBlocks.hh>
#include <Stntuple/obj/TStnJetBlock.hh>
#include <Stntuple/obj/TStnTrackBlock.hh>
#include <Stntuple/obj/TStnVertexBlock.hh>
#include <Stntuple/mod/StntupleUtilities.hh>
#include <Stntuple/obj/TStnErrorLogger.hh>
#include <Stntuple/obj/TStnEvent.hh>


namespace {
  static Double_t         gConeSize = -1;
  static TFolder*         fgFolder  = NULL;
  static TStnErrorLogger* fgLogger  = NULL;

};

//_____________________________________________________________________________
Int_t StntupleInitJet(TStnJet* Pasha, const CdfJet* Pierre) {
  // initialize STNTUPLE jet starting from CDF RunII jet


  Pasha->fCdfJet = (CdfJet*) Pierre;

  Pasha->SetNTrackWord(-1,-1,-1,-1);

  Pasha->fMomentum.SetXYZT(Pierre->px(),Pierre->py(),
			   Pierre->pz(),Pierre->totEnergy());

  Pasha->fDetEta    = Pierre->etaDetector();
  Pasha->fEmfr      = Pierre->emFraction();

					// Et/sum(Pt) - wait until defined in 4.2.x
  Pasha->fEp        = -1;
					// status code not yet defined

  Pasha->fWord      = -1; // jetWord(*Pierre, tracks);

					// comment out until proved that the
					// functions behave

  Pasha->fCorfm      = -1; // jetCorFM(*Pierre,4,1,1,gConeSize);
  Pasha->fCorfd      = -1; // jetCorFD(*Pierre,4,1,1,gConeSize);
  Pasha->fCorfa      = -1.;

  Pasha->fTag        = -1;
  Pasha->fEgrd       = Pierre->guardEnergy();
  Pasha->fVbpb       = -1.e6;
  Pasha->fPbpb       = -1.e6;
  Pasha->fScbpb      = -1.e6;
  Pasha->fScnpb      = -1.e6;
  Pasha->fPppb       = -1.e6;
  Pasha->fPnpb       = -1.e6;
  Pasha->fTrackMass  = -1.e6;
  Pasha->fMaxTrackPt = -1.e6;
//-----------------------------------------------------------------------------
// Fill z-vertex and sigma The HLO function fills doubles, and we need floats
//-----------------------------------------------------------------------------
  CdfTrackView_h tracks;
  double zvtx = 0.;
  double dz   = 0.;
  JetVariables::jetTracksOriginal(*Pierre, tracks);
  JetVariables::jetZVertexInfo(*Pierre, tracks, zvtx, dz);
  Pasha->fZvtx  = zvtx;
  Pasha->fDz    = dz; 
  Pasha->fPtout = -1.; //jetPtOut(*Pierre, tracks);
//-----------------------------------------------------------------------------
// loop over the jet towers to find a seed
//-----------------------------------------------------------------------------
  int ieta, iphi;
  const vector<const PhysicsTower*>& towlist = Pierre->towers().contents();
  int ntow = towlist.size();

  int ntlzp = 0; // jetNTracksPositiveLz(*Pierre, tracks);

  Pasha->SetNTowerWord(ntow,ntlzp,0,0);
  
  Pasha->fSeedEt = 0;
  for (int itow=0; itow<ntow; itow++) {
    const PhysicsTower* tow = towlist[itow];
    if (tow) {
      float et    = tow->totEt();
      if ( et > Pasha->fSeedEt) {
	Pasha->fSeedEt = et;
	ieta = tow->iEta();
	iphi = tow->iPhi();
      }
    }
    else {
      fgLogger->Report(-106,"StntupleInitJet: NULL tower pointer");
    }
  }                                   
					// finally pack fTypeWord

  Pasha->SetTypeWord(TStnTypes::kJet,0,ieta,iphi);

  // New variables 9/27/2001 Matt Reece

  Pasha->fEtaeta     = Pierre->etaEtaMoment();
  Pasha->fEtaphi     = Pierre->etaPhiMoment();
  Pasha->fPhiphi     = Pierre->phiPhiMoment();
  return 0;
}

//_____________________________________________________________________________
Int_t StntupleInitJetBlock(TStnDataBlock* block, 
                           AbsEvent*      event,
                           int            mode) 
{
  // initialize jet data block with the `event' data

  int   ev_number, rn_number;
  char  process[100], description[100];

  int   inum[200];

  if (! fgFolder) {
//-----------------------------------------------------------------------------
//  initialize local static variables
//-----------------------------------------------------------------------------
    fgFolder = (TFolder*) gROOT->GetRootFolder()->FindObject("StntupleFolder");
    fgLogger = (TStnErrorLogger*) fgFolder->FindObject("StntupleErrorLogger");
  }
//-----------------------------------------------------------------------------
//  save time: don't do initialization twice for the same event
//-----------------------------------------------------------------------------
  ev_number = AbsEnv::instance()->trigNumber();
  rn_number = AbsEnv::instance()->runNumber();

  if (block->Initialized(ev_number,rn_number)) return 0;

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

  StntupleGetProcessName(data->CollName()->Data(),process,description);
  CdfJetColl_ch handle(StntupleGetIterator(event,
					   "CdfJetColl",
					   process,
					   description));
  if (handle.is_null()) {
     return -1;
  } 
  else {
    const JetList& jlist    = handle->contents(); 
    const JetAlgParams& par = handle->jetParams();

    if (&par == 0) {
      fgLogger->Report(-105,"StntupleInitJetBlock: handle->jetParams() returns 0");
      data->SetConeSize(100.);
    }
    else {
      gConeSize = par.coneRadius();
      data->SetConeSize(gConeSize);
    }

    data->fZVertex = handle->zVertex();
					// order jets in Et
    const CdfJet *jet1, *jet2;

    int njets = jlist.size();

    for (int i=0; i<njets; i++) inum[i] = i;

    for (int i1=0; i1<njets-1; i1++) {
      jet1 = &jlist[i1];
      for (int i2=i1+1; i2<njets; i2++) {
	jet2 = &jlist[i2];
	if (jet1->et() < jet2->et()) {
					// swap the jets
	  int itmp = inum[i1];
	  inum[i1] = inum[i2];
	  inum[i2] = itmp;
          jet1     = jet2;
	}
      }
    }

    TStnJet*      pasha;
    const CdfJet* pierre;

    int nj, nj10, nj15, nj20, nj25;

    nj   = 0;
    nj10 = 0;
    nj15 = 0;
    nj20 = 0;
    nj25 = 0;

    for (int i=0; i<njets; i++) {
      pierre = &jlist[inum[i]];
      pasha  = data->NewJet();
      StntupleInitJet(pasha,pierre);
      nj++;
      if (pasha->Et() > 10.) nj10++;
      if (pasha->Et() > 15.) nj15++;
      if (pasha->Et() > 20.) nj20++;
      if (pasha->Et() > 25.) nj25++;
    }

    data->SetNJets  (nj  );
    data->SetNJets10(nj10);
    data->SetNJets15(nj15);
    data->SetNJets20(nj20);
    data->SetNJets25(nj25);
  }
//-----------------------------------------------------------------------------
//  on return mark block as initialized for given event/run
//-----------------------------------------------------------------------------
  data->f_RunNumber   = rn_number;
  data->f_EventNumber = ev_number;

  return 0;
}


//_____________________________________________________________________________
Int_t StntupleJetBlockLinks(TStnDataBlock* Block, AbsEvent* Event, int Mode) 
{
  // set links between the jets and the objects in the other blocks
  // return 0 if everything is OK
  // Initialize links and other variables that need information about other blocks.
  // Added November 2001 Matt Reece

  int                ev_number, rn_number;

  TStnEvent*       ev;
  TStnTrackBlock*  tblock;
  TStnVertexBlock* vblock;
  TLorentzVector   tmom;
					// assume that the block exists, it 
					// is supposed to be initialized

  ev_number = AbsEnv::instance()->trigNumber();
  rn_number = AbsEnv::instance()->runNumber();

  if (! Block->Initialized(ev_number,rn_number)) return -1;

					// do not do initialize links 2nd time

  if (Block->LinksInitialized()) return 0;

  TStnJetBlock*    jet_block = (TStnJetBlock*) Block;

  int njets = jet_block->NJets();
  jet_block->TrackLinkList()->Init(njets);

  ev     = Block->GetEvent();
  vblock = (TStnVertexBlock*) ev->GetDataBlock("VertexBlock");
  tblock = (TStnTrackBlock* ) ev->GetDataBlock("TrackBlock");

  int      ntrk;
  double   trkEta, trkPhi, r2, jetEta, jetPhi, sumpt, cone_size, pt, max_pt;
//-----------------------------------------------------------------------------
//  kludge
//-----------------------------------------------------------------------------
  cone_size = jet_block->ConeSize(); if (cone_size == 0) cone_size = 0.4;

  r2 = cone_size*cone_size;

  for(int i=0; i<njets; ++i) {
    TStnJet* jet = jet_block->Jet(i);
    jet_block->TrackLinkList()->InitLinks(i);

					// Get jet centroid eta, phi
    jetEta = jet->Momentum()->Eta();
    jetPhi = jet->Momentum()->Phi();
    ntrk   = 0;
    sumpt  = 0;
    max_pt = 0;
    tmom.SetXYZT(0,0,0,0);
    
    for(int j=0; j<tblock->NTracks(); ++j) {
      TStnTrack* trk = tblock->Track(j);

					// Extrapolate track eta, phi at detector

      TLorentzVector* p = trk->Momentum();
      pt = p->Pt();
      Hep3Vector momentum(p->X(),p->Y(),p->Z());
      //      int charge = trk->Charge();

      HepPoint3D position(trk->X0(),trk->Y0(),trk->Z0()); 

      trkEta = trk->Momentum()->Eta();
      trkPhi = trk->Momentum()->Phi();
					// doesnt exist in 4.2.0

//       JetVariables::extrapolatedEndEtaPhi(position,
// 					  momentum,
// 					  trk->Charge(),
// 					  trkEta,
// 					  trkPhi);

      double phiDiff = TVector2::Phi_mpi_pi(trkPhi - jetPhi);
      double etaDiff = trkEta - jetEta;

      if ((etaDiff*etaDiff + phiDiff*phiDiff) < r2) {

					// count this track

        jet_block->TrackLinkList()->Add(i,j);

	ntrk++;
	sumpt += p->Pt();
	tmom  += *p;
	if (pt > max_pt) max_pt = pt;
      }
    }
    jet->SetNTrackWord(ntrk,-1,-1,-1);
					// make this p/E - name confusing, but 
					// the result - OK
    jet->SetEp(sumpt/jet->Et());
    jet->fTrackMass  = tmom.M();
    jet->fMaxTrackPt = max_pt;
  }
//-----------------------------------------------------------------------------
// Added October 2001 Matt Reece:
// Fill jet/track link information & fill jet information that requires vertices
//-----------------------------------------------------------------------------
  int ntowers;
  for(int i=0; i<njets; ++i) {
    TStnJet* jet          = jet_block->Jet(i);
    const CdfJet*  pierre = jet->GetCdfJet();
//-----------------------------------------------------------------------------
// fill tower link block
//-----------------------------------------------------------------------------
    const vector<const PhysicsTower*>& towlist = pierre->towers().contents();
    ntowers = towlist.size();
    for (int itower=0; itower<ntowers; itower++) {
      const PhysicsTower* tower =  towlist[itower];

      if (tower) {

	int tower_key = (((tower->iPhi() & 0x3f  )      ) |
			 ((tower->iEta() & 0x3f  ) <<  8) |
			 ((tower->key().isValid()) << 16)    );

	jet_block->TowerLinkList()->Add(i,tower_key);
      }
      else {
	fgLogger->Report(-106,"StntupleInitJetBlockLinks: NULL tower pointer");
      }
    }
//-----------------------------------------------------------------------------
// Now, loop over the vertices and find the one closest to jet->fZvtx
//-----------------------------------------------------------------------------
    int best_v     = -1;
    double jetZv   = jet->Vz();
    if(jetZv>-999.0) {  // if the jet did indicate a z vertex
      double minDist = 99999.0;
      for(int k=0; k < vblock->NVertices(); ++k) {
	TStnVertex* v = vblock->Vertex(k);
	double dist = fabs(v->Z()-jetZv);
	if(dist < minDist) {
	  best_v  = k;
	  minDist = dist;
	}
      }
    }
    jet->SetVertexNumber(best_v);
					// Compute the sum Pt of tracks not 
					// pointing to a vertex, and of tracks
					// pointing at fNearvtx.
    double novpt, vpt;
    vpt   = 0.0;
    novpt = 0.0;
    int numTracks = jet_block->TrackLinkList()->NLinks(i);

					// Loop over tracks associated with 
					// the given jet to see which vertex 
					// they are coming from
    for(int j=0; j<numTracks; ++j) {
      int it = jet_block->TrackLinkList()->Index(i,j);
      TStnTrack* track = tblock->Track(it);
      
					// Get the index of the closest vertex 
					// to the track (filled above)

      if (track->VertexNumber() == jet->VertexNumber()) {
	vpt += track->Pt();
      }
      else {
	novpt += track->Pt();
      }
    }
    jet->SetVpt  (vpt);
    jet->SetNovpt(novpt);
  }  
//-----------------------------------------------------------------------------
// mark links as initialized
//-----------------------------------------------------------------------------
  Block->SetLinksInitialized();

  return 0;
}