// 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"
#include "StEmcUtil/projection/StEmcPosition.h"
#include "StContainers.h"
#include "StEmcUtil/database/StBemcTables.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
  //function calculates the proper scale factor given a hit type and an
  //eta so that energy can be calculated
    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){
  //This function finds the soft tower id of the highest tower a particular
  //track has.  It also sets class veriables eleEsum and posEsum, which are
  //the sum of the energies left in the BEMC by the electron and positron
  //daughters.  I would like to change this to the highest 3 towers, but then
  //I need to think about how to do acceptance.  It also records the maximum
  //ADC values and energy for the electron and positron
  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()
{
  mEmcDecoder = new StEmcDecoder;
  mBemcTables = new StBemcTables;



//Candidate upsilon file was something used for checking errors I had initiall
  //can be removed
  TString OutFile("CandidateUpsilon");
  OutFile.Append(mRunNumber);
  OutFile.Append("_");
  OutFile.Append(mIndex);
  OutFile.Append(".txt");
  ofs.open(OutFile.Data());
  // Create tree list
  mTreeList = new TList;

  // 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:eleADCreco:posADCreco:eleERc:posERc:upsFound:eleBEMChits:posBEMChits:eleEsum:posEsum:eleTPChits:posTPChits:eleTPChitspos:posTPChitspos:elePRc:posPRc:elePtRc:posPtRc:cosThetaRc:L0isfound:eleFlag:posFlag:nCandidate:elereco:posreco:L0Ecluster:L0eta:L0phi:L2Ecluster:L2eta:L2phi:eleRcEtaProj:eleRcPhiProj:posRcEtaProj:posRcPhiProj:upsMRc:eleTP:posTP:l2eleE:l2posE:l2invmass:l2costheta:UpsMRc";
  mTuple = new TNtuple("ups", "Upsion ntuple", varlist);
  mTreeList->Add(mTuple);

  return StMaker::Init();
}

int StUpsilonEmbedMaker::InitRun(int runNumber)
{
  mEmcDecoder->SetDateTime(GetDate(), GetTime());
  mBemcTables->loadTables(this);
  return StMaker::InitRun(runNumber);
}

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(117602) << endm;
  LOG_INFO << "Upsilon(200602) = " << emcTrig->isTrigger(137602) << endm;
  //LOG_INFO << "Upsilon(137603) = " << emcTrig->isTrigger(137603) << 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->primaryVertex()->numberOfDaughters(); ++i) {
    StMcTrack* mcTrack = mMcEvent->primaryVertex()->daughter(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;bemcData[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;
      //Int_t daqId;
      //Float_t pedestal, rms
      //mEmcDecoder->GetDaqIdFrom TowerId(softId,daqId);
      // mBemcTables->getPedestal(BTOW,softId,0,pedestal,rms);
      //bemcData[softId] = hit->adc() - pedestal;
      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);
      int status;
      float pedestal; float rms;
      //mEmcDecoder->GetDaqIdFrom TowerId(softId,daqId);
     
      mBemcTables->getStatus(1, softId, status);
      mBemcTables->getPedestal(1,softId,0,pedestal,rms);
      //cout<<"softid = "<<softId<<" adc = "<<hit->adc()<<" ped = "<<pedestal<<" energy = "<<hit->energy()<<endl;
      bemcData[softId] = hit->adc() - pedestal;
      //cout<<"bemcData["<<softId<<"] = "<<bemcData[softId]<<endl;
      bemcArray[softId] = hit->energy();
      //cout<<"softid "<<softId<<" has energy "<<hit->energy()<<endl;
    }}

  // 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) {
      StMcVertex* UpsStartVertex = upsMc->startVertex();
      cout<<"UPS VERTEX = "<<UpsStartVertex->position()<<endl;
      cout<<"PRIMARY VERTEX = "<<mMcEvent->primaryVertex()->position()<<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;
        }

      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 nL0Candidates = ups.L0candidates.size();
      int nCandidates1 = ups.candidates.size();

       cout<<"# candidates in embed = "<<nCandidates1<<endl;
      if (nCandidates1 > 0){
	cout<<"first = "<<ups.candidates[0].first<<" second = "<<ups.candidates[0].second<<endl;}
      bool isFound1 = false;
      if (nCandidates1 >0){
	if ((ups.candidates[0].first == eleSoftId) && (ups.candidates[0].second == posSoftId)){
	  isFound1 = true;
	  l2eleE = ups.mEnergyDaughter1;
	  l2posE = ups.mEnergyDaughter2;}
	if((ups.candidates[0].second == eleSoftId) && (ups.candidates[0].first == posSoftId)){
	  isFound1 = true;
	  l2eleE = ups.mEnergyDaughter2;
	  l2posE = ups.mEnergyDaughter1;
	}}



      bool isL0 = false;
      float Ecluster2 = -1000;float Ecluster1 = -1000;
      pair<float,float> etaphi1;pair<float,float> etaphi2;
      etaphi1.first = -1000;etaphi1.second = -1000;
      etaphi2.first = -1000;etaphi2.second = -1000;
      //      int TP_ADC = -1000; 

      for (int i = 0;i<nL0Candidates;i++){
	if ((eleSoftId == ups.L0candidates[i]) || (posSoftId == ups.L0candidates[i])){
	  //cout<<"L0 ID IS "<<ups.L0candidates[i]<<" with adc = "<<bemcData[ups.L0candidates[i]];
	  isL0 = true;
	  //TP_ADC = getTPADC(ups.L0candidates[i]);
	  //cout<<" TP ADC = "<<TP_ADC<<endl;
	}
      }
      int nCandidates = ups.candidates.size();
      cout<<"# candidates in embed = "<<nCandidates<<endl;
      //for (int i = 0; i <nCandidates;i++ ){
      //cout<<"first = "<<ups.candidates[i].first<<" second = "<<ups.candidates[i].second<<endl;}
      bool isFound = false;
      for (int i = 0;i <nCandidates; i++){
	if ((ups.candidates[i].first == eleSoftId) && (ups.candidates[i].second == posSoftId)){
	  isFound = true;
	  getClusterEtaPhiE(ups.candidates[i].first,etaphi1,Ecluster1);
	  getClusterEtaPhiE(ups.candidates[i].second,etaphi2,Ecluster2);
	  //  TP_ADC = getTPADC(ups.candidates[i].first);
	  //cout<<"candidate "<<ups.candidates[i].first<<" "<<ups.candidates[i].second<< "with adc = "<<bemcData[ups.candidates[i].first]<<" "<<bemcData[ups.candidates[i].second]<<" and TP = "<<TP_ADC<<endl;
	}
	if((ups.candidates[i].second == eleSoftId) && (ups.candidates[i].first == posSoftId)){
	  isFound = true;
	  getClusterEtaPhiE(ups.candidates[i].first,etaphi1,Ecluster1);
	  getClusterEtaPhiE(ups.candidates[i].second,etaphi2,Ecluster2);
	  //TP_ADC = getTPADC(ups.candidates[i].first);
	  //cout<<"candidate "<<ups.candidates[i].first<<" "<<ups.candidates[i].second<< "with adc = "<<bemcData[ups.candidates[i].first]<<" "<<bemcData[ups.candidates[i].second]<<" and TP = "<<TP_ADC<<endl;
	}}

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

      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;
      short eleFlag = -999;
      short posFlag = -999;
      float eleTPChits = -1000;
      float posTPChits = -1000;
      bool elereco = false;
      bool posreco = false;
      pair<float,float> eleetaphiRcTrack;
      pair<float,float> posetaphiRcTrack;
      eleetaphiRcTrack.first = -1000;
      eleetaphiRcTrack.second = -1000;
      posetaphiRcTrack.first = -1000;
      posetaphiRcTrack.second = -1000;
      float invMRc = -1000;

      if (eleRc)
	elereco = true;
      if (posRc)
	posreco = true;

      
      if (eleRc && posRc) {
	eleFlag = eleRc->flag();
	posFlag = posRc->flag();
	StThreeVectorF eleMom = eleRc->geometry()->momentum();
	StThreeVectorF posMom = posRc->geometry()->momentum();
	StLorentzVectorF p1(eleMom,eleMom.massHypothesis(emass));
        StLorentzVectorF p2(posMom,posMom.massHypothesis(emass));
        StLorentzVectorF p3 = p1 + p2;
	invMRc = p3.m();

	cosThetaRc = eleMom.unit() * posMom.unit();
	eleRcp = eleMom.mag();
	posRcp = posMom.mag();
	eleRcpt = eleMom.perp();
	posRcpt = posMom.perp();
	eleTPCpos = eleRc->numberOfPossiblePoints(kTpcId);
	posTPCpos = posRc->numberOfPossiblePoints(kTpcId);
	eleTPChits = eleRc->fitTraits().numberOfFitPoints(kTpcId);
	posTPChits = posRc->fitTraits().numberOfFitPoints(kTpcId);
	TrackEtaPhi(eleRc,eleetaphiRcTrack);
	TrackEtaPhi(posRc,posetaphiRcTrack);

      }


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

      // 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),
	//emcTrig->isTrigger(137603),
	mUpsTrig->isUpsilonTrigger(),
	//bbcTrig->isTrigger(),
	//rosi added stuff
	/*eleADC*/bemcData[eleSoftId],/*posADC*/bemcData[posSoftId],
	eleE,posE,isFound,
	eleMc->bemcHits().size(),
	posMc->bemcHits().size(),
	eleEsum,
	posEsum,
	eleTPChits,
	posTPChits,
	eleTPCpos,posTPCpos,
	eleRcp,posRcp,eleRcpt,posRcpt,cosThetaRc,
	isL0,eleFlag,posFlag,nCandidates,elereco,posreco,
	Ecluster1,etaphi1.first,etaphi1.second,
	Ecluster2,etaphi2.first,etaphi2.second,
	eleetaphiRcTrack.first,eleetaphiRcTrack.second,
	posetaphiRcTrack.first,posetaphiRcTrack.second,
	invMRc,getTPADC(eleSoftId),getTPADC(posSoftId),
        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;
}

