//
// Pibero Djawotho <pibero@indiana.edu>
// Indiana University
// 24 June 2008
//
// StUpsilonEmbedMaker runs over the Upsilon embedding and
// calculates acceptance, L0 and L2 trigger efficiencies
// for the Upsilon.
//
#include <iostream>
#include <fstream>

// Will not compile if moved under "local includes" ???
#include "StUpsilonTriggerMaker/StUpsilonTriggerMaker.h"
#include "StUpsilonTriggerMaker/Quarkonium.hh"
#include "StDaqLib/EMC/StEmcDecoder.h"

// ROOT includes
#include "TNtuple.h"
#include "TH1.h"
#include "TH2.h"

// STAR includes
#include "StEventTypes.h"
#include "StMcEvent/StMcEventTypes.hh"
#include "StAssociationMaker/StAssociationMaker.h"
#include "StAssociationMaker/StTrackPairInfo.hh"
#include "StEmcTriggerMaker/StEmcTriggerMaker.h"
#include "tables/St_g2t_event_Table.h"
#include "tables/St_g2t_pythia_Table.h"
#include "StEventUtilities/StuRefMult.hh"
#include "StEmcUtil/geometry/StEmcGeom.h"

// Local includes
//#include "StBbcTriggerMaker/StBbcTriggerMaker.h"
#include "StUpsilonEmbedMaker.h"

ClassImp(StUpsilonEmbedMaker);

float scaleFactor(double Eta, int hitType=0) {
//     hitType = 0 - towers
//     hitType = 1 - pre shower
//     hitType = 2 - shower max Eta
//     hitType = 3 - shower max Phi
    float P0[]={14.69,559.7,0.1185e6,0.1260e6};
    float P1[]={-0.1022,-109.9,-0.3292e5,-0.1395e5};
    float P2[]={0.7484,-97.81,0.3113e5,0.1971e5};
    
    float x=fabs(Eta);
    return P0[hitType]+P1[hitType]*x+P2[hitType]*x*x;
}

void StUpsilonEmbedMaker::SetHighTowerVar(StMcTrack* mcTrack, bool isele){
  StMcCalorimeterHit* MaxCalHit = 0;
  cout<<"Is Electron = "<<isele<<" and bemc size = "<<mcTrack->bemcHits().size()<<endl;
  //StPtrVecMcCalorimeterHit& theEmcHits =  mcTrack->bemcHits();
  double maxEnergy = 0;
  double sumE = 0;
  short unsigned int MaxRecoADC = 0;
  int MaxSoftId = 0;
  for (size_t i=0; i<mcTrack->bemcHits().size(); ++i) {
    StMcCalorimeterHit* calHit = mcTrack->bemcHits()[i];
    int modN = calHit->module();
    int etaN = calHit->eta();
    int subN = calHit->sub();
    float dE = calHit->dE();
    Float_t eta = 0;
    Float_t phi = 0;	
    StEmcGeom *geomBemc = StEmcGeom::getEmcGeom(1);
    geomBemc->getEta(modN,etaN,eta);
    geomBemc->getPhi(modN,subN,phi);
    // if (!inEmcGeometry(eta,phi)) continue;
    //sumEnergy += dE*scaleFactor(eta); // factor to convert adc -> GeV
    sumE = sumE + dE*scaleFactor(eta);
    if (dE*scaleFactor(eta) > maxEnergy){
      maxEnergy = dE*scaleFactor(eta);
      MaxCalHit = calHit;
    }
  }
  if (isele)
    eleEsum = sumE;
  else
    posEsum = sumE;

  if (MaxCalHit != 0){
    int MaxmodN = MaxCalHit->module();
    int MaxetaN = MaxCalHit->eta();
    int MaxsubN = MaxCalHit->sub();
    StEmcGeom *geomBemc = StEmcGeom::getEmcGeom(1);
    geomBemc->getId(MaxmodN, MaxetaN, MaxsubN, MaxSoftId);
    StEmcDecoder* EmcDecoder = new StEmcDecoder;
    int MaxDAQId;
    EmcDecoder->GetDaqIdFromTowerId(MaxSoftId,MaxDAQId);
    MaxRecoADC =  mUpsTrig->ReturnADCvalueFromDAQID(MaxDAQId);
    //cout<<"MaxRecoADC = "<<MaxRecoADC<<endl;
  }

  if (isele){
    eleADC = MaxRecoADC;eleE = maxEnergy;eleSoftId=MaxSoftId;}
  else{
    posADC = MaxRecoADC;posE = maxEnergy;posSoftId=MaxSoftId;}
    
  return;
}

