//
// Pibero Djawotho <pibero@indiana.edu>
// Indiana University
// 8 October 2006
//

// C++
#include <fstream>
#include <set>
using namespace std;

// ROOT
#include "TFile.h"
#include "TNtuple.h"
#include "TH1.h"
#include "TH2.h"
#include "TChain.h"

// STAR
#include "StEventTypes.h"
#include "StMuDSTMaker/COMMON/StMuTypes.hh"
//#include "StDaqLib/TRG/trgStructures2007.h"
#include "StDaqLib/TRG/trgStructures2008.h"
//#include "StGenericVertexMaker/StCtbMatcher.h"

// Local
#include "Event.h"
#include "Track.h"
#include "EmcHit.h"
#include "EmcCluster.h"
#include "CtbHit.h"
#include "Particle.h"
#include "SmdHit.h"
#include "SmdCluster.h"
#include "L2Result.h"
#include "StEventMixer.h"
#include "StUpsilonTreeReader.h"

ClassImp(StUpsilonTreeReader);

StUpsilonTreeReader::StUpsilonTreeReader(int nfiles ,const char* filelist)
{
  // Open output ROOT file
  //mFile = new TFile("ups.histo.root", "recreate");
  //assert(mFile && !mFile->IsZombie());

const char* name1 = "ups.histo.root";
 if (gSystem->Getenv("JOBID"))
 name1 = gSystem->ConcatFileName(gSystem->Getenv("SCRATCH"),name1);
 mFile = new TFile(name1, "recreate");
 assert(mFile && !mFile->IsZombie());

  // Create ntuples
  TString varlist = "run:event:refmult:etotal:level2:charge:vertexZ:pt:y:phi:m:cosTheta:L2_E1:L2_E2:L2_m:L2_cosTheta";
  varlist += ":charge1:p1:pt1:eta1:phi1:dEdx1:nHitsFit1:etrig1:ecluster1:ntracks1:ctb1:nctbtracks1:bsmde1:bsmdp1:tracks_with_hit1";
  varlist += ":charge2:p2:pt2:eta2:phi2:dEdx2:nHitsFit2:etrig2:ecluster2:ntracks2:ctb2:nctbtracks2:bsmde2:bsmdp2:tracks_with_hit2";
  mUps = new TNtuple("ups", "ups", varlist);
  //mIso = new TNtuple("iso", "iso", varlist);
  //mIso2 = new TNtuple("iso2", "iso2", varlist);
  //mMix = new TNtuple("mix", "mix", "refmult:vz1:vz2:level2:charge:pt:y:phi:m:cosTheta:charge1:charge2:p1:p2:pt1:pt2:eta1:eta2:phi1:phi2:dEdx1:dEdx2:nhits1:nhits2:etrig1:etrig2:ecluster1:ecluster2");

  // Create histograms
  mTriggerIds = new TH1F("triggerIds", "Trigger IDs", 5, 0, 5);
  mRefMult = new TH1F("refMult", "Reference multiplicity", 100, 0, 1000);
  mZvertex = new TH1F("zVertex", "z-vertex", 200, -200, 200);
  RefMult_vs_vz = new TH2F("RefMult_vs_vz", "RefMult_vs_vz", 200,-200,200,100,0,100);
  good_RefMult = new TH1F("good_RefMult", "Reference multiplicity with vz cut", 100, 0, 100);



  // Create chain and add files
  mChain = new TChain("upsTree");
  ifstream in(filelist);
  assert(in);
  string line;
  for (int i = 0; i < nfiles && getline(in, line); ++i) mChain->Add(line.c_str());
  in.close();

  // Event buffer
  mEvent = new Event;
  mChain->SetBranchAddress("eventBranch", &mEvent);
}

StUpsilonTreeReader::~StUpsilonTreeReader()
{
  mFile->Write();
  mFile->Close();
}

