csander/HEPTutorial/MyAnalysis.C

#define MyAnalysis_cxx
// The class definition in MyAnalysis.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("MyAnalysis.C")
// Root > T->Process("MyAnalysis.C","some options")
// Root > T->Process("MyAnalysis.C+")
//

#include "MyAnalysis.h"
#include <iostream>
#include <TH1F.h>
#include <TLatex.h>

using namespace std;

void MyAnalysis::BuildEvent() {
   Jets.clear();
   for (int i = 0; i < NJet; ++i) {
      MyJet jet(Jet_Px[i], Jet_Py[i], Jet_Pz[i], Jet_E[i]);
      jet.SetBTagDiscriminator(Jet_btag[i]);
      Jets.push_back(jet);
   }
   Muons.clear();
   for (int i = 0; i < NMuon; ++i) {
      MyMuon muon(Muon_Px[i], Muon_Py[i], Muon_Pz[i], Muon_E[i]);
      muon.SetIsolation(Muon_Iso[i]);
      muon.SetCharge(Muon_Charge[i]);
      Muons.push_back(muon);
   }
   Electrons.clear();
   for (int i = 0; i < NElectron; ++i) {
      MyElectron electron(Electron_Px[i], Electron_Py[i], Electron_Pz[i], Electron_E[i]);
      electron.SetIsolation(Electron_Iso[i]);
      electron.SetCharge(Electron_Charge[i]);
      Electrons.push_back(electron);
   }
   Photons.clear();
   for (int i = 0; i < NPhoton; ++i) {
      MyPhoton photon(Photon_Px[i], Photon_Py[i], Photon_Pz[i], Photon_E[i]);
      photon.SetIsolation(Photon_Iso[i]);
      Photons.push_back(photon);
   }
   hadB.SetXYZM(MChadronicBottom_px, MChadronicBottom_py, MChadronicBottom_pz, 4.8);
   lepB.SetXYZM(MCleptonicBottom_px, MCleptonicBottom_py, MCleptonicBottom_pz, 4.8);
   hadWq.SetXYZM(MChadronicWDecayQuark_px, MChadronicWDecayQuark_py, MChadronicWDecayQuark_pz, 0.0);
   hadWqb.SetXYZM(MChadronicWDecayQuarkBar_px, MChadronicWDecayQuarkBar_py, MChadronicWDecayQuarkBar_pz, 0.0);
   lepWl.SetXYZM(MClepton_px, MClepton_py, MClepton_pz, 0.0);
   lepWn.SetXYZM(MCneutrino_px, MCneutrino_py, MCneutrino_pz, 0.0);
   met.SetXYZM(MET_px, MET_py, 0., 0.);

   if (EventWeight < 1.e-5)
      EventWeight = 1.;
   EventWeight *= weight_factor;

}

void MyAnalysis::Begin(TTree * /*tree*/) {
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

}

