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

// C++ STL
#include <bitset>
#include <algorithm>
// ROOT
#include "TFile.h"
#include "TTree.h"

// STAR
#include "StEventTypes.h"
#include "StMuDSTMaker/COMMON/StMuTypes.hh"
#include "StEmcUtil/geometry/StEmcGeom.h"
#include "StEmcUtil/database/StBemcTables.h"
#include "StEmcUtil/projection/StEmcPosition.h"
#include "StDaqLib/EMC/StEmcDecoder.h"
#include "StEmcRawMaker/defines.h"
#include "StDaqLib/TRG/trgStructures2007.h"
#include "StEvent/StTrackDetectorInfo.cxx"//added for track study


// 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 "StUpsilonTreeMaker.h"
#include "L2Result.h"

const double StUpsilonTreeMaker::BEMC_RADIUS = StEmcGeom::instance("bemc")->Radius();
const double StUpsilonTreeMaker::BEMC_OUTER_RADIUS = 253.799;
const double StUpsilonTreeMaker::BSMDE_RADIUS = StEmcGeom::instance("bsmde")->Radius();
const double StUpsilonTreeMaker::BSMDP_RADIUS = StEmcGeom::instance("bsmdp")->Radius();
const double StUpsilonTreeMaker::CTB_RADIUS = 213.6;

using namespace std;



ClassImp(StUpsilonTreeMaker);


int StUpsilonTreeMaker::Init()
{
  const char* name = "ups.tree.root";
  if (gSystem->Getenv("JOBID"))
    name = gSystem->ConcatFileName(gSystem->Getenv("SCRATCH"), name);
  mFile = new TFile(name, "recreate");
  assert(mFile && !mFile->IsZombie());
  mTree = new TTree("upsTree", "Upsilon Tree");
  mEvent = new Event;
  mTree->Branch("eventBranch", "Event", &mEvent, 32000, 99);
  mTree->BranchRef();		// Save reference table
  mEmcDecoder = new StEmcDecoder;
  mEmcPosition = new StEmcPosition;
  mEmcGeom = StEmcGeom::instance("bemc");
  mBemcTables = new StBemcTables;
  return StMaker::Init();

  /*
  float L2_E1 = upsL2Result->energyOfClusterL0;
  float L2_E2 = upsL2Result->energyOfClusterL2;
  float L2_m = upsL2Result->invMass;
  int TId_L0 = upsL2Result->triggerTowerL0;
  int TId_L2 = upsL2Result->triggerTowerL2;
  cout<<L2_E1<<" "<<L2_E2<<" "<<L2_m<<" "<<TId_L0<<" "<<TId_L2<<endl;
  */
  


}

int StUpsilonTreeMaker::InitRun(int runNumber)
{
  mEmcDecoder->SetDateTime(GetDate(), GetTime());
  mBemcTables->loadTables(this);
/*
  // Load BEMC status tables from Ahmed Hamed
  string filename = getenv("HOME");
  if (7090046 <= runNumber && runNumber <= 7129067)
  filename += "/heavy/upsilon/run7/STATUS_pp2006/pp2006_TriId117602.status";
  if (7133052 <= runNumber && runNumber <= 7156040)
  filename += "/heavy/upsilon/run7/STATUS_pp2006/pp2006_TriId137603.status";

//earlier lines are commented as they are for pp, we need new tables for Aurun7


  ifstream in(filename.c_str());
  int id, status;
  while (in >> id >> status) {
    assert(1 <= id && id <= 4800);
    assert(status == 1 || status == 0);
    mStatus[id] = status;
  }
  in.close();
*/
  return StMaker::InitRun(runNumber);
}

void StUpsilonTreeMaker::Clear(Option_t* option)
{
  mEvent->Clear(option);
  StMaker::Clear(option);
  mVecSoftIdTrig.clear();
}