int StUpsilonTreeReader::EventLoop()
{
  // Event mixer
  //StEventMixer eventMixer;

  // Event loop
  int nentries = mChain->GetEntries();
  for (int ev = 0; ev < nentries; ++ev) {
    // Load entire event
    mChain->GetEntry(ev);

    set<int>& triggerIds = mEvent->triggerIds();

    for (set<int>::iterator i = triggerIds.begin(); i != triggerIds.end(); ++i)
      mTriggerIds->Fill(Form("%d", *i), 1);

    mRefMult->Fill(mEvent->refMult());
    mZvertex->Fill(mEvent->vertex().z());
	RefMult_vs_vz->Fill(mEvent->vertex().z(),mEvent->refMult());
	if(mEvent->vertex().z()>-30&&mEvent->vertex().z()<30)
		good_RefMult->Fill(mEvent->refMult());

	/*
#if 0
    // Event mixer
    eventMixer.addEvent(mEvent);

    int refmult = mEvent->refMult();

    if (eventMixer.ready(refmult)) {
      for (int k = 0; k < 10; ++k) {
	Event* event1 = eventMixer.getEvent(refmult);
	Event* event2 = eventMixer.getEvent(refmult);
	if (event1 && event2 && *event1 != *event2) {
	  // Get electrons from event1
	  vector<Track*> electrons;
	  for (int i = 0; i < event1->tracks().GetEntriesFast(); ++i) {
	    Track* track = (Track*)event1->tracks().At(i);
	    if (accept(track) && track->charge() < 0) electrons.push_back(track);
	  }

	  // Get positrons from event2
	  vector<Track*> positrons;
	  for (int i = 0; i < event2->tracks().GetEntriesFast(); ++i) {
	    Track* track = (Track*)event2->tracks().At(i);
	    if (accept(track) && track->charge() > 0) positrons.push_back(track);
	  }

	  // Mix electrons and positrons
	  for (size_t i = 0; i < electrons.size(); ++i) {
	    Track* track1 = electrons[i];
	    for (size_t j = 0; j < positrons.size(); ++j) {
	      Track* track2 = positrons[j];

	      Particle particle(track1, track2);
	      float tuple[] = {
		refmult,
		event1->vertex().z(),
		event2->vertex().z(),
                particle.level2(),
		particle.charge(),
		particle.pt(),
		particle.y(),
		particle.phi(),
		particle.m(),
		particle.cosTheta(),
		track1->charge(),
		track2->charge(),
		track1->momentum().Mag(),
		track2->momentum().Mag(),
		track1->momentum().Pt(),
		track2->momentum().Pt(),
		track1->momentum().Eta(),
		track2->momentum().Eta(),
		track1->momentum().Phi(),
		track2->momentum().Phi(),
		track1->dEdx(),
		track2->dEdx(),
		track1->numberOfHits(),
		track2->numberOfHits(),
		track1->emcCluster()->seed()->energy(),
		track2->emcCluster()->seed()->energy(),
		track1->emcCluster()->energy(),
		track2->emcCluster()->energy()
	      };

	      mMix->Fill(tuple);
	    }
	  }      
	}
      }
    }
#endif
*/


    // L2 result
//    L2Result* upsL2Result = (L2Result*)&mEvent->L2Result()[L2RESULTS_2007_OFFSET_UPS];
    L2Result* upsL2Result = (L2Result*)&mEvent->L2Result()[L2RESULTS_2008_OFFSET_UPS];
    float L2_E1 = upsL2Result->energyOfClusterL0;
    float L2_E2 = upsL2Result->energyOfClusterL2;
    float L2_m = upsL2Result->invMass;
    float L2_cosTheta = 1 - (L2_m * L2_m) / (2 * L2_E1 * L2_E2);

    //cout<<" L2_mass:"<<L2_m<<endl;

    // Loop over virtual particles
    for (int i = 0; i < mEvent->particles().GetEntriesFast(); ++i) {
      Particle* particle = (Particle*)mEvent->particles().At(i);
      Track* track1 = (Track*)particle->daughters().First();
      Track* track2 = (Track*)particle->daughters().Last();

      if (track1->pt() < track2->pt()) swap(track1, track2);
//this is to accept tracks with matched EMC hits//(needs check)

//     cout<<"Particles:"<<mEvent->particles().GetEntriesFast()<<endl;

     if (!accept(track1)) continue;
     if (!accept(track2)) continue;

      int ntracks1 = 0;
      int tracks_with_hit1=1;
      TRefArray& hits1 = track1->emcCluster()->hits();
      for (int j = 0; j < hits1.GetEntriesFast(); ++j) {
	EmcHit* hit = (EmcHit*)hits1.At(j);
	ntracks1 += hit->tracks().GetEntriesFast();
	for (int k = 0; k < hit->tracks().GetEntriesFast(); ++k) { 
		Track* temp_track = (Track*)hit->tracks().At(k);
		float delta_pt=track1->pt()-temp_track->pt();
		float delta_phi=track1->phi()-temp_track->phi();
		float delta_eta=track1->eta()-temp_track->eta();
		if(!(delta_pt>-1.e-5&&delta_pt<1.e-5//check same track
			   &&delta_phi>-1.e-5&&delta_phi<1.e-5
			   &&delta_eta>-1.e-5&&delta_eta<1.e-5
			))tracks_with_hit1++;
	}
 //       cout<<"Ntracks1:"<<ntracks1<<endl;
      }

      int ntracks2 = 0;
      int tracks_with_hit2=1;
      TRefArray& hits2 = track2->emcCluster()->hits();
      for (int j = 0; j < hits2.GetEntriesFast(); ++j) {
	EmcHit* hit = (EmcHit*)hits2.At(j);
	ntracks2 += hit->tracks().GetEntriesFast();

	for (int k = 0; k < hit->tracks().GetEntriesFast(); ++k) { 
		Track* temp_track = (Track*)hit->tracks().At(k);
		float delta_pt=track2->pt()-temp_track->pt();
		float delta_phi=track2->phi()-temp_track->phi();
		float delta_eta=track2->eta()-temp_track->eta();
		if(!(delta_pt>-1.e-5&&delta_pt<1.e-5//check same track
			   &&delta_phi>-1.e-5&&delta_phi<1.e-5
			   &&delta_eta>-1.e-5&&delta_eta<1.e-5
			))tracks_with_hit2++;
	}

	}

      int bsmde1 = track1->smdEtaCluster() ? track1->smdEtaCluster()->hits().GetEntriesFast() : 0;
      int bsmde2 = track2->smdEtaCluster() ? track2->smdEtaCluster()->hits().GetEntriesFast() : 0;
      int bsmdp1 = track1->smdPhiCluster() ? track1->smdPhiCluster()->hits().GetEntriesFast() : 0;
      int bsmdp2 = track2->smdPhiCluster() ? track2->smdPhiCluster()->hits().GetEntriesFast() : 0;



      // Fill ntuple
      float tuple[] = {
	mEvent->runId(),
	mEvent->eventId(),
	mEvent->refMult(),
	mEvent->energy(),
	particle->level2(),
	particle->charge(),
	mEvent->vertex().z(),
	particle->pt(),
	particle->y(),
	particle->phi(),
	particle->m(),
	particle->cosTheta(),
	L2_E1,
	L2_E2,
	L2_m,
	L2_cosTheta,

	track1->charge(),
	track1->momentum().Mag(),
	track1->momentum().Pt(),
	track1->momentum().Eta(),
	track1->momentum().Phi(),
	track1->dEdx(),
	track1->numberOfHits(),
	track1->emcCluster()->seed()->energy(),
	track1->emcCluster()->energy(),
	ntracks1,
	track1->ctbHit() ? track1->ctbHit()->adc() : 0,
	track1->ctbHit() ? track1->ctbHit()->tracks().GetEntriesFast() : 0,
	bsmde1,
	bsmdp1,
	tracks_with_hit1,

	track2->charge(),
	track2->momentum().Mag(),
	track2->momentum().Pt(),
	track2->momentum().Eta(),
	track2->momentum().Phi(),
	track2->dEdx(),
	track2->numberOfHits(),
	track2->emcCluster()->seed()->energy(),
	track2->emcCluster()->energy(),
	ntracks2,
	track2->ctbHit() ? track2->ctbHit()->adc() : 0,
	track2->ctbHit() ? track2->ctbHit()->tracks().GetEntriesFast() : 0,
	bsmde2,
	bsmdp2,
	tracks_with_hit2
      };

      mUps->Fill(tuple);

      //if (isolated(track1)) {
	//mIso->Fill(tuple);
	//if (isolated(track2)) 
         //mIso2->Fill(tuple);
      //}
    }

    mEvent->Clear();
  }

  return nentries;
}