void MyAnalysis::SlaveBegin(TTree * /*tree*/) {
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

   h_Mmumu = new TH1F("Mmumu", "Invariant di-muon mass", 50, 0, 150);
   h_Mmumu->SetXTitle("m_{#mu#mu}");
   h_Mmumu->Sumw2();

   h_Mbqqb_mc = new TH1F("Mbqqb_mc", "Invariant hadronic top mass", 50, 170, 175);
   h_Mbqqb_mc->SetXTitle("m_{hadronic top}");
   h_Mbqqb_mc->Sumw2();

   h_Mbln_mc = new TH1F("Mbln_mc", "Invariant leptonic top mass", 50, 170, 175);
   h_Mbln_mc->SetXTitle("m_{leptonic top}");
   h_Mbln_mc->Sumw2();

   h_Mbqqb_reco = new TH1F("Mbqqb_reco", "Invariant hadronic top mass", 50, 0, 500);
   h_Mbqqb_reco->SetXTitle("m_{hadronic top}");
   h_Mbqqb_reco->Sumw2();

   h_Mbln_reco = new TH1F("Mbln_reco", "Invariant transverse leptonic top mass", 50, 0, 500);
   h_Mbln_reco->SetXTitle("m_{T, leptonic top}");
   h_Mbln_reco->Sumw2();

   h_NJet = new TH1F("NJet", "Number of jets", 7, 0, 7);
   h_NJet->SetXTitle("NJet");
   h_NJet->Sumw2();

   h_NBJet = new TH1F("NBJet", "Number of b-jets", 7, 0, 7);
   h_NBJet->SetXTitle("NBJet");
   h_NBJet->Sumw2();

   h_NMuon = new TH1F("NMuon", "Number of muons", 7, 0, 7);
   h_NMuon->SetXTitle("NMuon");
   h_NMuon->Sumw2();

   h_Jet1_Pt = new TH1F("Jet1Pt", "Pt of jet1", 50, 0, 250);
   h_Jet1_Pt->SetXTitle("p_{T} [GeV]");
   h_Jet1_Pt->Sumw2();

   h_Jet2_Pt = new TH1F("Jet2Pt", "Pt of jet2", 50, 0, 250);
   h_Jet2_Pt->SetXTitle("p_{T} [GeV]");
   h_Jet2_Pt->Sumw2();

   h_Jet3_Pt = new TH1F("Jet3Pt", "Pt of jet3", 50, 0, 250);
   h_Jet3_Pt->SetXTitle("p_{T} [GeV]");
   h_Jet3_Pt->Sumw2();

   h_Jet1_Eta = new TH1F("Jet1Eta", "#eta of jet1", 50, -4, 4);
   h_Jet1_Eta->SetXTitle("#eta");
   h_Jet1_Eta->Sumw2();

   h_Jet2_Eta = new TH1F("Jet2Eta", "#eta of jet2", 50, -4, 4);
   h_Jet2_Eta->SetXTitle("#eta");
   h_Jet2_Eta->Sumw2();

   h_Jet3_Eta = new TH1F("Jet3Eta", "#eta of jet3", 50, -4, 4);
   h_Jet3_Eta->SetXTitle("#eta");
   h_Jet3_Eta->Sumw2();

   h_BJet1_Pt = new TH1F("BJet1Pt", "Pt of b-jet1", 50, 0, 250);
   h_BJet1_Pt->SetXTitle("p_{T} [GeV]");
   h_BJet1_Pt->Sumw2();

   h_BJet2_Pt = new TH1F("BJet2Pt", "Pt of b-jet2", 50, 0, 250);
   h_BJet2_Pt->SetXTitle("p_{T} [GeV]");
   h_BJet2_Pt->Sumw2();

   h_BJet1_Eta = new TH1F("BJet1Eta", "#eta of b-jet1", 50, -4, 4);
   h_BJet1_Eta->SetXTitle("#eta");
   h_BJet1_Eta->Sumw2();

   h_BJet2_Eta = new TH1F("BJet2Eta", "#eta of b-jet2", 50, -4, 4);
   h_BJet2_Eta->SetXTitle("#eta");
   h_BJet2_Eta->Sumw2();

   h_Muon1_Pt = new TH1F("Muon1Pt", "Pt of muon1", 50, 0, 250);
   h_Muon1_Pt->SetXTitle("p_{T} [GeV]");
   h_Muon1_Pt->Sumw2();

   h_Muon2_Pt = new TH1F("Muon2Pt", "Pt of muon2", 50, 0, 250);
   h_Muon2_Pt->SetXTitle("p_{T} [GeV]");
   h_Muon2_Pt->Sumw2();

   h_Muon1_Eta = new TH1F("Muon1Eta", "#eta of muon1", 50, -4, 4);
   h_Muon1_Eta->SetXTitle("#eta");
   h_Muon1_Eta->Sumw2();

   h_Muon2_Eta = new TH1F("Muon2Eta", "#eta of muon2", 50, -4, 4);
   h_Muon2_Eta->SetXTitle("#eta");
   h_Muon2_Eta->Sumw2();

   h_Muon1_Iso = new TH1F("Muon1Iso", "Isolation of muon1", 50, 0, 25);
   h_Muon1_Iso->SetXTitle("Isolation [GeV]");
   h_Muon1_Iso->Sumw2();

   h_Muon2_Iso = new TH1F("Muon2Iso", "Isolation of muon2", 50, 0, 25);
   h_Muon2_Iso->SetXTitle("Isolation [GeV]");
   h_Muon2_Iso->Sumw2();

   h_MET = new TH1F("MET", "MET", 30, 0, 300.);
   h_MET->SetXTitle("MET");
   h_MET->Sumw2();

   h_minDeltaPhi_MET_Muon = new TH1F("minDeltaPhi_MET_Muon", "min#Delta#phi(MET,muon)", 50, 0., TMath::Pi());
   h_minDeltaPhi_MET_Muon->SetXTitle("min#Delta#phi [rad]");
   h_minDeltaPhi_MET_Muon->Sumw2();

   h_minDeltaPhi_MET_BJet = new TH1F("minDeltaPhi_MET_BJet", "min#Delta#phi(MET,b-jet)", 50, 0., TMath::Pi());
   h_minDeltaPhi_MET_BJet->SetXTitle("min#Delta#phi [rad]");
   h_minDeltaPhi_MET_BJet->Sumw2();

   h_selectedEvents_Muon1_Pt = new TH1F("selectedEvents_Muon1_Pt", "Pt of muon1 (selected)", 50, 0, 250);
   h_selectedEvents_Muon1_Pt->SetXTitle("p_{T} [GeV]");
   h_selectedEvents_Muon1_Pt->Sumw2();

   h_selectedEvents_triggered_Muon1_Pt = new TH1F("selectedEvents_triggered_Muon1_Pt",
         "Pt of muon1 (selected and triggered)", 50, 0, 250);
   h_selectedEvents_triggered_Muon1_Pt->SetXTitle("p_{T} [GeV]");
   h_selectedEvents_triggered_Muon1_Pt->Sumw2();

}

