#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TF1.h>
#include <TProfile.h>
#include <TGraphAsymmErrors.h>

#include <PamCalo.h>

void posSelection(TString inputFile){

  // Create some histograms.
  // electrons
  TH1D *elesample = new TH1D("elesample","Electron sample;GeV;counts",100,0.,20.);
  TH1D *elesel = new TH1D("elesel","Selected electrons;GeV;counts",100,0.,20.);
  // positrons
  TH1D *possel = new TH1D("possel","Selected positrons;GeV;counts",100,0.,20.);
  // sum ep + em
  TH1D *leptosel = new TH1D("leptpsel","Selected positrons+electrons;GeV;counts",100,0.,20.);

  // scatter plots
  TH2D *qtotene = new TH2D("qtotene","qtot/energy;GeV;qtot/energy",1000,-20.,20.,1000,0.,500.); // energy momentum match
  TH2D *noint = new TH2D("noint","noint;GeV;noint",1000,-20.,20.,1000,0.,35.); // noint distribution
  TH2D *qncore = new TH2D("qncore","qcore/ncore;GeV;qcore/ncore",1000,-20.,20.,5000,0.,50.); // qcore/ncore distribution

  // Functions for selections
  TF1 *fnoint = new TF1("fnoint","abs(0.5+(-100./x))",-30,30);

  TF1 *fqncorelow = new TF1("fqncorelow","0.7345+0.5911*abs(x)",-30,30);
  TF1 *fqncoreup = new TF1("fqncoreup","4.7345+0.5911*abs(x)",-30,30);

  // Open the file
  TFile *file = TFile::Open(inputFile);
  TTree *tree = (TTree*)file->Get("pamcalotree");
  PamCalo *pc = new PamCalo();  
  tree->SetBranchAddress("PamCalo",&pc);

  // loop over events
  for ( Int_t i=0; i<tree->GetEntries(); i++) {
    //for ( Int_t i=0; i<20000; i++) {
    
    if ( i%10000 == 0 ) cout << " get entry " << i << "\n";

    tree->GetEntry(i);
    
    if ( pc->energy < 0. ) elesample->Fill(fabs(pc->energy));

    // selection condition (no selection at the moment)
    if ( pc->noint < fnoint->Eval(pc->energy) // noint selection
	 && pc->ncore > 0 && pc->qcore/pc->ncore > fqncorelow->Eval(pc->energy) &&  pc->qcore/pc->ncore < fqncoreup->Eval(pc->energy) // qncore selection
	 ){ 
      
      // fill scatter plots
      noint->Fill(pc->energy, pc->noint);
      if ( pc->ncore > 0 ) qncore->Fill(pc->energy,pc->qcore/(float)pc->ncore);
      if ( pc->energy != 0. ) qtotene->Fill(pc->energy,pc->qtot/fabs(pc->energy));
      // fill histrograms
      if ( pc->energy < 0. ){ // negatives particles
	elesel->Fill(fabs(pc->energy));
      } else { // positive particles
	possel->Fill(fabs(pc->energy));
      }
      leptosel->Fill(fabs(pc->energy));

    }

  }

  // draw histograms 
  TCanvas *cqtotenergy = new TCanvas("cqtotenergy","Energy momentum match",200,10,600,400);
  cqtotenergy->SetGrid();
  cqtotenergy->SetTicks();
  cqtotenergy->SetLogz();
  cqtotenergy->Draw();
  qtotene->Draw("colz");

  TCanvas *cnoint = new TCanvas("cnoint","Distribution: noint",200,10,600,400);
  cnoint->SetGrid();
  cnoint->SetTicks();
  cnoint->SetLogz();
  cnoint->Draw();
  noint->Draw("colz");
  fnoint->Draw("same");

  TCanvas *cqncore = new TCanvas("cqncore","Distribution: qcore/ncore",200,10,600,400);
  cqncore->SetGrid();
  cqncore->SetTicks();
  cqncore->SetLogz();
  cqncore->Draw();
  qncore->Draw("colz");
  fqncorelow->Draw("same");
  fqncoreup->Draw("same");


  // draw electron efficiency and positron fraction
  TCanvas *chisto = new TCanvas("chisto","Electron efficiency and positron fraction",200,10,800,800);
  chisto->Divide(1,2); // columns, raws
  chisto->GetPad(1)->SetLogx();
  chisto->GetPad(2)->SetLogx();
  chisto->GetPad(1)->SetGrid();
  chisto->GetPad(2)->SetGrid();
  chisto->GetPad(1)->SetTicks();
  chisto->GetPad(2)->SetTicks();
  chisto->Draw();

  // create electron efficiency, selected electrons / sample electrons
  TGraphAsymmErrors* effele = new TGraphAsymmErrors(elesel,elesample); 
  effele->SetMarkerStyle(20);

  // create a TH2 to be used as background
  TH2D *hbkgef = new TH2D("hbkgef",";Energy (GeV);e^{-} selection efficiency",1000,0.,30.,1000,0.,1.2);
  hbkgef->SetStats(0);

  // go to the first panel and draw
  chisto->cd(1);
  hbkgef->Draw();
  effele->Draw("Psame");

  // create positron fraction, selected positrons / (selected electrons + selected positrons)
  TGraphAsymmErrors* posfraction = new TGraphAsymmErrors(possel,leptosel); 
  posfraction->SetMarkerStyle(20);

  // create a TH2 to be used as background
  TH2D *hbkgpf = new TH2D("hbkgpf",";Energy (GeV);Positron fraction",1000,0.,30.,1000,0.,1.2);
  hbkgpf->SetStats(0);

  // go to the second panel and draw
  chisto->cd(2);
  hbkgpf->Draw();
  posfraction->Draw("Psame");


}