#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TLorentzVector.h"
#include <iostream>

using namespace std;

void read()
{
  TChain* fin = new TChain("h300");
  fin->Add("*.root/h300");


// This works as .x read.cc or .x read.cc+

  float px3,py3,pz3,E3;
  float px4,py4,pz4,E4;
  float px5,py5,pz5,E5;
  float px6,py6,pz6,E6;
  float px7,py7,pz7,E7;
  float wt,gg,gq,qq,qqb;

  h300->SetBranchAddress("px3",&px3);
  h300->SetBranchAddress("py3",&py3);
  h300->SetBranchAddress("pz3",&pz3);
  h300->SetBranchAddress("E3" ,&E3);
  h300->SetBranchAddress("px4",&px4);
  h300->SetBranchAddress("py4",&py4);
  h300->SetBranchAddress("pz4",&pz4);
  h300->SetBranchAddress("E4" ,&E4);
  h300->SetBranchAddress("px5",&px5);
  h300->SetBranchAddress("py5",&py5);
  h300->SetBranchAddress("pz5",&pz5);
  h300->SetBranchAddress("E5" ,&E5);
  h300->SetBranchAddress("px6",&px6);
  h300->SetBranchAddress("py6",&py6);
  h300->SetBranchAddress("pz6",&pz6);
  h300->SetBranchAddress("E6" ,&E6);
  h300->SetBranchAddress("px7",&px7);
  h300->SetBranchAddress("py7",&py7);
  h300->SetBranchAddress("pz7",&pz7);
  h300->SetBranchAddress("E7" ,&E7);
  h300->SetBranchAddress("wt" ,&wt);
  h300->SetBranchAddress("gg" ,&gg);
  h300->SetBranchAddress("gq" ,&gq);
  h300->SetBranchAddress("qq" ,&qq);
  h300->SetBranchAddress("qqb",&qqb);

  TFile* fout = new TFile("mcfm_histograms.root","recreate");

  TH1F* h_mtau_tau = new TH1F("h_mtau_tau","tau-tau transverse mass",100,0,500);
  TH1F* h_y_j1   = new TH1F("h_y_j1","1st jet rapidity",40,-5,5);
  TH1F* h_y_j2   = new TH1F("h_y_j2","2nd jet rapidity",40,-5,5);
  TH1F* h_y_j3   = new TH1F("h_y_j3","3rd jet rapidity",40,-5,5);

  TH1F* h_pt_j1  = new TH1F("h_pt_j1","1st jet pt",40,0,200);
  TH1F* h_pt_j2  = new TH1F("h_pt_j2","2nd jet pt",40,0,200);
  TH1F* h_pt_j3  = new TH1F("h_pt_j3","3rd jet pt",40,0,200);

  TH1F* h_DeltaR_j1j2 = new TH1F("h_DeltaR_j1_j2","deltaR j1 j2",50,0,10);

  TH1F* h_Delta_y1y2 = new TH1F("h_Delta_y1y2","delta y1-y2",100,0,10);

  TH1F* h_DEta_j1j2 = new TH1F("h_DEta_j1j2","deta eta1-eta2",100,0,10);

  TH1F* h_DeltaR_j1e = new TH1F("h_DeltaR_j1_e","deltaR j1 e",50,0,10);
  TH1F* h_DeltaR_j2e = new TH1F("h_DeltaR_j2_e","deltaR j2 e",50,0,10);
  TH1F* h_DeltaR_j1j3 = new TH1F("h_DeltaR_j1_j3","deltaR j1 j3",50,0.01,10.01);
  //  TH1F* h_DEta_j1j2 = new TH1F("h_DEta_j1j2","delta eta j1 j2",50,0,5);
  TH1F* h_pt_j3_2 = new TH1F("h_pt_j3_2","pt j3 Dy > 2",40,0,200);
  TH1F* h_pt_j3_4 = new TH1F("h_pt_j3_4","pt j3 Dy > 4",40,0,200);
  TH1F* h_pt_j3_6 = new TH1F("h_pt_j3_6","pt j3 Dy > 6",40,0,200);
  TH1F* h_eta3s_2 = new TH1F("h_eta3s_2","eta3 Deta > 2",40,-5,5);
  TH1F* h_eta3s_4 = new TH1F("h_eta3s_4","eta3 Deta > 4",40,-5,5);
  TH1F* h_eta3s_6 = new TH1F("h_eta3s_6","eta3 Deta > 6",40,-5,5);

  int nEntries = (int)h300->GetEntries();
  cout << "nEntries = " << nEntries << endl;

  for(int iEntry = 0; iEntry < nEntries; iEntry++){

    if(!(iEntry%50000))
      cout << "... reading event " << iEntry << endl;
    h300->GetEntry(iEntry);

    TLorentzVector ve(px3,py3,pz3,E3);   // electron 4-vector
    TLorentzVector vnu(px4,py4,pz4,E4);  // neutrino 4-vector
    TLorentzVector vj1(px5,py5,pz5,E5);  // 1st jet
    TLorentzVector vj2(px6,py6,pz6,E6);  // 2nd jet
    TLorentzVector vj3(px7,py7,pz7,E7);  // 3rd jet
	
    float mtau_tau = (ve + vnu).M();      // e-nu transverse mass
    
    float y_e = ve.Rapidity();           // electron rapidity
    float eta_e = ve.PseudoRapidity();   // electron pseudorapidity
    float phi_e = ve.Phi();              // electron phi
    float pt_e = ve.Perp();              // electron pT

    float y_nu = vnu.Rapidity();         // neutrino rapidity
    float eta_nu = vnu.PseudoRapidity();  // neutrino pseudorapidity
    float phi_nu = vnu.Phi();            // neutrino phi
    float pt_nu = vnu.Perp();            // neutrino pT

    float y_j1   = vj1.Rapidity();       // 1st jet rapidity
    float eta_j1 = vj1.PseudoRapidity(); // 1st jet pseudorapidity
    float phi_j1 = vj1.Phi();            // 1st jet phi
    float pt_j1 = vj1.Perp();            // 1st  jet  pT

    float y_j2   = vj2.Rapidity();       // 2nd jet ...
    float eta_j2 = vj2.PseudoRapidity();
    float phi_j2 = vj2.Phi();
    float pt_j2 = vj2.Perp();            // 2nd jet pT

    float DEta_j1j2 = fabs(eta_j1-eta_j2);

    float DeltaR_j1j2 = vj1.DeltaR(vj2);
    float Delta_y1y2 = fabs(y_j1-y_j2);

    // Note: This is DeltaR = (DeltaEta*DeltaEta + DeltaPhi*DeltaPhi). PseudoRapidity!!!

        if(pt_j2 < 15)
        cout << "pt_j2 = " << pt_j2 << ", "
        << "px6 = " << px6 << ", py6 = " << py6 << ", pz6 " << pz6 << ", "
        << "E6 = " << E6 << ", wt = " << wt  << ", pt_j1 = " << pt_j1 << ", "
        << "E5 = " << E5;
    
    float y_j3 = -10, eta_j3 = 10, phi_j3 = -10, DeltaR_j1j3 = -10, DeltaR_j2j3 = -10;
    float pt_j3 = 1, eta3s = -10;
    if(vj3.E() != 0){
    y_j3   = vj3.Rapidity();       // 3rd jet ...
    eta_j3 = vj3.PseudoRapidity();
    phi_j3 = vj3.Phi();
    pt_j3 = vj3.Perp();

    eta3s = eta_j3 - (eta_j1 + eta_j2)/2;

    DeltaR_j1j3 = vj1.DeltaR(vj3);
    DeltaR_j2j3 = vj2.DeltaR(vj3);

      //      if(iEntry < 10)
      //	cout << ", y_j3 = " << y_j3 << ", eta_j3 = " << eta_j3 << ", phi_j3 = " << phi_j3 << ", "
      //	     << ", DeltaR_j1j3 = " << DeltaR_j1j3 << ", DeltaR_j2j3 = " << DeltaR_j2j3;
          }

    //   if(iEntry < 10)
    //     cout << endl;

    if(pt_j1 > 50){
    h_mtau_tau->Fill(mtau_tau,wt);
    h_y_j1->Fill(y_j1,wt);
    h_y_j2->Fill(y_j2,wt);
    h_y_j3->Fill(y_j3,wt);
    
    h_Delta_y1y2->Fill(Delta_y1y2,wt);
    h_DEta_j1j2->Fill(DEta_j1j2,wt);

    h_pt_j1->Fill(pt_j1,wt);
    h_pt_j2->Fill(pt_j2,wt);
    h_pt_j3->Fill(pt_j3,wt);
    if(DEta_j1j2 > 2){
      h_pt_j3_2->Fill(pt_j3,wt); 
      h_eta3s_2->Fill(eta3s,wt);}
    if(DEta_j1j2 > 4){
      h_pt_j3_4->Fill(pt_j3,wt);
      h_eta3s_4->Fill(eta3s,wt);}
    if(DEta_j1j2 > 6){
      h_pt_j3_6->Fill(pt_j3,wt); 
      h_eta3s_6->Fill(eta3s,wt);}
    }
  }

  fout->Write();
  fout->Close();
}