#include "StHiJet.h"
#include "TSystem.h"
#include "TClonesArray.h"
#include "TChain.h"
#include <TStopwatch.h>
#include <TFile.h>
#include <TTree.h>
#include <sstream>
#include "TF2.h"
#include "TProfile.h"
#include "TVector2.h"
#include "TString.h"
#include "TROOT.h"
#include "TObjString.h"
#include "../StHiMicroEvent/StHiMicroEvent.h"
#include "StarClassLibrary/StThreeVectorF.hh"


//define binning
const float StHiJet::mMinSeed=5.0;
const float StHiJet::mJetR=0.3;
const float StHiJet::mEtaCut=1;

const bool includeFtpc=false;

// Help Functions

int getHighPtCentBin(int centrality)
{

  Int_t centhighpt=-1;

  // analog zu highpt ?
  if (centrality==9) centhighpt=6;
  else if (centrality==8) centhighpt=5;
  else if (centrality==7) centhighpt=4;
  else if (centrality==6) centhighpt=3;
  else if (centrality==5) centhighpt=2;
  else if (centrality==4 || centrality==3) centhighpt=1;
  else if (centrality==2 || centrality==1) centhighpt=0;
  else centhighpt=-1;

  return centhighpt;
}

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


StHiJet::StHiJet()
{
  mMinAssoc = 4.0;
  //out_root_dir="/eliza2/rnc/bhaag/dAudata/";
  out_root_dir="~/AuTest/AuAu/centFix/sortTriggers/assoc40/";

}

StHiJet::~StHiJet()
{
}

Float_t StHiJet::dPhi(Float_t phi1, Float_t phi2) {
  static float twopi=2*TMath::Pi();

  Float_t d_phi=phi2-phi1;

  if (d_phi < -0.5*TMath::Pi()) {d_phi = twopi + d_phi;} 
  if (d_phi > 1.5*TMath::Pi()) {d_phi = d_phi - twopi;}


  // DEBUG:
  //if (fabs(d_phi)>TMath::Pi())
  //cout<<d_phi<<endl;

  return d_phi;
}

Float_t StHiJet::GetR(StHiMicroTrack *tr1, StHiMicroTrack *tr2) {
  Float_t dphi=tr2->PhiPr()-tr1->PhiPr();
  if (dphi < -TMath::Pi())
    dphi += 2*TMath::Pi();
  if (dphi > TMath::Pi())
    dphi -= 2*TMath::Pi();
  Float_t deta=tr2->EtaPr()-tr1->EtaPr();
  return sqrt(dphi*dphi+deta*deta);
}