void StUpsilonEmbedMaker::getClusterEtaPhiE(int& id, pair<float,float>& etaphi, float& Ecluster)
{
  cout<<"In getclusteretaphi ";
  if (id == 0){
    cout<<"soft id is zero!!! exit etaphi"<<endl;
    return;}
  Float_t E1 = bemcArray[id];
  cout<<"got hightower energy = "<<E1<<" with soft id = "<<id<<endl;
  vector<int> neighbors;
  StEmcPosition* mEmcPosition = new StEmcPosition;
  for (int deta = -1; deta <= 1; ++deta) {
    for (int dphi = -1; dphi <= 1; ++dphi) {
      if (deta || dphi) {
        int id2 = mEmcPosition->getNextTowerId(id, deta, dphi);
        //cout<<"id2 = "<<id2<<endl;
        if (id2) {
          neighbors.push_back(id2);
        }
      }
    }
  }
  delete mEmcPosition;
  //cout<<"found neighbors"<<endl;
  float E2 = 0;
  float E3 = 0;
  int ID2 = 0;
  int ID3 = 0;

  for (size_t i = 0;i<neighbors.size();i++){
    int idcur = neighbors[i];
    float Ecur = bemcArray[idcur];
    if (Ecur > E2){
      E3 = E2;
      ID3 = ID2;
      E2 = Ecur;
      ID2 = idcur;}
    else if (Ecur > E3){
      E3 = Ecur;
      ID3 = idcur;}}
  StEmcGeom *geomBemc = StEmcGeom::getEmcGeom(1);
  Ecluster = E1+E2+E3;
  float eta1;float eta2;float eta3;
  float phi1;float phi2;float phi3;
  geomBemc->getEta(id,eta1);
  geomBemc->getPhi(id,phi1);
  if (ID2 > 0){
    geomBemc->getEta(ID2,eta2);
    geomBemc->getPhi(ID2,phi2);}
  if (ID3 > 0){
    geomBemc->getEta(ID3,eta3);
    geomBemc->getPhi(ID3,phi3);}
  if ((E1+E2+E3)>0){
    etaphi.first = (eta1*E1+eta2*E2+eta3*E3)/Ecluster;
    etaphi.second = (phi1*E1+phi2*E2+phi3*E3)/Ecluster;}
  else{
    etaphi.first = -1000;etaphi.second = -1000;}
  //cout<<"Finished etaphi"<<endl;
  return;
}

