#include "SimpleConeJetFinder.h"
#include <iostream>
#include "d0root_analysis/objects/D0Track.h"
#include "d0root_analysis/objects/D0TrackJet.h"
#include<vector>
using namespace std;
using std::vector;

ClassImp(SimpleConeJetFinder)

///
 SimpleConeJetFinder::SimpleConeJetFinder(float dr, float dz, float ptMin, 
					 float ptJet, float ptSeed, int hits, 
					 float dca, float zdca) :
  TrackJetFinder(dr,dz,ptMin,ptJet,ptSeed,hits,dca,zdca) 
{
  _masterList.SetOwner(kTRUE);
}

///
 SimpleConeJetFinder::SimpleConeJetFinder(float* cuts) : 
    TrackJetFinder(cuts[0],cuts[1],cuts[2],cuts[3],cuts[4],
		   static_cast<int> (cuts[5]), cuts[6], cuts[7]) 
{
    _masterList.SetOwner(kTRUE);
}


///
 SimpleConeJetFinder::~SimpleConeJetFinder() {}



///
 void SimpleConeJetFinder::FindJets(TObjArray& jets)
{  
  if(_debug>=1) {
    cout << endl << "+++ Simple Cone D0Jet Algorithm +++" << endl;
    cout         << "---------------------------------" << endl;
    cout << endl;
    cout << "- Cone size (R,Z)      = " << _parameter << ", " << _dZ << endl;
    cout << "- minimum multiplicity = " << _mult << endl;
    cout << "- minimum seed pT      = " << _ptSeed << endl;
    cout << "- minimum jet pT       = " << _ptJet << endl;
    cout << "- minimum track pT     = " << _ptMin << endl;
    cout << "- SMT hits            >= " << _hits << endl;
    if(_pv) {
      cout << "Using Primary D0Vertex: " << endl;
      _pv->print();
    }
    cout << endl;
    cout << "Number of input tracks = " << _tracks.GetEntriesFast() << endl;
    cout << endl;
  }

  vector<D0Track*>::iterator i;
  vector<D0Track*>::iterator j;

  bool found = false;
  TVectorD dca(2); TVectorD dcaerr(2);
  float ip, ipz;

  //fill vector of selected tracks
  vector<D0Track*> lvtracks;
  for(int h=0; h<_tracks.GetEntriesFast(); h++) {
    D0Track* trk = (D0Track*) _tracks[h];

    ip = 0.0; ipz = 0.0;
    if(_pv) {
      trk->Get_impactParameters(_pv,dca,dcaerr);
      ip = TMath::Abs(dca(0));
      ipz = TMath::Abs(dca(1));
    }


    if(_debug) {
      cout << trk->Pt() << " "<< trk->GetHitMask().num_smt_ladders() +
       trk->GetHitMask().num_smt_f_wedges() << " ";
      cout << ip << " " << ipz << " " << trk->Z() << endl;
    }

    if(trk->Pt()>=_ptMin && trk->GetHitMask().num_smt_ladders() +
       trk->GetHitMask().num_smt_f_wedges() >=_hits && 
       ip<=_dca && ipz<=_zdca && trk->Z()!=0) lvtracks.push_back(trk);
  }
  

  if(_debug>=2) {
    cout << "Number of selected tracks = "<< lvtracks.size() << endl;
  }


  jets.Clear();


  int id = 0;
  
  // outer loop
  i=lvtracks.begin();
  while(!found && i !=lvtracks.end()){
    D0Track* trk = (*i);            // get pointer to particle
    if(trk->Pt()<_ptSeed) {
      found = true;

      if(_debug>=2) {
	cout << "track #id = " << trk->Id() << ", pT = " << trk->Pt();
	cout << endl << "No more tracks to cluster" << endl;
      }
    } else {
      TObjArray particles; particles.SetOwner(kFALSE);
      
      //replaze trk by trk_i
      particles.Add(trk);
      float etaSeed = trk->Eta();
      float phiSeed = trk->Phi();
      float zSeed = trk->Z();
      TLorentzVector sum(trk->Px(),trk->Py(),trk->Pz(),trk->E());

      if(_debug>=2) {
	cout << "+ Starting new jet" << endl;
	cout << "track #id = " << trk->Id() << ", pT = " << trk->Pt(); 
	cout << endl << "Eta = " << etaSeed << ", Phi = " << phiSeed;
	cout << ", Z = " << zSeed << endl << endl;
      }


      //inner loop
      for(j=i; j != lvtracks.end(); ++j) {
	if(j!=i){  //  do not use seed again
	  trk = (*j);                    // get pointer to Vtrack

	  float dr = trk->DeltaR(sum);
	  float dz = trk->Z()-zSeed;
	  
	  if(_debug>=2) {
	    cout << "  track id #" << trk->Id() << ", pT = "<<trk->Pt(); 
	    cout << endl << "  Delta R = " << dr << ", Delta Z = " << dz;
	    cout << endl << endl;
	  }
	  
	  if(dr<_parameter && TMath::Abs(dz)<_dZ) {
	    
	    sum += (*trk);
	    particles.Add(trk);
	    lvtracks.erase(j--);  

	    //re-calculate center
	    etaSeed = sum.Eta();
	    phiSeed = sum.Phi();
	    zSeed += (1./particles.GetEntriesFast())*dz;
	    
	    if(_debug>=2) {
	      cout << "   ----- track within cone " << endl;
	      cout << "         New center = " << etaSeed << ", " << phiSeed;
	      cout << "; " << zSeed << endl;
	      cout << endl << endl;
	    }	      
	  }
	  if(j == lvtracks.end()) break; // quit if list ends prematurely
	}
      } //end inner loop
      
      id ++;
      D0TrackJet* jet = new D0TrackJet(particles,zSeed,id);
      if(_pv) jet->PrimaryVertex(_pv);
      _masterList.Add(jet);


      if(jet->Pt() > _ptJet && jet->Multiplicity()>=_mult) {
	jets.Add(jet);
	
	if(_debug>=2) {
	  cout << endl;
	  cout << "D0Jet found " << endl;
	  jet->print();
	}
      }
      if(i == lvtracks.end()) break;  // quit if list ends prematurely
      ++i;
    }
    
    
  } // end outer loop
  

  jets.Sort();

  if(_debug >= 1) {
    cout << endl << "************************************" << endl;
    cout << jets.GetEntriesFast() << " track-jets found" << endl << endl;
    cout << "List of jets found: " << endl << endl;
    for(int i=0; i<jets.GetEntriesFast(); i++) {
      D0TrackJet* jet = (D0TrackJet*) jets[i];
      jet->print();
    }
    cout << endl << "************************************" << endl;
  }
}