int StUpsilonTreeMaker::Make()
{
  if (!GetDataSet("MuDst")) {
    LOG_WARN << "No MuDst" << endm;
    return kStWarn;
  }
  L2Result* upsL2Result = (L2Result*)&mEvent->L2Result()[L2RESULTS_2007_OFFSET_UPS];
  float L2_E1 = upsL2Result->energyOfClusterL0;
  float L2_E2 = upsL2Result->energyOfClusterL2;
  float L2_m = upsL2Result->invMass;
  int TId_L0 = upsL2Result->triggerTowerL0;
  int TId_L2 = upsL2Result->triggerTowerL2;
  cout<<L2_E1<<" "<<L2_E2<<" "<<L2_m<<" "<<TId_L0<<" "<<TId_L2<<endl;
  
  

  StEmcPosition emcPos;
  for (int deta1 = -1; deta1 <= 1; ++deta1) 
    {
    for (int dphi1 = -1; dphi1 <= 1; ++dphi1) 
      {
      int id11 = emcPos.getNextTowerId(TId_L2, deta1, dphi1);
      int id22 = emcPos.getNextTowerId(TId_L0, deta1, dphi1);
    //  cout<<"id1s= "<<id11<<" "<<endl;
    //  cout<<"id2s= "<<id22<<" "<<endl;
      if (id11!=0) mVecSoftIdTrig.push_back(id11);
      if (id22!=0) mVecSoftIdTrig.push_back(id22);
    }
  }
  /*
  cout << "Vector mVecSoftIdTrig = ( " ;
  for (vector<int>::iterator Iter1 = mVecSoftIdTrig.begin( ) ; Iter1 != mVecSoftIdTrig.end( ) ; Iter1++ )
    cout << *Iter1 << " ";
  cout << ")" << endl;
  */
  
   if (!accept(StMuDst::event())) return kStOk;

  mEvent->setEventInfo(StMuDst::event());
  makeEmcHits();
  makeSmdClusters();
  makeCtbHits();
  makeTracks();
  makeEmcClusters();
  makeParticles();

  mTree->Fill();

  return kStOk;
}


int StUpsilonTreeMaker::Finish()
{
  mFile->Write();
  mFile->Close();
  return kStOk;
}

void StUpsilonTreeMaker::makeEmcHits()
{

  TClonesArray& hits = mEvent->emcHits();
  for (int id = 1; id <= 4800; ++id) {
    //if (!mStatus[id]) continue;	// Use Ahmed's status tables
//now temporarily using default tables


/*
#if 1
    // Use status tables from BEMC database
    int status;
    mBemcTables->getStatus(BTOW, id, status);
    if (status != 1) continue;
#endif
*///this is what I had initally used (first pass in Oct 2007)

// Use online status tables from BEMC database
    int crate;
    int sequence;
    mEmcDecoder->GetCrateFromTowerId(id, crate, sequence);
    int status;
    mBemcTables->getTriggerTowerStatus(crate, sequence, status);
    if (status != 1) continue;
    //this is from Pibero//

//this is not used earlier as Ahmed's status tables used 
   int adc = StMuDst::muEmcCollection()->getTowerADC(id);
    if (!adc) continue;
    float pedestal, rms;
    mBemcTables->getPedestal(BTOW, id, 0, pedestal, rms);
    if (adc < pedestal + 2 * rms) continue;
    float calib;
    mBemcTables->getCalib(BTOW, id, 1, calib);
    float energy = calib * (adc - pedestal);
    if (energy <= 0) continue;

    // Make hit
    EmcHit* hit = new (hits[hits.GetEntriesFast()]) EmcHit;
    hit->setId(id);
    hit->setAdc(adc);

    // Set position
    float x, y, z;
    StEmcGeom::instance("bemc")->getXYZ(id, x, y, z);
    hit->setPosition(TVector3(x, y, z));

    // Set momentum
    TVector3 momentum = hit->position() + mEvent->vertex();
    momentum.SetMag(energy);
    hit->setMomentum(momentum);
  }
}

EmcHit* StUpsilonTreeMaker::getEmcHit(int id) const
{
  for (int i = 0; i < mEvent->emcHits().GetEntriesFast(); ++i) {
    EmcHit* hit = (EmcHit*)mEvent->emcHits().At(i);
    if (hit->id() == id) return hit;
  }
  return 0;
}