Bool_t MyAnalysis::Process(Long64_t entry) {
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // It can be passed to either MyAnalysis::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. When processing
   // keyed objects with PROOF, the object is already loaded and is available
   // via the fObject pointer.
   //
   // This function should contain the "body" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

   ++TotalEvents;

   GetEntry(entry);

   if (TotalEvents % 10000 == 0)
      cout << "Next event -----> " << TotalEvents << endl;

   BuildEvent();

   //   cout << "Jets: " << endl;
   //   for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
   //      cout << "pt, eta, phi, btag: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->IsBTagged()
   //            << endl;
   //   }
   //   cout << "Muons: " << endl;
   //   for (vector<MyMuon>::iterator it = Muons.begin(); it != Muons.end(); ++it) {
   //      cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
   //            << it->GetIsolation() << ", " << it->GetCharge() << endl;
   //   }
   //   cout << "Electrons: " << endl;
   //   for (vector<MyElectron>::iterator it = Electrons.begin(); it != Electrons.end(); ++it) {
   //      cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
   //            << it->GetIsolation() << ", " << it->GetCharge() << endl;
   //   }
   //   cout << "Photons: " << endl;
   //   for (vector<MyPhoton>::iterator it = Photons.begin(); it != Photons.end(); ++it) {
   //      cout << "pt, eta, phi, iso: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->GetIsolation()
   //            << endl;
   //   }

   //////////////////////////////
   // Exercise XX: Invariant Di-Muon mass
   if (NMuon > 1) {
      if (Muons.at(0).GetCharge() * Muons.at(1).GetCharge() < 1 && Muons.at(0).IsIsolated() && Muons.at(1).IsIsolated()) {
         h_Mmumu->Fill((Muons.at(0) + Muons.at(1)).M(), EventWeight);
      }
   }
   //////////////////////////////


   //////////////////////////////
   // Exercise XX: MC / data comparisons

   // require at least one isolapted muon
   // +trigger bit
   if (triggerIsoMu24 && NMuon > 0) {
      if (Muons.at(0).Pt() > 24. && Muons.at(0).IsIsolated()) {

         int N_BJet = 0;
         int N_Jet = 0;
         for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
            ++N_Jet;
            if (N_Jet == 1) {
               h_Jet1_Pt->Fill(it->Pt(), EventWeight);
               h_Jet1_Eta->Fill(it->Eta(), EventWeight);
            } else if (N_Jet == 2) {
               h_Jet2_Pt->Fill(it->Pt(), EventWeight);
               h_Jet2_Eta->Fill(it->Eta(), EventWeight);
            } else if (N_Jet == 3) {
               h_Jet3_Pt->Fill(it->Pt(), EventWeight);
               h_Jet3_Eta->Fill(it->Eta(), EventWeight);
            }
            if (it->IsBTagged()) {
               ++N_BJet;
               if (N_BJet == 1) {
                  h_BJet1_Pt->Fill(it->Pt(), EventWeight);
                  h_BJet1_Eta->Fill(it->Eta(), EventWeight);
               } else if (N_BJet == 2) {
                  h_BJet2_Pt->Fill(it->Pt(), EventWeight);
                  h_BJet2_Eta->Fill(it->Eta(), EventWeight);
               }
            }
         }
         h_NBJet->Fill(N_BJet, EventWeight);
         h_NJet->Fill(N_Jet, EventWeight);

         int N_Muon = 0;
         for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
            if (jt->IsIsolated()) {
               ++N_Muon;
               if (N_Muon == 1) {
                  h_Muon1_Pt->Fill(jt->Pt(), EventWeight);
                  h_Muon1_Eta->Fill(jt->Eta(), EventWeight);
                  h_Muon1_Iso->Fill(jt->GetIsolation(), EventWeight);
               } else if (N_Muon == 2) {
                  h_Muon2_Pt->Fill(jt->Pt(), EventWeight);
                  h_Muon2_Eta->Fill(jt->Eta(), EventWeight);
                  h_Muon2_Iso->Fill(jt->GetIsolation(), EventWeight);
               }
            }
         }
         h_NMuon->Fill(N_Muon, EventWeight);

         h_MET->Fill(met.Pt(), EventWeight);

         float dPhiMuon = 999;
         for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
            float tmp;
            if (jt->Pt() > 24.) {
               tmp = fabs(jt->DeltaPhi(met));
               if (tmp < dPhiMuon)
                  dPhiMuon = tmp;
            }
         }
         if (dPhiMuon < 10.)
            h_minDeltaPhi_MET_Muon->Fill(dPhiMuon, EventWeight);

         float dPhiBJet = 999;
         for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
            float tmp;
            if (it->IsBTagged()) {
               tmp = fabs(it->DeltaPhi(met));
               if (tmp < dPhiBJet)
                  dPhiBJet = tmp;
            }
         }
         if (dPhiBJet < 10.)
            h_minDeltaPhi_MET_BJet->Fill(dPhiBJet, EventWeight);

      } // end of muon requirement
   } //end of trigger requirement
   //////////////////////////////

   //////////////////////////////
   // Exercise XXa: Trigger turn om
   if (Muons.size() == 1) {
      if (Muons.at(0).IsIsolated()) {
         h_selectedEvents_Muon1_Pt->Fill(Muons.at(0).Pt(), EventWeight);
         if (triggerIsoMu24)
            h_selectedEvents_triggered_Muon1_Pt->Fill(Muons.at(0).Pt(), EventWeight);
      }
   }
   //////////////////////////////

   //////////////////////////////
   // Exercise XXb: Acceptance and Trigger Efficiencies
   GeneratedEvents += EventWeight;
   bool IsSelected = false;
   if (NMuon > 0) {
      if (Muons.at(0).Pt() > 24. && Muons.at(0).IsIsolated()) {

         int N_BJet = 0;
         for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
            if (it->IsBTagged()) {
               ++N_BJet;
            }
         }
         if (N_BJet > 1) {
            SelectedEvents += EventWeight;
            if (triggerIsoMu24) {
               SelectedEvents_triggered += EventWeight;
               IsSelected = true;
            }
         }

      }
   }
   //////////////////////////////

   //////////////////////////////
   // Exercixe XXa: Top mass (MC Truth)
   h_Mbqqb_mc->Fill((hadB + hadWq + hadWqb).M(), EventWeight);
   h_Mbln_mc->Fill((lepB + lepWl + lepWn).M(), EventWeight);
   //////////////////////////////

   //////////////////////////////
   // Exercixe XXb: Top mass (hadroninc top decay) (reco)
   // Exercixe XXc: Top mass (leptonic top decay) (reco)
   if (IsSelected) {

      for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
         if (it->IsBTagged()) {
            for (vector<MyJet>::iterator jt = Jets.begin(); jt != Jets.end(); ++jt) {
               if (jt != it && !(jt->IsBTagged())) {
                  for (vector<MyJet>::iterator kt = jt; kt != Jets.end(); ++kt) {
                     if (kt != it && kt != jt && !(kt->IsBTagged())) {
                        if (((*jt) + (*kt)).M() > 65 && ((*jt) + (*kt)).M() < 95)
                           h_Mbqqb_reco->Fill(((*it) + (*jt) + (*kt)).M(), EventWeight);
                     }
                  }
               }
            }
         }

         for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
            if (it->IsBTagged()) {
               for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
                  if (jt->IsIsolated() && jt->Pt() > 24.) {
                     double px = met.Px();
                     double py = met.Py();
                     double mW = 80.4;
                     double A = pow(jt->E(), 2) - pow(jt->Pz(), 2);
                     double B = mW * mW / 2. + (jt->Px() * px + jt->Py() * py);
                     double D = pow(jt->E(), 2) * (pow(B, 2) - pow(met.Pt(), 2) * A);
                     cout << D << endl;
                     // first solution
                     if (D >= 0) {
                        double pz = jt->Pz() * B / A + sqrt(D) / A;
                        double E = sqrt(px * px + py * py + pz * pz);
                        TLorentzVector neutrino1(px, py, pz, E);
                        h_Mbln_reco->Fill(((*it) + (*jt) + neutrino1).M(), EventWeight);
                     }
                     // second solution
                     if (D > 0) {
                        double pz = jt->Pz() * B / A - sqrt(D) / A;
                        double E = sqrt(px * px + py * py + pz * pz);
                        TLorentzVector neutrino2(px, py, pz, E);
                        h_Mbln_reco->Fill(((*it) + (*jt) + neutrino2).M(), EventWeight);
                     }
                  }
               }
            }
         }
      }
   }
   //////////////////////////////

   return kTRUE;
}

