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

// C++ STL
#include <bitset>
#include <algorithm>
// ROOT
#include "TFile.h"
#include "TTree.h"
#include "TNtuple.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 "StDaqLib/TRG/trgStructures2008.h"
#include "StEvent/StTrackDetectorInfo.cxx"//added for track study

// Local
//#include "struct_Event.h"
#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"
#include "SurroundingTowers.h"//added to access the surrounding towers

#include "StHeavyTrigger/Upsilon.hh"// this is added to test
  #include "read_in_geog.h"

extern Upsilon ups;


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());

  his_dedx_vs_p_all_primary=new TH2F("his_dedx_vs_p_all_primary","his_dedx_vs_p_all_primary",24*100,-12,12,4*200,1,5);
  his_dedx_vs_p_track_with_emc=new TH2F("his_dedx_vs_p_track_with_emc","his_dedx_vs_p_track_with_emc",24*100,-12,12,4*200,1,5);
  his_eta_vs_nHitFit_all_primary=new TH2F("his_eta_vs_nHitFit_all_primary","his_eta_vs_nHitFit_all_primary",49,0.5,49.5,1000,-5,5);
  his_eta_vs_nHitFit_track_with_emc=new TH2F("his_eta_vs_nHitFit_track_with_emc","his_eta_vs_nHitFit_track_with_emc",49,0.5,49.5,1000,-5,5);
  delta_eta_vs_vz_highpt=new TH2F("delta_eta_vs_vz_highpt","delta_eta_vs_vz_highpt",200,-200,200,200,-1,1);
  delta_eta_vs_vz_highpt_after_correction=new TH2F("delta_eta_vs_vz_highpt_after_correction","delta_eta_vs_vz_highpt_after_correction",200,-200,200,200,-1,1);

tower_eta_tower_phi_run=new TH3F("tower_eta_tower_phi_run","tower_eta_tower_phi_run",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);
track_eta_track_phi_run=new TH3F("track_eta_track_phi_run","track_eta_track_phi_run",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);
track_proj_eta_track_phi_run=new TH3F("track_proj_eta_track_phi_run","track_proj_eta_track_phi_run",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);
match_eta_match_phi_run=new TH3F("match_eta_match_phi_run","match_eta_match_phi_run",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);
match_eta_match_phi_run_check=new TH3F("match_eta_match_phi_run_check","match_eta_match_phi_run_check",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);
my_own_match_eta_match_phi_run=new TH3F("my_own_match_eta_match_phi_run","my_own_match_eta_match_phi_run",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);
//my_own_match_eta_match_phi_run2=new TH3F("my_own_match_eta_match_phi_run2","my_own_match_eta_match_phi_run2",44,-1.1,1.1,124,-3.255,3.255,700,8330000,9030000);

  tower_fired_softid = new TH1F("tower_fired_softid", "tower_fired_softid", 4801, -0.5, 4800.5);
  track_proj_softid = new TH1F("track_proj_softid", "track_proj_softid", 4801, -0.5, 4800.5);
  matched_softid = new TH1F("matched_softid", "matched_softid", 4801, -0.5, 4800.5);
 

//check_matched_entries=new TH2F("check_matched_entries","check_matched_entries",44,-1.1,1.1,124,-3.255,3.255);
//match_id_correlation_geomID_vs_towerID=new TH2F("match_id_correlation_geomID_vs_towerID","match_id_correlation_geomID_vs_towerID",4801,-0.5,4800.5,4801,-0.5,4800.5);

  // 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);
  // Create histograms
//  mTriggerIds = new TH1F("triggerIds", "Trigger IDs", 5, 0, 5);
  mRefMult = new TH1F("refMult", "Reference multiplicity", 100, 0, 100);
  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);

   mEvent = new Event;
//  mTree = new TTree("upsTree", "Upsilon Tree");
//  mTree->Branch("eventBranch", "Event", &mEvent, 32000, 99);
//  mTree->BranchRef();		// Save reference table

//   myevent=new struct_myevent;
  myTree = new TTree("event_tree", "event_tree");
  definetree();
//	tree->Branch("trkId",&m_trkId,"trkId/I");
  //myTree->Branch("eventBranch", "struct_Event", &myevent, 32000, 99);
 // myTree->BranchRef();		// Save reference table

  mEmcDecoder = new StEmcDecoder;
  mEmcPosition = new StEmcPosition;
  mEmcGeom = StEmcGeom::instance("bemc");
  mBemcTables = new StBemcTables;
  return StMaker::Init();
  
}

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();
  mTrigtowers.clear();
  mStMuTrackTowerMap.clear();
  mTrackTowerMap.clear();
}

