//
// Pibero Djawotho <pibero@indiana.edu>
// Indiana University
// Feb 10, 2008
//

// ROOT
#include "TChain.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TClonesArray.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH1.h"
#include "TLine.h"

// STAR
#include "SystemOfUnits.h"
#include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
#include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
#include "StEEmcUtil/EEmcSmdMap/EEmcSmdMap.h"
#include "StarCallf77.h"

// Local
#include "StGammaEvent.h"
#include "StGammaCandidate.h"
#include "StEndcapSmdCluster.h"
#include "StEndcapEmcCluster.h"
#include "StEndcapEmcPoint.h"
#include "StEtaFinder.h"

// Declaration of CERNLIB routine (mathlib, H301) assndx to be used for matching clusters
#define assndx F77_NAME(assndx,ASSNDX)

extern "C" void type_of_call assndx(int&, float*, int&, int&, int&, int*, float&, int*, int&);

struct StripCompare {
  bool operator()(StGammaStrip* strip1, StGammaStrip* strip2) const
  {
    return strip1->energy > strip2->energy;
  }
};

struct TowerCompare {
  bool operator()(StGammaTower* tower1, StGammaTower* tower2) const
  {
    return tower1->pt() > tower2->pt();
  }
};

ClassImp(StEtaFinder);

void StEtaFinder::clear()
{
  memset(mTowers, 0, sizeof(mTowers));
//   memset(mPreshower1, 0, sizeof(mPreshower1));
//   memset(mPreshower2, 0, sizeof(mPreshower2));
//   memset(mPostshower, 0, sizeof(mPostshower));
  memset(mStrips, 0, sizeof(mStrips));

  for (int sector = 0; sector < 12; ++sector)
    for (int plane = 0; plane < 2; ++plane)
      mClusters[sector][plane]->Clear();

  mTracks.clear();
  mEmcClusters->Clear();
  mMatchedClusters.clear();
  mPoints->Clear();
  mCandidates.clear();
  mGammaPoints.clear();
  mTowerPoints.clear();
}

void StEtaFinder::init(const char* filelist, const char* ofile)
{
  // Load input files
  mChain = new TChain("gammas");

  ifstream in(filelist);
  string line;
  while (getline(in, line)) mChain->AddFile(line.c_str());
  in.close();

  // Event buffer
  mEvent = 0;
  mChain->SetBranchAddress("GammaEvent", &mEvent);

  // Create output file
  mFile = new TFile(ofile, "recreate");
  assert(mFile && !mFile->IsZombie());

  // Create ntuple
  TString varlist = "vx:vy:vz:cosTheta:m:p:pt:eta:phi";
  varlist += ":";
  varlist += "p1:pt1:eta1:phi1";
  varlist += ":";
  varlist += "p2:pt2:eta2:phi2";
  varlist += ":";
  varlist += "dca";
  varlist += ":";
  varlist += "ntrk1:ntrk2";
  varlist += ":";
  varlist += "eid1:eid2";
  varlist += ":";
  varlist += "iso1:iso2";
  varlist += ":";
  varlist += "npoints";
  varlist += ":";
  varlist += "umax1:umax2:usum1:usum2";
  varlist += ":";
  varlist += "vmax1:vmax2:vsum1:vsum2";

  mTuple = new TNtuple("etas", "etas", varlist);

  // Create TClonesArray for SMD clusters
  for (int sector = 0; sector < 12; ++sector)
    for (int plane = 0; plane < 2; ++plane)
      mClusters[sector][plane] = new TClonesArray("StEndcapSmdCluster");

  // Create TClonesArray for EMC clusters
  mEmcClusters = new TClonesArray("StEndcapEmcCluster");
  mPoints = new TClonesArray("StEndcapEmcPoint");

  // Create postscript file
  gStyle->SetOptStat(0);
  gStyle->SetPalette(1);
  gStyle->SetPadGridX(false);
  gStyle->SetPadGridY(false);
  mCanvas = new TCanvas;
  mPsFile = ofile;
  mPsFile.ReplaceAll(".root", ".ps");
  mCanvas->Print(mPsFile + "[");

  // Keep tracks of number of points after successive cuts
  mPointCounter = new TH1F("mPointCounter", "mPointCounter", 5, 0, 5);
}

void StEtaFinder::finish()
{
  mPointCounter->LabelsDeflate();

  mFile->Write();
  mFile->Close();

  mCanvas->Print(mPsFile + "]");
}

void StEtaFinder::eventLoop()
{
  // Event loop
  for (int iEntry = 0;; ++iEntry) {
    clear();
    if (!mChain->GetEntry(iEntry)) break;
    if (acceptEvent()) {
//       getCandidates();
      getTowers();
//       getPreshower1();
//       getPreshower2();
//       getPostshower();
      getStrips();
      getTracks();
      findTowerPoints();
      processEvent();
    } // Accept event
  } // Loop over events
}

