// Top Analysis
// Rosi Reed
// Phy 252b Final Project
// Fall '07
//
// Sample Analysis
// Bob McElrath (bob@mcelrath.org) 11/17/2005
// 
//   Slightly adapted for round 2 of the LHC Olympics from Gordon Watt's
//   original code:
//      http://d0.phys.washington.edu/~gwatts/research/Conferences/2005-LHCOlympics/default.htm

#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TMath.h"

#include <iostream>
using std::cout;
using std::endl;

#include "OlympicEvent.h"

void TopAnalysis() {

  //TFile *f = new TFile ("~erbacher/252/data/W-only.root", "READ");
  TFile *f = new TFile ("~erbacher/252/data/top-only.root", "READ");
  TTree *t = static_cast<TTree*> (f->Get("LHCTree"));
  OlympicEvent *e = new OlympicEvent();
  t->SetBranchAddress ("Event", &e);

  TFile *output = new TFile("plots.root", "RECREATE");

  const double PI = TMath::Pi();

  TH1F *h_nele = new TH1F ("ele_n", "Number of electrons", 10, 0.0, 10.0);
  TH1F *h_njet = new TH1F ("jet_n", "Number of jets", 30, 0.0, 30.0);
  TH1F *h_npho = new TH1F ("pho_n", "Number of photons", 10, 0.0, 10.0);
  TH1F *h_ntau = new TH1F ("tau_n", "Number of hadronic taus", 10, 0.0, 10.0);
  TH1F *h_nmuo = new TH1F ("muo_n", "Number of muons", 10, 0.0, 10.0);

  TH1F *h_ptele = new TH1F ("ele_pt", "Pt of electrons", 50, 0.0, 250.0);
  TH1F *h_ptjet = new TH1F ("jet_pt", "Pt of jets", 50, 0.0, 250.0);
  TH1F *h_ptpho = new TH1F ("pho_pt", "Pt of photons", 50, 0.0, 200.0);
  TH1F *h_pttau = new TH1F ("tau_pt", "Pt of hadronic taus", 50, 0.0, 200.0);
  TH1F *h_ptmuo = new TH1F ("muo_pt", "Pt of muons", 50, 0.0, 250.0);
  TH1F *h_ptmet = new TH1F ("met_pt", "Pt of missing energy", 50, 0.0, 200.0);

  TH1F *h_etaele = new TH1F ("ele_eta", "Eta of electrons", 50, -2.5, 2.5);
  TH1F *h_etajet = new TH1F ("jet_eta", "Eta of jets", 50, -2.5, 2.5);
  TH1F *h_etapho = new TH1F ("pho_eta", "Eta of photons", 50, -2.5, 2.5);
  TH1F *h_etatau = new TH1F ("tau_eta", "Eta of hadronic taus", 50, -2.5, 2.5);
  TH1F *h_etamuo = new TH1F ("muo_eta", "Eta of muons", 50, -2.5, 2.5);
  TH1F *h_etamet = new TH1F ("met_eta", "Eta of missing energy", 50, -2.5, 2.5);

  TH1F *h_phiele = new TH1F ("ele_phi", "Phi of electrons", 50, -PI, PI);
  TH1F *h_phijet = new TH1F ("jet_phi", "Phi of jets", 50, -PI, PI);
  TH1F *h_phipho = new TH1F ("pho_phi", "Phi of photons", 50, -PI, PI);
  TH1F *h_phitau = new TH1F ("tau_phi", "Phi of hadronic taus", 50, -PI, PI);
  TH1F *h_phimuo = new TH1F ("muo_phi", "Phi of muons", 50, -PI, PI);
  TH1F *h_phimet = new TH1F ("met_phi", "Phi of missing energy", 50, -PI, PI);

  TH1F *h_dijetmass = new TH1F("dijet_mass", "Invariant mass of jet pairs", 50, 0, 200);
  TH1F *h_dileptonmass = new TH1F("dilep_mass", "Invariant mass of lepton pairs", 50, 0, 200);
  TH1F *h_Ht = new TH1F("Ht","Event Transverse Energy",50,0,1000);
  TH1F *h_lepsize = new TH1F("lepsize","Number of leptons(e,muo)",50,0.0,10.0);


  double ptlepjets;
  int IsGood;
  double elemax;
  double muomax;
  double jetmax;

  int event_number = 1;
   while (t->GetEntry(event_number) > 0) {

     //The following loops give values that I then use for selection criteria
     elemax = 0;muomax = 0;jetmax = 0;
     for (int i = 0; i < e->_electrons.size(); i++){
       if (e->_electrons[i].Pt() > elemax)
	 elemax = e->_electrons[i].Pt();}
     for(int i = 0; i < e->_muons.size(); i++){
       if (e->_muons[i].Pt() > muomax)
	 muomax = e->_muons[i].Pt();}
     for(int i = 0; i < e->_jets.size(); i++){
       if (e->_jets[i].Pt() > jetmax)
	 jetmax = e->_jets[i].Pt();}

     //Uncomment the following if statement to selection on dilepton events  
     //if ((e->_met.Pt() > 20) and (((e->_electrons.size() > 1) and (elemax > 20)) or ((e->_muons.size() > 1) and (muomax > 20)) or (((e->_electrons.size()+e->_muons.size()) > 1) and ((muomax > 20) or (elemax > 20)))))
     //
     //The following if statement selections lepton+jets events
     //if ((e->_jets.size()>2) and (e->_met.Pt() > 20) and ((muomax > 20) or (elemax > 20)))
     //
     //The following if statement selects hadronic decays
     // if ((e->_jets.size() > 3) and ((e->_electrons.size() + e->_muons.size()) < 1) and (jetmax > 20))
     //
     //If statement for extra questions 10 and 11
     //if ((e->_jets.size() == 3) and ((e->_electrons.size() + e->_muons.size()) > 0))
       {
   
	IsGood = 0;
	
	ptlepjets = 0;
	h_nele->Fill (e->_electrons.size());
	h_nmuo->Fill (e->_muons.size());
	h_ntau->Fill (e->_hadronic_taus.size());
	h_njet->Fill (e->_jets.size());
	h_npho->Fill (e->_photons.size());
	h_lepsize->Fill((e->_electrons.size() + e->_muons.size()));

	h_ptmet->Fill (e->_met.Pt());
	h_etamet->Fill (e->_met.Eta());
	h_phimet->Fill (e->_met.Phi());

	for (int i = 0; i < e->_electrons.size(); i++) {
	  ptlepjets = ptlepjets + e->_electrons[i].Pt();
	  h_ptele->Fill(e->_electrons[i].Pt());
	  h_etaele->Fill(e->_electrons[i].Eta());
	  h_phiele->Fill(e->_electrons[i].Phi());
          for(int j = i+1; j < e->_electrons.size(); j++) {
	    if (e->_electrons[i].Pt() > 20 & e->_electrons[j].Pt() > 20){
	    TLorentzVector leptonpair = e->_electrons[i] + e->_electrons[j];
	    h_dileptonmass->Fill(leptonpair.M());}
          }
	}
	for (int i = 0; i < e->_muons.size(); i++) {
	  ptlepjets = ptlepjets + e->_muons[i].Pt();
	  h_ptmuo->Fill(e->_muons[i].Pt());
	  h_etamuo->Fill(e->_muons[i].Eta());
	  h_phimuo->Fill(e->_muons[i].Phi());
          for(int j = i+1; j < e->_muons.size(); j++) {
	    if ((e->_muons[i].Pt() > 20 & e->_muons[j].Pt() > 20)) 
	    {
		TLorentzVector leptonpair = e->_muons[i] + e->_muons[j];
		h_dileptonmass->Fill(leptonpair.M());}
          }
	}
	for (int i = 0; i < e->_hadronic_taus.size(); i++) {
	  h_pttau->Fill(e->_hadronic_taus[i].Pt());
	  h_etatau->Fill(e->_hadronic_taus[i].Eta());
	  h_phitau->Fill(e->_hadronic_taus[i].Phi());
	}
	for (int i = 0; i < e->_jets.size(); i++) {
	  ptlepjets = ptlepjets + e->_jets[i].Pt();
	  h_ptjet->Fill(e->_jets[i].Pt());
	  h_etajet->Fill(e->_jets[i].Eta());
	  h_phijet->Fill(e->_jets[i].Phi());
          for(int j = i+1; j < e->_jets.size(); j++) {
	    if (e->_jets[i].Pt()>25 & e->_jets[j].Pt()>25){
	    TLorentzVector jetpair = e->_jets[i] + e->_jets[j];
	    h_dijetmass->Fill(jetpair.M());}
          }
	}
	for (int i = 0; i < e->_photons.size(); i++) {
	  h_ptpho->Fill(e->_photons[i].Pt());
	  h_etapho->Fill(e->_photons[i].Eta());
	  h_phipho->Fill(e->_photons[i].Phi());
	}
	//transverse energy is the sum of transverse energies of the leptons,
	//jets and missing Et.
	h_Ht->Fill((e->_met.Pt() + ptlepjets));


	///
	/// Event number
	///
	
	}//end selection if statements
	event_number++;
	if (event_number % 5000 == 0) {
	  cout << "Just did event " << event_number << endl;
	}
	e->Clear();
  }

  f->Close();
  output->Write();
  output->Close();
}

int main(void) {
    SampleAnalysis();
}