int StUpsilonTreeMaker::Make()
{
//	cout<<"I'm in StUpsilonTreeMaker::Make()"<<endl;
  if (!GetDataSet("MuDst")) {
    LOG_WARN << "No MuDst" << endm;
    return kStWarn;
  }
//  cout<<"StUpsilonTreeMaker: Make 0"<<endl;
// cout<<"===============NEW EVENT(from data) ; StUpsilonTreeMaker: Make============="<<endl;
//  L2Result* upsL2Result = (L2Result*)&(StMuDst::event()->L2Result().GetArray()[L2RESULTS_2007_OFFSET_UPS]);
//if(!(StMuDst::event()->L2Result())) {cout<<"---wrong !!  StMuDst::event()->L2Result()"<<endl; return kStWarn;  }

/*  -- hdliu
  L2Result* upsL2Result = (L2Result*)&(StMuDst::event()->L2Result().GetArray()[L2RESULTS_2008_OFFSET_UPS]);
  float L2_E1 = upsL2Result->energyOfClusterL0;
  float L2_E2 = upsL2Result->energyOfClusterL2;
  float L2_m = upsL2Result->invMass;

  int TId_L0_l2arr = upsL2Result->triggerTowerL0;
  int TId_L2_l2arr = upsL2Result->triggerTowerL2;
  cout<<"L0 towerID, L2 TowerID (from MuDst):"<<" "<<TId_L0_l2arr<<" "<<TId_L2_l2arr<<endl;
*/


  StEmcPosition emcPos;
  // Obtain candidate softIds
  for (vector<pair<int,int> >::iterator cIt=ups.mL2Candidates.begin();
       cIt!=ups.mL2Candidates.end(); ++cIt) {
    int TId_L0 = (*cIt).first;
    int TId_L2 = (*cIt).second;
 //   cout<<"ups L2 Pairs: "<<TId_L0<<", "<<TId_L2<<endl;
    SurroundingTowers L0SurrTow;
    L0SurrTow.setTower(0,TId_L0);
    SurroundingTowers L2SurrTow;
    L2SurrTow.setTower(0,TId_L2);
    size_t counter(0);
    // obtain the 8 towers surrounding the 
    // candidate softId's (3x3 area)
    //
    for (int deta1 = -1; deta1 <= 1; ++deta1)  {
	for (int dphi1 = -1; dphi1 <= 1; ++dphi1) {
	  if (deta1==dphi1 && deta1==0) continue;
	  counter++;
	  int id11 = emcPos.getNextTowerId(TId_L0, deta1, dphi1);
	  int id22 = emcPos.getNextTowerId(TId_L2, deta1, dphi1);
	  if (id11!=0) {
	    L0SurrTow.setTower(counter,id11);
	  }
	  if (id22!=0) {
	    L2SurrTow.setTower(counter,id22);
	  }	    
	    
	} //dphi loop
    } //deta loop

//    cout << "Surr Towers for pair: " << endl;
//    for (int is=0; is<9; ++is) {
//      cout << L0SurrTow.tower(is) << " ";
//    }
//    cout << endl;

//    for (int is2=0; is2<9; ++is2) {
//      cout << L2SurrTow.tower(is2) << " ";
//    }
//    cout << endl;

    // push each of the 3x3 tower Ids
    // associated with a trigger tower
    // into a vector containing the soft Ids
    // which will be used for selecting tracks
    mVecSoftIdTrig.push_back(L0SurrTow);
    mVecSoftIdTrig.push_back(L2SurrTow);

  } // candidate pair loop
//  cout<<"StUpsilonTreeMaker: Make 1"<<endl;

/*
  cout << "Vector mVecSoftId0Trig = ( " ;
for (vector<int>::iterator Iter1 = mVecSoftId0Trig.begin( ) ; Iter1 != mVecSoftId0Trig.end( ) ; Iter1++ )
   cout << *Iter1 << " ";
   cout << ")" << endl;
  
  cout << "Vector mVecSoftId2Trig = ( " ;
for (vector<int>::iterator Iter2 = mVecSoftId2Trig.begin( ) ; Iter2 != mVecSoftId2Trig.end( ) ; Iter2++ )
   cout << *Iter2 << " ";
   cout << ")" << endl;  
*/

//cout << "StUpsilonPPMuDstMaker::Make -- Upsilon.fNumberOfL2Candiates = "<< ups.fNumberOfL2Candidates << endl;


    //cout << "StUpsilonPPMuDstMaker::Make -- Upsilon.fNumberOfL2Candiates = "
        // << ups.fNumberOfL2Candidates << endl;
    

/*
    pair<int,int> anL2Pairnew;
    size_t candCounternew=0;
    for (vector<pair<int,int> >::iterator it1 = mTrigtowers.begin();
         it1 != mTrigtowers.end();
         ++it1) {
        anL2Pairnew = *it1;
         
        cout << "Pair " << candCounternew++ << " : " << anL2Pairnew.first << ", " << anL2Pairnew.second << endl;
    }

*/

   if (!accept(StMuDst::event())) return kStOk;

  mEvent->setEventInfo(StMuDst::event());
  fill_event_histograms();
  makeEmcHits();
  makeSmdClusters();
  makeCtbHits();
  makeTracks();
  makeEmcClusters();
  makeParticles();
	fill_histograms();
    project_to_tower_fill_his();
	fill_ntuples_by_TreeReader();
	scan_tracks();
	cout<<" to clear_tree_array()"<<endl;
	clear_tree_array();
	cout<<" to fill_tree_par()"<<endl;
	fill_tree_par();
	fill_his_tower_eff();
	cout<<" to myTree->Fill()"<<endl;
	myTree->Fill();
	clear_vector();
	cout<<" done ------ StUpsilonTreeMaker"<<endl;

//  mTree->Fill();
  //	cout<<"I'm leaving StUpsilonTreeMaker::Make()"<<endl;


  return kStOk;
}


int StUpsilonTreeMaker::Finish()
{
	cout<<"  Finish!!!    StUpsilonTreeMaker::Finish()"<<endl;
  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);
   // cout<<"Crate:Sequence:Status::"<<" "<<crate<<" "<<sequence<<" "<<status<<endl;
    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;
// cout<<"11 hits.GetEntriesFast()="<<hits.GetEntriesFast()<<endl;
   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();
    TVector3 momentum = hit->position() - mEvent->vertex();
//	cout<<"momentum.x()="<<momentum.x()<<"  momentum.y()="<<momentum.y()<<"  momentum.z()="<<momentum.z()<<endl;
    momentum.SetMag(energy);
    hit->setMomentum(momentum);
//	cout<<"x="<<x<<"  y="<<y<<"  z="<<z<<endl;
//	cout<<"hit->position().x()="<<hit->position().x()<<"  hit->position().y()="<<hit->position().y()<<"  hit->position().z()="<<hit->position().z()<<endl;
//	cout<<"mEvent->vertex().x()="<<mEvent->vertex().x()<<"  mEvent->vertex().y()="<<mEvent->vertex().y()<<"  mEvent->vertex().z()="<<mEvent->vertex().z()<<endl;
//	cout<<"momentum.x()="<<momentum.x()<<"  momentum.y()="<<momentum.y()<<"  momentum.z()="<<momentum.z()<<endl;
//	cout<<"energy="<<energy<<endl;
  }
//cout<<"22 hits.GetEntriesFast()="<<hits.GetEntriesFast()<<endl;
//cout<<"hits[0].position().x()="<<hits[0].position().x()<<endl;
}

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:
int hdliu_jump_flag=0;
    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();
	position*=energy; 
    // 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);
					hdliu_jump_flag=1;
					  break;  //modified by hdliu
				//	  ++i;
				//	  goto nextCluster;
					}
					energy += hit2->energy();
					TVector3 position2 = hit2->position();
					position += position2 * hit2->energy();
					cluster->hits().Add(hit2);
					used.set(id2);
			  }
			  if(hdliu_jump_flag==1) break;
			}
    if(hdliu_jump_flag==1) continue;
    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()) {
//    double pMagnitude = muTrack->momentum().mag();

    if ( accept(muTrack) && acceptSoftId(muTrack) ) {
      Track* track = new (tracks[tracks.GetEntriesFast()]) Track(muTrack);
      transferMuTrackToTrackMultimap(muTrack,track);
      track->setEvent(mEvent);
      getFirstEmcHit(track);
      getLastEmcHit(track);
      getEmcCluster(track);
      getSmdEtaCluster(track);
      getSmdPhiCluster(track);
      getCtbHit(track);
    }
  }