int StUpsilonEmbedMaker::Init()
{
  // TString OutFile("CandidateUpsilon");
  //OutFile.Append(mRunNumber);
  //OutFile.Append("_");
  //OutFile.Append(mIndex);
  //OutFile.Append(".txt");
  //ofs.open(OutFile.Data());
  // Create tree list
  //float bemcArray[4801];

  mTreeList = new TList;
 
 PythiaVertex = new TH1F("PythiaVertex","Vertex of pythia event",2000,-200,200);
 PythiaVertex->GetXaxis()->SetTitle("Vz (cm)");
 PythiaVertex->GetYaxis()->SetTitle("# of counts");
 
 UpsilonVertex = new TH1F("UpsilonVertex","Start Vertex of Upsilon",2000,-200,200);
 UpsilonVertex->GetXaxis()->SetTitle("Vz (cm)");
 UpsilonVertex->GetYaxis()->SetTitle("# of counts");
 
 DeltaVertex = new TH1F("DeltaVertex","#Delta Vz between Pythia and Upsilon",2000,-200,200);
 DeltaVertex->GetXaxis()->SetTitle("Vz (cm) Pythia - #Upsilon Start");
 DeltaVertex->GetYaxis()->SetTitle("# of counts");
 
 CosinethetaRc = new TH1F("CosinethetaRc","reco Cosine Theta",200,-2,2);
 CosinethetaRc->GetXaxis()->SetTitle("Cosine Theta Reco");
 CosinethetaRc->GetYaxis()->SetTitle("# of counts");
 
 
 // Create ntuples
  TString varlist = "vx:vy:vz:pid:ntracks:upsGeantId:upsM:upsCosTheta:upsRap:upsP:upsPt:upsPz:upsEta:upsPhi:eleP:elePt:elePz:eleEta:elePhi:posP:posPt:posPz:posEta:posPhi:isTrig200601:isTrig200602:isTrigUps:isTrigBbc:eleADCreco:posADCreco:eleERc:posERc:upsFound:eleBEMChits:posBEMChits:eleEsum:posEsum:eleTPChits:posTPChits:eleTPChitspos:posTPChitspos:elePRc:posPRc:elePtRc:posPtRc:cosThetaRc:l2eleE:l2posE:l2invmass:l2costheta:UpsMRc";
  mTuple = new TNtuple("ups", "Upsion ntuple", varlist);
  mTreeList->Add(mTuple);
  mTreeList->Add(PythiaVertex);
  mTreeList->Add(UpsilonVertex);
  mTreeList->Add(DeltaVertex);
  mTreeList->Add(CosinethetaRc);
  return StMaker::Init();
}

