#include <algorithm>
#include <sstream>
#include <iostream>
#include <map>
#include <cmath>
#include <vector>
#include <algorithm>

#include <TFile.h>
#include <TTree.h>
#include <TMath.h>
#include <TGraphErrors.h>
#include <TCut.h>
#include <TEventList.h>

using namespace std;

const double degree = TMath::Pi() / 180.;

struct Event {
  double mjd;
  double zen;
  double azi;
};

struct by_mjd {
  bool operator()(Event const &a, Event const &b) {
    return a.mjd < b.mjd;
  }
};

double phiDiff(double phi1, double phi2) {
  double diff = fabs(phi1 - phi2);

  if (diff > 360 * degree)
    return diff - 360 * degree;

  return diff;
}

bool isEventInWindow(double eventAzi, double moonAzi, double eventZen, double moonZen, double aziWindow, double zenWindow) {
  return phiDiff(eventAzi,moonAzi) < aziWindow && fabs(eventZen-moonZen) < zenWindow;
}

int main(int argc, char * argv[]) {

  TFile * ifile = new TFile("./testData.root","READ");
  TTree * tree = (TTree*) ifile->Get("aTree");

  double deltaAzi = 10 * degree;
  double deltaZen = 10 * degree;
  int NChanCut = 12;
  int NStringCut = 3;

  ostringstream geoCutStr;
  geoCutStr << "NCh >= " << NChanCut << " && NString >= " << NStringCut;

  TCut geoCut(geoCutStr.str().c_str()); 

  double mZen, mAzi, mMJD, TrueMJD;
  double LLHAzi, LLHZen;

  tree->SetBranchAddress("CorsikaMoonZenith",&mZen);
  tree->SetBranchAddress("CorsikaMoonAzimuth",&mAzi);
  tree->SetBranchAddress("CorsikaMoonMJD",&mMJD);
  tree->SetBranchAddress("MJDTime",&TrueMJD);

  int NEvents = tree->GetEntries();
  vector<Event> moonEvents;

  tree->GetEntry(0);
  double trueMJDStart = TrueMJD;
  tree->GetEntry(NEvents-1);
  double trueMJDEnd = TrueMJD;
  double DeltaT = (trueMJDEnd - trueMJDStart) * 24 * 3600;

  for (int i = 0; i < NEvents; i++) {
    tree->GetEntry(i);
    if (mMJD == mMJD) {
      Event myEvent;
      myEvent.mjd = mMJD;
      myEvent.zen = mZen;
      myEvent.azi = mAzi;

      moonEvents.push_back(myEvent);
    }
  }

  sort(moonEvents.begin(), moonEvents.end(), by_mjd());  

  double MJDStart = moonEvents[0].mjd;
  double MJDEnd = moonEvents[moonEvents.size()-1].mjd;

  int Nperiods = 30;
  cout << "Events in file: " << moonEvents.size() << endl;
  cout << "Real time length (s): " << DeltaT << endl;

  tree->SetBranchAddress("PoleMuonLlhFit_Azimuth",&LLHAzi);
  tree->SetBranchAddress("PoleMuonLlhFit_Zenith",&LLHZen);

  for (int i = 0; i < Nperiods; i++) {
    double dt = moonEvents.size() / double(Nperiods);
    int index = int(i * dt);

    double currentMZen = moonEvents[index].zen;
    double currentMAzi = moonEvents[index].azi;

    ostringstream moonCutStr;

    moonCutStr << "isEventInWindow(PoleMuonLlhFit_Azimuth," << currentMAzi << ",";
    moonCutStr << "PoleMuonLlhFit_Zenith, " << currentMZen;
    moonCutStr << ", " << deltaAzi;
    moonCutStr << ", " << deltaZen << ")";

    TCut moonCut(moonCutStr.str().c_str());
    cout << moonCutStr.str() << endl;

    tree->Draw(">>elist", moonCut + geoCut);
    TEventList *elist = (TEventList*)gDirectory->Get("elist");
    int Npass = elist->GetN(); 
    cout << moonEvents[index].mjd << "\t" << Npass/DeltaT << endl;
  }

  return 0;

}