//  cout<<"The No of track-tower matches to L2 TowerId:"<<mStMuTrackTowerMap.size()<<endl;
}

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();
 /*
 for (int i = 0; i < mEvent->particles().GetEntriesFast(); ++i) {
      Particle* particle = (Particle*)mEvent->particles().At(i);
      //cout<<"Mass from Particles: "<< particle->m()<<endl;
 }
 */
 // cout<<" The number of tracks in Loop :"<< tracks.GetEntriesFast() <<endl;
  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);
//	cout << "i,j " << i << ", " << j << endl;
	bool l2SimRes= checkForValidTrackPair(track1,track2);
	if (l2SimRes) {
	  // Only fill particles if the tracks match tower Ids given by the L2 simulator
  //        cout<<"Now filling particles"<<endl;
	  size_t indx = particles.GetEntriesFast();
	  new (particles[particles.GetEntriesFast()]) Particle(track1, track2);
	  Particle* part = (Particle*)particles[indx];
	  part->setL2Sim(l2SimRes);
//	  cout << "Particle cosTheta: " << part->cosTheta() << endl;
//	  cout << "Particle invMass : " << part->m() << endl;
//	  if (part->m()>6.0) {
//	    cout << "Found a large INV MASS candidate!!!" << endl;
//	  }
	} //if check for Valid Track Pair
      }// track j loop
    }// track i loop
  }// if there are more than 2 tracks
}

bool StUpsilonTreeMaker::checkForValidTrackPair(Track* tr1, Track* tr2) {
  // get id's for tr1
  //use the multimap
//  cout << "Check valid track pair for : " << endl;
//  cout << "Track1 id " << tr1->id() << ", mom: " << tr1->momentum().Mag() << ", chg: " << tr1->charge() << endl;
//  cout << "Track2 id " << tr2->id() << ", mom: " << tr2->momentum().Mag() << ", chg: " << tr1->charge() << endl;
  pair<multimap<Track*,int>::iterator, multimap<Track*,int>::iterator> bounds = mTrackTowerMap.equal_range(tr1);
  // Loop over id's of tr1
  for (multimap<Track*,int>::iterator it = bounds.first; it!=bounds.second; ++it) {
    int TowId1 = (*it).second; // returns the id matched to tr1
    pair<multimap<Track*,int>::iterator, multimap<Track*,int>::iterator> bounds2 = mTrackTowerMap.equal_range(tr2);
    // Loop over id's of tr2
    for (multimap<Track*,int>::iterator it2 = bounds2.first; it2!=bounds2.second; ++it2) {
      int TowId2 = (*it2).second; // returns the id matched to tr2      
      
      // Create a pair to test against the candidate pair vector
      pair<int,int> towPair1(TowId1,TowId2);
//      cout<<" checking the combination of "<< towPair1.first << ", " << towPair1.second<<endl;
      // Check if this combination was found in the L2 simulator
      if (find(ups.mL2Candidates.begin(),ups.mL2Candidates.end(),towPair1)!=ups.mL2Candidates.end()) {
	// Found a match to the L2 Simulator!
	// Debugging:
//        cout<<" I now should get pairs"<<endl;
//	cout << "From Multimap " << towPair1.first << ", " << towPair1.second << endl;
	return true;
      }
      // Check swapped combination
      pair<int,int> towPair2(TowId2,TowId1);
      // Check if this combination was found in the L2 simulator
      if (find(ups.mL2Candidates.begin(),ups.mL2Candidates.end(),towPair2)!=ups.mL2Candidates.end()) {
	// Found a match to the L2 Simulator!
	// Debugging:
//	cout << "From Multimap " << towPair2.first << ", " << towPair2.second << endl;
	return true;
      }
    }// loop over id's of tr2
  } // loop over id's of tr1
  return false;
}


// 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;
}


void StUpsilonTreeMaker::transferMuTrackToTrackMultimap(StMuTrack* trk1, Track* trk2)
{
  pair<multimap<StMuTrack*,int>::iterator, multimap<StMuTrack*,int>::iterator> bounds = mStMuTrackTowerMap.equal_range(trk1);
  // Loop over id's of trk1
  for (multimap<StMuTrack*,int>::iterator it = bounds.first; it!=bounds.second; ++it) {
    int TowId = (*it).second;
    mTrackTowerMap.insert(make_pair(trk2,TowId));
    
  }
  return;
}
  
  

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

bool StUpsilonTreeMaker::accept(StMuTrack* track) const
{
  return (track->flag() > 0 &&
	  track->flag() != 701 &&
	  track->flag() != 801 &&
	  track->flag() != 901 &&
//	  track->pt() > 0.2 &&
//	  0.0 < track->p().mag()  &&  track->p().mag() < 100 &&  
//	  track->nHitsFit() > 20 &&
//	  track->nHitsFit() / (float)track->nHitsPoss() > 0.52 &&
	  3.0 * keV < track->dEdx() && track->dEdx() < 6.0 * keV);
}

bool StUpsilonTreeMaker::accept(Track* track) const
{
  return track->emcCluster();
}