bool StEtaFinder::acceptEvent() const
{
  // L2gamma
  return (mEvent->isTrigger(117641) ||
	  mEvent->isTrigger(127641) ||
	  mEvent->isTrigger(137641));

  // ----------- DEAD CODE -----------

  // bjp, ejp, dijet
  return (mEvent->isTrigger(137273) ||
	  mEvent->isTrigger(137652) ||
	  mEvent->isTrigger(137551) ||
	  mEvent->isTrigger(137222) ||
	  mEvent->isTrigger(137622));

  // eemc-jp1-mb || eemc-jp0-mb
  return (mEvent->isTrigger(117271) ||
	  mEvent->isTrigger(127271) ||
	  mEvent->isTrigger(137271) ||
	  mEvent->isTrigger(137272) ||
	  mEvent->isTrigger(137273) ||
	  mEvent->isTrigger(117551) ||
	  mEvent->isTrigger(127551) ||
	  mEvent->isTrigger(137551));
  // ---------------------------------
}

void StEtaFinder::processEvent()
{
  cout << "----------------- run = " << mEvent->runNumber() << ", event = " << mEvent->eventNumber() << " ----------------" << endl;
#if 0
  // Find SMD clusters
  findClusters();

  // Match SMD clusters
  matchClusters();

  // Find EMC clusters
  findEmcClusters();

  // Combine EMC+SMD clusters
  makePoints();

  // Match gamma candidates to SMD clusters
  matchCandidates();
#endif

  // Calculate invariant mass
  for (size_t i = 0; i < mTowerPoints.size(); ++i) {
    StTowerPoint& point1 = mTowerPoints[i];
    TVector3 p1 = point1.momentum(mEvent->vertex());

    for (size_t j = i + 1; j < mTowerPoints.size(); ++j) {
      StTowerPoint& point2 = mTowerPoints[j];
      TVector3 p2 = point2.momentum(mEvent->vertex());

      float cosTheta = p1.Unit() * p2.Unit();
      float m = sqrt(2 * p1.Mag() * p2.Mag() * (1 - cosTheta));
      TVector3 p = p1 + p2;
      TVector3 dca = point1.position - point2.position;

      float tuple[] = {
	mEvent->vertex().x(),
	mEvent->vertex().y(),
	mEvent->vertex().z(),
	cosTheta,
	m,
	p.Mag(),
	p.Pt(),
	p.Eta(),
	p.Phi(),
	p1.Mag(),
	p1.Pt(),
	p1.Eta(),
	p1.Phi(),
	p2.Mag(),
	p2.Pt(),
	p2.Eta(),
	p2.Phi(),
	dca.Mag(),
	0,
	0,
	0,
	0,
	0,
	0,
	mTowerPoints.size(),
	0,
	0,
	0,
	0,
	0,
	0,
	0,
	0
      };

      mTuple->Fill(tuple);

      if (0.4 < m && m < 0.6) {
	print(point1);
	print(point2);
      }
    } // Loop over 2nd point
  } // Loop over 1st point

#if 0
  // Calculate invariant mass
  int npoints = mPoints->GetEntriesFast();
  cout << "Npoints = " << npoints << endl;

  for (size_t i = 0; i < mGammaPoints.size(); ++i) {
    StGammaPoint& point1 = mGammaPoints[i];
    TVector3 p1 = point1.momentum(mEvent->vertex());

    for (size_t j = i + 1; j < mGammaPoints.size(); ++j) {
      StGammaPoint& point2 = mGammaPoints[j];
      TVector3 p2 = point2.momentum(mEvent->vertex());

      float cosTheta = p1.Unit() * p2.Unit();
      float m = sqrt(2 * p1.Mag() * p2.Mag() * (1 - cosTheta));
      TVector3 p = p1 + p2;
      TVector3 dca = point1.position - point2.position;

      float tuple[] = {
	mEvent->vertex().x(),
	mEvent->vertex().y(),
	mEvent->vertex().z(),
	cosTheta,
	m,
	p.Mag(),
	p.Pt(),
	p.Eta(),
	p.Phi(),
	p1.Mag(),
	p1.Pt(),
	p1.Eta(),
	p1.Phi(),
	p2.Mag(),
	p2.Pt(),
	p2.Eta(),
	p2.Phi(),
	dca.Mag(),
	point1.candidate->numberOfTracks(),
	point2.candidate->numberOfTracks(),
	0,
	0,
	0,
	0,
	mGammaPoints.size(),
	point1.smdu->maxStrip()->energy,
	point2.smdu->maxStrip()->energy,
	point1.smdu->energy(),
	point2.smdu->energy(),
	point1.smdv->maxStrip()->energy,
	point2.smdv->maxStrip()->energy,
	point1.smdv->energy(),
	point2.smdv->energy()
      };

      mTuple->Fill(tuple);      
    } // Loop over 2nd point
  } // Loop over 1st point
#endif

#if 0
  for (int i = 0; i < npoints; ++i) {
    bool etaCandidate = false;

    StEndcapEmcPoint* point1 = (StEndcapEmcPoint*)mPoints->UncheckedAt(i);
    TVector3 p1 = point1->momentum(mEvent->vertex());
    bool eid1 = electron(point1);
    float iso1 = 0.5 * (point1->smdCluster(0)->isolation() + point1->smdCluster(1)->isolation());

    for (int j = i + 1; j < npoints; ++j) {
      StEndcapEmcPoint* point2 = (StEndcapEmcPoint*)mPoints->UncheckedAt(j);
      TVector3 p2 = point2->momentum(mEvent->vertex());
      bool eid2 = electron(point2);
      float iso2 = 0.5 * (point2->smdCluster(0)->isolation() + point2->smdCluster(1)->isolation());

      float cosTheta = p1.Unit() * p2.Unit();
      float invMass = sqrt(2 * p1.Mag() * p2.Mag() * (1 - cosTheta));
      TVector3 p = p1 + p2;
      TVector3 dca = point1->position() - point2->position();

      float tuple[] = {
	mEvent->vertex().x(),
	mEvent->vertex().y(),
	mEvent->vertex().z(),
	cosTheta,
	invMass,
	p.Mag(),
	p.Pt(),
	p.Eta(),
	p.Phi(),
	p1.Mag(),
	p1.Pt(),
	p1.Eta(),
	p1.Phi(),
	p2.Mag(),
	p2.Pt(),
	p2.Eta(),
	p2.Phi(),
	dca.Mag(),
	point1->numberOfTracks(),
	point2->numberOfTracks(),
	eid1,
	eid2,
	iso1,
	iso2,
	npoints,
	point1->smdCluster(0)->maxStrip()->energy,
	point2->smdCluster(0)->maxStrip()->energy,
	point1->smdCluster(0)->energy(),
	point2->smdCluster(0)->energy(),
	point1->smdCluster(1)->maxStrip()->energy,
	point2->smdCluster(1)->maxStrip()->energy,
	point1->smdCluster(1)->energy(),
	point2->smdCluster(1)->energy()
      };

      mTuple->Fill(tuple);

      // Test for eta candidate
      if (p1.Pt() > 4.0 && 0.4 < invMass && invMass < 0.7 && dca.Mag() > 10) etaCandidate = true;
    }

    // Print to postscript file
    if (etaCandidate) {
      // Create canvas
      if (mCanvas) delete mCanvas;
      mCanvas = new TCanvas("c1", "c1", 700, 900);
      mCanvas->Divide(1, 2);

      // 1st photon
      int u1 = point1->smdCluster(0)->maxStrip()->index;
      int v1 = point1->smdCluster(1)->maxStrip()->index;

      TString tU1 = Form("run=%d, event=%d, sector=%d", mEvent->runNumber(), mEvent->eventNumber(), point1->smdCluster(0)->sector());
      TString tV1 = Form("p_{T}=%.1lf GeV, #eta_{det}=%.1lf", (double)point1->momentum(mEvent->vertex()).Pt(), (double)point1->position().Eta());

      TH1F hU1("hU1", ";u-strip", 40, u1 - 20, u1 + 20);
      TH1F hV1("hV1", ";v-strip", 40, v1 - 20, v1 + 20);

      hU1.SetTitle(tU1);
      hV1.SetTitle(tV1);

      for (int id = 0; id < 288; ++id) {
	StGammaStrip* strip = mStrips[point1->smdCluster(0)->sector()][point1->smdCluster(0)->plane()][id];
	if (strip) hU1.Fill(id, strip->energy);
      }

      for (int id = 0; id < 288; ++id) {
	StGammaStrip* strip = mStrips[point1->smdCluster(1)->sector()][point1->smdCluster(1)->plane()][id];
	if (strip) hV1.Fill(id, strip->energy);
      }

      float ymax1 = 1.1 * max(hU1.GetMaximum(), hV1.GetMaximum());
      hU1.SetMaximum(ymax1);
      hV1.SetMaximum(ymax1);

      TLine lU1(u1, 0, u1, ymax1);
      TLine lV1(v1, 0, v1, ymax1);

      lU1.SetLineColor(kRed);
      lV1.SetLineColor(kRed);

      mCanvas->cd(1); hU1.Draw(); lU1.Draw("same");
      mCanvas->cd(2); hV1.Draw(); lV1.Draw("same");

      // Print to postscript file
      mCanvas->Print(mPsFile);
    }
  }
#endif
}