void StUpsilonTreeMaker::makeEmcClusters()
{
  TClonesArray& clusters = mEvent->emcClusters();
  bitset<4801> used;
  StEmcPosition emcPosition;
  TVector3 vertex(StMuDst::event()->primaryVertexPosition().xyz());
  for (int i = 0; i < mEvent->emcHits().GetEntriesFast(); ++i) {
  nextCluster:
    EmcHit* hit = (EmcHit*)mEvent->emcHits().At(i);
    int id = hit->id();
    if (used.test(id)) continue;
    if (hit->energy() < 1.0) continue;
    EmcCluster* cluster = new (clusters[clusters.GetEntriesFast()]) EmcCluster;
    cluster->hits().Add(hit);
    used.set(id);
    float energy = hit->energy();
    TVector3 position = hit->position();
    // Get neighbors
    for (int deta = -1; deta <= 1; ++deta) {
      for (int dphi = -1; dphi <= 1; ++dphi) {
	if (deta == 0 && dphi == 0) continue;
	int id2 = emcPosition.getNextTowerId(id, deta, dphi);
	if (!id2) continue;
	if (used.test(id2)) continue;
	EmcHit* hit2 = getEmcHit(id2);
	if (!hit2) continue;
	if (hit->energy() < hit2->energy()) {
	  // Seed tower is not highest. Remove cluster.
	  for (int j = 0; j < cluster->hits().GetEntriesFast(); ++j) {
	    EmcHit* hit3 = (EmcHit*)cluster->hits().At(j);
	    int id3 = hit3->id();
	    used.set(id3, 0);
	  }
	  clusters.Remove(cluster);
	  ++i;
	  goto nextCluster;
	}
	energy += hit2->energy();
	TVector3 position2 = hit2->position();
	position += position2 * hit2->energy();
	cluster->hits().Add(hit2);
	used.set(id2);
      }
    }

    // Cluster energy
    cluster->setEnergy(energy);

    // Cluster position
    position *= 1 / energy;
    cluster->position() = position;

    // Cluster momentum
    TVector3 momentum = position - vertex;
    momentum.SetMag(energy);
    cluster->momentum() = momentum;

    // Set cluster for tracks
    for (int j = 0; j < cluster->hits().GetEntriesFast(); ++j) {
      EmcHit* hit = (EmcHit*)cluster->hits().At(j);
      for (int k = 0; k < hit->tracks().GetEntriesFast(); ++k) {
	Track* track = (Track*)hit->tracks().At(k);
	track->setEmcCluster(cluster);
      }
    }

    // Remove empty entries
    clusters.Compress();
  }
}

void StUpsilonTreeMaker::makeCtbHits()
{
  TClonesArray& hits = mEvent->ctbHits();
  StCtbTriggerDetector& ctb = StMuDst::event()->ctbTriggerDetector();
  for (unsigned int tray = 0; tray < ctb.numberOfTrays(); ++tray) {
    for (unsigned int slat = 0; slat < ctb.numberOfSlats(); ++slat) {
      int adc = static_cast<int>(ctb.mips(tray, slat));
      if (adc > 0) {
	new (hits[hits.GetEntriesFast()]) CtbHit(tray, slat, adc);
      }
    }
  }
}

void StUpsilonTreeMaker::makeTracks()
{
  TClonesArray& tracks = mEvent->tracks();
  TIter nextTrack(StMuDst::primaryTracks());
  while (StMuTrack* muTrack = (StMuTrack*)nextTrack()) {
    if ( accept(muTrack) ) {
      Track* track = new (tracks[tracks.GetEntriesFast()]) Track(muTrack);
      track->setEvent(mEvent);
      getFirstEmcHit(track);
      getLastEmcHit(track);
      getEmcCluster(track);
      getSmdEtaCluster(track);
      getSmdPhiCluster(track);
      getCtbHit(track);
    }
  }
}

void StUpsilonTreeMaker::getFirstEmcHit(Track* track)
{
  try {
    TVector3 position = track->positionAtRadius(BEMC_RADIUS);
    int id;
    if (StEmcGeom::instance("bemc")->getId(position.Phi(), position.Eta(), id) == 0) {
      if (EmcHit* hit = getEmcHit(id)) {
	track->emcHits().Add(hit);
	hit->tracks().Add(track);
      }
    }
  }
  catch (Track::Exception& e) {}
}

void StUpsilonTreeMaker::getLastEmcHit(Track* track)
{
  try {
    TVector3 position = track->positionAtRadius(BEMC_OUTER_RADIUS);
    int id;
    if (StEmcGeom::instance("bemc")->getId(position.Phi(), position.Eta(), id) == 0) {
      if (EmcHit* hit = getEmcHit(id)) {
	track->emcHits().Add(hit);
	if (!hit->tracks().FindObject(track)) hit->tracks().Add(track);
      }
    }
  }
  catch (Track::Exception& e) {}
}

void StUpsilonTreeMaker::getEmcCluster(Track* track)
{
  // Done in makeEmcClusters() ???
}

