#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TLorentzVector.h"
#include <iostream>
#include <math.h>
using namespace std;

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


// This works as .x read.cc or .x read.cc+
// adapted from routine  for WW fusion

  float wt_ALL,wt_gq,wt_gg,wt_qq,wt_qqb;
  float PDF01,PDF02,PDF03,PDF04,PDF05,PDF06,PDF07,PDF08;
  float PDF09,PDF10,PDF11,PDF12,PDF13,PDF14,PDF15,PDF16;
  float PDF17,PDF18,PDF19,PDF20,PDF21,PDF22,PDF23,PDF24;
  float PDF25,PDF26,PDF27,PDF28,PDF29,PDF30,PDF31,PDF32;
  float PDF33,PDF34,PDF35,PDF36,PDF37,PDF38,PDF39,PDF40;
  float PDF41,PDF42,PDF43,PDF44;

  float px3,px4,px5,px6,px7;
  float py3,py4,py5,py6,py7;
  float pz3,pz4,pz5,pz6,pz7;
  float E_3,E_4,E_5,E_6,E_7;

  h10->SetBranchAddress("px3",&px3);
  h10->SetBranchAddress("px4",&px4);
  h10->SetBranchAddress("px5",&px5);
  h10->SetBranchAddress("px6",&px6);
  //  h10->SetBranchAddress("px7",&px7);
  h10->SetBranchAddress("py3",&py3);
  h10->SetBranchAddress("py4",&py4);
  h10->SetBranchAddress("py5",&py5);
  h10->SetBranchAddress("py6",&py6);
  //  h10->SetBranchAddress("py7",&py7);
  h10->SetBranchAddress("pz3",&pz3);
  h10->SetBranchAddress("pz4",&pz4);
  h10->SetBranchAddress("pz5",&pz5);
  h10->SetBranchAddress("pz6",&pz6);
  //  h10->SetBranchAddress("pz7",&pz7);
  h10->SetBranchAddress("E_3",&E_3);
  h10->SetBranchAddress("E_4",&E_4);
  h10->SetBranchAddress("E_5",&E_5);
  h10->SetBranchAddress("E_6",&E_6);
  //  h10->SetBranchAddress("E_7",&E_7);
  h10->SetBranchAddress("wt_ALL",&wt_ALL);
  h10->SetBranchAddress("wt_gg",&wt_gg);
  h10->SetBranchAddress("wt_gq",&wt_gq);
  h10->SetBranchAddress("wt_qq",&wt_qq);
  h10->SetBranchAddress("wt_qqb",&wt_qqb);
  h10->SetBranchAddress("PDF01",&PDF01);
  h10->SetBranchAddress("PDF02",&PDF02);
  h10->SetBranchAddress("PDF03",&PDF03);
  h10->SetBranchAddress("PDF04",&PDF04);
  h10->SetBranchAddress("PDF05",&PDF05);
  h10->SetBranchAddress("PDF06",&PDF06);
  h10->SetBranchAddress("PDF07",&PDF07);
  h10->SetBranchAddress("PDF08",&PDF08);
  h10->SetBranchAddress("PDF09",&PDF09);
  h10->SetBranchAddress("PDF10",&PDF10);
  h10->SetBranchAddress("PDF11",&PDF11);
  h10->SetBranchAddress("PDF12",&PDF12);
  h10->SetBranchAddress("PDF13",&PDF13);
  h10->SetBranchAddress("PDF14",&PDF14);
  h10->SetBranchAddress("PDF15",&PDF15);
  h10->SetBranchAddress("PDF16",&PDF16);
  h10->SetBranchAddress("PDF17",&PDF17);
  h10->SetBranchAddress("PDF18",&PDF18);
  h10->SetBranchAddress("PDF19",&PDF19);
  h10->SetBranchAddress("PDF20",&PDF20);
  h10->SetBranchAddress("PDF21",&PDF21);
  h10->SetBranchAddress("PDF22",&PDF22);
  h10->SetBranchAddress("PDF23",&PDF23);
  h10->SetBranchAddress("PDF24",&PDF24);
  h10->SetBranchAddress("PDF25",&PDF25);
  h10->SetBranchAddress("PDF26",&PDF26);
  h10->SetBranchAddress("PDF27",&PDF27);
  h10->SetBranchAddress("PDF28",&PDF28); 
  h10->SetBranchAddress("PDF29",&PDF29);
  h10->SetBranchAddress("PDF30",&PDF30);
  h10->SetBranchAddress("PDF31",&PDF31);
  h10->SetBranchAddress("PDF32",&PDF32);
  h10->SetBranchAddress("PDF33",&PDF33);
  h10->SetBranchAddress("PDF34",&PDF34);
  h10->SetBranchAddress("PDF35",&PDF35);
  h10->SetBranchAddress("PDF36",&PDF36);
  h10->SetBranchAddress("PDF37",&PDF37);
  h10->SetBranchAddress("PDF38",&PDF38);
  h10->SetBranchAddress("PDF39",&PDF39);
  h10->SetBranchAddress("PDF40",&PDF40);
  h10->SetBranchAddress("PDF41",&PDF41);
  h10->SetBranchAddress("PDF42",&PDF42);
  h10->SetBranchAddress("PDF43",&PDF43);
  h10->SetBranchAddress("PDF44",&PDF44);

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

  TH1F* h_m_h = new TH1F("h_m_h","higgs mass",100,0,500);
  TH1F* h_y_tau1   = new TH1F("h_y_tau1","tau1 rapidity",50,-5,5);
  TH1F* h_y_tau2   = new TH1F("h_y_tau2","tau2 rapidity",50,-5,5);
  TH1F* h_y_h = new TH1F("h_y_h","h rapidity",20,-5,5);

  TH1F* h_pt_tau1 = new TH1F("h_pt_tau1","tau1 pT",50,0,500);
  TH1F* h_pt_tau2 = new TH1F("h_pt_tau2","tau2 pT",50,0,500);
  TH1F* h_pt_h = new TH1F("h_pt_h","higgs pT",50,0,500);

  TH1F* h_pt_j1 = new TH1F("h_pt_j1","jet 1 pT",50,0,500);
  TH1F* h_et_j1 = new TH1F("h_et_j1","jet 1 ET",50,0,500);
 
  TH1F* h_y_j1 = new TH1F("h_y_j1","y jet 1",50,-5,5);

  TH1F* h_pt_j2 = new TH1F("h_pt_j2","jet 2 pT",50,0,500);

  TH1F* h_et_j2 = new TH1F("h_et_j2","jet 2 ET",50,0,500);

  TH1F* h_y_j2 = new TH1F("h_y_j2","y jet 2",50,-5,5);
  
  TH1F* h_y_dif = new TH1F("h_y_dif","y dif between two tag jets",20,0,10);
  TH1F* h_y_dif_20 = new TH1F("h_y_dif_20","y dif between two tag jets 20",20,0,10);
  TH1F* h_y_dif_30 = new TH1F("h_y_dif_30","y dif between two tag jets 30",20,0,10);
  TH1F* h_y_dif_40 = new TH1F("h_y_dif_40","y dif between two tag jets 40",20,0,10);

  TH1F* h_pt_j3 = new TH1F("h_pt_j3","jet 3 pT",50,0,500);
  TH1F* h_y_j3 = new TH1F("y_j3","jet 3 y",50,0,500);

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

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

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

    TLorentzVector vtau1(px3,py3,pz3,E_3);   // neutrino 4-vector
    TLorentzVector vtau2(px4,py4,pz4,E_4);  // electron 4-vector
    TLorentzVector vj1(px5,py5,pz5,E_5);  // 1st jet
    TLorentzVector vj2(px6,py6,pz6,E_6);  // 2nd jet
    //    TLorentzVector vj3(px7,py7,pz7,E_7);  // 3rd jet

    TLorentzVector vh = vtau1 + vtau2;      // higgs 4-vector
	
    float m_h = (vtau1 + vtau2).M();      // higgs  mass
    
    float m_j1 = vj1.M();            // mass of jet 1

    float y_tau1 = vtau1.Rapidity();           // y_tau1 rapidity