void StEtaFinder::getTowers()
{
  for (int i = 0; i < mEvent->numberOfTowers(); ++i) {
    StGammaTower* tower = mEvent->tower(i);
    if (tower->layer == kEEmcTower)
      mTowers[tower->etabin()][tower->phibin()] = tower;
  }
}

void StEtaFinder::getPreshower1()
{
  for (int i = 0; i < mEvent->numberOfPreshower1(); ++i) {
    StGammaTower* tower = mEvent->preshower1(i);
    if (tower->layer == kEEmcPre1)
      mPreshower1[tower->etabin()][tower->phibin()] = tower;
  }
}

void StEtaFinder::getPreshower2()
{
  for (int i = 0; i < mEvent->numberOfPreshower2(); ++i) {
    StGammaTower* tower = mEvent->preshower2(i);
    if (tower->layer == kEEmcPre2)
      mPreshower2[tower->etabin()][tower->phibin()] = tower;
  }
}

void StEtaFinder::getPostshower()
{
  for (int i = 0; i < mEvent->numberOfTowers(); ++i) {
    StGammaTower* tower = mEvent->tower(i);
    if (tower->layer == kEEmcPost)
      mTowers[tower->etabin()][tower->phibin()] = tower;
  }
}

void StEtaFinder::getStrips()
{
  for (int i = 0; i < mEvent->numberOfStrips(); ++i) {
    StGammaStrip* strip = mEvent->strip(i);
    if (strip->plane == kEEmcSmdu || strip->plane == kEEmcSmdv)
      mStrips[strip->sector][strip->plane][strip->index] = strip;
  }
}