bool StUpsilonTreeMaker::acceptSoftId(StMuTrack *track)
{

   StEmcPosition *pos = new StEmcPosition();
 
  float bFld1 = StMuDst::event()->magneticField();
  float bFld = bFld1/10.;//in Tesla
  //cout <<" magneticField = "<<bFld1<<" "<<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();
//  double delta_eta=3.82339e-03*mEvent->vertex().z();
//	 eta=eta-delta_eta;
  //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);

/*
  //only see tower, not cluster
//  double track_p=track->p().mag();
//  if(track_p<2.5) return false;
	for (vector<pair<int,int> >::iterator cIt=ups.mL2Candidates.begin();cIt!=ups.mL2Candidates.end(); ++cIt) {
			int TId_L0 = (*cIt).first;
			int TId_L2 = (*cIt).second;
			if(id1==TId_L2){mStMuTrackTowerMap.insert(make_pair(track,id1)); return true;}
			if(id1==TId_L0){mStMuTrackTowerMap.insert(make_pair(track,id1)); return true;}
//			if(id1==TId_L0&&track_p>4){mStMuTrackTowerMap.insert(make_pair(track,id1)); return true;}
	   }
	   return false;
  //end of only see tower, not cluster
// my own match
			for(int i=0;i<int(ups.myevent.all_tower_list.size());i++){
				float tower_eta=ups.myevent.all_tower_list[i].eta;
				float tower_phi=ups.myevent.all_tower_list[i].phi;
				int tower_softid=ups.myevent.all_tower_list[i].softid;
				float tower_E=ups.myevent.all_tower_list[i].E;
				if(abs(tower_eta-eta)<0.025&&abs(tower_phi-phi)<0.02625){
					id1=tower_softid;
					if(tower_E>3&&track->p().perp()>3)check_matched_entries->Fill(eta,phi);
					break;
				}
			}
// end of my own match
*/


  bool returnvalue = false;
  for(vector<SurroundingTowers>::iterator It1 = mVecSoftIdTrig.begin(); It1 != mVecSoftIdTrig.end(); ++It1)
    { 
      SurroundingTowers surrTowers = *It1;
      for(size_t i = 0; i<9; i++) { 
	int towerId = surrTowers.tower(i);
 
	
	//cout<< "Getting TowerId from surrounding class:"<<towerId<<endl;
	if (id1==towerId)
	  { 
	    int Htow = surrTowers.highTower();
            mStMuTrackTowerMap.insert(make_pair(track,Htow));
            returnvalue = true;
	    if (0&&track->momentum().mag()>2) {
	      cout << "Found track with p>2.0 GeV, and matching a Trigger tower, dumping info:" << endl;
	      cout<<" Track: " <<track->id() << " is paired with the High Tower:"<<Htow<<endl;	    
	      cout << "Flag : " << track->flag() << endl;
	      cout << "Pt   : " << track->pt() << " should be >0.2 " << endl; 
	      cout << "Mom  : " << track->p().mag() << endl;
	      cout << "HitFit:" << track->nHitsFit() << " should be > 20 " << endl;
	      cout << "F/Pos :" << track->nHitsFit() / (float)track->nHitsPoss() << " should be > 0.52 " << endl;
	      cout << "dE/dx :" << track->dEdx() << " should be > " << 3.0 * keV << endl;
	      
	      
	      cout << "Checking extrapolation " << endl;
	      cout << "Track eta, phi : " << track->momentum().pseudoRapidity() << ", " << track->momentum().phi() << endl;
	      cout << "Tower eta, phi : " << eta << ", " << phi << endl;
	      cout << "THE extrapolated TOWER ID ===="<<id1<<endl;

	    }          

	  }            
	
      } // closing 3x3 loop
    }// loop over surrounding tower vector
//  if (!returnvalue) cout << "Track is NOT paired with L2 towers" << endl;
  return returnvalue; 
}

void StUpsilonTreeMaker::fill_event_histograms(){
//	cout<<"in StUpsilonTreeMaker::fill_event_histograms()"<<endl;
	float vz=mEvent->vertex().z();
	float refmult=	mEvent->refMult();
	mRefMult->Fill(refmult);
	mZvertex->Fill(vz);
	RefMult_vs_vz->Fill(vz,refmult);
	if(vz>-30&&vz<30)good_RefMult->Fill(refmult);
}


void StUpsilonTreeMaker::fill_histograms(){
//	cout<<"in StUpsilonTreeMaker::fill_histograms()"<<endl;
//	cout<<"-- hdliu I'm in fill_histograms()  1"<<endl;
  TIter nextTrack(StMuDst::primaryTracks());
//	cout<<"-- hdliu I'm in fill_histograms()  2"<<endl;

  for (;StMuTrack* muTrack = (StMuTrack*)nextTrack();) {
//	cout<<"-- hdliu I'm in fill_histograms() 3"<<endl;
	 if(!(muTrack->flag() > 0)) continue;
	 double vz=mEvent->vertex().z();                     
     int track_nHitsFit = muTrack->nHitsFit();          
     double track_dedx = muTrack->dEdx();          
     double track_eta = muTrack->eta();   
	 int track_charge=muTrack->charge();   
	 double track_p=muTrack->p().mag();   
	 track_p*=track_charge;
//cout<<"dedx="<<track_dedx<<"   vz="<<vz<<"   track_eta="<<track_eta<<"  track_nHitsFit="<<track_nHitsFit<<endl;
	 if(!(vz>-30&&vz<30)) continue;

	his_eta_vs_nHitFit_all_primary->Fill(track_nHitsFit,track_eta);
	 if(acceptSoftId(muTrack))  his_eta_vs_nHitFit_track_with_emc->Fill(track_nHitsFit,track_eta);
	 
	 if(!(track_eta>-1&&track_eta<1)) continue;
	 if(!(track_nHitsFit>=23)) continue;
	 if(!(track_dedx>1.e-8)) continue;
	 track_dedx*=1.e6;

	 his_dedx_vs_p_all_primary->Fill(track_p,track_dedx);
	 if(acceptSoftId(muTrack)) his_dedx_vs_p_track_with_emc->Fill(track_p,track_dedx);
//	 cout<<"dedx="<<track_dedx<<endl;
  }
}

void StUpsilonTreeMaker::fill_ntuples_by_TreeReader(){
//	cout<<"in StUpsilonTreeMaker::fill_ntuples_by_TreeReader()"<<endl;
    // L2 result
 //   L2Result* upsL2Result = (L2Result*)&mEvent->L2Result()[L2RESULTS_2007_OFFSET_UPS];
    L2Result* upsL2Result = (L2Result*)&mEvent->L2Result()[L2RESULTS_2008_OFFSET_UPS];
//cout<<"read good--L2Result* upsL2Result "<<endl; 
	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);
 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);
 }

}

