//
// Pibero Djawotho <pibero@indiana.edu>
// Indiana University
// 12 July 2007
//

class StGammaEvent;

TH1F* triggerIds;
TNtuple* photons;
TCanvas* c1;

void ReadGammaTree2(const char* filename = "@gamma2.list")
{
  // Load shared libraries
  gROOT->Macro("StRoot/StGammaMaker/macros/loadGammaLibs.C");

  // Create output ROOT file
  TFile* ofile = new TFile("gamma.root", "recreate");
  assert(ofile);

  // Create histograms
  triggerIds = new TH1F("triggerIds", "triggerIds", 5, 0, 5);
  photons = new TNtuple("photons", "Photons Tuple", "pt:E:eta:phi:pre1:pre2:post:x:y:ntowers:unhits:vnhits:uyield:vyield:uresidual:vresidual:umean:vmean:usigma:vsigma:smdx:smdy");

  // Create canvas
  c1 = new TCanvas("c1", "c1", 800, 600);
  c1->Divide(2, 1);

  c1->Print("gamma.ps[");

  // Read input ROOT file(s)
  if (*filename == '@') {
    // File list
    TChain* chain = new TChain("gammas");
    ++filename;
    ifstream in(filename);
    assert(in);
    string line;
    while (getline(in, line)) chain->AddFile(line.c_str());
    in.close();

    StGammaEvent* event = new StGammaEvent;
    chain->SetBranchAddress("GammaEvent", &event);

    int nentries = chain->GetEntries();
    cout << "nentries = " << nentries << endl;

    for (int ev = 0; ev < nentries; ++ev) {
      chain->GetEntry(ev);
      processEvent(event);
    }
  }
  else {
    // Single file
    TFile* file = new TFile(filename);
    assert(file && !file->IsZombie());

    TTree* gammas = (TTree*)file->Get("gammas");
    assert(gammas);

    StGammaEvent* event = new StGammaEvent;
    gammas->SetBranchAddress("GammaEvent", &event);

    int nentries = gammas->GetEntriesFast();
    cout << "nentries = " << nentries << endl;

    for (int ev = 0; ev < gammas->GetEntriesFast(); ++ev) {
      gammas->GetEntry(ev);
      processEvent(event);
    }
  }

  c1->Print("gamma.ps]");

  triggerIds->LabelsDeflate();

  // Write and close output ROOT file
  ofile->Write();
  ofile->Close();
}