void StEtaFinder::getTracks()
{
  for (int i = 0; i < mEvent->numberOfTracks(); ++i) {
    StGammaTrack* track = mEvent->track(i);
    try {
      static const EEmcGeomSimple& geo = EEmcGeomSimple::Instance();
      TVector3 position = track->positionAtZ(geo.getZSMD());
      int sector, subsector, etabin;
      if (geo.getTower(position, sector, subsector, etabin)) {
	int phibin = sector * 5 + subsector;
	StGammaTower* tower = mTowers[etabin][phibin];
	if (tower) {
	  mTracks.insert(make_pair(tower, track));
	}
      }
    }
    catch (StGammaTrack::Exception& e) {}
  }
}

void StEtaFinder::findClusters()
{
  const int HALF_WIDTH = 3;

  for (int sector = 0; sector < 12; ++sector) {
    for (int plane = 0; plane < 2; ++plane) {
      // Find seed candidates (> 2 MeV)
      vector<StGammaStrip*> seeds;
      for (int i = 0; i < 288; ++i) {
	StGammaStrip* strip = mStrips[sector][plane][i];
	if (strip && strip->energy > 2 * MeV) seeds.push_back(strip);
      } // Loop over strips

      // Sort seeds in descending order of energy
      sort(seeds.begin(), seeds.end(), StripCompare());

      // Build (2*HALF_WIDTH+1)-strip clusters
      for (size_t i = 0; i < seeds.size(); ++i) {
	StGammaStrip* seed = mStrips[sector][plane][seeds[i]->index];
	if (seed) {
#if 0
	  const int RANGE = 10;
	  float energy = 0;
	  float residual = 0;
	  for (int id = (int)seed->index - RANGE; id <= (int)seed->index + RANGE; ++id) {
	    if (0 <= id && id < 288) {
	      StGammaStrip* strip = mStrips[sector][plane][id];
	      if (strip) {
		if ((int)seed->index - HALF_WIDTH <= id && id <= (int)seed->index + HALF_WIDTH)
		  energy += strip->energy;
		else
		  residual += strip->energy;
	      }
	    }
	  }

	  float fraction = residual / (energy + residual);
	  if (fraction > 0.25) continue;
#endif
	  TClonesArray& clusters = *mClusters[sector][plane];
	  int nentries = clusters.GetEntriesFast();
	  StEndcapSmdCluster* cluster = new (clusters[nentries]) StEndcapSmdCluster;
	  for (int dx = -HALF_WIDTH; dx <= HALF_WIDTH; ++dx) {
	    int id = seed->index + dx;
	    if (0 <= id && id < 288) {
	      StGammaStrip* strip = mStrips[sector][plane][id];
	      if (strip) {
		cluster->addStrip(strip);
		mStrips[sector][plane][id] = 0;	// Remove strip from further consideration
	      }
	    }
	  } // Loop over strips

	  // Sort strips in cluster in order of descending energy
	  cluster->sort();

	  // Isolation cut
	  //if (cluster->isolation() > 0.5) clusters.RemoveLast();
	}
      } // Loop over seeds
    } // Loop over planes
  } // Loop over sectors

  // Reload strips from event
  getStrips();
}

