#include <ROOT/TDataFrame.hxx>
#include <ROOT/TResultProxy.hxx>
#include <TApplication.h>
#include <TStopwatch.h>
#include <TFile.h>
#include <TH1F.h>
#include <fstream>
#include <iostream>
#include <vector>
#include <TChain.h>
#include <TPaveText.h>
#include <TPaveLabel.h>
#include <TLegend.h>
#include <string>
#include <TVector3.h>
#include <TVectorD.h>
#include <TMath.h>
#include <TStyle.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TGraph2D.h>
#include <TGraph.h>
#include <TGaxis.h>
#include <TLatex.h>
#include <TLine.h>
#include <TProfile.h>

using namespace std;

#define DFHisto(name,nbins,lower,upper,var,weight) Histo1D<double>(TH1D(name, "", nbins, lower, upper), var, weight).GetValue()

class Analysis {

public:
  Analysis (string title, 
	string pid, string mom, string id, 
	double theta_min, double theta_max, string file, 
	const bool cut_weights = false, const double gem_track_cut = 1.0E30,
	const double vz = 300.0,
  const double scint_threshold = 1.0,
  const string rescatter = "frame");

public:
	string title_;
  string particle_;
  string momentum_;
  string id_;
  string input_file_;
  double theta_min_;
  double theta_max_;
  string rescatter_;
  double scint_threshold_;

  int number_of_events_;
  int number_of_thrown_events_;
  bool is_normalized_;
  bool run_ok_;

  TH1D h_theta_;
  TH1D h_theta_tgt_;
  TH1D h_theta_no_tgt_;
  TH1D h_vertex_dz_;
  TH1D h_vertex_dx_;
  TH1D h_vertex_dy_;
  TH1D h_dtheta_;

};