//    float eta_e = ve.PseudoRapidity();   //  pseudorapidity
    float phi_tau1 = vtau1.Phi();              // tau1 phi
    float pt_tau1 = vtau1.Perp();              // tau1 pT

    float y_tau2 = vtau2.Rapidity();         // tau2 rapidity
    //    float eta_nu = vnu.PseudoRapidity();  // neutrino pseudorapidity
    float phi_tau2 = vtau2.Phi();            // tau2 phi
    float pt_tau2 = vtau2.Perp();            // tau2 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 et_j1 = sqrt(m_j1*m_j1 + pt_j1*pt_j1);

    float y_j2   = vj2.Rapidity();       // 2nd jet rapidity
//    float eta_j2 = vj2.PseudoRapidity(); // 2nd jet pseudorapidity
    float phi_j2 = vj2.Phi();            // 2nd jet phi
    float pt_j2 = vj2.Perp();            // 2nd  jet  pT
	
	float m_j2 = vj2.M();            // mass of jet 2
	float et_j2 = sqrt(m_j2*m_j2 + pt_j2*pt_j2);

	//    float y_j3   = vj3.Rapidity();       // 3rd jet rapidity
	//    float eta_j3 = vj3.PseudoRapidity(); // 2nd jet pseudorapidity
	//    float phi_j3 = vj3.Phi();            // 3rd jet phi
	//    float pt_j3 = vj3.Perp();            // 3rd  jet  pT



    float y_h   = vh.Rapidity();       // W rapidity
    //    float eta_w = vw.PseudoRapidity(); // W pseudorapidity
    float phi_h = vh.Phi();            // W phi
    float pt_h = vh.Perp();            // W pT