void StEtaFinder::matchClusters()
{
  // Hard limit on number of clusters per plane
  const int MAX_CLUSTERS = 20;

  for (int sector = 0; sector < 12; ++sector) {
    // Skip sector if at least one of the plane has no cluster
    if (mClusters[sector][0]->IsEmpty() || mClusters[sector][1]->IsEmpty()) continue;

    // Clear array of energy asymmetries between SMDU and SMDV clusters
    float asym[MAX_CLUSTERS][MAX_CLUSTERS];
    memset(asym, 0, sizeof(asym));

    // Should probably sort clusters in order of descending energy

    // Compute energy asymmetries between SMDU and SMDV clusters
    // Loop over SMDU clusters
    for (int i = 0; i < mClusters[sector][0]->GetEntriesFast() && i < MAX_CLUSTERS; ++i) {
      StEndcapSmdCluster* cluster1 = (StEndcapSmdCluster*)mClusters[sector][0]->UncheckedAt(i);
      float energy1 = cluster1->energy();
      // Loop over SMDV clusters
      for (int j = 0; j < mClusters[sector][1]->GetEntriesFast() && j < MAX_CLUSTERS; ++j) {
	StEndcapSmdCluster* cluster2 = (StEndcapSmdCluster*)mClusters[sector][1]->UncheckedAt(j);
	float energy2 = cluster2->energy();
	asym[i][j] = fabs(energy1 - energy2) / (energy1 + energy2);
      } // Loop over SMDV clusters
    } // Loop over SMDU clusters

    int mode = 1;
    int na = mClusters[sector][0]->GetEntriesFast();
    int ma = mClusters[sector][1]->GetEntriesFast();
    int ida = MAX_CLUSTERS;
    int k[MAX_CLUSTERS];
    float smin;
    int iw[MAX_CLUSTERS][MAX_CLUSTERS];
    int idw = MAX_CLUSTERS;

    // Clear arrays
    memset(k,  0, sizeof(k ));
    memset(iw, 0, sizeof(iw));

    assndx(mode, asym[0], na, ma, ida, k, smin, iw[0], idw);

    for (int i = 0; i < na; ++i) {
      if (k[i]) {
	int j = k[i]-1;
	StEndcapSmdCluster* cluster1 = (StEndcapSmdCluster*)mClusters[sector][0]->At(i);
	StEndcapSmdCluster* cluster2 = (StEndcapSmdCluster*)mClusters[sector][1]->At(j);
	float e1 = cluster1->energy();
	float e2 = cluster2->energy();
	float asym = (e1 - e2) / (e1 + e2);
	if (asym > 0.75) continue;
	TVector3 position = getIntersection(cluster1, cluster2);
	if (position.z() != -999) {
	  StGammaTower* tower = getTower(position);
	  if (tower) {
	    mMatchedClusters.push_back(make_pair(cluster1, cluster2));
	  }
	}
      }
    }

  } // loop over sectors
}

float StEtaFinder::conePt(StGammaTower* tower, float radius) const
{
  float pt = 0;

  for (int etabin = 0; etabin < 12; ++etabin) {
    for (int phibin = 0; phibin < 60; ++phibin) {
      StGammaTower* tower2 = mTowers[etabin][phibin];
      if (tower2) {
	float deta = tower->eta - tower2->eta;
	float dphi = tower->phi - tower2->phi;
	dphi = TVector2::Phi_mpi_pi(dphi);
	float r = hypot(deta, dphi);
	if (r <= radius) pt += tower2->pt();
      }
    }
  }

  return pt;
}

void StEtaFinder::findEmcClusters()
{
  // Find seed towers (pT > 0.5 GeV)
  vector<StGammaTower*> seeds;

  for (int etabin = 0; etabin < 12; ++etabin) {
    for (int phibin = 0; phibin < 60; ++phibin) {
      StGammaTower* tower = mTowers[etabin][phibin];
      if (tower && tower->pt() > 0.5) seeds.push_back(tower);
    } // Loop over phibin
  } // Loop over etabin

  // Sort seeds in order of descending transverse energy
  sort(seeds.begin(), seeds.end(), TowerCompare());

  // Create 3x3 EMC clusters by adding neighboring towers
  for (size_t i = 0; i < seeds.size(); ++i) {
    int etabin = seeds[i]->etabin();
    int phibin = seeds[i]->phibin();
    StGammaTower* seed = mTowers[etabin][phibin];
    if (seed) {
      TClonesArray& clusters = *mEmcClusters;
      int nentries = mEmcClusters->GetEntriesFast();
      StEndcapEmcCluster* cluster = new (clusters[nentries]) StEndcapEmcCluster;
      cluster->addTower(seed);
      mTowers[seed->etabin()][seed->phibin()] = 0; // Remove seed from further consideration
      for (int deta = -1; deta <= 1; ++deta) {
	for (int dphi = -1; dphi <= 1; ++dphi) {
	  if (deta || dphi) {
	    StGammaTower* tower = getNextTower(seed, deta, dphi);
	    if (tower) {
	      cluster->addTower(tower);
	      mTowers[tower->etabin()][tower->phibin()] = 0; // Remove tower from further consideration
	    }
	  }
	} // Loop over phibin
      }	// Loop over etabin
    }
  } // Loop over seeds

  // Reload towers from event
  getTowers();
}