void processEvent(StGammaEvent* event)
{
  // Select trigger id
  if (!(isTrigger(event, 117641) || // eemc-http-mb-l2gamma
	isTrigger(event, 127641) || // eemc-http-mb-l2gamma
	isTrigger(event, 137641)))  // eemc-http-mb-l2gamma
    return;

//   if (isTrigger(event, 127652) || // eemc-jjp0-etot-mb-L2jet
//       isTrigger(event, 127271) || // eemc-jp1-mb
//       isTrigger(event, 127551))	  // eemc-jp0-mb
//     return;

  for (int i = 0; i < event->triggerIds().GetSize(); ++i)
    triggerIds->Fill(Form("%d", event->triggerIds().At(i)), 1);

  for (int i = 0; i < event->numberOfCandidates(); ++i) {
    StGammaCandidate* candidate = event->candidate(i);

    // EEMC candidates only
    if (candidate->detectorId() != StGammaCandidate::kEEmc) continue;

    // Require a minimum of 3 SMD hits in each plane
    if (candidate->numberOfSmdu() < 3 || candidate->numberOfSmdv() < 3) continue;

    // Require *exactly* no track in candidate
    if (candidate->numberOfMyTracks()) continue;

    // Isolation cut
    cout << "sum track pt = " << candidate->sumTrackPt(0.3) << endl;
    cout << "cluster pt = " << candidate->momentum().Pt() << endl;

    if (candidate->sumTrackPt(0.3) > 0.15 * candidate->momentum().Pt()) continue;

    cout << "Passed isolation cut" << endl;

    // Photon selection criteria
    if (candidate->pre1Energy() < 0.5e-3) continue;
    if (candidate->pre2Energy() < 2.0e-3) continue;
    if (candidate->postEnergy() > 0.5e-3) continue;

    // Fit
    StGammaFitterResult u;
    StGammaFitterResult v;

    if (!StGammaFitter::instance()->fitSector(candidate, &u, &v)) continue;

#if 1
    // Fill SMD response histograms
    TH1F* hU = new TH1F("hSmdu", "hSmdu", 288, 0, 288);
    TH1F* hV = new TH1F("hSmdv", "hSmdv", 288, 0, 288);

    for (int k = 0; k < candidate->numberOfSmdu(); ++k) {
      StGammaStrip* strip = candidate->smdu(k);
      hU->SetBinContent(strip->index + 1, 1000 * strip->energy);
    }

    for (int k = 0; k < candidate->numberOfSmdv(); ++k) {
      StGammaStrip* strip = candidate->smdv(k);
      hV->SetBinContent(strip->index + 1, 1000 * strip->energy);
    }

    const char* formula = "[0]*(0.69*exp(-0.5*((x-[1])/0.87)**2)/(sqrt(2*pi)*0.87)+0.31*exp(-0.5*((x-[1])/3.3)**2)/(sqrt(2*pi)*3.3))";
    TF1* fU = new TF1("fU", formula, 0, 288);
    TF1* fV = new TF1("fV", formula, 0, 288);

    fU->SetParameters(u.yield(), u.mean() + 1);
    fV->SetParameters(v.yield(), v.mean() + 1);

    // Draw
    const char* title1 = Form("Run=%d, Event=%d", event->runNumber(), event->eventNumber());
    const char* title2 = Form("p_{T}=%.1lf GeV, E=%.1lf GeV", (double)candidate->momentum().Pt(), (double)candidate->energy());
    float ymax = 1.1 * max(hU->GetMaximum(), hV->GetMaximum());

    hU->SetAxisRange(u.mean() - 50, u.mean() + 50);
    hV->SetAxisRange(v.mean() - 50, v.mean() + 50);

    double x1 = max(u.mean() - 45, 5.);
    double y1 = 0.60 * ymax;
    double x2 = x1 + 20;
    double y2 = 0.91 * ymax;

    TPaveText* pave = new TPaveText(x1, y1, x2, y2);
    pave->AddText("Triggers")->SetTextColor(kBlue);

    for (int k = 0; k < event->triggerIds().GetSize(); ++k) {
      int id = event->triggerIds().At(k);
      const char* name = Form("%d", id);
      TText* text = pave->AddText(name);
      if (id == 117641 || id == 127641 || id == 137641) text->SetTextColor(kRed);
    }

    c1->cd(1);
    hU->SetTitle(title1);
    hU->SetMaximum(ymax);
    hU->Draw();
    fU->Draw("same");
    pave->Draw();

    c1->cd(2);
    hV->SetTitle(title2);
    hV->SetMaximum(ymax);
    hV->Draw();
    fV->Draw("same");

    c1->Print("gamma.ps");

//     char c;
//     do {
//       c1->Update();
//       cout << "Continue? (y/n) ";
//       cin >> c;
//     } while (c == 'n');

    // Clean up
    delete hU;
    delete hV;
    delete fU;
    delete fV;
    delete pave;
#endif

    // Get interssection of u- and v-strip in SMD plane
    TVector3 position = EEmcSmdGeom::instance()->getIntersection(candidate->smdu(0)->sector, u.mean(), v.mean());
    if (position.z() == -999) position.SetXYZ(0,0,0);

    float tuple[] = {
      candidate->momentum().Pt(),
      candidate->energy(),
      candidate->momentum().Eta(),
      candidate->momentum().Phi(),
      candidate->pre1Energy(),
      candidate->pre2Energy(),
      candidate->postEnergy(),
      candidate->position().x(),
      candidate->position().y(),
      candidate->numberOfMyTowers(),
      candidate->numberOfSmdu(),
      candidate->numberOfSmdv(),
      u.yield(),
      v.yield(),
      u.residual(),
      v.residual(),
      u.mean(),
      v.mean(),
      u.sigma(),
      v.sigma(),
      position.x(),
      position.y()
    };

    photons->Fill(tuple);
  }
}

bool isTrigger(StGammaEvent* event, int id)
{
  for (int i = 0; i < event->triggerIds().GetSize(); ++i)
    if (event->triggerIds().At(i) == id) return true;
  return false;
}