inline bool StUpsilonTreeReader::accept(Track* track) const
{
  return track->emcCluster();
}
/*
#include "StEmcUtil/projection/StEmcPosition.h"

bool StUpsilonTreeReader::isolated(Track* track) const
{
  if (track->emcHits().IsEmpty()) return false;
  EmcHit* hit = (EmcHit*)track->emcHits().First();
  int id = hit->id();
  StEmcPosition emcPosition;
  vector<int> towers3x3;
  vector<int> towers5x5;
  float energy3x3 = 0;
  float energy5x5 = 0;
  for (int deta = -2; deta <= 2; ++deta) {
    for (int dphi = -2; dphi <= 2; ++dphi) {
      int id2 = emcPosition.getNextTowerId(id, deta, dphi);
      if (!id2) continue;
      EmcHit* hit2 = getTower(id2);
      if (!hit2) continue;
      towers5x5.push_back(id2);
      energy5x5 += hit2->energy();
      if (abs(deta) < 2 && abs(dphi) < 2) {
	towers3x3.push_back(id2);
	energy3x3 += hit2->energy();
      }
    }
  }
  return energy3x3 > 0.85 * energy5x5;
}

EmcHit* StUpsilonTreeReader::getTower(int id) const
{
  TClonesArray& emcHits = mEvent->emcHits();
  for (int i = 0; i < emcHits.GetEntriesFast(); ++i) {
    EmcHit* hit = (EmcHit*)emcHits.At(i);
    if (hit->id() == id) return hit;
  }
  return 0;
}
*/