int StUpsilonEmbedMaker::Make()
{
  const double emass = 0.511e-3;
  float l2eleE= -999;
  float l2posE= -999;
  // Get StMcEvent
  mMcEvent = (StMcEvent*)GetDataSet("StMcEvent");
  if (!mMcEvent) {
    LOG_WARN << "No StMcEvent" << endm;
    return kStWarn;
  }
 
  // Print StMcEvent
  LOG_INFO << "eventGeneratorEventLabel = " << mMcEvent->eventGeneratorEventLabel() << endm;
  LOG_INFO << "eventNumber = " << mMcEvent->eventNumber() << endm;
  LOG_INFO << "runNumber = " << mMcEvent->runNumber() << endm;
  LOG_INFO << "type = " << mMcEvent->type() << endm;
  LOG_INFO << "zWest = " << mMcEvent->zWest() << endm;
  LOG_INFO << "nWest = " << mMcEvent->nWest() << endm;
  LOG_INFO << "zEast = " << mMcEvent->zEast() << endm;
  LOG_INFO << "nEast = " << mMcEvent->nEast() << endm;
  LOG_INFO << "eventGeneratorFinalStateTracks = " << mMcEvent->eventGeneratorFinalStateTracks() << endm;
  LOG_INFO << "numberOfPrimaryTracks = " << mMcEvent->numberOfPrimaryTracks() << endm;
  LOG_INFO << "subProcessId = " << mMcEvent->subProcessId() << endm;
  LOG_INFO << "impactParameter = " << mMcEvent->impactParameter() << endm;
  LOG_INFO << "phiReactionPlane = " << mMcEvent->phiReactionPlane() << endm;
  LOG_INFO << "triggerTimeOffset = " << mMcEvent->triggerTimeOffset() << endm;
  LOG_INFO << "nBinary = " << mMcEvent->nBinary() << endm;
  LOG_INFO << "nWoundedEast = " << mMcEvent->nWoundedEast() << endm;
  LOG_INFO << "nWoundedWest = " << mMcEvent->nWoundedWest() << endm;
  LOG_INFO << "nJets = " << mMcEvent->nJets() << endm;

  // Get Pythia record
  TDataSet* geant = GetDataSet("geant");
  TDataSetIter geantIter(geant);

  St_g2t_pythia* pythia = (St_g2t_pythia*)geantIter("g2t_pythia");

  if (pythia) {
    g2t_pythia_st* g2t_pythia = pythia->GetTable();
    assert(g2t_pythia);
    LOG_INFO << "Pythia subprocess_id = " << g2t_pythia->subprocess_id << endm;
  }
  else {
    LOG_WARN << "No Pythia record" << endm;
  }

  // Get StEvent
  mEvent = (StEvent*)GetDataSet("StEvent");
  if (!mEvent) {
    LOG_WARN << "No StEvent" << endm;
    return kStWarn;
  }

  // Get BBC trigger maker
  //StBbcTriggerMaker* bbcTrig = (StBbcTriggerMaker*)GetMakerInheritsFrom("StBbcTriggerMaker");
  //assert(bbcTrig);

  //LOG_INFO << "BBC = " << bbcTrig->isTrigger() << endm;

  // Get EMC trigger maker
  StEmcTriggerMaker* emcTrig = (StEmcTriggerMaker*)GetMakerInheritsFrom("StEmcTriggerMaker");
  assert(emcTrig);

  LOG_INFO << "Upsilon(200601) = " << emcTrig->isTrigger(200601) << endm;
  LOG_INFO << "Upsilon(200602) = " << emcTrig->isTrigger(200602) << endm;

  // Get Upsilon trigger
  mUpsTrig = (StUpsilonTriggerMaker*)GetMakerInheritsFrom("StUpsilonTriggerMaker");
  assert(mUpsTrig);

  LOG_INFO << "Upsilon(L2) = " << mUpsTrig->isUpsilonTrigger() << endm;

  // Print particles in event
  LOG_INFO << "Subprocess Id = " << mMcEvent->subProcessId() << endm;
  LOG_INFO << "Ntracks = " << mMcEvent->numberOfPrimaryTracks() << endm;
  cout << "Particles: ";
  for (size_t i = 0; i < mMcEvent->tracks().size(); ++i) {
    StMcTrack* mcTrack = mMcEvent->tracks()[i];
    if (mcTrack->particleDefinition())
      cout << mcTrack->particleDefinition()->name() << " ";
    else
      cout << mcTrack->geantId() << " ";
  }
  cout << endl;

  StEmcDetector* bemcDet = mEvent->emcCollection()->detector(kBarrelEmcTowerId);
  for (int i = 0;i<4801;i++){    bemcArray[i] = 0;}

  for (unsigned int m = 1; m<=bemcDet->numberOfModules(); ++m){
    StSPtrVecEmcRawHit& hits = bemcDet->module(m)->hits();
    for (StSPtrVecEmcRawHitIterator i = hits.begin(); i != hits.end(); ++i){
      StEmcRawHit* hit = *i;
      if (hit->energy()<=0) continue;
      int softId;
      int modN = hit->module();
      int etaN = hit->eta();
      int subN = hit->sub();
      StEmcGeom *geomBemc = StEmcGeom::getEmcGeom(1);
      geomBemc->getId(modN, etaN, subN, softId);
      bemcArray[softId] = hit->energy();
      cout<<"softid "<<softId<<" has energy "<<hit->energy()<<endl;
    }}



  // Loop over primary tracks CHANGED FOR ANDREWS EMBEDDING 4/14/09
  //Upsilon does not appear to be defined as a primary vertex
  //is the upsilon vertex the same as the primary vertex?
  //for (size_t i = 0; i < mMcEvent->tracks().size(); ++i) {
  //StMcTrack* upsMc = mMcEvent->tracks()[i];

  // Loop over primary tracks
  for (size_t i = 0; i < mMcEvent->primaryVertex()->numberOfDaughters(); ++i) {
    StMcTrack* upsMc = mMcEvent->primaryVertex()->daughter(i);
    
    
    // Get Upsilon
    if (upsMc->geantId() >= 160) {
      cout<<"Found Upsilon with geant ID"<<upsMc->geantId()<<endl;
      // Loop over daughters and asign electron and positron
      StMcTrack* eleMc = 0;
      StMcTrack* posMc = 0;
      if (!upsMc->stopVertex())
        {
          cout<<"Found upsilon without stop vertex, bailing out"<<endl;
          return kStOk;
        }
      
      StMcVertex* UpsStartVertex = upsMc->startVertex();
    
      PythiaVertex->Fill(mMcEvent->primaryVertex()->position().z());
      
      UpsilonVertex->Fill(UpsStartVertex->position().z());
      
      DeltaVertex->Fill(mMcEvent->primaryVertex()->position().z() - UpsStartVertex->position().z());
      
      cout<<"UPS VERTEX = "<<UpsStartVertex->position()<<endl;
      cout<<"PRIMARY VERTEX = "<<mMcEvent->primaryVertex()->position()<<endl;
      
      for (size_t j = 0; j < upsMc->stopVertex()->numberOfDaughters(); ++j) {
	StMcTrack* mcTrack = upsMc->stopVertex()->daughter(j);
	
	switch (mcTrack->geantId()) {
	case 2:
	  posMc = mcTrack;
	  break;
	case 3:
	  eleMc = mcTrack;
	  break;
	default:
	  LOG_WARN << "Not an electron nor a positron - geantId = " << mcTrack->geantId() << endm;
	  return kStWarn;
	}
      }
      cout<<"USING THIS ONE!!!"<<endl;

      // Opening angle between daughters
      float cosThetaMc = eleMc->momentum().unit() * posMc->momentum().unit();
      SetHighTowerVar(eleMc,true);
      SetHighTowerVar(posMc,false);
      int nCandidates = ups.candidates.size();
      cout<<"# candidates in embed = "<<nCandidates<<endl;
      if (nCandidates > 0){
	cout<<"first = "<<ups.candidates[0].first<<" second = "<<ups.candidates[0].second<<endl;}
      bool isFound = false;
      if (nCandidates >0){
	if ((ups.candidates[0].first == eleSoftId) && (ups.candidates[0].second == posSoftId)){
	  isFound = true;
	  l2eleE = ups.mEnergyDaughter1;
	  l2posE = ups.mEnergyDaughter2;}
	if((ups.candidates[0].second == eleSoftId) && (ups.candidates[0].first == posSoftId)){
	  isFound = true;
	  l2eleE = ups.mEnergyDaughter2;
	  l2posE = ups.mEnergyDaughter1;
	}}

      //**********

      StTrack* eleRc = findPartner(eleMc);
      StTrack* posRc = findPartner(posMc);
    
      /*
      float cosThetaRc = -999;
      float eleRcp = -999;
      float eleRcpt = -999;
      float posRcp = -999;
      float posRcpt = -999;
      int eleTPCpos = -999;
      int posTPCpos = -999;
      float invMRc = -1000;
      */

      float cosThetaRc;
      float eleRcp;
      float eleRcpt;
      float posRcp;
      float posRcpt;
      int eleTPCpos;
      int posTPCpos;
      float invMRc;
      
      if (eleRc && posRc) {

	StThreeVectorF eleMom = eleRc->geometry()->momentum();
	StThreeVectorF posMom = posRc->geometry()->momentum();


        cout<<"EleMoM.unit"<<eleMom.unit()<<endl;

        //float energyelectron = eleRc->energy();

	StLorentzVectorF p1(eleMom,eleMom.massHypothesis(emass));
        StLorentzVectorF p2(posMom,posMom.massHypothesis(emass));
        StLorentzVectorF p3 = p1 + p2;
        invMRc = p3.m();



	cosThetaRc = eleMom.unit()*posMom.unit();
        cout<<"reco angle"<<cosThetaRc<<endl;  

	CosinethetaRc->Fill(cosThetaRc);
	
	eleRcp = eleMom.mag();
	posRcp = posMom.mag();
	eleRcpt = eleMom.perp();
	posRcpt = posMom.perp();
	eleTPCpos = eleRc->numberOfPossiblePoints(kTpcId);
	posTPCpos = posRc->numberOfPossiblePoints(kTpcId);


      }


      //**********

      // Fill tuple
      float tuple[] = {
	mMcEvent->primaryVertex()->position().x(),
	mMcEvent->primaryVertex()->position().y(),
	mMcEvent->primaryVertex()->position().z(),
	mMcEvent->subProcessId(),
	mMcEvent->numberOfPrimaryTracks(),

	upsMc->geantId(),
	upsMc->fourMomentum().m(),
	cosThetaMc,
	upsMc->rapidity(),
	upsMc->momentum().mag(),
	upsMc->pt(),
	upsMc->momentum().z(),
	upsMc->pseudoRapidity(),
	upsMc->momentum().phi(),

	eleMc->momentum().mag(),
	eleMc->pt(),
	eleMc->momentum().z(),
	eleMc->pseudoRapidity(),
	eleMc->momentum().phi(),

	posMc->momentum().mag(),
	posMc->pt(),
	posMc->momentum().z(),
	posMc->pseudoRapidity(),
	posMc->momentum().phi(),

	emcTrig->isTrigger(200601),
	emcTrig->isTrigger(200602),
	mUpsTrig->isUpsilonTrigger(),
	-1,//taking bbcTrig out because it appears not to be working with
	//new embeding in pdsf 4/14/09
	//bbcTrig->isTrigger(),
	//rosi added stuff
	eleADC,posADC,eleE,posE,isFound,
	eleMc->bemcHits().size(),
	posMc->bemcHits().size(),
	eleEsum,
	posEsum,
	eleMc->tpcHits().size(),
	posMc->tpcHits().size(),
	eleTPCpos,posTPCpos,
	eleRcp,posRcp,eleRcpt,posRcpt,
	cosThetaRc,
	l2eleE,l2posE,ups.mInvMass,ups.mCosTheta,
	invMRc
	
      };
      //for (int i = 0;i<(sizeof(tuple)/sizeof(*tuple)); i++){
      //cout<<"tuple["<<i<<"] = "<<tuple[i]<<endl;
      //}
      //char temp;
      //cin>>temp;

      mTuple->Fill(tuple);
      //if (nCandidates>0){
     
	//  ofs<<"Run#, Event#, EmbedJob = "<<mEvent->runId()<<", "<<mEvent->id()<<", "<<mIndex<<endl;
      //ofs<<"# candidates = "<<nCandidates<<endl;
      //for (int i = 0;i<nCandidates;i++){
      //ofs<<ups.candidates[i].first<<","<<ups.candidates[i].second<<endl;}
      //ofs<<"eleSoftId = "<<eleSoftId<<" posSoftId = "<<posSoftId<<endl;
      //ofs<<"l2eleE = "<<l2eleE<<" l2posE = "<<l2posE<<endl;
      //ofs<<"l2mEnergyDaughter1 = "<<ups.mEnergyDaughter1<<"l2mEnergyDaughter2 = "<<ups.mEnergyDaughter2<<endl;
      //ofs<<"*"<<endl;}
    }
  }

  return kStOk;
}

StTrack* StUpsilonEmbedMaker::findPartner(StMcTrack* mcTrack)
{
  StAssociationMaker* assoc = (StAssociationMaker*)GetMakerInheritsFrom("StAssociationMaker");
  pair<mcTrackMapIter, mcTrackMapIter> p = assoc->mcTrackMap()->equal_range(mcTrack);
  StTrack* maxTrack = 0;
  int maxCommonTpcHits = 0;
  for (mcTrackMapIter k = p.first; k != p.second; ++k) {
    int commonTpcHits = k->second->commonTpcHits();
    StTrack* track = k->second->partnerTrack()->node()->track(primary);
    if (track && commonTpcHits > maxCommonTpcHits) {
      maxTrack = track;
      maxCommonTpcHits = commonTpcHits;
    }
  }
  return maxTrack;
}