void StUpsilonTreeMaker::getSmdEtaCluster(Track* track)
{
  try {
    TVector3 position = track->positionAtRadius(BSMDE_RADIUS);
    double minDistance = 999;
    SmdCluster* minCluster = 0;
    TIter nextCluster(&mEvent->smdEtaClusters());
    while (SmdCluster* cluster = (SmdCluster*)nextCluster()) {
      double deta = position.Eta() - cluster->position().Eta();
      double dphi = TVector2::Phi_mpi_pi(position.Phi() - cluster->position().Phi());
      double distance = hypot(deta, dphi);
      if (distance < 0.05 && distance < minDistance) {
	minDistance = distance;
	minCluster = cluster;
      }
    }
    if (minCluster) {
      track->setSmdEtaCluster(minCluster);
      minCluster->tracks().Add(track);
    }
  }
  catch (Track::Exception& e) {}
}

void StUpsilonTreeMaker::getSmdPhiCluster(Track* track)
{
  try {
    TVector3 position = track->positionAtRadius(BSMDP_RADIUS);
    double minDistance = 999;
    SmdCluster* minCluster = 0;
    TIter nextCluster(&mEvent->smdPhiClusters());
    while (SmdCluster* cluster = (SmdCluster*)nextCluster()) {
      double deta = position.Eta() - cluster->position().Eta();
      double dphi = TVector2::Phi_mpi_pi(position.Phi() - cluster->position().Phi());
      double distance = hypot(deta, dphi);
      if (distance < 0.05 && distance < minDistance) {
	minDistance = distance;
	minCluster = cluster;
      }
    }
    if (minCluster) {
      track->setSmdPhiCluster(minCluster);
      minCluster->tracks().Add(track);
    }
  }
  catch (Track::Exception& e) {}
}

void StUpsilonTreeMaker::getCtbHit(Track* track)
{
  try {
    TVector3 position = track->positionAtRadius(CTB_RADIUS);
    TIter nextCtbHit(&mEvent->ctbHits());
    while (CtbHit* ctbHit = (CtbHit*)nextCtbHit()) {
      double eta;
      double phi;
      ctbHit->getEtaPhi(eta, phi);
      double deta = position.Eta() - eta;
      double dphi = TVector2::Phi_mpi_pi(position.Phi() - phi);
      if (fabs(deta) < 0.25 && fabs(dphi) < 3 * degree) {
	track->setCtbHit(ctbHit);
	ctbHit->tracks().Add(track);
	break;
      }
    }
  }
  catch (Track::Exception& e) {}
}

void StUpsilonTreeMaker::makeParticles()
{
  TClonesArray& tracks = mEvent->tracks();
  TClonesArray& particles = mEvent->particles();
  if (tracks.GetEntriesFast() >= 2) {
    for (int i = 0; i < tracks.GetEntriesFast(); ++i) {
      Track* track1 = (Track*)tracks.At(i);
      for (int j = i + 1; j < tracks.GetEntriesFast(); ++j) {
	Track* track2 = (Track*)tracks.At(j);
	new (particles[particles.GetEntriesFast()]) Particle(track1, track2);
      }
    }
  }
}

int StUpsilonTreeMaker::makeSmdClusters()
{
  StEvent* event = (StEvent*)GetDataSet("StEvent");

  if (!event) {
    LOG_WARN << "No StEvent" << endm;
    return kStWarn;
  }

  StEmcCollection* emc = event->emcCollection();
  if (!emc) {
    LOG_WARN << "No EMC collection" << endm;
    return kStWarn;
  }

  // Get SMD-eta clusters
  if (StEmcDetector* det = emc->detector(kBarrelSmdEtaStripId)) {
    if (det->cluster())
      makeSmdEtaClusters(det->cluster()->clusters());
  }

  // Get SMD-phi clusters
  if (StEmcDetector* det = emc->detector(kBarrelSmdPhiStripId)) {
    if (det->cluster())
      makeSmdPhiClusters(det->cluster()->clusters());
  }

  return kStOk;
}

int StUpsilonTreeMaker::makeSmdEtaClusters(const StSPtrVecEmcCluster& clusters)
{
  TClonesArray& smdClusters = mEvent->smdEtaClusters();
  for (size_t i = 0; i < clusters.size(); ++i) {
    SmdCluster* cluster = new (smdClusters[smdClusters.GetEntriesFast()]) SmdCluster;
    cluster->setDetectorId(kBarrelSmdEtaStripId);
    cluster->setEnergy(clusters[i]->energy());
    TVector3 position;
    position.SetPtEtaPhi(BSMDE_RADIUS, clusters[i]->eta(), clusters[i]->phi());
    cluster->setPosition(position);
    TClonesArray& smdEtaHits = mEvent->smdEtaHits();
    StPtrVecEmcRawHit& hits = clusters[i]->hit();
    for (size_t j = 0; j < hits.size(); ++j) {
      SmdHit* hit = new (smdEtaHits[smdEtaHits.GetEntriesFast()]) SmdHit;
      hit->setDetectorId(kBarrelSmdEtaStripId);
      hit->setId(hits[j]->softId(BSMDE));
      hit->setEnergy(hits[j]->energy());
      cluster->hits().Add(hit);
    }
  }
  return kStOk;
}