//void StUpsilonTreeMaker::scan_towers(){
////	cout<<"ups.all_tower_list.size()="<<ups.myevent.all_tower_list.size()<<endl;
//	for(unsigned i=0;i<ups.myevent.all_tower_list.size();i++){
//		struct_BemcTower tower=ups.myevent.all_tower_list[i];
////		 myevent->all_tower_list.push_back(tower);
//		cout<<"i="<<i<<"  softid="<<tower.softid<<"  x"<<tower.x<<"  y="<<tower.y<<"z=  "<<tower.z<<"  E="<<tower.E<<endl;
//	}
////	myevent.all_tower_list=ups.myevent.all_tower_list;
//}

void StUpsilonTreeMaker::scan_tracks(){
  TIter nextTrack(StMuDst::primaryTracks());
//  ups.myevent.vz=mEvent->vertex().z();        
  for (;StMuTrack* muTrack = (StMuTrack*)nextTrack();) {
	  struct_matched_track  matched_track;
	 if(!(muTrack->flag() > 0)) continue;
	 if(track_project_to_tower(muTrack)){ 
		 matched_track.nHitsFit = muTrack->nHitsFit();          
		 matched_track.dedx = muTrack->dEdx();
		 matched_track.nSigE = muTrack->nSigmaElectron();
		 matched_track.nSigPi = muTrack->nSigmaPion();
		 matched_track.nSigP = muTrack->nSigmaProton();
		 matched_track.nSigK = muTrack->nSigmaKaon();
		 matched_track.eta = muTrack->eta();   
		 matched_track.phi = muTrack->phi();   
		 matched_track.charge=muTrack->charge();   
		 matched_track.p=muTrack->p();
		 matched_track.project_eta = project_eta;
		 matched_track.project_phi = project_phi;
		 for(int i=0;i<3;i++)for(int j=0;j<3;j++)matched_track.project_softid[i][j] = project_softid[i][j];
		 ups.myevent.matched_track_list.push_back(matched_track);
	 }
  }
}

bool StUpsilonTreeMaker::track_project_to_tower(StMuTrack* track){
	StEmcPosition *pos = new StEmcPosition();
	float bFld = (StMuDst::event()->magneticField())/10.;//in Tesla
	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();
	int project_softid_center;
	StEmcGeom* emcgeom = StEmcGeom::getEmcGeom("bemc");
	emcgeom->getId(phi,eta,project_softid_center);
	if(!(project_softid_center>0&&project_softid_center<=4800)) return false;
    for (int deta = -1; deta <= 1; ++deta){
		for (int dphi = -1; dphi <= 1; ++dphi) {
			if (dphi==0 && deta==0) continue;
			project_softid[deta+1][dphi+1]=pos->getNextTowerId(project_softid_center, deta, dphi);
		}
	}
	project_softid[1][1]=project_softid_center;
	project_eta=eta;
	project_phi=phi;

	return true;
}

void StUpsilonTreeMaker::clear_vector(){
	ups.myevent.matched_track_list.clear();
	ups.myevent.all_tower_list.clear();
}

void StUpsilonTreeMaker::clear_tree_array(){
	m_Run=0;
	m_evtID=0;
	m_refmult=0;
	m_vz=0;
	m_mag_field=0;
	m_vpdVertex=0;
	m_zdcVertex=0;
	m_bbcVertex=0;
	m_NPrimaryTracks=0;
	m_NTowers=0;

	m_refmult_EFTPC=0;
	m_refmult_WFTPC=0;

	for(int i=0;i<5000;i++){
		m_px[i]=0;
		m_py[i]=0;
		m_pz[i]=0;
		m_project_distance[i]=0;
		m_project_cluster_E[i]=0;
		m_project_distance_top3[i]=0;
		m_project_cluster_E_top3[i]=0;
		m_project_tower_E[i]=0;
		m_project_softid_center[i]=0;
		
		m_pt[i]=0;
		m_p[i]=0;
		m_eta[i]=0;
		m_phi[i]=0;
		m_project_eta[i]=0;
		m_project_phi[i]=0;
		m_dedx[i]=0;
		m_q[i]=0;
		m_nFitPoints[i]=0;
		for(int m=0;m<3;m++)for(int n=0;n<3;n++)m_project_softid[i][m][n]=0;

		m_adc[i]=0;
		m_tower_eta[i]=0;
		m_tower_phi[i]=0;
		m_tower_E[i]=0;
		m_softid[i]=0;
		m_NNeighbors[i]=0;
		for(int m=0;m<3;m++)m_neighbor[i][m]=0;
	}

}

