// largely based on Mike Hance's code

#include <sstream>
#include <iostream>
#include <stdlib.h>
#include <unistd.h>

// Input file (uncomment only one)

// My files:
//TFile* fin = new TFile("Conversions37.aan.root"); const double erMin=-0.06; const double erMax=0.06;
//TFile* fin = new TFile("Conversions35.aan.root"); const double erMin=-0.06; const double erMax=0.06;
//TFile* fin = new TFile("Conversions33.aan.root");  const double erMin=-0.06; const double erMax=0.16;
//TFile *fin = new TFile("ConversionsTOPO.aan.root"); const double erMin=-0.3; const double erMax=0.3;

// Mike's files:
//TFile* fin = new TFile("Conversions.finalHiggs4.root"); const double erMin=-0.06; const double erMax=0.06;
TFile* fin = new TFile("Conversions.singlPhotn.root"); const double erMin=-0.06; const double erMax=0.06;

TTree* ntuple = (TTree*)(fin->Get("CollectionTree"));

//----------------------------------------------------------------------
// Looking at Etruth AND E histograms for three different eta-cuts and R-cuts
void z0(){
  std::cout << "Ethruth and E histograms for 3x3 cuts in R and in eta" << std::endl;
  // Canvas
  c1 = new TCanvas("c1","Ethruth and E histograms for 3x3 cuts in R and in eta",1280,768);
  c1->SetBorderMode(0);
  TCut Cut1("eConvEnergyInCalo1==eConvEnergyInCalo2");
  TCut Cut2("eConvEnergyInCalo1>-5");
  TCut Cut3_1("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=400");
  TCut Cut3_2("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))>400 && sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=800");
  TCut Cut3_3("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))>800 && sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=1000");
  TCut Cut3_a("eConvVerEta>0 && eConvVerEta<0.5");
  TCut Cut3_b("eConvVerEta>0.5 && eConvVerEta<1");
  TCut Cut3_c("eConvVerEta>1 && eConvVerEta<1.5");
  
  // Truth histos  
  ntuple->Draw("eConvGammaEnergy >> h1a(30,0,150e3",Cut1 && Cut2 && Cut3_1 && Cut3_a);
  ntuple->Draw("eConvGammaEnergy >> h1b(30,0,150e3)",Cut1 && Cut2 && Cut3_1 && Cut3_b);
  ntuple->Draw("eConvGammaEnergy >> h1c(30,0,150e3)",Cut1 && Cut2 && Cut3_1 && Cut3_c);

  ntuple->Draw("eConvGammaEnergy >> h2a(30,0,150e3)",Cut1 && Cut2 && Cut3_2 && Cut3_a);
  ntuple->Draw("eConvGammaEnergy >> h2b(30,0,150e3)",Cut1 && Cut2 && Cut3_2 && Cut3_b);
  ntuple->Draw("eConvGammaEnergy >> h2c(30,0,150e3)",Cut1 && Cut2 && Cut3_2 && Cut3_c);

  ntuple->Draw("eConvGammaEnergy >> h3a(30,0,150e3)",Cut1 && Cut2 && Cut3_3 && Cut3_a);
  ntuple->Draw("eConvGammaEnergy >> h3b(30,0,150e3)",Cut1 && Cut2 && Cut3_3 && Cut3_b);
  ntuple->Draw("eConvGammaEnergy >> h3c(30,0,150e3)",Cut1 && Cut2 && Cut3_3 && Cut3_c);

  TH1F *h1a = (TH1F*)gDirectory->Get("h1a");
  TH1F *h1b = (TH1F*)gDirectory->Get("h1b");
  TH1F *h1c = (TH1F*)gDirectory->Get("h1c");
  TH1F *h2a = (TH1F*)gDirectory->Get("h2a");
  TH1F *h2b = (TH1F*)gDirectory->Get("h2b");
  TH1F *h2c = (TH1F*)gDirectory->Get("h2c");
  TH1F *h3a = (TH1F*)gDirectory->Get("h3a");
  TH1F *h3b = (TH1F*)gDirectory->Get("h3b");
  TH1F *h3c = (TH1F*)gDirectory->Get("h3c");

  // Calo energy histos
  ntuple->Draw("eConvEnergyInCalo1 >> f1a(30,0,150e3",Cut1 && Cut2 && Cut3_1 && Cut3_a);
  ntuple->Draw("eConvEnergyInCalo1 >> f1b(30,0,150e3)",Cut1 && Cut2 && Cut3_1 && Cut3_b);
  ntuple->Draw("eConvEnergyInCalo1 >> f1c(30,0,150e3)",Cut1 && Cut2 && Cut3_1 && Cut3_c);

  ntuple->Draw("eConvEnergyInCalo1 >> f2a(30,0,150e3)",Cut1 && Cut2 && Cut3_2 && Cut3_a);
  ntuple->Draw("eConvEnergyInCalo1 >> f2b(30,0,150e3)",Cut1 && Cut2 && Cut3_2 && Cut3_b);
  ntuple->Draw("eConvEnergyInCalo1 >> f2c(30,0,150e3)",Cut1 && Cut2 && Cut3_2 && Cut3_c);

  ntuple->Draw("eConvEnergyInCalo1 >> f3a(30,0,150e3)",Cut1 && Cut2 && Cut3_3 && Cut3_a);
  ntuple->Draw("eConvEnergyInCalo1 >> f3b(30,0,150e3)",Cut1 && Cut2 && Cut3_3 && Cut3_b);
  ntuple->Draw("eConvEnergyInCalo1 >> f3c(30,0,150e3)",Cut1 && Cut2 && Cut3_3 && Cut3_c);

  TH1F *f1a = (TH1F*)gDirectory->Get("f1a");
  TH1F *f1b = (TH1F*)gDirectory->Get("f1b");
  TH1F *f1c = (TH1F*)gDirectory->Get("f1c");
  TH1F *f2a = (TH1F*)gDirectory->Get("f2a");
  TH1F *f2b = (TH1F*)gDirectory->Get("f2b");
  TH1F *f2c = (TH1F*)gDirectory->Get("f2c");
  TH1F *f3a = (TH1F*)gDirectory->Get("f3a");
  TH1F *f3b = (TH1F*)gDirectory->Get("f3b");
  TH1F *f3c = (TH1F*)gDirectory->Get("f3c");
  
  // 1 2 3
  // 4 5 6
  // 7 8 9
  c1->Clear();
  c1->Divide(3,3);

  // R<400
  c1->cd(1);
  h1a->SetFillColor(8);
  h1a->Draw();
  f1a->SetLineColor(kRed);
  f1a->Draw("histsame");
  h1a->SetTitle("R<400, 0<|#eta|<0.5");
  h1a->GetXaxis()->SetTitle("Energy");

  c1->cd(2);
  h1b->SetFillColor(8);
  h1b->Draw();
  f1b->SetLineColor(kRed);
  f1b->Draw("histsame");
  h1b->SetTitle("R<400, 0.5<|#eta|<1");
  h1b->GetXaxis()->SetTitle("Energy");

  c1->cd(3);
  h1c->SetFillColor(8);
  h1c->Draw();
  f1c->SetLineColor(kRed);
  f1c->Draw("histsame");
  h1c->SetTitle("R<400, 1<|#eta|<1.5");
  h1c->GetXaxis()->SetTitle("Energy");


  // 400<R<800
  c1->cd(4);
  h2a->SetFillColor(8);
  h2a->Draw();
  f2a->SetLineColor(kRed);
  f2a->Draw("histsame");
  h2a->SetTitle("400<R<800, 0<|#eta|<0.5");
  h2a->GetXaxis()->SetTitle("Energy");

  c1->cd(5);
  h2b->SetFillColor(8);
  h2b->Draw();
  f2b->SetLineColor(kRed);
  f2b->Draw("histsame");
  h2b->SetTitle("400<R<800, 0.5<|#eta|<1");
  h2b->GetXaxis()->SetTitle("Energy");

  c1->cd(6);
  h2c->SetFillColor(8);
  h2c->Draw();
  f2c->SetLineColor(kRed);
  f2c->Draw("histsame");
  h2c->SetTitle("400<R<800, 1<|#eta|<1.5");
  h2c->GetXaxis()->SetTitle("Energy");


  // 800<R<1000
  c1->cd(7);
  h3a->SetFillColor(8);
  h3a->Draw();
  f3a->SetLineColor(kRed);
  f3a->Draw("histsame");
  h3a->SetTitle("800<R<1000, 0<|#eta|<0.5");
  h3a->GetXaxis()->SetTitle("Energy");

  c1->cd(8);
  h3b->SetFillColor(8);
  h3b->Draw();
  f3b->SetLineColor(kRed);
  f3b->Draw("histsame");
  h3b->SetTitle("800<R<1000, 0.5<|#eta|<1");
  h3b->GetXaxis()->SetTitle("Energy");

  c1->cd(9);
  h3c->SetFillColor(8);
  h3c->Draw();
  f3c->SetLineColor(kRed);
  f3c->Draw("histsame");
  h3c->SetTitle("800<R<1000, 1<|#eta|<1.5");
  h3c->GetXaxis()->SetTitle("Energy");
  
  c1->Modified();
  c1->Update();
  printHisto(c1,"energy_eta_histos");
  return;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// Looking at <(Etruth-E)/Etruth> as a function of eta for three different R-cuts
void z1(){
  std::cout << "Energy response as a function of eta for three different cuts in R" << std::endl;
  // Canvas
  c2 = new TCanvas("c2","<(Et-E)/Et> as a function of eta for three different R-cuts",1000,700);
  c2->SetBorderMode(0);
  TCut Cut1("eConvEnergyInCalo1==eConvEnergyInCalo2");
  TCut Cut2("eConvEnergyInCalo1>-5");
  TCut Cut3_1("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=400");
  TCut Cut3_2("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))>400 && sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=800");
  TCut Cut3_3("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))>800 && sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=1000");

  const Double_t eta_step = 0.05, eta_max=2.35, eta_min=0;
  const Int_t n_bins = 52; // cannot dynamically calculate in ROOT macro! Only in CINT!
  Double_t etas[n_bins+1];
  Double_t means[3][n_bins+1]; // three r-regions
  Double_t entries[3][n_bins+1];
  int i=0;
  for (Double_t eta=eta_min; eta<eta_max; eta+=eta_step) {
    ostringstream cutStr;
    cutStr << "eConvVerEta>" << eta << " && eConvVerEta<" << (eta+eta_step);
    string cutS = cutStr.str();
    TCut Cut(cutS.c_str());
    ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h1(100,-1,1)",Cut1 && Cut2 && Cut && Cut3_1);
    ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h2(100,-1,1)",Cut1 && Cut2 && Cut && Cut3_2);
    ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h3(100,-1,1)",Cut1 && Cut2 && Cut && Cut3_3);
    TH1F *h1 = (TH1F*)gDirectory->Get("h1"); // this avoids segfaults due to recylcing h1!
    TH1F *h2 = (TH1F*)gDirectory->Get("h2"); // this avoids segfaults due to recylcing h1!
    TH1F *h3 = (TH1F*)gDirectory->Get("h3"); // this avoids segfaults due to recylcing h1!
    etas[i]=eta;
    means[0][i]=h1->GetMean(); means[1][i]=h2->GetMean(); means[2][i]=h3->GetMean();
    entries[0][i]=h1->GetEntries(); entries[1][i]=h2->GetEntries(); entries[2][i]=h3->GetEntries();
    i++;
  }

  // 1 2 3
  // 4 5 6
  c2->Clear();
  c2->Divide(3,2);

  c2->cd(1);
  TGraph *g1_1 = new TGraph(i,etas,means[0]);
  g1_1->SetTitle("R<400");
  g1_1->GetXaxis()->SetTitle("|#eta|");
  g1_1->GetYaxis()->SetTitle("<#frac{E_{t}-E_{c}}{E_{t}}>");
  g1_1->GetYaxis()->SetRangeUser(erMin,erMax);
  g1_1->Draw("AC*");

  c2->cd(2);
  TGraph *g2_1 = new TGraph(i,etas,means[1]);
  g2_1->SetTitle("400<R<800 ");
  g2_1->GetXaxis()->SetTitle("|#eta|");
  g2_1->GetYaxis()->SetTitle("<#frac{E_{t}-E_{c}}{E_{t}}>");
  g2_1->GetYaxis()->SetRangeUser(erMin,erMax);
  g2_1->Draw("AC*");

  c2->cd(3);
  TGraph *g3_1 = new TGraph(i,etas,means[2]);
  g3_1->SetTitle("800<R<1000 ");
  g3_1->GetXaxis()->SetTitle("|#eta|");
  g3_1->GetYaxis()->SetTitle("<#frac{E_{t}-E_{c}}{E_{t}}>");
  g3_1->GetYaxis()->SetRangeUser(erMin,erMax);
  g3_1->Draw("AC*");

  c2->cd(4);
  TGraph *g1_2 = new TGraph(i,etas,entries[0]);
  g1_2->SetTitle("R<400");
  g1_2->GetXaxis()->SetTitle("|#eta|");
  g1_2->GetYaxis()->SetTitle("# pts in avg");
  g1_2->Draw("AC*");

  c2->cd(5);
  TGraph *g2_2 = new TGraph(i,etas,entries[1]);
  g2_2->SetTitle("400<R<800 ");
  g2_2->GetXaxis()->SetTitle("|#eta|");
  g2_2->GetYaxis()->SetTitle("# pts in avg");
  g2_2->Draw("AC*");

  c2->cd(6);
  TGraph *g3_2 = new TGraph(i,etas,entries[2]);
  g3_2->SetTitle("800<R<1000 ");
  g3_2->GetXaxis()->SetTitle("|#eta|");
  g3_2->GetYaxis()->SetTitle("# pts in avg");
  g3_2->Draw("AC*");

  c2->Modified();
  c2->Update();
  printHisto(c2,"energy_eta");
  return;
}
//----------------------------------------------------------------------

//----------------------------------------------------------------------
// Looking at (Etruth-E)/Etruth HISTOS for 6 different eta regions
void z2(){
  std::cout << "Energy response HISTOS for 6 different cuts in eta" << std::endl;
  // Canvas
  c3 = new TCanvas("c3","Looking at (Et-E)/Et HISTOS for 6 different eta regions",1000,700);
  c3->SetBorderMode(0);
  TCut Cut1("eConvEnergyInCalo1==eConvEnergyInCalo2");
  TCut Cut2("eConvEnergyInCalo1>-5");
  TCut Cut3_a("eConvVerEta>0 && eConvVerEta<0.5");
  TCut Cut3_b("eConvVerEta>0.5 && eConvVerEta<1");
  TCut Cut3_c("eConvVerEta>1 && eConvVerEta<1.5");
  TCut Cut3_d("eConvVerEta>1.5 && eConvVerEta<2.0");
  TCut Cut3_e("eConvVerEta>2.0 && eConvVerEta<2.5");
  TCut Cut3_f("eConvVerEta>2.5 && eConvVerEta<3.0");

  ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h1(40,-2,2)",Cut1 && Cut2 && Cut3_a);
  ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h2(40,-2,2)",Cut1 && Cut2 && Cut3_b);
  ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h3(40,-2,2)",Cut1 && Cut2 && Cut3_c);
  ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h4(40,-2,2)",Cut1 && Cut2 && Cut3_d);
  ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h5(40,-2,2)",Cut1 && Cut2 && Cut3_e);
  ntuple->Draw("(eConvGammaEnergy-eConvEnergyInCalo1)/(eConvGammaEnergy) >> h6(40,-2,2)",Cut1 && Cut2 && Cut3_f);
  TH1F *h1 = (TH1F*)gDirectory->Get("h1"); // this avoids segfaults due to recylcing h1!
  h1->SetTitle("0.0<|#eta|<0.5");
  TH1F *h2 = (TH1F*)gDirectory->Get("h2"); // this avoids segfaults due to recylcing h1!
  h2->SetTitle("0.5<|#eta|<1.0");
  TH1F *h3 = (TH1F*)gDirectory->Get("h3"); // this avoids segfaults due to recylcing h1!
  h3->SetTitle("1.0<|#eta|<1.5");
  TH1F *h4 = (TH1F*)gDirectory->Get("h4"); // this avoids segfaults due to recylcing h1!
  h4->SetTitle("1.5<|#eta|<2.0");
  TH1F *h5 = (TH1F*)gDirectory->Get("h5"); // this avoids segfaults due to recylcing h1!
  h5->SetTitle("2.0<|#eta|<2.5");
  TH1F *h6 = (TH1F*)gDirectory->Get("h6"); // this avoids segfaults due to recylcing h1!
  h6->SetTitle("2.5<|#eta|<3.0");
  h1->GetXaxis()->SetTitle("(Et-E)/Et");
  h2->GetXaxis()->SetTitle("(Et-E)/Et");
  h3->GetXaxis()->SetTitle("(Et-E)/Et");
  h4->GetXaxis()->SetTitle("(Et-E)/Et");
  h5->GetXaxis()->SetTitle("(Et-E)/Et");
  h6->GetXaxis()->SetTitle("(Et-E)/Et");

  // 1 2 3
  // 4 5 6
  c3->Clear();
  c3->Divide(3,2);

  c3->cd(1);
  gPad->SetLogy();
  h1->Draw();

  c3->cd(2);
  gPad->SetLogy();
  h2->Draw();

  c3->cd(3);
  gPad->SetLogy();
  h3->Draw();

  c3->cd(4);
  gPad->SetLogy();
  h4->Draw();

  c3->cd(5);
  gPad->SetLogy();
  h5->Draw();

  c3->cd(6);
  gPad->SetLogy();
  h6->Draw();

  c3->Modified();
  c3->Update();
  printHisto(c3,"energy_eta_diff_histo");
  return;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// Looking at E/Et HISTOS for 6 different eta regions
void z3(){
  std::cout << "E/Et HISTOS for 6 different cuts in eta" << std::endl;
  // Canvas
  c3a = new TCanvas("c3a","Looking at E/Et HISTOS for 6 different eta regions",1000,700);
  c3a->SetBorderMode(0);
  c3a->cd();
  TCut Cut1("eConvEnergyInCalo1==eConvEnergyInCalo2");
  TCut Cut2("eConvEnergyInCalo1>-5");
  TCut Cut3_a("eConvVerEta>0 && eConvVerEta<0.5");
  TCut Cut3_b("eConvVerEta>0.5 && eConvVerEta<1");
  TCut Cut3_c("eConvVerEta>1 && eConvVerEta<1.5");
  TCut Cut3_d("eConvVerEta>1.5 && eConvVerEta<2.0");
  TCut Cut3_e("eConvVerEta>2.0 && eConvVerEta<2.5");
  TCut Cut3_f("eConvVerEta>2.5 && eConvVerEta<3.0");

  ntuple->Draw("(eConvEnergyInCalo1)/(eConvGammaEnergy) >> h1A(40,0,4)",Cut1 && Cut2 && Cut3_a);
  ntuple->Draw("(eConvEnergyInCalo1)/(eConvGammaEnergy) >> h2A(40,0,4)",Cut1 && Cut2 && Cut3_b);
  ntuple->Draw("(eConvEnergyInCalo1)/(eConvGammaEnergy) >> h3A(40,0,4)",Cut1 && Cut2 && Cut3_c);
  ntuple->Draw("(eConvEnergyInCalo1)/(eConvGammaEnergy) >> h4a(40,0,4)",Cut1 && Cut2 && Cut3_d);
  ntuple->Draw("(eConvEnergyInCalo1)/(eConvGammaEnergy) >> h5a(40,0,4)",Cut1 && Cut2 && Cut3_e);
  ntuple->Draw("(eConvEnergyInCalo1)/(eConvGammaEnergy) >> h6a(40,0,4)",Cut1 && Cut2 && Cut3_f);
  TH1F *h1A = (TH1F*)gDirectory->Get("h1A"); // this avoids segfaults due to recylcing h1!
  h1A->SetTitle("0.0<|#eta|<0.5");
  TH1F *h2A = (TH1F*)gDirectory->Get("h2A"); // this avoids segfaults due to recylcing h1!
  h2A->SetTitle("0.5<|#eta|<1.0");
  TH1F *h3A = (TH1F*)gDirectory->Get("h3A"); // this avoids segfaults due to recylcing h1!
  h3A->SetTitle("1.0<|#eta|<1.5");
  TH1F *h4a = (TH1F*)gDirectory->Get("h4a"); // this avoids segfaults due to recylcing h1!
  h4a->SetTitle("1.5<|#eta|<2.0");
  TH1F *h5a = (TH1F*)gDirectory->Get("h5a"); // this avoids segfaults due to recylcing h1!
  h5a->SetTitle("2.0<|#eta|<2.5");
  TH1F *h6a = (TH1F*)gDirectory->Get("h6a"); // this avoids segfaults due to recylcing h1!
  h6a->SetTitle("2.5<|#eta|<3.0");
  h1A->GetXaxis()->SetTitle("E/Et");
  h2A->GetXaxis()->SetTitle("E/Et");
  h3A->GetXaxis()->SetTitle("E/Et");
  h4a->GetXaxis()->SetTitle("E/Et");
  h5a->GetXaxis()->SetTitle("E/Et");
  h6a->GetXaxis()->SetTitle("E/Et");

  // 1 2 3
  // 4 5 6
  c3a->Clear();
  c3a->Divide(3,2);

  c3a->cd(1);
  gPad->SetLogy();
  h1A->Draw();

  c3a->cd(2);
  gPad->SetLogy();
  h2A->Draw();

  c3a->cd(3);
  gPad->SetLogy();
  h3A->Draw();

  c3a->cd(4);
  gPad->SetLogy();
  h4a->Draw();

  c3a->cd(5);
  gPad->SetLogy();
  h5a->Draw();

  c3a->cd(6);
  gPad->SetLogy();
  h6a->Draw();

  c3a->Modified();
  c3a->Update();
  printHisto(c3a,"energy_eta_ratio_histo");
  return;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// Conversion Rate as a function of Eta for three regions in R
void z4() {
  std::cout << "Conversion rate as a function of eta for 3 different cuts in R" << std::endl;
  // Canvas
  c4 = new TCanvas("c4","Conversion Rate as a function of Eta for three regions in R",400,700);
  c4->SetBorderMode(0);
  TCut Cut3_1("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=500");
  TCut Cut3_2("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2)) && sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=650");
  TCut Cut3_3("sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))>650 && sqrt(pow(eConvVerX,2)+pow(eConvVerY,2))<=1000");

  TH1F *photCnvRate[3]; //  for three regions in R!
  ntuple->Draw("abs(TruthPhotonEta) >> ph1a(25,0,4)","abs(TruthPhotonEta)<10");
  TH1F *ph1a = (TH1F*)gDirectory->Get("ph1a");
  ntuple->Draw("abs(eConvVerEta) >> ph1b(25,0,4)","abs(TruthPhotonEta)<10 && NConversions==1" && Cut3_1);
  TH1F *ph1b = (TH1F*)gDirectory->Get("ph1b");

  ntuple->Draw("abs(TruthPhotonEta) >> ph2a(25,0,4)","abs(TruthPhotonEta)<10");
  TH1F *ph2a = (TH1F*)gDirectory->Get("ph2a");
  ntuple->Draw("abs(eConvVerEta) >> ph2b(25,0,4)","abs(TruthPhotonEta)<10 && NConversions==1"&& Cut3_2);
  TH1F *ph2b = (TH1F*)gDirectory->Get("ph2b");

  ntuple->Draw("abs(TruthPhotonEta) >> ph3a(25,0,4)","abs(TruthPhotonEta)<10");
  TH1F *ph3a = (TH1F*)gDirectory->Get("ph3a");
  ntuple->Draw("abs(eConvVerEta) >> ph3b(25,0,4)","abs(TruthPhotonEta)<10 && NConversions==1"&& Cut3_3);
  TH1F *ph3b = (TH1F*)gDirectory->Get("ph3b");

  c4->Clear();
  c4->Divide(1,3);
  c4->cd(1); photCnvRate[0] = new TH1F("photCnvRate_0","Photon Conversion Rate",25,0,4);
  c4->cd(2); photCnvRate[1] = new TH1F("photCnvRate_1","Photon Conversion Rate",25,0,4);
  c4->cd(3); photCnvRate[2] = new TH1F("photCnvRate_2","Photon Conversion Rate",25,0,4);

  for(int i=1; i<25; i++){
    float totalentries = ph1a->GetBinContent(i);
    float conventries = ph1b->GetBinContent(i);
    if(totalentries!=0)
      photCnvRate[0]->SetBinContent(i,conventries/totalentries);
    else
      photCnvRate[0]->SetBinContent(i,0);
  }
 
  for(int i=1; i<25; i++){
    float totalentries = ph2a->GetBinContent(i);
    float conventries = ph2b->GetBinContent(i);
    if(totalentries!=0)
      photCnvRate[1]->SetBinContent(i,conventries/totalentries);
    else
      photCnvRate[1]->SetBinContent(i,0);
  }
 
  for(int i=1; i<25; i++){
    float totalentries = ph3a->GetBinContent(i);
    float conventries = ph3b->GetBinContent(i);
    if(totalentries!=0)
      photCnvRate[2]->SetBinContent(i,conventries/totalentries);
    else
      photCnvRate[2]->SetBinContent(i,0);
  }

  photCnvRate[0]->SetTitle("R<500");
  photCnvRate[1]->SetTitle("500>R>650");
  photCnvRate[2]->SetTitle("650>R>1000");

  for (int i=0; i<3; i++) {
    //    photCnvRate[i]->GetXaxis()->SetRangeUser(0,3);
    //    photCnvRate[i]->GetYaxis()->SetRangeUser(.4,.8);
    photCnvRate[i]->GetXaxis()->SetTitle("|#eta|");
    photCnvRate[i]->GetYaxis()->SetTitle("rate");
    photCnvRate[i]->SetStats(kFALSE);
    photCnvRate[i]->SetFillColor(40+i); //42
    c4->cd(i+1);
    photCnvRate[i]->Draw();
  }
  
  c4->Modified();
  c4->Update();
  printHisto(c4,"cnvRate");
  return;
}
//----------------------------------------------------------------------


//----------------------------------------------------------------------
// Conversion Rate as a function of Eta - overall!
void z5() {
  std::cout << "Conversion rate as a function of Eta, summed over all R" << std::endl;
  // Canvas
  c5 = new TCanvas("c5","Conversion rate as a function of Eta, summed over all R",800,600);
  c5->SetBorderMode(0);

  TH1F *photCnvRate;
  ntuple->Draw("abs(TruthPhotonEta) >> ph1a(25,0,4)","abs(TruthPhotonEta)<10");
  TH1F *ph1a = (TH1F*)gDirectory->Get("ph1a");
  ntuple->Draw("abs(eConvVerEta) >> ph1b(25,0,4)"," NConversions>0");
  TH1F *ph1b = (TH1F*)gDirectory->Get("ph1b");

  c5->Clear();
  c5->cd();
  photCnvRate = new TH1F("photCnvRate","Photon Conversion Rate",25,0,4);

  for(int i=1; i<25; i++){
    float totalentries = ph1a->GetBinContent(i);
    float conventries = ph1b->GetBinContent(i);
    if(totalentries!=0)
      photCnvRate->SetBinContent(i,conventries/totalentries);
    else
      photCnvRate->SetBinContent(i,0);
  }

  photCnvRate->SetTitle("Conv. rate");

  //    photCnvRate->GetXaxis()->SetRangeUser(0,3);
  //    photCnvRate->GetYaxis()->SetRangeUser(.4,.8);
  photCnvRate->GetXaxis()->SetTitle("|#eta|");
  photCnvRate->GetYaxis()->SetTitle("rate");
  photCnvRate->SetStats(kFALSE);
  photCnvRate->SetFillColor(40);
  photCnvRate->Draw();
    
  c5->Modified();
  c5->Update();
  printHisto(c5,"cnvRate_all_R");
  return;
}
//----------------------------------------------------------------------

void zhelp() {
  std::cout << "z0(): Et and E histograms for 3 cuts in R and 3 cuts in eta" << std::endl;
  std::cout << "z1(): <(Et-E)/Et> as a function of eta for three different cuts in R. VERY SLOW!" << std::endl;
  std::cout << "z2(): (Et-E)/Et HISTOS for 6 different cuts in eta" << std::endl;
  std::cout << "z3(): E/Et HISTOS for 6 different cuts in eta" << std::endl;
  std::cout << "z4(): Conversion rate as a function of eta for 3 different cuts in R" << std::endl;
  std::cout << "z5(): Conversion rate as a function of Eta, summed over all R" << std::endl;
  std::cout << "zz(): Make all plots. VERY SLOW!" << std::endl;
  return;
}

void zz() {
  z0();
  z1();
  z2();
  z3();
  z4();
  z5();
  return;
}

void printHisto(TCanvas* c1, string name){
  system("mkdir -p plots-png");
  //  system("mkdir -p plots-eps");
  string pngdir("plots-png/");
  string epsdir("plots-eps/");
  string dot(1,'.');
  string eps("eps");
  string png("png");
  //  c1->Print((epsdir+name+dot+eps).c_str(),eps.c_str());
  c1->Print((pngdir+name+dot+png).c_str(),png.c_str());
  return;
}