#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.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_gg,wt_gq,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;
  float py3,py4,py5;
  float pz3,pz4,pz5;
  float E_3,E_4,E_5;

  h10->SetBranchAddress("E_3" ,&E_3);
  h10->SetBranchAddress("E_4" ,&E_4);
  h10->SetBranchAddress("E_5" ,&E_5);

//  t3->SetBranchAddress("nparticle", &nparticle);
  h10->SetBranchAddress("px3",&px3);
  h10->SetBranchAddress("px4",&px4);
  h10->SetBranchAddress("px5",&px5);

  h10->SetBranchAddress("py3",&py3);
  h10->SetBranchAddress("py4",&py4);
  h10->SetBranchAddress("py5",&py5);

  h10->SetBranchAddress("pz3",&pz3);
  h10->SetBranchAddress("pz4",&pz4);
  h10->SetBranchAddress("pz5",&pz5);
 
  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_y_t = new TH1F("h_y_t","top rapidity",20,-5,5);
  TH1F* h_y_tb   = new TH1F("h_y_tb","tbar rapidity",20,-5,5);
  TH1F* h_y_j   = new TH1F("h_y_j","jet rapidity",20,-5,5);
  TH1F* h_y_j_pt_30 = new TH1F("h_y_j_pt_30","jet rapidity pt>30",20,-5,5);
  TH1F* h_pt_j = new TH1F("h_pt_j","jet pt",100,0,500);

  TH1F* h_m_tT = new TH1F("h_m_tT","tT mass",34,300,2000);
  TH1F* h_m_tT_1 = new TH1F("h_m_tT_1","tT mass 1",34,300,2000);
  TH1F* h_m_tT_2 = new TH1F("h_m_tT_2","tT mass 2",34,300,2000);
  TH1F* h_m_tT_3 = new TH1F("h_m_tT_3","tT mass 3",34,300,2000);
  TH1F* h_m_tT_4 = new TH1F("h_m_tT_4","tT mass 4",34,300,2000);
  TH1F* h_m_tT_5 = new TH1F("h_m_tT_5","tT mass 5",34,300,2000);
  TH1F* h_m_tT_6 = new TH1F("h_m_tT_6","tT mass 6",34,300,2000);
  TH1F* h_m_tT_7 = new TH1F("h_m_tT_7","tT mass 7",34,300,2000);
  TH1F* h_m_tT_8 = new TH1F("h_m_tT_8","tT mass 8",34,300,2000);
  TH1F* h_m_tT_9 = new TH1F("h_m_tT_9","tT mass 9",34,300,2000);
  TH1F* h_m_tT_10 = new TH1F("h_m_tT_10","tT mass 10",34,300,2000);
  TH1F* h_m_tT_11 = new TH1F("h_m_tT_11","tT mass 11",34,300,2000);
  TH1F* h_m_tT_12 = new TH1F("h_m_tT_12","tT mass 12",34,300,2000);
  TH1F* h_m_tT_13 = new TH1F("h_m_tT_13","tT mass 13",34,300,2000);
  TH1F* h_m_tT_14 = new TH1F("h_m_tT_14","tT mass 14",34,300,2000);
  TH1F* h_m_tT_15 = new TH1F("h_m_tT_15","tT mass 15",34,300,2000);
  TH1F* h_m_tT_16 = new TH1F("h_m_tT_16","tT mass 16",34,300,2000);
  TH1F* h_m_tT_17 = new TH1F("h_m_tT_17","tT mass 17",34,300,2000);
  TH1F* h_m_tT_18 = new TH1F("h_m_tT_18","tT mass 18",34,300,2000);
  TH1F* h_m_tT_19 = new TH1F("h_m_tT_19","tT mass 19",34,300,2000);
  TH1F* h_m_tT_20 = new TH1F("h_m_tT_20","tT mass 20",34,300,2000);
  TH1F* h_m_tT_21 = new TH1F("h_m_tT_21","tT mass 21",34,300,2000);
  TH1F* h_m_tT_22 = new TH1F("h_m_tT_22","tT mass 22",34,300,2000);
  TH1F* h_m_tT_23 = new TH1F("h_m_tT_23","tT mass 23",34,300,2000);
  TH1F* h_m_tT_24 = new TH1F("h_m_tT_24","tT mass 24",34,300,2000);
  TH1F* h_m_tT_25 = new TH1F("h_m_tT_25","tT mass 25",34,300,2000);
  TH1F* h_m_tT_26 = new TH1F("h_m_tT_26","tT mass 26",34,300,2000);
  TH1F* h_m_tT_27 = new TH1F("h_m_tT_27","tT mass 27",34,300,2000);
  TH1F* h_m_tT_28 = new TH1F("h_m_tT_28","tT mass 28",34,300,2000);
  TH1F* h_m_tT_29 = new TH1F("h_m_tT_29","tT mass 29",34,300,2000);
  TH1F* h_m_tT_30 = new TH1F("h_m_tT_30","tT mass 30",34,300,2000);
  TH1F* h_m_tT_31 = new TH1F("h_m_tT_31","tT mass 31",34,300,2000);
  TH1F* h_m_tT_32 = new TH1F("h_m_tT_32","tT mass 32",34,300,2000);
  TH1F* h_m_tT_33 = new TH1F("h_m_tT_33","tT mass 33",34,300,2000);
  TH1F* h_m_tT_34 = new TH1F("h_m_tT_34","tT mass 34",34,300,2000);
  TH1F* h_m_tT_35 = new TH1F("h_m_tT_35","tT mass 35",34,300,2000);
  TH1F* h_m_tT_36 = new TH1F("h_m_tT_36","tT mass 36",34,300,2000);
  TH1F* h_m_tT_37 = new TH1F("h_m_tT_37","tT mass 37",34,300,2000);
  TH1F* h_m_tT_38 = new TH1F("h_m_tT_38","tT mass 38",34,300,2000);
  TH1F* h_m_tT_39 = new TH1F("h_m_tT_39","tT mass 39",34,300,2000);
  TH1F* h_m_tT_40 = new TH1F("h_m_tT_40","tT mass 40",34,300,2000);
  TH1F* h_m_tT_41 = new TH1F("h_m_tT_41","tT mass 41",34,300,2000);
  TH1F* h_m_tT_42 = new TH1F("h_m_tT_42","tT mass 42",34,300,2000);
  TH1F* h_m_tT_43 = new TH1F("h_m_tT_43","tT mass 43",34,300,2000);
  TH1F* h_m_tT_44 = new TH1F("h_m_tT_44","tT mass 44",34,300,2000);

  TH1F* h_y_t_1 = new TH1F("h_y_t_1","top rapidity pdf 1",20,-5,5);
  TH1F* h_y_t_2 = new TH1F("h_y_t_2","top rapidity pdf 2",20,-5,5);
  TH1F* h_y_t_3 = new TH1F("h_y_t_3","top rapidity pdf 3",20,-5,5);
  TH1F* h_y_t_4 = new TH1F("h_y_t_4","top rapidity pdf 4",20,-5,5);
  TH1F* h_y_t_5 = new TH1F("h_y_t_5","top rapidity pdf 5",20,-5,5);
  TH1F* h_y_t_6 = new TH1F("h_y_t_6","top rapidity pdf 6",20,-5,5);
  TH1F* h_y_t_7 = new TH1F("h_y_t_7","top rapidity pdf 7",20,-5,5);
  TH1F* h_y_t_8 = new TH1F("h_y_t_8","top rapidity pdf 8",20,-5,5);
  TH1F* h_y_t_9 = new TH1F("h_y_t_9","top rapidity pdf 9",20,-5,5);
  TH1F* h_y_t_10 = new TH1F("h_y_t_10","top rapidity pdf 10",20,-5,5);
  TH1F* h_y_t_11 = new TH1F("h_y_t_11","top rapidity pdf 11",20,-5,5);
  TH1F* h_y_t_12 = new TH1F("h_y_t_12","top rapidity pdf 12",20,-5,5);
  TH1F* h_y_t_13 = new TH1F("h_y_t_13","top rapidity pdf 13",20,-5,5);
  TH1F* h_y_t_14 = new TH1F("h_y_t_14","top rapidity pdf 14",20,-5,5);
  TH1F* h_y_t_15 = new TH1F("h_y_t_15","top rapidity pdf 15",20,-5,5);
  TH1F* h_y_t_16 = new TH1F("h_y_t_16","top rapidity pdf 16",20,-5,5);
  TH1F* h_y_t_17 = new TH1F("h_y_t_17","top rapidity pdf 17",20,-5,5);
  TH1F* h_y_t_18 = new TH1F("h_y_t_18","top rapidity pdf 18",20,-5,5);
  TH1F* h_y_t_19 = new TH1F("h_y_t_19","top rapidity pdf 19",20,-5,5);
  TH1F* h_y_t_20 = new TH1F("h_y_t_20","top rapidity pdf 20",20,-5,5);
  TH1F* h_y_t_21 = new TH1F("h_y_t_21","top rapidity pdf 21",20,-5,5);
  TH1F* h_y_t_22 = new TH1F("h_y_t_22","top rapidity pdf 22",20,-5,5);
  TH1F* h_y_t_23 = new TH1F("h_y_t_23","top rapidity pdf 23",20,-5,5);
  TH1F* h_y_t_24 = new TH1F("h_y_t_24","top rapidity pdf 24",20,-5,5);
  TH1F* h_y_t_25 = new TH1F("h_y_t_25","top rapidity pdf 25",20,-5,5);
  TH1F* h_y_t_26 = new TH1F("h_y_t_26","top rapidity pdf 26",20,-5,5);
  TH1F* h_y_t_27 = new TH1F("h_y_t_27","top rapidity pdf 27",20,-5,5);
  TH1F* h_y_t_28 = new TH1F("h_y_t_28","top rapidity pdf 28",20,-5,5);
  TH1F* h_y_t_29 = new TH1F("h_y_t_29","top rapidity pdf 29",20,-5,5);
  TH1F* h_y_t_30 = new TH1F("h_y_t_30","top rapidity pdf 30",20,-5,5);
  TH1F* h_y_t_31 = new TH1F("h_y_t_31","top rapidity pdf 31",20,-5,5);
  TH1F* h_y_t_32 = new TH1F("h_y_t_32","top rapidity pdf 32",20,-5,5);
  TH1F* h_y_t_33 = new TH1F("h_y_t_33","top rapidity pdf 33",20,-5,5);
  TH1F* h_y_t_34 = new TH1F("h_y_t_34","top rapidity pdf 34",20,-5,5);
  TH1F* h_y_t_35 = new TH1F("h_y_t_35","top rapidity pdf 35",20,-5,5);
  TH1F* h_y_t_36 = new TH1F("h_y_t_36","top rapidity pdf 36",20,-5,5);
  TH1F* h_y_t_37 = new TH1F("h_y_t_37","top rapidity pdf 37",20,-5,5);
  TH1F* h_y_t_38 = new TH1F("h_y_t_38","top rapidity pdf 38",20,-5,5);
  TH1F* h_y_t_39 = new TH1F("h_y_t_39","top rapidity pdf 39",20,-5,5);
  TH1F* h_y_t_40 = new TH1F("h_y_t_40","top rapidity pdf 40",20,-5,5);
  TH1F* h_y_t_41 = new TH1F("h_y_t_41","top rapidity pdf 41",20,-5,5);
  TH1F* h_y_t_42 = new TH1F("h_y_t_42","top rapidity pdf 42",20,-5,5);
  TH1F* h_y_t_43 = new TH1F("h_y_t_43","top rapidity pdf 43",20,-5,5);
  TH1F* h_y_t_44 = new TH1F("h_y_t_44","top rapidity pdf 44",20,-5,5);

  TH1F* h_y_j_29   = new TH1F("h_y_j_29","jet rapidity pdf 29",20,-5,5);
  TH1F* h_y_j_30   = new TH1F("h_y_j_30","jet rapidity pdf 30",20,-5,5);
  

  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 vt(px3,py3,pz3,E_3);   // top 4-vector
    TLorentzVector vtb(px4,py4,pz4,E_4);  // tbar 4-vector
    TLorentzVector vj(px5,py5,pz5,E_5);  //  jet 4-vector

    float m_tT = (vt + vtb).M();        // tT mass

    float y_t = vt.Rapidity();           // top rapidity
    float eta_t = vt.PseudoRapidity();   // top pseudorapidity
    float phi_t = vt.Phi();              // top phi
    float pt_t = vt.Perp();              // top pT

    float y_tb = vtb.Rapidity();         // tbar rapidity
	float eta_tb = vtb.PseudoRapidity();  // tbar pseudorapidity
    float phi_tb = vtb.Phi();            // tbar phi
    float pt_tb = vtb.Perp();            // tbar pT

    if(vj.E() !=0){

    float y_j   = vj.Rapidity();       // jet rapidity
	float eta_j = vj.PseudoRapidity(); // jet pseudorapidity
    float phi_j = vj.Phi();            // jet phi
    float pt_j = vj.Perp();            // jet  pT
    }

    h_y_t->Fill(y_t,wt_ALL);
    h_y_tb->Fill(y_tb,wt_ALL);
    h_m_tT->Fill(m_tT,wt_ALL);
    h_y_j->Fill(y_j,wt_ALL);
    h_pt_j->Fill(pt_j,wt_ALL);
       
        h_m_tT_1->Fill(m_tT,PDF01);
	h_m_tT_2->Fill(m_tT,PDF02);
	h_m_tT_3->Fill(m_tT,PDF03);
	h_m_tT_4->Fill(m_tT,PDF04);
	h_m_tT_5->Fill(m_tT,PDF05);
	h_m_tT_6->Fill(m_tT,PDF06);
	h_m_tT_7->Fill(m_tT,PDF07);
	h_m_tT_8->Fill(m_tT,PDF08);
	h_m_tT_9->Fill(m_tT,PDF09);
        h_m_tT_10->Fill(m_tT,PDF10); 
	h_m_tT_11->Fill(m_tT,PDF11);
	h_m_tT_12->Fill(m_tT,PDF12);
	h_m_tT_13->Fill(m_tT,PDF13);
	h_m_tT_14->Fill(m_tT,PDF14);
	h_m_tT_15->Fill(m_tT,PDF15);
	h_m_tT_16->Fill(m_tT,PDF16);
	h_m_tT_17->Fill(m_tT,PDF17);
	h_m_tT_18->Fill(m_tT,PDF18);
	h_m_tT_19->Fill(m_tT,PDF19);
	h_m_tT_20->Fill(m_tT,PDF20);
	h_m_tT_21->Fill(m_tT,PDF21);
	h_m_tT_22->Fill(m_tT,PDF22);
	h_m_tT_23->Fill(m_tT,PDF23);
	h_m_tT_24->Fill(m_tT,PDF24);
	h_m_tT_25->Fill(m_tT,PDF25); 
	h_m_tT_26->Fill(m_tT,PDF26); 
	h_m_tT_27->Fill(m_tT,PDF27);
	h_m_tT_28->Fill(m_tT,PDF28);
	h_m_tT_29->Fill(m_tT,PDF29);
	h_m_tT_30->Fill(m_tT,PDF30);
	h_m_tT_31->Fill(m_tT,PDF31);
 	h_m_tT_32->Fill(m_tT,PDF32);
 	h_m_tT_33->Fill(m_tT,PDF33);
 	h_m_tT_34->Fill(m_tT,PDF34);
 	h_m_tT_35->Fill(m_tT,PDF35);
 	h_m_tT_36->Fill(m_tT,PDF36);
 	h_m_tT_37->Fill(m_tT,PDF37);
 	h_m_tT_38->Fill(m_tT,PDF38);
 	h_m_tT_39->Fill(m_tT,PDF39);
 	h_m_tT_40->Fill(m_tT,PDF40);
        h_m_tT_41->Fill(m_tT,PDF41);
        h_m_tT_42->Fill(m_tT,PDF42);
        h_m_tT_43->Fill(m_tT,PDF43);
        h_m_tT_44->Fill(m_tT,PDF44);
	

	h_y_t_1->Fill(y_t,PDF01);
	h_y_t_2->Fill(y_t,PDF02);
	h_y_t_3->Fill(y_t,PDF03);
	h_y_t_4->Fill(y_t,PDF04);
	h_y_t_5->Fill(y_t,PDF05);
	h_y_t_6->Fill(y_t,PDF06);
	h_y_t_7->Fill(y_t,PDF07);
	h_y_t_8->Fill(y_t,PDF08);
	h_y_t_9->Fill(y_t,PDF09);
	h_y_t_10->Fill(y_t,PDF10);
	h_y_t_11->Fill(y_t,PDF11);
	h_y_t_12->Fill(y_t,PDF12);
	h_y_t_13->Fill(y_t,PDF13); 
	h_y_t_14->Fill(y_t,PDF14);
	h_y_t_15->Fill(y_t,PDF15);
	h_y_t_16->Fill(y_t,PDF16);
	h_y_t_17->Fill(y_t,PDF17);
	h_y_t_18->Fill(y_t,PDF18);
	h_y_t_19->Fill(y_t,PDF19);
	h_y_t_20->Fill(y_t,PDF20);
	h_y_t_21->Fill(y_t,PDF21);
	h_y_t_22->Fill(y_t,PDF22);
	h_y_t_23->Fill(y_t,PDF23);
	h_y_t_24->Fill(y_t,PDF24);
	h_y_t_25->Fill(y_t,PDF25);
	h_y_t_26->Fill(y_t,PDF26);
	h_y_t_27->Fill(y_t,PDF27);
	h_y_t_28->Fill(y_t,PDF28);
	h_y_t_29->Fill(y_t,PDF29);
	h_y_t_30->Fill(y_t,PDF30);
	h_y_t_31->Fill(y_t,PDF31);
	h_y_t_32->Fill(y_t,PDF32);
	h_y_t_33->Fill(y_t,PDF33);
	h_y_t_34->Fill(y_t,PDF34);
	h_y_t_35->Fill(y_t,PDF35);
	h_y_t_36->Fill(y_t,PDF36);
	h_y_t_37->Fill(y_t,PDF37);
	h_y_t_38->Fill(y_t,PDF38);
	h_y_t_39->Fill(y_t,PDF39);
	h_y_t_40->Fill(y_t,PDF40);
        h_y_t_41->Fill(y_t,PDF41);
        h_y_t_42->Fill(y_t,PDF42);
        h_y_t_43->Fill(y_t,PDF43);
        h_y_t_44->Fill(y_t,PDF44);
	
  }
  fout->Write();
  fout->Close();
}