void StUpsilonTreeMaker::definetree(){
	cout<<"define tree"<<endl;
	myTree->Branch("m_Run",&m_Run,"m_Run/I");
	myTree->Branch("m_evtID",&m_evtID,"m_evtID/I");
	myTree->Branch("m_refmult",&m_refmult,"m_refmult/I");
	myTree->Branch("m_mag_field",&m_mag_field,"m_mag_field/F");
	myTree->Branch("m_vz",&m_vz,"m_vz/F");
	myTree->Branch("m_vpdVertex",&m_vpdVertex,"m_vpdVertex/F");
	myTree->Branch("m_zdcVertex",&m_zdcVertex,"m_zdcVertex/F");
	myTree->Branch("m_bbcVertex",&m_bbcVertex,"m_bbcVertex/F");

	myTree->Branch("m_NPrimaryTracks",&m_NPrimaryTracks,"m_NPrimaryTracks/I");
	myTree->Branch("m_q",m_q,"m_q[5000]/I");
	myTree->Branch("m_nFitPoints",m_nFitPoints,"m_nFitPoints[5000]/I");
	myTree->Branch("m_pt",m_pt,"m_pt[5000]/F");
//	myTree->Branch("m_p",m_p,"m_p[5000]/F");
	myTree->Branch("m_eta",m_eta,"m_eta[5000]/F");
	myTree->Branch("m_phi",m_phi,"m_phi[5000]/F");
	myTree->Branch("m_dedx",m_dedx,"m_dedx[5000]/F");
	myTree->Branch("m_nSigE",m_nSigE,"m_nSigE[5000]/F");
	myTree->Branch("m_nSigPi",m_nSigPi,"m_nSigPi[5000]/F");
	myTree->Branch("m_nSigP",m_nSigP,"m_nSigP[5000]/F");
	myTree->Branch("m_nSigK",m_nSigK,"m_nSigK[5000]/F");
//	myTree->Branch("m_project_eta",m_project_eta,"m_project_eta[5000]/F");
//	myTree->Branch("m_project_phi",m_project_phi,"m_project_phi[5000]/F");
//	myTree->Branch("m_project_softid",m_project_softid,"m_project_softid[5000][3][3]/I");

	myTree->Branch("m_NTowers",&m_NTowers,"m_NTowers/I");
//	myTree->Branch("m_adc",m_adc,"m_adc[5000]/F");
//	myTree->Branch("m_tower_eta",m_tower_eta,"m_tower_eta[5000]/F");
//	myTree->Branch("m_tower_phi",m_tower_phi,"m_tower_phi[5000]/F");
//	myTree->Branch("m_tower_E",m_tower_E,"m_tower_E[5000]/F");
//	myTree->Branch("m_softid",m_softid,"m_softid[5000]/I");
//	myTree->Branch("m_NNeighbors",m_NNeighbors,"m_NNeighbors[5000]/I");
//	myTree->Branch("m_neighbor",m_neighbor,"m_neighbor[5000][8]/I");

	myTree->Branch("m_refmult_EFTPC",&m_refmult_EFTPC,"m_refmult_EFTPC/I");
	myTree->Branch("m_refmult_WFTPC",&m_refmult_WFTPC,"m_refmult_WFTPC/I");
	myTree->Branch("m_project_softid_center",m_project_softid_center,"m_project_softid_center[5000]/I");
	myTree->Branch("m_project_tower_E",m_project_tower_E,"m_project_tower_E[5000]/F");
	myTree->Branch("m_project_cluster_E",m_project_cluster_E,"m_project_cluster_E[5000]/F");
	myTree->Branch("m_project_cluster_E_top3",m_project_cluster_E_top3,"m_project_cluster_E_top3[5000]/F");
	myTree->Branch("m_project_distance",m_project_distance,"m_project_distance[5000]/F");
	myTree->Branch("m_project_distance_top3",m_project_distance_top3,"m_project_distance_top3[5000]/F");
	myTree->Branch("m_px",m_px,"m_px[5000]/F");
	myTree->Branch("m_py",m_py,"m_py[5000]/F");
	myTree->Branch("m_pz",m_pz,"m_pz[5000]/F");

}