int StUpsilonTreeMaker::makeSmdPhiClusters(const StSPtrVecEmcCluster& clusters)
{
  TClonesArray& smdClusters = mEvent->smdPhiClusters();
  for (size_t i = 0; i < clusters.size(); ++i) {
    SmdCluster* cluster = new (smdClusters[smdClusters.GetEntriesFast()]) SmdCluster;
    cluster->setDetectorId(kBarrelSmdPhiStripId);
    cluster->setEnergy(clusters[i]->energy());
    TVector3 position;
    position.SetPtEtaPhi(BSMDP_RADIUS, clusters[i]->eta(), clusters[i]->phi());
    cluster->setPosition(position);
    TClonesArray& smdPhiHits = mEvent->smdPhiHits();
    StPtrVecEmcRawHit& hits = clusters[i]->hit();
    for (size_t j = 0; j < hits.size(); ++j) {
      SmdHit* hit = new (smdPhiHits[smdPhiHits.GetEntriesFast()]) SmdHit;
      hit->setDetectorId(kBarrelSmdPhiStripId);
      hit->setId(hits[j]->softId(BSMDP));
      hit->setEnergy(hits[j]->energy());
      cluster->hits().Add(hit);
    }
  }
  return kStOk;
}

bool StUpsilonTreeMaker::accept(StMuEvent* event) const
{
  return (event->triggerIdCollection().nominal().isTrigger(200601) ||
	  event->triggerIdCollection().nominal().isTrigger(200602));
}

bool StUpsilonTreeMaker::accept(StMuTrack* track) const
{
  return (track->flag() > 0 &&
	  track->flag() != 701 &&
	  track->flag() != 801 &&
	  track->flag() != 901 &&
	  track->pt() > 0.2 &&
	  2.0 < track->p().mag() && track->p().mag() < 100 && //added from 3 to 100 GeV/c
	  track->nHitsFit() > 20 &&
	  track->nHitsFit() / (float)track->nHitsPoss() > 0.52 &&
	  3.5 * keV < track->dEdx() && track->dEdx() < 5.0 * keV);
}


bool StUpsilonTreeMaker::acceptSoftId(StMuTrack *track) const
{

  // StEvent *event = NULL;//initialize Stevent
  //event = (StEvent*)GetInputDS("StEvent");
    //StEvent* event = (StEvent*)GetDataSet("StEvent");
  /*  
  if (!event) { 
    LOG_WARN << "No StEvent" << endm;
    return kStWarn;
  }
*/
 
  StEmcPosition *pos = new StEmcPosition();
  //StMuEvent* event = new StMuEvent(); 
  float bFld1 = StMuDst::event()->magneticField();
  float bFld = bFld1/10.;//in Tesla
  //cout <<" magneticField = "<<bFld1<<" "<<bFld<<" (Tesla)"<< endl;
  /*
  Double_t bFld;
  StEventSummary* summary = event->summary();
  if(summary) {
    bFld = summary->magneticField(); // bFld in Tesla
    cout<<" StEvent::summary()->magneticField() = "<<bFld<<"(Tesla)"<< endl;
  }
 */
  StThreeVectorD position1, momentum1;
  Bool_t ok = pos->trackOnEmc(&position1, &momentum1, track, bFld); // bFld in Tesla        
  if(!ok) return false;
  double eta=position1.pseudoRapidity();
  double phi=position1.phi();
  //cout<<" eta position "<< eta <<" phi position "<< phi <<endl;
	 
  //tower id and module number///
  //Int_t id1;
  StEmcGeom* emcgeom = StEmcGeom::getEmcGeom("bemc");
  //emcgeom->getId(phi,eta,id1);

  Int_t module,eta1,sub1,id1;
  emcgeom->getBin(phi,eta,module,eta1,sub1);
  emcgeom->getId(module,eta1,sub1,id1);



  //cout<<"THE extrapolated TOWER ID ===="<<id1<<endl;

  //cout << "Function returned " << (bool) (find(mVecSoftIdTrig.begin(),mVecSoftIdTrig.end(),id1)!=mVecSoftIdTrig.end()) << endl;
  return (find(mVecSoftIdTrig.begin(),mVecSoftIdTrig.end(),id1)!=mVecSoftIdTrig.end());
}