void StUpsilonEmbedMaker :: TrackEtaPhi(StTrack* mcTrack,pair<float,float>& etaphi)
{
  //cout<<"in track etaphi StTrack ";
  StEmcPosition* projPos = new StEmcPosition;
  StThreeVectorD* PosVec = new StThreeVectorD;
  StThreeVectorD* MomVec = new StThreeVectorD;
  Double_t mMagneticField = mEvent->summary()->magneticField()*0.1;//kilogause/tesla
  //cout<<"Defined variables"<<endl;
  projPos->projTrack(PosVec,MomVec,mcTrack,mMagneticField);
  //cout<<"called projection"<<endl;
  etaphi.first = PosVec->pseudoRapidity();
  etaphi.second = PosVec->phi();
  delete projPos;
  delete PosVec;
  delete MomVec;
  return;
}

int StUpsilonEmbedMaker :: getTPADC(int id){
  int TPid;int TPADC = 0;
  StEmcDecoder* EmcDecoder = new StEmcDecoder;
  EmcDecoder->GetTriggerPatchFromTowerId(id,TPid);
  for (int softId = 1;softId<=4800;++softId){
    int TPidcur;
    EmcDecoder->GetTriggerPatchFromTowerId(softId,TPidcur);
    if (TPidcur == TPid) TPADC+=bemcData[softId];}
  
  return TPADC;
}