File contents

#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("h300");
  fin->Add("./*.root");


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

  float wt,gg,gq,qq,qqb;
  float w_1,w_2,w_3,w_4,w_5,w_6,w_7,w_8,w_9,w10;
  float w11,w12,w13,w14,w15,w16,w17,w18,w19,w20;
  float w21,w22,w23,w24,w25,w26,w27,w28,w29,w30;
  float w31,w32,w33,w34,w35,w36,w37,w38,w39,w40;

  float px3,px4,px5;
  float py3,py4,py5;
  float pz3,pz4,pz5;
  float E3,E4,E5;

  h300->SetBranchAddress("E3" ,&E3);
  h300->SetBranchAddress("E4" ,&E4);
  h300->SetBranchAddress("E5" ,&E5);

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

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

  h300->SetBranchAddress("pz3",&pz3);
  h300->SetBranchAddress("pz4",&pz4);
  h300->SetBranchAddress("pz5",&pz5);
 
  h300->SetBranchAddress("wt",&wt);
  h300->SetBranchAddress("gg",&gg);
  h300->SetBranchAddress("gq",&gq);
  h300->SetBranchAddress("qq",&qq);
  h300->SetBranchAddress("qqb",&qqb);

  h300->SetBranchAddress("w_1",&w_1);
  h300->SetBranchAddress("w_2",&w_2);
  h300->SetBranchAddress("w_3",&w_3);
  h300->SetBranchAddress("w_4",&w_4);
  h300->SetBranchAddress("w_5",&w_5);
  h300->SetBranchAddress("w_6",&w_6);
  h300->SetBranchAddress("w_7",&w_7);
  h300->SetBranchAddress("w_8",&w_8);
  h300->SetBranchAddress("w_9",&w_9);
  h300->SetBranchAddress("w10",&w10);
  h300->SetBranchAddress("w11",&w11);
  h300->SetBranchAddress("w12",&w12);
  h300->SetBranchAddress("w13",&w13);
  h300->SetBranchAddress("w14",&w15);
  h300->SetBranchAddress("w15",&w15);
  h300->SetBranchAddress("w16",&w16);
  h300->SetBranchAddress("w17",&w17);
  h300->SetBranchAddress("w18",&w18);
  h300->SetBranchAddress("w19",&w19);
  h300->SetBranchAddress("w20",&w20);
  h300->SetBranchAddress("w21",&w21);
  h300->SetBranchAddress("w22",&w22);
  h300->SetBranchAddress("w23",&w23);
  h300->SetBranchAddress("w24",&w24);
  h300->SetBranchAddress("w25",&w25);
  h300->SetBranchAddress("w26",&w26);
  h300->SetBranchAddress("w27",&w27);
  h300->SetBranchAddress("w28",&w28);
  h300->SetBranchAddress("w29",&w29);
  h300->SetBranchAddress("w30",&w30);
  h300->SetBranchAddress("w31",&w31);
  h300->SetBranchAddress("w32",&w32);
  h300->SetBranchAddress("w33",&w33);
  h300->SetBranchAddress("w34",&w34);
  h300->SetBranchAddress("w35",&w35);
  h300->SetBranchAddress("w36",&w36);
  h300->SetBranchAddress("w37",&w37);
  h300->SetBranchAddress("w38",&w38);
  h300->SetBranchAddress("w39",&w39);
  h300->SetBranchAddress("w40",&w40);
 
  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_y_t_20 = new TH1F("h_y_t_20","top y pt_j> 20",20,-5,5);
  TH1F* h_y_t_30 = new TH1F("h_y_t_30","top y pt_j> 30",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)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 vt(px3,py3,pz3,E3);   // top 4-vector
    TLorentzVector vtb(px4,py4,pz4,E4);  // tbar 4-vector
    TLorentzVector vj(px5,py5,pz5,E5);  //  jet 4-vector

    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

    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);
    h_y_tb->Fill(y_tb,wt);
    h_y_j->Fill(y_j,wt);
    h_y_j_29->Fill(y_j,w29);
    h_y_j_30->Fill(y_j,w30);
    h_pt_j->Fill(pt_j,wt);
	if(pt_j>20){
		h_y_t_20->Fill(y_t,wt);
		}
	if(pt_j>30){
		h_y_t_30->Fill(y_t,wt);
		h_y_j_pt_30->Fill(y_j,wt);
		}

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