StGammaTower* StEtaFinder::getNextTower(StGammaTower* tower, int deta, int dphi) const
{
  int etabin = tower->etabin() + deta;
  if (etabin < 0 || etabin >= 12) return 0;
  int phibin = tower->phibin() + dphi;
  if (phibin <   0) phibin += 60;
  if (phibin >= 60) phibin -= 60;
  return mTowers[etabin][phibin];
}

void StEtaFinder::makePoints()
{
  for (int i = 0; i < mEmcClusters->GetEntriesFast(); ++i) {
    StEndcapEmcCluster* emcCluster = (StEndcapEmcCluster*)mEmcClusters->UncheckedAt(i);
    float smdEnergy = 0;
    vector<StEndcapEmcPoint*> points;
    for (size_t j = 0; j < mMatchedClusters.size(); ++j) {
      StEndcapSmdCluster* cluster1 = mMatchedClusters[j].first;
      StEndcapSmdCluster* cluster2 = mMatchedClusters[j].second;
      TVector3 position = getIntersection(cluster1, cluster2);
      if (position.z() != -999) {
	StGammaTower* tower = getTower(position);
	if (tower && emcCluster->contains(tower)) {
	  StEndcapEmcPoint* point = new ((*mPoints)[mPoints->GetEntriesFast()]) StEndcapEmcPoint;
	  point->setEmcCluster(emcCluster);
	  point->setSmdCluster(cluster1, 0);
	  point->setSmdCluster(cluster2, 1);
	  point->setPosition(position);
	  float energy = 0.5 * (cluster1->energy() + cluster2->energy());
	  point->setEnergy(energy);
	  smdEnergy += energy;
	  points.push_back(point);

	  // Find tracks that are within 5 cm of the cluster position
	  for (int k = 0; k < mEvent->numberOfTracks(); ++k) {
	    StGammaTrack* track = mEvent->track(k);
	    try {
	      static const EEmcGeomSimple& geo = EEmcGeomSimple::Instance();
	      TVector3 position2 = track->positionAtZ(geo.getZSMD());
	      TVector3 dca = position - position2;
	      if (dca.Mag() < 5) point->addTrack(track);
	    }
	    catch (StGammaTrack::Exception& e) {}
	  }
	}
      }
    } // Loop over matched clusters

    // Reassign energies to points associated with this EMC cluster
    // by sharing the EMC cluster energy in proportion to the individual
    // SMD clusters energies.
    float emcEnergy = emcCluster->energy();
    float scale = emcEnergy / smdEnergy;
    for (size_t j = 0; j < points.size(); ++j) {
      StEndcapEmcPoint* point = points[j];
      point->setEnergy(scale * point->energy());
    }

  } // Loop over EMC clusters
}

int StEtaFinder::numberOfTracks(StEndcapEmcCluster* cluster) const
{
  int ntracks = 0;
  for (int i = 0; i < cluster->numberOfTowers(); ++i) {
    StGammaTower* tower = cluster->tower(i);
    ntracks += mTracks.count(tower);
  }
  return ntracks;
}

TVector3 StEtaFinder::getIntersection(StEndcapSmdCluster* cluster1, StEndcapSmdCluster* cluster2) const
{
  if (cluster1->sector() != cluster2->sector()) return TVector3(0,0,-999);
  if (cluster1->plane () == cluster2->plane ()) return TVector3(0,0,-999);
  int sector = cluster1->sector();
  int uid = cluster1->maxStrip()->index;
  int vid = cluster2->maxStrip()->index;
  return EEmcSmdGeom::instance()->getIntersection(sector, uid, vid);
}

StGammaTower* StEtaFinder::getTower(const TVector3& position) const
{
  int sector, subsector, etabin;
  float deta, dphi;
  if (EEmcGeomSimple::Instance().getTower(position, sector, subsector, etabin, deta, dphi)) {
    if (fabs(deta) < 0.7 && fabs(dphi) < 0.7) {
      int phibin = sector * 5 + subsector;
      return mTowers[etabin][phibin];
    }
  }
  return 0;
}

float StEtaFinder::isolation(StGammaTower* tower) const
{
  float pt = 0;
  for (int deta = -1; deta <= 1; ++deta) {
    for (int dphi = -1; dphi <= 1; ++dphi) {
      if (deta || dphi) {
	StGammaTower* tower2 = getNextTower(tower, deta, dphi);
	if (tower2) pt += tower2->pt();
      }
    }
  }
  return pt / tower->pt();
}

bool StEtaFinder::electron(StEndcapEmcPoint* point) const
{
  if (point->numberOfTracks() == 1) {
    StGammaTrack* track = point->track(0);
    float E = point->energy();
    float p = track->momentum().Mag();
    float EoverP = E / p;
    EoverP = 1.0;
    float dEdx = track->dEdx();
    return ((0.7 < EoverP && EoverP < 1.3) && (3.3 < dEdx && dEdx < 4.1));
  }
  return false;
}