Analysis::Analysis(string title, 
	string pid, string mom, string id, 
	double theta_min, double theta_max, string file, 
	const bool cut_weights, 
	const double gem_track_cut,
	const double vz,
  const double scint_threshold,
  const string rescatter) {

	title_ = title;
	particle_ = pid;
	momentum_ = mom + " MeV/c";
	theta_min_ = theta_min;
	theta_max_ = theta_max;
	scint_threshold_ = scint_threshold;
	rescatter_ = rescatter;
	run_ok_ = false;

	id_ = Form("%s_%.0f_%.0f", id.c_str(), theta_min, theta_max);
  input_file_ = file;
  is_normalized_ = false;

  TFile f(input_file_.c_str());
  if (!f.IsOpen()) return;

  TVectorD* v1 = (TVectorD*) f.Get("RC");
  number_of_thrown_events_ = (*v1)[0];
  number_of_events_ = (*v1)[1];
  f.Close();

  std::cout << "We are looping over " << number_of_events_ << " events." << std::endl;

  TStopwatch watch;

  //ROOT::EnableImplicitMT(1);

  watch.Reset();
  watch.Start();
  ROOT::Experimental::TDataFrame df("T", file.c_str());//Create TDataFrame
  //If we're cutting weights
  /*auto cut_weight = [&](double w) {if (w >=1.0) return 0.0; else return w;};
  df.Define("wgt", cut_weight, {"weight"});*/

  //Which rescattering were interesed in
  /*auto rescatter_opt = [&](bool r) {return r;};
  if (rescatter_ == "chamber") {
    df.Define("rescatter_option", rescatter_opt, {"CHAMBER_InScatter"});
  }
  if (rescatter_ == "frame") {
    df.Define("rescatter_option", rescatter_opt, {"CHAMBER_hit_ds_frame"});
  }*/

  //Filters for Reconstructed events
  auto Select_Reco = [&](ROOT::Experimental::TDataFrame &dataFrame) {
    auto ret = 
      dataFrame.Filter([](bool veto_hit) {return veto_hit == false;}, {"VSC_hit"})
      .Filter([](bool frame_hit) {return frame_hit == false;}, {"frame_hit"})
      .Filter([](bool hit_blsc) {return hit_blsc == false;}, {"BLSC_hit"})
      .Filter([](bool MuonDecay) {return MuonDecay == false;}, {"MuonDecay"})
      .Filter([&](double theta) {return theta < theta_max * TMath::DegToRad();}, {"recon_theta"})
      .Filter([&](double theta) {return theta > theta_min * TMath::DegToRad();}, {"recon_theta"})
      .Filter([](double doca) {return doca < 25;}, {"recon_doca"})
      .Filter([](TVector3 &vertex) {return fabs(vertex.z()) < 90.0;}, {"recon_vertex"})
      .Filter([](TVector3 &vertex) {return fabs(vertex.y()) < 100.0;}, {"recon_vertex"})
      .Filter([](TVector3 &vertex) {return fabs(vertex.x()) < 45.0;}, {"recon_vertex"})
      .Filter([&](TVector3 &recon_tgt_pos) {return recon_tgt_pos.Perp() < gem_track_cut;}, {"recon_tgt_hit_pos"});
    return ret;
  };

  //Filters for Reconstructed Signal Events
  auto Select_Reco_Signal = [&](ROOT::Experimental::TDataFrame &dataFrame) {
    auto ret = 
      dataFrame.Filter([](bool veto_hit) {return veto_hit == false;}, {"VSC_hit"})
        .Filter([](bool frame_hit) {return frame_hit == false;}, {"frame_hit"})
        .Filter([](bool hit_blsc) {return hit_blsc == false;}, {"BLSC_hit"})
        .Filter([](bool MuonDecay) {return MuonDecay == false;}, {"MuonDecay"})
        .Filter([&](double theta) {return theta < (theta_max*TMath::DegToRad() && theta > theta_min*TMath::DegToRad());}, {"recon_theta"})
        .Filter([](double doca) {return doca < 25;}, {"recon_doca"})
        .Filter([](TVector3 &vertex) {return fabs(vertex.z()) < 90.0;}, {"recon_vertex"})
        .Filter([](TVector3 &vertex) {return fabs(vertex.y()) < 100.0;}, {"recon_vertex"})
        .Filter([](TVector3 &vertex) {return fabs(vertex.x()) < 45.0;}, {"recon_vertex"})
        .Filter([&](TVector3 &recon_tgt_pos) {return recon_tgt_pos.Perp() < gem_track_cut;}, {"recon_tgt_hit_pos"})
        .Filter([](bool tgt_event) {return tgt_event == true;}, {"CHAMBER_target_event"})
        .Filter([](bool is_coulomb) {return is_coulomb == true;}, {"CHAMBER_is_coulomb"})
        .Filter([](bool is_lh2) {return is_lh2 == true;}, {"CHAMBER_is_lh2"})
        .Filter([&](double tgt_theta) {return (tgt_theta < theta_max*TMath::DegToRad() && tgt_theta > theta_min*TMath::DegToRad()); }, {"CHAMBER_Theta"});
      return ret;
  };

  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
  //Reconstructed Events
  auto selected_reco = Select_Reco(df);

  auto coulomb = selected_reco.Filter([](bool c) {return c == true;}, {"CHAMBER_is_coulomb"}).Count();
  auto eioni = selected_reco.Filter([](bool e) {return e == true;}, {"CHAMBER_is_eioni"}).Count();
  auto ebrem = selected_reco.Filter([](bool b) {return b == true;}, {"CHAMBER_is_ebrem"}).Count();

  h_theta_ = (TH1D)selected_reco.DFHisto(("h_theta"+id).c_str(), 80, 0, 2, "recon_theta", "weight");
  h_theta_tgt_ = (TH1D)selected_reco.Filter([](bool is_lh2) {return is_lh2 == true;}, {"CHAMBER_is_lh2"}).DFHisto(("h_theta_tgt_"+id).c_str(), 80, 0, 2, "recon_theta", "weight");
  h_theta_no_tgt_ = (TH1D)selected_reco.Filter([](bool is_lh2) {return is_lh2 == false;}, {"CHAMBER_is_lh2"}).DFHisto(("h_theta_no_tgt_"+id).c_str(), 80, 0, 2, "recon_theta", "weight");

  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
  //Reconstructed Signal Events
  auto selected_signal_reco = Select_Reco_Signal(df);

  auto dtheta = [](double x, double y) { return (x - y) * 1000.0; };
  h_dtheta_ = (TH1D)selected_signal_reco.Define("dtheta", dtheta, {"recon_theta", "CHAMBER_Theta"}).DFHisto(("h_dtheta_"+id).c_str(), 100, -80, 80, "dtheta", "weight");


  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
  cout << "  Norm: tgt yield = "
   << h_theta_tgt_.Integral("width") / number_of_thrown_events_ * 1.0E6
   << " per 10^6 beam particles\n";
  cout << "  TGT processes: Coulomb = " << *coulomb << " eIoni = " << *eioni << " eBrem = " << *ebrem << "\n";
  cout << ""<< "\n";  

  run_ok_ = true;

  watch.Stop();

  cout << "The code took " << watch.RealTime() << " seconds." << endl;

}