//    float DeltaR_j1e = vj1.DeltaR(ve);
//    float DeltaR_j2e = vj2.DeltaR(ve);
//    float DeltaR_j3e = 10;
	//	if(vj3.E()!0){
	//
	//      float DeltaR_j3e = vj3.DeltaR(ve);
	//    }

    h_m_h->Fill(m_h,wt_ALL);
    h_y_tau1->Fill(y_tau1,wt_ALL);
    h_y_tau2->Fill(y_tau2,wt_ALL);
    h_y_h->Fill(y_h,wt_ALL);
    //   h_eta_w->Fill(eta_w,wgt);
//    h_eta_e->Fill(eta_e,wt_ALL);

    h_pt_tau1->Fill(pt_tau1,wt_ALL);
    h_pt_tau2->Fill(pt_tau2,wt_ALL);
    h_pt_h->Fill(pt_h,wt_ALL);

    h_y_j1->Fill(y_j1,wt_ALL);
//    h_eta_j1->Fill(eta_j1,wt_ALL);
    h_pt_j1->Fill(pt_j1,wt_ALL);
    h_et_j1->Fill(et_j1,wt_ALL);
	
    h_y_j2->Fill(y_j2,wt_ALL);
//    h_eta_j2->Fill(eta_j2,wt_ALL);
    h_pt_j2->Fill(pt_j2,wt_ALL);
    h_et_j2->Fill(et_j2,wt_ALL);

    float y_dif = fabs(y_j1-y_j2);
//    float eta_dif = fabs(eta_j1-eta_j2);

    h_y_dif->Fill(y_dif,wt_ALL);
//    h_eta_dif->Fill(eta_dif,wt_ALL);

//    h_pt_j3->Fill(pt_j3,wt_ALL);
//    h_y_j3->Fill(y_j3,wt_ALL);

    if(pt_j1 > 20){
      if(pt_j2 > 20){
	  float y_dif_20 = fabs(y_j1-y_j2);
	  h_y_dif_20->Fill(y_dif_20,wt_ALL);
	  if(pt_j1 > 30){
	    if(pt_j2 > 30){
		float y_dif_30 = y_dif_20;
		h_y_dif_30->Fill(y_dif_30,wt_ALL);
		if(pt_j1 > 40){
		  if(pt_j2 > 40){
		      float y_dif_40 = y_dif_20;
		      h_y_dif_40->Fill(y_dif_40,wt_ALL);
		  }
		}
	    }
	  }
      }
    }
  }
 
  fout->Write();
  fout->Close();
}