int StEtaFinder::getCandidates()
{
  for (int iCandidate = 0; iCandidate < mEvent->numberOfCandidates(); ++iCandidate) {
    StGammaCandidate* candidate = mEvent->candidate(iCandidate);
    // 0=EEMC, 1=BEMC
    if (candidate->detectorId() == 0) {
      // No tracks must point to candidate
      if (candidate->numberOfTracks()) continue;

      // 90% of the EM energy over a cone of radius 0.7
      // in eta-phi space must be contained in the candidate
      // EM cluster.
      const float radius = 0.7;
      float ratio = candidate->momentum().Pt() / candidate->sumTowerPt(radius);
      if (ratio < 0.9 || ratio > 1.0) continue;

      mCandidates.push_back(candidate);
    }
  }
  return mCandidates.size();
}

void StEtaFinder::matchCandidates()
{
  for (size_t i = 0; i < mCandidates.size(); ++i) {
    StGammaCandidate* candidate = mCandidates[i];
    vector<StGammaPoint> points;
    float sumEnergy = 0;
    for (size_t j = 0; j < mMatchedClusters.size(); ++j) {
      StEndcapSmdCluster* cluster1 = mMatchedClusters[j].first;
      StEndcapSmdCluster* cluster2 = mMatchedClusters[j].second;
      TVector3 position = getIntersection(cluster1, cluster2);
      if (position.z() != -999) {
	TVector3 separation = candidate->position() - position;
	if (separation.Mag() < 30) {
	  StGammaPoint point;
	  point.candidate = candidate;
	  point.smdu = cluster1;
	  point.smdv = cluster2;
	  point.energy = 0.5 * (cluster1->energy() + cluster2->energy());
	  point.position = position;
	  sumEnergy += point.energy;
	  points.push_back(point);
	}
      }
    } // Loop over matched clusters

    float scale = candidate->energy() / sumEnergy;
    for (size_t j = 0; j < points.size(); ++j) {
      StGammaPoint& point = points[j];
      point.energy *= scale;
      mGammaPoints.push_back(point);
    }
  } // Loop over candidates
}