void StHiJet::initHistos() {
  
  HistB = new TH1D("HistB","dNdphi;#Delta#phi;#frac{dN}{d(#delta#phi)}",24,-0.5*TMath::Pi(),1.5*TMath::Pi());	
  HistBbig = new TH2D("dPhiJetPtAssoc","dPhiJetPtAssoc;p_{T,assoc} (GeV);#Delta#phi",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi());	
  HistBbigJetPt = new TH3D("dPhiJetPtAssocJetPt","dPhiJetPtAssoc;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),24,6.0,30.0);	
  mdPhiSeedAssocSingleJetPt = new TH3D("dPhiSeedAssocSingleJetPt","dPhiJetPtAssoc;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),26,4.0,30.0);//Single Particle JetPt, associated Pt, delta phi	
  mdPhiSeedAssocSingleClusterJetPt = new TH3D("dPhiSeedAssocSingleClusterJetPt","dPhiJetPtAssocSingleCluster;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),26,4.0,30.0);//Clusters with only 1 particle  JetPt, associated Pt, delta phi	
  mdPhiSeedAssocMultipleClusterJetPt = new TH3D("dPhiSeedAssocMultipleClusterJetPt","dPhiJetPtAssocMultipleCluster;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),26,4.0,30.0);//Clusters with more than 1 particle JetPt, associated Pt, delta phi	
  mPtSeeds = new TH1D("PtSeeds","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//Single Particle Triggers	
  mPtSingleSeeds = new TH1D("PtSingleSeeds","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//Single Particle Cluster Triggers	
  mPtMultiSeeds = new TH1D("PtMultiSeeds","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//Multiple Particle Cluster Triggers	
  Triggers = new TH1D("dNdPtJet","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//All Cluster Triggers	
  
  mPt = new TH1D("Pt","Pt;p_{T} (GeV)",25,0,25);
  mPtJet = new TH1D("PtJet","PtJet;p_{T,jet} (GeV)",25,0,25);
  mNTrackPtJet = new TH2D("NTrackPtJet","NTrackPtJet;p_{T,jet} (GeV);n_{trk}",25,0,25,10,0.5,10.5);
  mNTrackPtSeed = new TH2D("NTrackPtSeed","NTrackPtSeed;p_{T,jet} (GeV);n_{trk}",25,0,25,10,0.5,10.5);
  mNSeed = new TH1D("NSeed","NSeed;n_{seed}",10,-0.5,9.5);
  mdPhiPtMax  = new TH2D("dPhiPtMax","dPhiPtMax;#Delta#phi;p_{T,max}",24,-0.5*TMath::Pi(),1.5*TMath::Pi(),25,0,25);
  mdPhiSeedAll  = new TH1D("dPhiSeedAll","dPhiSeedAll;#Delta#phi",24,-0.5*TMath::Pi(),1.5*TMath::Pi());
  mdPhiSeedAssoc  = new TH1D("dPhiSeedAssoc","dPhiSeedAssoc;#Delta#phi",24,-0.5*TMath::Pi(),1.5*TMath::Pi());
  mdPhidEta = new TH2D("dPhidEta","dPhidEta;#Delta#phi;#DeltaEta",24,-0.5*TMath::Pi(),1.5*TMath::Pi(),20,-2,2);
  mPtPt = new TH2D("PtPt","PtPt;p_{T} (GeV);p_{T} (GeV)",25,0,25,25,0,25);
  mdPhiPtPt = new TH3D("dPhiPtPt","dPhiPtPt;#Delta#phi;p_{T} (GeV);p_{T} (GeV)",24,-0.5*TMath::Pi(),1.5*TMath::Pi(),25,0,25,25,0,25);
  mNJetPtMax = new TH2D("NJetPtMax","NJetPtMax;p_{T,max} (GeV);n_{jet}",25,0,25,5,0.5,5.5);
}

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

void StHiJet::histofill(Int_t maxevents, Int_t cent)
{
  TStopwatch timer;
  timer.Start();
  Float_t mbytes =0;

  StHiMicroEvent* event = new StHiMicroEvent();

  TChain* pPicoChain = new TChain("StHiMicroTree");
  pPicoChain->SetBranchAddress("StHiMicroEvent", &event);

  for(int ii=0; ii<mInputFileList->GetLast()+1;ii++){
    //cout << "Adding file " << ((TObjString*) mInputFileList->At(ii))->GetString() << endl;
    pPicoChain->Add(((TObjString*) mInputFileList->At(ii))->GetString());
  }

  cout<<" ..."<<endl;
  //cout << "total events: " << pPicoChain->GetEntries() << endl;
  TString mLastFile=((TObjString*) mInputFileList->At(mInputFileList->GetLast()))->GetString();
  // DEBUG:
  //mLastFile="/eliza2/rnc/putschke/himicro_mid/st_physics_5033127_raw_4010003.himicro.root";
  cout << "Last file in list : "<<mLastFile<<endl;

  cout << "requested number of events to process: " << maxevents << endl;

  // open output file
  TFile *hfile = (TFile*) gROOT->FindObject("histos.root");
  if (hfile) hfile->Close();
  char fileName[80];
  if (cent<99)
    sprintf(fileName, "%sjethist_pt%2d_Cent%1d_.root",out_root_dir.Data(),Int_t(10*mMinAssoc),cent);
  else if (cent==99)
    sprintf(fileName, "%sjethist_pt%2d_MinBias.root",out_root_dir.Data(),Int_t(10*mMinAssoc));

  //sprintf(fileName, "root.flow/ana.central.%2d.root", cent+10);
  hfile = new TFile(fileName,"RECREATE","file with histograms");
  if (!hfile->IsOpen())
	  cout << "Problem opening file " << fileName << endl;
  initHistos();
  
  Int_t run_no, evt_no, i_seed, n_trk;
  Float_t pt_jet, pt_seed, eta_seed, phi_seed;
  TTree *jet_tree = new TTree("jet_tree","jet_tree"); 
  jet_tree->Branch("run",&run_no,"run/I");  
  jet_tree->Branch("event",&evt_no,"event/I");
  jet_tree->Branch("i_seed",&i_seed,"i_seed/I");  
  jet_tree->Branch("n_trk",&n_trk,"n_trk/I");
  jet_tree->Branch("pt",&pt_jet,"pt");
  jet_tree->Branch("pt_seed",&pt_seed,"pt_seed");
  jet_tree->Branch("eta",&eta_seed,"eta");
  jet_tree->Branch("phi",&phi_seed,"phi");

  Int_t n_trk1, n_trk2;
  Float_t dphi_jets;
  Float_t pt1, pt_seed1, eta1, phi1;  
  Float_t pt2, pt_seed2, eta2, phi2;
  TTree *pair_tree = new TTree("pair_tree","pair_tree"); 
  pair_tree->Branch("run",&run_no,"run/I");  
  pair_tree->Branch("event",&evt_no,"event/I");
  pair_tree->Branch("dphi",&dphi_jets,"dphi");
  pair_tree->Branch("pt1",&pt1,"pt1");
  pair_tree->Branch("n_trk1",&n_trk1,"n_trk1/I");
  pair_tree->Branch("pt_seed1",&pt_seed1,"pt_seed1");
  pair_tree->Branch("eta1",&eta1,"eta1");
  pair_tree->Branch("phi1",&phi1,"phi1");
  pair_tree->Branch("pt2",&pt2,"pt2");
  pair_tree->Branch("n_trk2",&n_trk2,"n_trk2/I");
  pair_tree->Branch("pt_seed2",&pt_seed2,"pt_seed2");
  pair_tree->Branch("eta2",&eta2,"eta2");
  pair_tree->Branch("phi2",&phi2,"phi2");

  // fill correlation histos
  Int_t nev=0;
  Int_t nevt_seed=0;
  Int_t n_seed_tot=0;

  TString oldfile;
  TString newfile;

  for(Int_t i=0; i<maxevents;i++){
    pPicoChain->Clear();
  
    if(nev%10000 == 0) cout << "event " << nev << endl;
    
    Float_t nbyte_evt;
    if ((nbyte_evt = pPicoChain->GetEntry(i)) <= 0) // may happen once due to bad file
      continue;
    mbytes += nbyte_evt*1.e-6;

    newfile=pPicoChain->GetFile()->GetName();

    
    if(newfile != oldfile)
      cout << "analyzing file " <<  pPicoChain->GetFile()->GetName() << endl;
    
    oldfile=newfile;
    
    nev++;

    run_no = event->RunId();
    evt_no = event->EventId();

    if (event->RunId() < 5000000) { // assume dAu
       if (!event->IsTrigger(2001) && !event->IsTrigger(2003))
          continue;
       if (fabs(event->VertexZ()) > 60)
          continue;
    }
    else if (event->RunId() < 6000000) { // AuAu year-4
       if (cent == 10) { //central trigger
         if (!event->IsTrigger(15105))
            continue;
         // No vtx-z cut needed
       }
       else { // Minbias trigger
         if (!event->IsTrigger(15007))
            continue;
         if (fabs(event->VertexZ()) > 30)
            continue;
         if (cent<99)
           if(getHighPtCentBin(event->Centrality())!=cent)
             continue;
       }
    }
    else {
      cout << "Runnumber " << event->RunId() << " out of bounds" << endl;
      continue;
    }


     //centrality check - only running over minbias though and centrality bin of 7 achieves 10-20% cut or so - fraction of most central

/*    if (!event->IsTrigger(15007))
    {
	continue;
    }
    else
	{
    	if (cent<99)
        //if(getHighPtCentBin(event->Centrality())!=cent) 
        if(event->Centrality()!=cent) //flow centrality definition 
	continue;
	}
*/
  
  //THIS IS ONLY COMMENTED OUT FOR DAU
    
/*if (!event->IsTrigger(2001) || !event->IsTrigger(2003))
	{
	    if(event -> EventId() % 1000 == 0)
		{
			int trigger_counter=0;
			while(event->TriggerId(trigger_counter) !=0)
			{
			cout << "Bad trigger was:" << 
			event->TriggerId(trigger_counter)<<endl;
			trigger_counter++;
			}
		}
	                continue;
	}
*/
    // check Vertex-z position

    if (fabs(event->VertexZ())>30) continue;
    
    TObjArray seedTracks;
    TObjArray assocTracks;

    for (Int_t i_tr=0; i_tr < event->NTrack(); i_tr++) {
      
      StHiMicroTrack* track = (StHiMicroTrack*) event->tracks()->At(i_tr);
      
      if(track->FitPts()<20)continue;
      if(fabs(track->DcaGl())>1.0)continue;
      if(fabs(track->EtaPr())>mEtaCut)continue;
      mPt->Fill(track->PtPr());
      if(track->PtPr() > mMinAssoc)
	assocTracks.AddLast(track);

      if(track->PtPr() > mMinSeed)
	seedTracks.AddLast(track);
    }
    
    mNSeed->Fill(seedTracks.GetEntries());

    if(seedTracks.GetEntries()==0)continue;

    nevt_seed++;
    
    // Fill Vertex Z histogram

    //hVertZ->Fill(event->VertexZ());
    //hCentrality->Fill(getHighPtCentBin(event->Centrality())); // Flow centrality definition -> HighPt def.

    Int_t n_seed = seedTracks.GetEntries();
    Int_t n_assoc = assocTracks.GetEntries();
    pt1=0;
    StHiMicroTrack *Seed1=0;
    pt2=0;
    StHiMicroTrack *Seed2=0;
    Float_t singlept_jet; //get pt_jet for single particle comparison
    
    for (i_seed=0; i_seed < n_seed; i_seed++) {
      n_seed_tot++;
      StHiMicroTrack *seed_tr = (StHiMicroTrack *)seedTracks[i_seed];
      cout << "Seed " << i_seed << " pt " << seed_tr->PtPr() << ", phi " << seed_tr->PhiPr() << endl;
      for (Int_t i_seed2=i_seed+1; i_seed2 < n_seed; i_seed2++) {
	StHiMicroTrack *seed2_tr = (StHiMicroTrack *)seedTracks[i_seed2];
	mdPhiSeedAll->Fill(dPhi(seed_tr->PhiPr(),seed2_tr->PhiPr()));
       }
      
      singlept_jet=seed_tr->PtPr();
      pt_jet=seed_tr->PtPr();
      n_trk=1;

      for (Int_t i_assoc=0; i_assoc < n_assoc; i_assoc++) {
	StHiMicroTrack *assoc_tr = (StHiMicroTrack *)assocTracks[i_assoc];
	if (assoc_tr==seed_tr)
	  continue;
	cout << "assoc " << i_assoc << " pt " << assoc_tr->PtPr() << ",phi " << assoc_tr->PhiPr() << ", dphi " << dPhi(seed_tr->PhiPr(),assoc_tr->PhiPr()) << endl;
	mdPhiSeedAssoc->Fill(dPhi(seed_tr->PhiPr(),assoc_tr->PhiPr()));
	mdPhiSeedAssocSingleJetPt->Fill(assoc_tr->PtPr(),dPhi(seed_tr->PhiPr(),assoc_tr->PhiPr()),singlept_jet);//get comparison histogram for single trigger comparison
	if (GetR(seed_tr,assoc_tr) < mJetR) {
	  pt_jet += assoc_tr->PtPr();
	  n_trk++;
	}
      }
      mPtSeeds->Fill(singlept_jet);
      mPtJet->Fill(pt_jet);
      mNTrackPtJet->Fill(pt_jet,n_trk);
      mNTrackPtSeed->Fill(seed_tr->PtPr(),n_trk);
      pt_seed = seed_tr->PtPr();
      phi_seed = seed_tr->PhiPr();      
      eta_seed = seed_tr->EtaPr();
      jet_tree->Fill();
      if (pt_jet > pt1) {
	if (pt1 > 0 && GetR(seed_tr,Seed1) > mJetR) {
	  pt2 = pt1;
	  n_trk2 = n_trk1;
	  Seed2 = Seed1;
	}
	pt1 = pt_jet;
	n_trk1 = n_trk;
	Seed1 = seed_tr;
      }
      else if (pt_jet > pt2 && GetR(seed_tr,Seed1) > mJetR) {
	pt2 = pt_jet;
	n_trk2 = n_trk;
	Seed2 = seed_tr;
      }
    }
    if (Seed1) {
      phi1 = Seed1->PhiPr();      
      eta1 = Seed1->EtaPr();      
    }
    if (Seed2) {
      phi2 = Seed2->PhiPr();      
      eta2 = Seed2->EtaPr();      
    }
    if (pt1 > 0 && pt2 > 0) {
      mPtPt->Fill(pt1,pt2);
      dphi_jets = dPhi(phi1,phi2);

      pair_tree->Fill();

      mdPhiPtMax->Fill(dphi_jets,pt1);
      mdPhiPtPt->Fill(dphi_jets,pt1,pt2);
      mdPhidEta->Fill(dphi_jets, Seed2->EtaPr() - Seed1->EtaPr());      
    }
    mNJetPtMax->Fill(pt1, seedTracks.GetEntries());
 
    cout << "pt1 " << pt1 << ", phi1 " << phi1 << endl;
    if (pt1 > 0) { 
	Triggers->Fill(pt1);
	for (Int_t i_assoc=0; i_assoc < n_assoc; i_assoc++){
	StHiMicroTrack *assoc_tr = (StHiMicroTrack *)assocTracks[i_assoc];
    	      	//if(assoc_tr->FitPts()<20)continue;
	 	//if(fabs(assoc_tr->DcaGl())>1.0)continue;
		//if(fabs(assoc_tr->EtaPr())>mEtaCut)continue;
	if(assoc_tr == Seed1)continue;			  
	cout << "assoc " << i_assoc << " pt " << assoc_tr->PtPr() << ",phi " << assoc_tr->PhiPr() << ", dphi " << dPhi(phi1,assoc_tr->PhiPr()) << endl;
//		HistB->Fill(dPhi(seed_tr->PhiPr(),assoc_tr->PhiPr()));
	HistB->Fill(dPhi(phi1,assoc_tr->PhiPr()));
	HistBbig->Fill(assoc_tr->PtPr(),dPhi(phi1,assoc_tr->PhiPr()));
	HistBbigJetPt->Fill(assoc_tr->PtPr(),dPhi(phi1,assoc_tr->PhiPr()),pt_jet);
	if(n_trk == 1){
	mdPhiSeedAssocSingleClusterJetPt->Fill(assoc_tr->PtPr(),dPhi(phi1,assoc_tr->PhiPr()),pt_jet);
	mPtSingleSeeds->Fill(pt_jet);
	}
	if(n_trk > 1){
	mdPhiSeedAssocMultipleClusterJetPt->Fill(assoc_tr->PtPr(),dPhi(phi1,assoc_tr->PhiPr()),pt_jet);	
	mPtMultiSeeds->Fill(pt_jet);
	}
	}	
     }
  }

  hfile->Write();
  cout << "Written to file" << hfile->GetName() << endl;
  hfile->ls();
  hfile->Close();

  cout << "seed threshold " << mMinSeed << endl;
  cout << "assoc threshold " << mMinAssoc << endl;
  cout << "radius cut " << mJetR << endl;
  cout << "tot evts: " << nev << endl;
  cout << "with seeds " << nevt_seed << ", tot seeds " << n_seed_tot << endl;

  //  Stop timer and print results
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  printf("You have scanned %f Mbytes/Realtime seconds\n",mbytes/rtime);
  printf("You have scanned %f Mbytes/Cputime seconds\n",mbytes/ctime);
}