void StUpsilonTreeMaker::fill_tree_par(){
	m_Run=mEvent->runId();
	m_evtID=mEvent->eventId();
	m_refmult=mEvent->refMult();
	m_vz=StMuDst::event()->primaryVertexPosition().z();
	m_mag_field=StMuDst::event()->magneticField();
	m_vpdVertex=ups.myevent.VpdVertex;
	m_zdcVertex=ups.myevent.ZdcVertex;
	m_bbcVertex=ups.myevent.BbcVertex;
	m_refmult_EFTPC=StMuDst::event()->refMultFtpcEast();
	m_refmult_WFTPC=StMuDst::event()->refMultFtpcWest ();

	m_NPrimaryTracks=ups.myevent.matched_track_list.size();
	for(int i=0;i<m_NPrimaryTracks;i++){
		m_q[i]=ups.myevent.matched_track_list[i].charge;
		m_dedx[i]=ups.myevent.matched_track_list[i].dedx;
		m_nSigE[i]=ups.myevent.matched_track_list[i].nSigE;
		m_nSigPi[i]=ups.myevent.matched_track_list[i].nSigPi;
		m_nSigP[i]=ups.myevent.matched_track_list[i].nSigP;
		m_nSigK[i]=ups.myevent.matched_track_list[i].nSigK;
		m_eta[i]=ups.myevent.matched_track_list[i].eta;
		m_phi[i]=ups.myevent.matched_track_list[i].phi;
		m_p[i]=ups.myevent.matched_track_list[i].p.mag();
		m_pt[i]=ups.myevent.matched_track_list[i].p.perp();
		m_nFitPoints[i]=ups.myevent.matched_track_list[i].nHitsFit;
		m_project_eta[i]=ups.myevent.matched_track_list[i].project_eta;
		m_project_phi[i]=ups.myevent.matched_track_list[i].project_phi;
		for(int k=0;k<3;k++)for(int t=0;t<3;t++)m_project_softid[i][k][t]=ups.myevent.matched_track_list[i].project_softid[k][t];
		for(int k=0;k<int(ups.myevent.all_tower_list.size());k++)if(m_project_softid[i][1][1]==ups.myevent.all_tower_list[k].softid&&ups.myevent.all_tower_list[k].E>3){
			match_eta_match_phi_run->Fill(m_project_eta[i],m_project_phi[i],mEvent->runId());
			match_eta_match_phi_run_check->Fill(ups.myevent.all_tower_list[k].eta,ups.myevent.all_tower_list[k].phi,mEvent->runId());
		}
	}

	m_NTowers=ups.myevent.all_tower_list.size();
	for(int i=0;i<m_NTowers;i++){
		m_adc[i]=ups.myevent.all_tower_list[i].adc;
		m_tower_eta[i]=ups.myevent.all_tower_list[i].eta;
		m_tower_phi[i]=ups.myevent.all_tower_list[i].phi;
		m_tower_E[i]=ups.myevent.all_tower_list[i].E;
		m_softid[i]=ups.myevent.all_tower_list[i].softid;
		m_NNeighbors[i]=ups.myevent.all_tower_list[i].numberOfNeighbors;
		for(int j=0;j<m_NNeighbors[i];j++)	m_neighbor[i][j]=ups.myevent.all_tower_list[i].neighbor[j];
//		if(m_tower_E[i]>3) tower_eta_tower_phi_run->Fill(m_tower_eta[i],m_tower_phi[i],mEvent->runId());
	}

	float tower_E_from_softid[5000];
	for(int i=0;i<5000;i++){
		tower_E_from_softid[i]=0;
	}

	for(int i=0;i<m_NTowers;i++){
		tower_E_from_softid[m_softid[i]]=m_tower_E[i];
	}
	for(int i=0;i<m_NPrimaryTracks;i++){
		m_px[i]=ups.myevent.matched_track_list[i].p.x();
		m_py[i]=ups.myevent.matched_track_list[i].p.y();
		m_pz[i]=ups.myevent.matched_track_list[i].p.z();

		int softid=ups.myevent.matched_track_list[i].project_softid[1][1];
		if(!(softid>=1&&softid<=4800)) continue;
		m_project_softid_center[i]=softid;
		m_project_tower_E[i]=tower_E_from_softid[softid];
		m_project_cluster_E[i]=m_project_tower_E[i];
		m_project_cluster_E_top3[i]=m_project_tower_E[i];

//-------------------
		bool gap_phi=false;
		for(int j=0;j<geog_Nneighbors[softid]-1&&(!gap_phi);j++){
			int nei_id1=geog_neighbor[softid][j];
			for(int t=j+1;t<geog_Nneighbors[softid];t++){
				int nei_id2=geog_neighbor[softid][t];
				if(geog_phi[nei_id1]-geog_phi[nei_id2]>6.2||geog_phi[nei_id1]-geog_phi[nei_id2]<-6.2){
					gap_phi=true;
					break;
				}
			}
		}
//-------------------
		if(!gap_phi){
			float temp_phi=tower_E_from_softid[softid]*geog_phi[softid];
			float temp_eta=tower_E_from_softid[softid]*geog_eta[softid];
			float temp_phi_top3=tower_E_from_softid[softid]*geog_phi[softid];
			float temp_eta_top3=tower_E_from_softid[softid]*geog_eta[softid];
			for(int j=0;j<geog_Nneighbors[softid];j++){
				int nei_id=geog_neighbor[softid][j];
				m_project_cluster_E[i]+=tower_E_from_softid[nei_id];
				temp_phi+=tower_E_from_softid[nei_id]*geog_phi[nei_id];
				temp_eta+=tower_E_from_softid[nei_id]*geog_eta[nei_id];
				if(IsTop2(tower_E_from_softid,softid,j)){
					m_project_cluster_E_top3[i]+=tower_E_from_softid[nei_id];
					temp_phi_top3+=tower_E_from_softid[nei_id]*geog_phi[nei_id];
					temp_eta_top3+=tower_E_from_softid[nei_id]*geog_eta[nei_id];
				}
			}
			if(m_project_cluster_E[i]<0.01) continue;
			if(m_project_cluster_E_top3[i]<0.01) continue;
			temp_phi=temp_phi/m_project_cluster_E[i];
			temp_eta=temp_eta/m_project_cluster_E[i];
			temp_phi_top3=temp_phi_top3/m_project_cluster_E_top3[i];
			temp_eta_top3=temp_eta_top3/m_project_cluster_E_top3[i];
			m_project_distance[i]=sqrt(pow(temp_phi-m_project_phi[i],2)+pow(temp_eta-m_project_eta[i],2));
			m_project_distance_top3[i]=sqrt(pow(temp_phi_top3-m_project_phi[i],2)+pow(temp_eta_top3-m_project_eta[i],2));
			if(m_project_distance[i]>0.3&&m_p[i]>2){
				cout<<"p="<<m_p[i]<<" eta="<<m_eta[i]<<" phi="<<m_phi[i]<<"  r="<< m_project_distance[i]<<"  id="<<softid<<"  tE="<<m_project_tower_E[i]
				<<"  cE="<<m_project_cluster_E[i]<<"  proj_eta="<< m_project_eta[i]<<"  proj_phi="<<m_project_phi[i]<<
				"  temp_eta="<<temp_eta<<"  temp_phi="<<temp_phi<<endl;
				cout<<"  geog_Nneighbors="<<geog_Nneighbors[softid]<<"  nei0_id="<<geog_neighbor[softid][0]<<"  nei1_id="<<geog_neighbor[softid][1]<<"  nei2_id="<<geog_neighbor[softid][2]<<endl;
				cout<<"                                                     "<<"  nei0_E="<<tower_E_from_softid[geog_neighbor[softid][0]]<<"  nei1_E="<<tower_E_from_softid[geog_neighbor[softid][1]]<<"  nei2_E="<<tower_E_from_softid[geog_neighbor[softid][2]]<<endl;
			}
		}
		else{
			float gap=6.2832;
			float temp_phi,temp_phi_top3;
			if(geog_phi[softid]>0){
				temp_phi=tower_E_from_softid[softid]*geog_phi[softid];
				temp_phi_top3=tower_E_from_softid[softid]*geog_phi[softid];
			}
			else {
				temp_phi=tower_E_from_softid[softid]*(gap+geog_phi[softid]); 
				temp_phi_top3=tower_E_from_softid[softid]*(gap+geog_phi[softid]); 
			}
			float temp_eta=tower_E_from_softid[softid]*geog_eta[softid];
			float temp_eta_top3=tower_E_from_softid[softid]*geog_eta[softid];
			for(int j=0;j<geog_Nneighbors[softid];j++){
				int nei_id=geog_neighbor[softid][j];
				m_project_cluster_E[i]+=tower_E_from_softid[nei_id];
				if(geog_phi[nei_id]>0)temp_phi+=tower_E_from_softid[nei_id]*geog_phi[nei_id];
				else temp_phi+=tower_E_from_softid[nei_id]*(gap+geog_phi[nei_id]);
				temp_eta+=tower_E_from_softid[nei_id]*geog_eta[nei_id];
				if(IsTop2(tower_E_from_softid,softid,j)){
					m_project_cluster_E_top3[i]+=tower_E_from_softid[nei_id];
					if(geog_phi[nei_id]>0)temp_phi_top3+=tower_E_from_softid[nei_id]*geog_phi[nei_id];
					else temp_phi_top3+=tower_E_from_softid[nei_id]*(gap+geog_phi[nei_id]);
					temp_eta_top3+=tower_E_from_softid[nei_id]*geog_eta[nei_id];
				}
			}
			if(m_project_cluster_E[i]<0.01) continue;
			if(m_project_cluster_E_top3[i]<0.01) continue;
			temp_phi=temp_phi/m_project_cluster_E[i];
			temp_eta=temp_eta/m_project_cluster_E[i];
			temp_phi_top3=temp_phi_top3/m_project_cluster_E_top3[i];
			temp_eta_top3=temp_eta_top3/m_project_cluster_E_top3[i];
			if(m_project_phi[i]>0){
				m_project_distance[i]=sqrt(pow(temp_phi-m_project_phi[i],2)+pow(temp_eta-m_project_eta[i],2));
				m_project_distance_top3[i]=sqrt(pow(temp_phi_top3-m_project_phi[i],2)+pow(temp_eta_top3-m_project_eta[i],2));
			}
			else{
				m_project_distance[i]=sqrt(pow(temp_phi-(m_project_phi[i]+gap),2)+pow(temp_eta-m_project_eta[i],2));
				m_project_distance_top3[i]=sqrt(pow(temp_phi_top3-(m_project_phi[i]+gap),2)+pow(temp_eta_top3-m_project_eta[i],2));
			}
			if(m_project_distance[i]>0.3&&m_p[i]>2){
				cout<<"p="<<m_p[i]<<" eta="<<m_eta[i]<<" phi="<<m_phi[i]<<"  r="<< m_project_distance[i]<<"  id="<<softid<<"  tE="<<m_project_tower_E[i]
				<<"  cE="<<m_project_cluster_E[i]<<"  proj_eta="<< m_project_eta[i]<<"  proj_phi="<<m_project_phi[i]<<
				"  temp_eta="<<temp_eta<<"  temp_phi="<<temp_phi<<endl;
				cout<<"  geog_Nneighbors="<<geog_Nneighbors[softid]<<"  nei0_id="<<geog_neighbor[softid][0]<<"  nei1_id="<<geog_neighbor[softid][1]<<"  nei2_id="<<geog_neighbor[softid][2]<<endl;
				cout<<"                                                     "<<"  nei0_E="<<tower_E_from_softid[geog_neighbor[softid][0]]<<"  nei1_E="<<tower_E_from_softid[geog_neighbor[softid][1]]<<"  nei2_E="<<tower_E_from_softid[geog_neighbor[softid][2]]<<endl;
			}
		}
		
	}

}