void StEtaFinder::findTowerPoints()
{
  // Find seed towers with pT > 0.8 GeV and no track
  vector<StGammaTower*> seedTowers;

  for (int etabin = 0; etabin < 12; ++etabin) {
    for (int phibin = 0; phibin < 60; ++phibin) {
      StGammaTower* tower = mTowers[etabin][phibin];
      if (tower && tower->pt() > 0.8) seedTowers.push_back(tower);
    }
  }

  // Loop over each seed tower and look for a SMD cluster
  for (size_t iSeed = 0; iSeed < seedTowers.size(); ++iSeed) {
    StGammaTower* tower = seedTowers[iSeed];
    mPointCounter->Fill("Seed tower pT > 0.8 GeV", 1);

    // Require no track to tower
    if (mTracks.count(tower)) continue;
    mPointCounter->Fill("Track isolation", 1);    

    // Get ranges of SMD strips under tower
    int sector    = tower->sector();
    int subsector = tower->subsector();
    int etabin    = tower->etabin();
    int xmin[2], xmax[2];
    EEmcSmdMap::instance()->getRangeU(sector, subsector, etabin, xmin[0], xmax[0]);
    EEmcSmdMap::instance()->getRangeV(sector, subsector, etabin, xmin[1], xmax[1]);

    // Find max strips under ranges
    StGammaStrip* maxStrips[2] = { 0, 0 };
    for (int plane = 0; plane < 2; ++plane) {
      for (int i = xmin[plane]; i <= xmax[plane]; ++i) {
	StGammaStrip* strip = mStrips[sector][plane][i];
	if (strip) {
	  if (!maxStrips[plane] || strip->energy > maxStrips[plane]->energy) maxStrips[plane] = strip;
	}
      }
    }

    // Skip if no max strip in one plane
    if (!(maxStrips[0] && maxStrips[1])) continue;
    mPointCounter->Fill("Max strips in 2 planes", 1);

    // 2 MeV minimum per max strip (no MIP)
    if (maxStrips[0]->energy < 2 * MeV || maxStrips[1]->energy < 2 * MeV) continue;
    mPointCounter->Fill("Max strips energy > 2 MeV", 1);

    // Get intersection of max strips and check that it lies within the inner 70%
    // of the fiducial volume of the seed tower.
    int uid = maxStrips[0]->index;
    int vid = maxStrips[1]->index;
    TVector3 position = EEmcSmdGeom::instance()->getIntersection(sector, uid, vid);

    // Skip if no valid intersection
    if (position.z() == -999) continue;
    mPointCounter->Fill("Max strips intersect", 1);

    // Skip if no tower at intersection
    int sector2, subsector2, etabin2;
    float dphi, deta;
    if (!EEmcGeomSimple::Instance().getTower(position, sector2, subsector2, etabin2, dphi, deta)) continue;
    mPointCounter->Fill("Tower at intersection", 1);
    if (fabs(dphi) > 0.7 || fabs(deta) > 0.7) continue;
    mPointCounter->Fill("Inside fiducial volume", 1);

    // Skip if tower at intersection doesn't match seed tower
    if (sector != sector2 || subsector != subsector2 || etabin != etabin2) continue;
    mPointCounter->Fill("Intersection not in seed tower", 1);

    // Enter & exit the same tower
    TVector3 smd = position - mEvent->vertex();
    float cosTheta = smd.CosTheta();
    TVector3 v1 = smd;
    v1.SetMag((EEmcGeomSimple::Instance().getZ1() - mEvent->vertex().z()) / cosTheta);
    v1 += mEvent->vertex();
    TVector3 v2 = smd;
    v2.SetMag((EEmcGeomSimple::Instance().getZ2() - mEvent->vertex().z()) / cosTheta);
    v2 += mEvent->vertex();
    int sector3, subsector3, etabin3;
    if (!EEmcGeomSimple::Instance().getTower(v1, sector3, subsector3, etabin3)) continue;
    mPointCounter->Fill("Tower at z1", 1);
    if (sector != sector3 || subsector != subsector3 || etabin != etabin3) continue;
    mPointCounter->Fill("Tower at z1 differs from seed tower", 1);
    int sector4, subsector4, etabin4;
    if (!EEmcGeomSimple::Instance().getTower(v2, sector4, subsector4, etabin4)) continue;
    mPointCounter->Fill("Tower at z2", 1);
    if (sector != sector4 || subsector != subsector4 || etabin != etabin4) continue;
    mPointCounter->Fill("Tower at z2 differs from seed tower", 1);

    // Get 7-strip and 11-strip energy sums
    float E7 [2] = { 0, 0 };
    float E11[2] = { 0, 0 };
    for (int plane = 0; plane < 2; ++plane) {
      int xmin7 = maxStrips[plane]->index - 5;
      int xmax7 = maxStrips[plane]->index + 5;
      int xmin11 = maxStrips[plane]->index - 20;
      int xmax11 = maxStrips[plane]->index + 20;
      for (int i = 0; i < 288; ++i) {
	if (StGammaStrip* strip = mStrips[sector][plane][i]) {
	  if (xmin11 <= i && i <= xmax11) {
	    E11[plane] += strip->energy;
	    if (xmin7 <= i && i <= xmax7) E7[plane] += strip->energy;
	  }
	}
      }
    }

    // Require 70% of energy to be in 7-strip cluster
    if (E7[0] / E11[0] < 0.7) continue;
    if (E7[1] / E11[1] < 0.7) continue;
    mPointCounter->Fill("SMD isolation", 1);

    // Require energy asymmetry between planes < 20%
    if (fabs(E7[0] - E7[1]) / (E7[0] + E7[1]) > 0.2) continue;
    mPointCounter->Fill("SMD energy asymmetry", 1);

    // Create point
    StTowerPoint point;
    point.energy = tower->energy;
    point.position = position;
    point.sector = sector;
    point.u_strip = maxStrips[0]->index;
    point.v_strip = maxStrips[1]->index;
    mTowerPoints.push_back(point);
  }
}

void StEtaFinder::print(const StTowerPoint& point)
{
  TH1F hU("hU", ";u-strip", 288, 0, 288);
  TString title1 = Form("run=%d, event=%d, sector=%d", mEvent->runNumber(), mEvent->eventNumber(), point.sector);
  hU.SetTitle(title1);
  for (int i = 0; i < 288; ++i) {
    StGammaStrip* strip = mStrips[point.sector][0][i];
    if (strip) hU.Fill(strip->index, strip->energy);
  }

  hU.SetAxisRange(point.u_strip - 40, point.u_strip + 40);
  TH1F hV("hV", ";v-strip", 288, 0, 288);
  double pt = point.momentum(mEvent->vertex()).Pt();
  TString title2 = Form("p_{T}=%.1lf GeV", pt);
  hV.SetTitle(title2);
  for (int i = 0; i < 288; ++i) {
    StGammaStrip* strip = mStrips[point.sector][1][i];
    if (strip) hV.Fill(strip->index, strip->energy);
  }

  hV.SetAxisRange(point.v_strip - 40, point.v_strip + 40);

  float ymax = 1.1 * max(hU.GetMaximum(), hV.GetMaximum());

  hU.SetMaximum(ymax);
  hV.SetMaximum(ymax);

  TLine line1(point.u_strip, 0, point.u_strip, ymax);
  TLine line2(point.v_strip, 0, point.v_strip, ymax);

  line1.SetLineColor(kRed);
  line2.SetLineColor(kRed);

  if (mCanvas) delete mCanvas;
  mCanvas = new TCanvas;

  mCanvas->Divide(2, 1);

  mCanvas->cd(1);
  hU.Draw();
  line1.Draw();

  mCanvas->cd(2);
  hV.Draw();
  line2.Draw();

  mCanvas->Print(mPsFile);
}