void MyAnalysis::SlaveTerminate() {
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void MyAnalysis::Terminate() {
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

}
Revision:	1.12
Committed:	Thu Mar 1 19:52:31 2012 UTC (13 years, 2 months ago) by csander
Content type:	text/plain
Branch:	MAIN
Changes since 1.11:	+38 -15 lines
Log Message:	now with top mass determination
#	User	Rev	Content
1	csander	1.1	#define MyAnalysis_cxx
2			// The class definition in MyAnalysis.h has been generated automatically
3			// by the ROOT utility TTree::MakeSelector(). This class is derived
4			// from the ROOT class TSelector. For more information on the TSelector
5			// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
6
7			// The following methods are defined in this file:
8			// Begin(): called every time a loop on the tree starts,
9			// a convenient place to create your histograms.
10			// SlaveBegin(): called after Begin(), when on PROOF called only on the
11			// slave servers.
12			// Process(): called for each event, in this function you decide what
13			// to read and fill your histograms.
14			// SlaveTerminate: called at the end of the loop on the tree, when on PROOF
15			// called only on the slave servers.
16			// Terminate(): called at the end of the loop on the tree,
17			// a convenient place to draw/fit your histograms.
18			//
19			// To use this file, try the following session on your Tree T:
20			//
21			// Root > T->Process("MyAnalysis.C")
22			// Root > T->Process("MyAnalysis.C","some options")
23			// Root > T->Process("MyAnalysis.C+")
24			//
25
26			#include "MyAnalysis.h"
27	csander	1.3	#include <iostream>
28	csander	1.2	#include <TH1F.h>
29	csander	1.9	#include <TLatex.h>
30	csander	1.1
31	csander	1.3	using namespace std;
32
33			void MyAnalysis::BuildEvent() {
34			Jets.clear();
35			for (int i = 0; i < NJet; ++i) {
36			MyJet jet(Jet_Px[i], Jet_Py[i], Jet_Pz[i], Jet_E[i]);
37			jet.SetBTagDiscriminator(Jet_btag[i]);
38			Jets.push_back(jet);
39			}
40			Muons.clear();
41			for (int i = 0; i < NMuon; ++i) {
42			MyMuon muon(Muon_Px[i], Muon_Py[i], Muon_Pz[i], Muon_E[i]);
43			muon.SetIsolation(Muon_Iso[i]);
44	csander	1.4	muon.SetCharge(Muon_Charge[i]);
45	csander	1.3	Muons.push_back(muon);
46			}
47	csander	1.4	Electrons.clear();
48			for (int i = 0; i < NElectron; ++i) {
49			MyElectron electron(Electron_Px[i], Electron_Py[i], Electron_Pz[i], Electron_E[i]);
50			electron.SetIsolation(Electron_Iso[i]);
51			electron.SetCharge(Electron_Charge[i]);
52			Electrons.push_back(electron);
53			}
54			Photons.clear();
55			for (int i = 0; i < NPhoton; ++i) {
56			MyPhoton photon(Photon_Px[i], Photon_Py[i], Photon_Pz[i], Photon_E[i]);
57			photon.SetIsolation(Photon_Iso[i]);
58			Photons.push_back(photon);
59			}
60	csander	1.6	hadB.SetXYZM(MChadronicBottom_px, MChadronicBottom_py, MChadronicBottom_pz, 4.8);
61			lepB.SetXYZM(MCleptonicBottom_px, MCleptonicBottom_py, MCleptonicBottom_pz, 4.8);
62			hadWq.SetXYZM(MChadronicWDecayQuark_px, MChadronicWDecayQuark_py, MChadronicWDecayQuark_pz, 0.0);
63			hadWqb.SetXYZM(MChadronicWDecayQuarkBar_px, MChadronicWDecayQuarkBar_py, MChadronicWDecayQuarkBar_pz, 0.0);
64			lepWl.SetXYZM(MClepton_px, MClepton_py, MClepton_pz, 0.0);
65			lepWn.SetXYZM(MCneutrino_px, MCneutrino_py, MCneutrino_pz, 0.0);
66	csander	1.8	met.SetXYZM(MET_px, MET_py, 0., 0.);
67	csander	1.9
68	csander	1.10	if (EventWeight < 1.e-5)
69	csander	1.9	EventWeight = 1.;
70			EventWeight *= weight_factor;
71
72	csander	1.3	}
73
74	csander	1.2	void MyAnalysis::Begin(TTree * /tree/) {
75	csander	1.1	// The Begin() function is called at the start of the query.
76			// When running with PROOF Begin() is only called on the client.
77			// The tree argument is deprecated (on PROOF 0 is passed).
78
79			TString option = GetOption();
80
81			}
82
83	csander	1.2	void MyAnalysis::SlaveBegin(TTree * /tree/) {
84	csander	1.1	// The SlaveBegin() function is called after the Begin() function.
85			// When running with PROOF SlaveBegin() is called on each slave server.
86			// The tree argument is deprecated (on PROOF 0 is passed).
87
88			TString option = GetOption();
89
90	csander	1.7	h_Mmumu = new TH1F("Mmumu", "Invariant di-muon mass", 50, 0, 150);
91	csander	1.4	h_Mmumu->SetXTitle("m_{#mu#mu}");
92			h_Mmumu->Sumw2();
93
94	csander	1.10	h_Mbqqb_mc = new TH1F("Mbqqb_mc", "Invariant hadronic top mass", 50, 170, 175);
95	csander	1.7	h_Mbqqb_mc->SetXTitle("m_{hadronic top}");
96			h_Mbqqb_mc->Sumw2();
97
98	csander	1.10	h_Mbln_mc = new TH1F("Mbln_mc", "Invariant leptonic top mass", 50, 170, 175);
99	csander	1.7	h_Mbln_mc->SetXTitle("m_{leptonic top}");
100			h_Mbln_mc->Sumw2();
101
102	csander	1.10	h_Mbqqb_reco = new TH1F("Mbqqb_reco", "Invariant hadronic top mass", 50, 0, 500);
103	csander	1.7	h_Mbqqb_reco->SetXTitle("m_{hadronic top}");
104			h_Mbqqb_reco->Sumw2();
105
106	csander	1.10	h_Mbln_reco = new TH1F("Mbln_reco", "Invariant transverse leptonic top mass", 50, 0, 500);
107	csander	1.8	h_Mbln_reco->SetXTitle("m_{T, leptonic top}");
108			h_Mbln_reco->Sumw2();
109
110	csander	1.9	h_NJet = new TH1F("NJet", "Number of jets", 7, 0, 7);
111			h_NJet->SetXTitle("NJet");
112			h_NJet->Sumw2();
113
114			h_NBJet = new TH1F("NBJet", "Number of b-jets", 7, 0, 7);
115			h_NBJet->SetXTitle("NBJet");
116			h_NBJet->Sumw2();
117
118	csander	1.11	h_NMuon = new TH1F("NMuon", "Number of muons", 7, 0, 7);
119			h_NMuon->SetXTitle("NMuon");
120			h_NMuon->Sumw2();
121
122	csander	1.10	h_Jet1_Pt = new TH1F("Jet1Pt", "Pt of jet1", 50, 0, 250);
123			h_Jet1_Pt->SetXTitle("p_{T} [GeV]");
124	csander	1.9	h_Jet1_Pt->Sumw2();
125
126	csander	1.10	h_Jet2_Pt = new TH1F("Jet2Pt", "Pt of jet2", 50, 0, 250);
127			h_Jet2_Pt->SetXTitle("p_{T} [GeV]");
128	csander	1.9	h_Jet2_Pt->Sumw2();
129
130	csander	1.10	h_Jet3_Pt = new TH1F("Jet3Pt", "Pt of jet3", 50, 0, 250);
131			h_Jet3_Pt->SetXTitle("p_{T} [GeV]");
132	csander	1.9	h_Jet3_Pt->Sumw2();
133
134	csander	1.10	h_Jet1_Eta = new TH1F("Jet1Eta", "#eta of jet1", 50, -4, 4);
135	csander	1.9	h_Jet1_Eta->SetXTitle("#eta");
136			h_Jet1_Eta->Sumw2();
137
138	csander	1.10	h_Jet2_Eta = new TH1F("Jet2Eta", "#eta of jet2", 50, -4, 4);
139	csander	1.9	h_Jet2_Eta->SetXTitle("#eta");
140			h_Jet2_Eta->Sumw2();
141
142	csander	1.10	h_Jet3_Eta = new TH1F("Jet3Eta", "#eta of jet3", 50, -4, 4);
143	csander	1.9	h_Jet3_Eta->SetXTitle("#eta");
144			h_Jet3_Eta->Sumw2();
145
146	csander	1.10	h_BJet1_Pt = new TH1F("BJet1Pt", "Pt of b-jet1", 50, 0, 250);
147			h_BJet1_Pt->SetXTitle("p_{T} [GeV]");
148	csander	1.9	h_BJet1_Pt->Sumw2();
149
150	csander	1.10	h_BJet2_Pt = new TH1F("BJet2Pt", "Pt of b-jet2", 50, 0, 250);
151			h_BJet2_Pt->SetXTitle("p_{T} [GeV]");
152	csander	1.9	h_BJet2_Pt->Sumw2();
153
154	csander	1.10	h_BJet1_Eta = new TH1F("BJet1Eta", "#eta of b-jet1", 50, -4, 4);
155	csander	1.9	h_BJet1_Eta->SetXTitle("#eta");
156			h_BJet1_Eta->Sumw2();
157
158	csander	1.10	h_BJet2_Eta = new TH1F("BJet2Eta", "#eta of b-jet2", 50, -4, 4);
159	csander	1.9	h_BJet2_Eta->SetXTitle("#eta");
160			h_BJet2_Eta->Sumw2();
161
162	csander	1.10	h_Muon1_Pt = new TH1F("Muon1Pt", "Pt of muon1", 50, 0, 250);
163			h_Muon1_Pt->SetXTitle("p_{T} [GeV]");
164	csander	1.9	h_Muon1_Pt->Sumw2();
165
166	csander	1.10	h_Muon2_Pt = new TH1F("Muon2Pt", "Pt of muon2", 50, 0, 250);
167			h_Muon2_Pt->SetXTitle("p_{T} [GeV]");
168	csander	1.9	h_Muon2_Pt->Sumw2();
169
170	csander	1.10	h_Muon1_Eta = new TH1F("Muon1Eta", "#eta of muon1", 50, -4, 4);
171	csander	1.9	h_Muon1_Eta->SetXTitle("#eta");
172			h_Muon1_Eta->Sumw2();
173
174	csander	1.10	h_Muon2_Eta = new TH1F("Muon2Eta", "#eta of muon2", 50, -4, 4);
175	csander	1.9	h_Muon2_Eta->SetXTitle("#eta");
176			h_Muon2_Eta->Sumw2();
177
178			h_Muon1_Iso = new TH1F("Muon1Iso", "Isolation of muon1", 50, 0, 25);
179	csander	1.10	h_Muon1_Iso->SetXTitle("Isolation [GeV]");
180	csander	1.9	h_Muon1_Iso->Sumw2();
181
182			h_Muon2_Iso = new TH1F("Muon2Iso", "Isolation of muon2", 50, 0, 25);
183	csander	1.10	h_Muon2_Iso->SetXTitle("Isolation [GeV]");
184	csander	1.9	h_Muon2_Iso->Sumw2();
185
186			h_MET = new TH1F("MET", "MET", 30, 0, 300.);
187			h_MET->SetXTitle("MET");
188			h_MET->Sumw2();
189
190	csander	1.10	h_minDeltaPhi_MET_Muon = new TH1F("minDeltaPhi_MET_Muon", "min#Delta#phi(MET,muon)", 50, 0., TMath::Pi());
191	csander	1.9	h_minDeltaPhi_MET_Muon->SetXTitle("min#Delta#phi [rad]");
192			h_minDeltaPhi_MET_Muon->Sumw2();
193
194	csander	1.10	h_minDeltaPhi_MET_BJet = new TH1F("minDeltaPhi_MET_BJet", "min#Delta#phi(MET,b-jet)", 50, 0., TMath::Pi());
195	csander	1.9	h_minDeltaPhi_MET_BJet->SetXTitle("min#Delta#phi [rad]");
196			h_minDeltaPhi_MET_BJet->Sumw2();
197
198	csander	1.10	h_selectedEvents_Muon1_Pt = new TH1F("selectedEvents_Muon1_Pt", "Pt of muon1 (selected)", 50, 0, 250);
199			h_selectedEvents_Muon1_Pt->SetXTitle("p_{T} [GeV]");
200	csander	1.9	h_selectedEvents_Muon1_Pt->Sumw2();
201
202			h_selectedEvents_triggered_Muon1_Pt = new TH1F("selectedEvents_triggered_Muon1_Pt",
203	csander	1.10	"Pt of muon1 (selected and triggered)", 50, 0, 250);
204			h_selectedEvents_triggered_Muon1_Pt->SetXTitle("p_{T} [GeV]");
205	csander	1.9	h_selectedEvents_triggered_Muon1_Pt->Sumw2();
206
207	csander	1.1	}
208
209	csander	1.2	Bool_t MyAnalysis::Process(Long64_t entry) {
210	csander	1.1	// The Process() function is called for each entry in the tree (or possibly
211			// keyed object in the case of PROOF) to be processed. The entry argument
212			// specifies which entry in the currently loaded tree is to be processed.
213			// It can be passed to either MyAnalysis::GetEntry() or TBranch::GetEntry()
214			// to read either all or the required parts of the data. When processing
215			// keyed objects with PROOF, the object is already loaded and is available
216			// via the fObject pointer.
217			//
218			// This function should contain the "body" of the analysis. It can contain
219			// simple or elaborate selection criteria, run algorithms on the data
220			// of the event and typically fill histograms.
221			//
222			// The processing can be stopped by calling Abort().
223			//
224			// Use fStatus to set the return value of TTree::Process().
225			//
226			// The return value is currently not used.
227
228	csander	1.3	++TotalEvents;
229
230			GetEntry(entry);
231
232	csander	1.4	if (TotalEvents % 10000 == 0)
233			cout << "Next event -----> " << TotalEvents << endl;
234
235	csander	1.3	BuildEvent();
236
237	csander	1.7	// cout << "Jets: " << endl;
238			// for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
239			// cout << "pt, eta, phi, btag: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->IsBTagged()
240			// << endl;
241			// }
242			// cout << "Muons: " << endl;
243			// for (vector<MyMuon>::iterator it = Muons.begin(); it != Muons.end(); ++it) {
244			// cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
245			// << it->GetIsolation() << ", " << it->GetCharge() << endl;
246			// }
247			// cout << "Electrons: " << endl;
248			// for (vector<MyElectron>::iterator it = Electrons.begin(); it != Electrons.end(); ++it) {
249			// cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
250			// << it->GetIsolation() << ", " << it->GetCharge() << endl;
251			// }
252			// cout << "Photons: " << endl;
253			// for (vector<MyPhoton>::iterator it = Photons.begin(); it != Photons.end(); ++it) {
254			// cout << "pt, eta, phi, iso: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->GetIsolation()
255			// << endl;
256			// }
257
258	csander	1.11	//////////////////////////////
259			// Exercise XX: Invariant Di-Muon mass
260	csander	1.10	if (NMuon > 1) {
261			if (Muons.at(0).GetCharge() * Muons.at(1).GetCharge() < 1 && Muons.at(0).IsIsolated() && Muons.at(1).IsIsolated()) {
262			h_Mmumu->Fill((Muons.at(0) + Muons.at(1)).M(), EventWeight);
263			}
264	csander	1.5	}
265	csander	1.11	//////////////////////////////
266
267	csander	1.4
268	csander	1.11	//////////////////////////////
269			// Exercise XX: MC / data comparisons
270	csander	1.7
271	csander	1.11	// require at least one isolapted muon
272			// +trigger bit
273	csander	1.10	if (triggerIsoMu24 && NMuon > 0) {
274			if (Muons.at(0).Pt() > 24. && Muons.at(0).IsIsolated()) {
275
276			int N_BJet = 0;
277			int N_Jet = 0;
278			for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
279			++N_Jet;
280			if (N_Jet == 1) {
281			h_Jet1_Pt->Fill(it->Pt(), EventWeight);
282			h_Jet1_Eta->Fill(it->Eta(), EventWeight);
283			} else if (N_Jet == 2) {
284			h_Jet2_Pt->Fill(it->Pt(), EventWeight);
285			h_Jet2_Eta->Fill(it->Eta(), EventWeight);
286			} else if (N_Jet == 3) {
287			h_Jet3_Pt->Fill(it->Pt(), EventWeight);
288			h_Jet3_Eta->Fill(it->Eta(), EventWeight);
289			}
290			if (it->IsBTagged()) {
291			++N_BJet;
292			if (N_BJet == 1) {
293			h_BJet1_Pt->Fill(it->Pt(), EventWeight);
294			h_BJet1_Eta->Fill(it->Eta(), EventWeight);
295			} else if (N_BJet == 2) {
296			h_BJet2_Pt->Fill(it->Pt(), EventWeight);
297			h_BJet2_Eta->Fill(it->Eta(), EventWeight);
298			}
299			}
300			}
301			h_NBJet->Fill(N_BJet, EventWeight);
302	csander	1.11	h_NJet->Fill(N_Jet, EventWeight);
303	csander	1.10
304	csander	1.11	int N_Muon = 0;
305	csander	1.10	for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
306			if (jt->IsIsolated()) {
307	csander	1.11	++N_Muon;
308			if (N_Muon == 1) {
309	csander	1.10	h_Muon1_Pt->Fill(jt->Pt(), EventWeight);
310			h_Muon1_Eta->Fill(jt->Eta(), EventWeight);
311			h_Muon1_Iso->Fill(jt->GetIsolation(), EventWeight);
312	csander	1.11	} else if (N_Muon == 2) {
313	csander	1.10	h_Muon2_Pt->Fill(jt->Pt(), EventWeight);
314			h_Muon2_Eta->Fill(jt->Eta(), EventWeight);
315			h_Muon2_Iso->Fill(jt->GetIsolation(), EventWeight);
316			}
317			}
318	csander	1.7	}
319	csander	1.11	h_NMuon->Fill(N_Muon, EventWeight);
320	csander	1.9
321	csander	1.10	h_MET->Fill(met.Pt(), EventWeight);
322	csander	1.9
323	csander	1.10	float dPhiMuon = 999;
324			for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
325			float tmp;
326			if (jt->Pt() > 24.) {
327			tmp = fabs(jt->DeltaPhi(met));
328			if (tmp < dPhiMuon)
329			dPhiMuon = tmp;
330			}
331			}
332			if (dPhiMuon < 10.)
333			h_minDeltaPhi_MET_Muon->Fill(dPhiMuon, EventWeight);
334	csander	1.9
335	csander	1.10	float dPhiBJet = 999;
336			for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
337			float tmp;
338			if (it->IsBTagged()) {
339			tmp = fabs(it->DeltaPhi(met));
340			if (tmp < dPhiBJet)
341			dPhiBJet = tmp;
342			}
343			}
344			if (dPhiBJet < 10.)
345			h_minDeltaPhi_MET_BJet->Fill(dPhiBJet, EventWeight);
346	csander	1.1
347	csander	1.10	} // end of muon requirement
348			} //end of trigger requirement
349	csander	1.11	//////////////////////////////
350	csander	1.9
351	csander	1.11	//////////////////////////////
352			// Exercise XXa: Trigger turn om
353	csander	1.9	if (Muons.size() == 1) {
354	csander	1.11	if (Muons.at(0).IsIsolated()) {
355			h_selectedEvents_Muon1_Pt->Fill(Muons.at(0).Pt(), EventWeight);
356			if (triggerIsoMu24)
357			h_selectedEvents_triggered_Muon1_Pt->Fill(Muons.at(0).Pt(), EventWeight);
358			}
359			}
360			//////////////////////////////
361
362			//////////////////////////////
363			// Exercise XXb: Acceptance and Trigger Efficiencies
364			GeneratedEvents += EventWeight;
365	csander	1.12	bool IsSelected = false;
366	csander	1.11	if (NMuon > 0) {
367			if (Muons.at(0).Pt() > 24. && Muons.at(0).IsIsolated()) {
368
369			int N_BJet = 0;
370			for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
371			if (it->IsBTagged()) {
372			++N_BJet;
373			}
374			}
375			if (N_BJet > 1) {
376			SelectedEvents += EventWeight;
377	csander	1.12	if (triggerIsoMu24) {
378	csander	1.11	SelectedEvents_triggered += EventWeight;
379	csander	1.12	IsSelected = true;
380			}
381	csander	1.11	}
382
383			}
384	csander	1.9	}
385	csander	1.11	//////////////////////////////
386	csander	1.9
387	csander	1.11	//////////////////////////////
388			// Exercixe XXa: Top mass (MC Truth)
389			h_Mbqqb_mc->Fill((hadB + hadWq + hadWqb).M(), EventWeight);
390			h_Mbln_mc->Fill((lepB + lepWl + lepWn).M(), EventWeight);
391			//////////////////////////////
392
393			//////////////////////////////
394			// Exercixe XXb: Top mass (hadroninc top decay) (reco)
395			// Exercixe XXc: Top mass (leptonic top decay) (reco)
396	csander	1.12	if (IsSelected) {
397	csander	1.11
398	csander	1.12	for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
399			if (it->IsBTagged()) {
400			for (vector<MyJet>::iterator jt = Jets.begin(); jt != Jets.end(); ++jt) {
401			if (jt != it && !(jt->IsBTagged())) {
402			for (vector<MyJet>::iterator kt = jt; kt != Jets.end(); ++kt) {
403			if (kt != it && kt != jt && !(kt->IsBTagged())) {
404			if (((jt) + (kt)).M() > 65 && ((jt) + (kt)).M() < 95)
405	csander	1.11	h_Mbqqb_reco->Fill(((it) + (jt) + (*kt)).M(), EventWeight);
406			}
407			}
408			}
409			}
410	csander	1.12	}
411	csander	1.11
412	csander	1.12	for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
413			if (it->IsBTagged()) {
414			for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
415			if (jt->IsIsolated() && jt->Pt() > 24.) {
416			double px = met.Px();
417			double py = met.Py();
418			double mW = 80.4;
419			double A = pow(jt->E(), 2) - pow(jt->Pz(), 2);
420			double B = mW * mW / 2. + (jt->Px() * px + jt->Py() * py);
421			double D = pow(jt->E(), 2) * (pow(B, 2) - pow(met.Pt(), 2) * A);
422			cout << D << endl;
423			// first solution
424			if (D >= 0) {
425			double pz = jt->Pz() * B / A + sqrt(D) / A;
426			double E = sqrt(px * px + py * py + pz * pz);
427			TLorentzVector neutrino1(px, py, pz, E);
428			h_Mbln_reco->Fill(((it) + (jt) + neutrino1).M(), EventWeight);
429			}
430			// second solution
431			if (D > 0) {
432			double pz = jt->Pz() * B / A - sqrt(D) / A;
433			double E = sqrt(px * px + py * py + pz * pz);
434			TLorentzVector neutrino2(px, py, pz, E);
435			h_Mbln_reco->Fill(((it) + (jt) + neutrino2).M(), EventWeight);
436			}
437	csander	1.11	}
438			}
439			}
440			}
441			}
442			}
443			//////////////////////////////
444	csander	1.8
445	csander	1.1	return kTRUE;
446			}
447
448	csander	1.2	void MyAnalysis::SlaveTerminate() {
449	csander	1.1	// The SlaveTerminate() function is called after all entries or objects
450			// have been processed. When running with PROOF SlaveTerminate() is called
451			// on each slave server.
452
453			}
454
455	csander	1.2	void MyAnalysis::Terminate() {
456	csander	1.1	// The Terminate() function is the last function to be called during
457			// a query. It always runs on the client, it can be used to present
458			// the results graphically or save the results to file.
459
460			}