void StUpsilonTreeMaker::project_to_tower_fill_his(){
	  TIter nextTrack(StMuDst::primaryTracks());
	  for (;StMuTrack* track = (StMuTrack*)nextTrack();) {
			StEmcPosition *pos = new StEmcPosition();
			float bFld = (StMuDst::event()->magneticField())/10.;//in Tesla
			StThreeVectorD position1, momentum1;
			Bool_t ok = pos->trackOnEmc(&position1, &momentum1, track, bFld); // bFld in Tesla        
			if(!ok)  continue;
			double eta=position1.pseudoRapidity();
//			double phi=position1.phi();

		  double track_pt=track->p().perp();
		  double track_eta=track->eta();   
		  double vz=mEvent->vertex().z();
//		 if(track_pt>3) track_eta_track_phi_run->Fill(eta,phi,mEvent->runId());
		  if(track_pt>2.5)delta_eta_vs_vz_highpt->Fill(vz,eta-track_eta);
		  //cout<<" eta position "<< eta <<" phi position "<< phi <<endl;
//			double delta_eta=3.82339e-03*mEvent->vertex().z();
//			eta=eta-delta_eta;
		   if(track_pt>2.5)delta_eta_vs_vz_highpt_after_correction->Fill(vz,eta-track_eta);
	  }

}

void StUpsilonTreeMaker::fill_his_tower_eff(){
	float ENERGY_CUT=3;

//	for(int i=0;i<int(ups.myevent.all_tower_list.size());i++){
//		if(ups.myevent.all_tower_list[i].E>ENERGY_CUT) tower_eta_tower_phi_run->Fill(ups.myevent.all_tower_list[i].eta,ups.myevent.all_tower_list[i].phi,mEvent->runId());
//	}

	  TIter nextTrack(StMuDst::primaryTracks());
	  for (;StMuTrack* track = (StMuTrack*)nextTrack();) {
			if(!(track->flag() > 0)) continue;
			if(!(track->nHitsFit() > 23)) continue;
//			if(!(track->p().perp()>ENERGY_CUT)) continue;
			if(!(track->p().mag()>ENERGY_CUT)) continue;
		  
			double track_pt=track->p().perp();
			double track_eta=track->eta();   
			double track_phi=track->phi();   
			track_eta_track_phi_run->Fill(track_eta,track_phi,mEvent->runId());

			StEmcPosition *pos = new StEmcPosition();
			float bFld = (StMuDst::event()->magneticField())/10.;//in Tesla
			StThreeVectorD position1, momentum1;
			Bool_t ok = pos->trackOnEmc(&position1, &momentum1, track, bFld); // bFld in Tesla        
			if(!ok)  continue;
			double eta=position1.pseudoRapidity();
			double phi=position1.phi();

			int project_softid_center;
			StEmcGeom* emcgeom = StEmcGeom::getEmcGeom("bemc");
			emcgeom->getId(phi,eta,project_softid_center);
			if(!(project_softid_center>0&&project_softid_center<=4800)) continue;
			track_proj_eta_track_phi_run->Fill(eta,phi,mEvent->runId());
			track_proj_softid->Fill(project_softid_center);

			for(int i=0;i<int(ups.myevent.all_tower_list.size());i++){
				float tower_eta=ups.myevent.all_tower_list[i].eta;
				float tower_phi=ups.myevent.all_tower_list[i].phi;
				int tower_softid=ups.myevent.all_tower_list[i].softid;
				float tower_E=ups.myevent.all_tower_list[i].E;
//				if(abs(tower_eta-eta)<0.025&&abs(tower_phi-phi)<0.02625&&tower_E>ENERGY_CUT){
				if(tower_E>ENERGY_CUT){
					tower_eta_tower_phi_run->Fill(tower_eta,tower_phi,mEvent->runId());
					tower_fired_softid->Fill(tower_softid);
					if(tower_softid==project_softid_center){
						my_own_match_eta_match_phi_run->Fill(tower_eta,tower_phi,mEvent->runId());
						matched_softid->Fill(tower_softid);
//						my_own_match_eta_match_phi_run2->Fill(tower_eta,tower_phi,mEvent->runId());
//						match_id_correlation_geomID_vs_towerID->Fill(tower_softid,project_softid_center);
					}
				}  //tower.E>3
			}// loop towers


	  }// MuTracks loop


}
bool StUpsilonTreeMaker::IsTop2(float tower_E_from_softid[5000],int softid, int index){
	int top=1;
	for(int j=0;j<geog_Nneighbors[softid];j++){
		if(j==index) continue;
		float E=tower_E_from_softid[geog_neighbor[softid][index]];
		int nei_id=geog_neighbor[softid][j];
		if(E<tower_E_from_softid[nei_id]) top++;
		if(top>2) break;
	}
	if(top==1||top==2) return true;
	return false;
}



