trigger_ef.cpp

#include <cstdio> 
#include <cstdlib>
#include <iostream>
#include <memory> 
#include <cmath> 
#include <vector>
#include <fstream>
#include <string> 

#include "TFile.h"     
#include "TTree.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "Math/Vector4D.h"
#include "TF1.h"
#include "TChain.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TLine.h"
#include "THStack.h"
#include "TH2D.h"


using ROOT::Math::PtEtaPhiMVector;

using namespace std;

// const double electron_mass = 0.000511; // GeV
const float isoCut = 0.15;

void fillHistogramFromTree(TTree* tree, const string& sampleName){
    
    TH1D hDenLeading("hDenLeading", "Denominator Leading", 24, -2.5, 2.5);
    TH1D hNumLeading("hNumLeading", "Numerator Leading", 24, -2.5, 2.5);

    TH1D hDenSub("hDenSub", "Denominator Subleading", 24, -2.5, 2.5);
    TH1D hNumSub("hNumSub", "Numerator Subleading", 24, -2.5, 2.5);

    TH2D hDen2D("hDen2D", "Denominator 2D;Leading #eta;Subleading #eta", 
                12, -2.5, 2.5, 12, -2.5, 2.5);
    TH2D hNum2D("hNum2D", "Numerator 2D;Leading #eta;Subleading #eta", 
                12, -2.5, 2.5, 12, -2.5, 2.5);

    tree->SetBranchStatus("*", false);
    tree->SetBranchStatus("Electron_charge", true);
    tree->SetBranchStatus("Electron_cutBased", true);
    tree->SetBranchStatus("Electron_pfRelIso03_all", true);
    tree->SetBranchStatus("nElectron", true);
    tree->SetBranchStatus("Electron_phi", true);
    tree->SetBranchStatus("Electron_eta", true);
    tree->SetBranchStatus("Electron_pt", true);
    tree->SetBranchStatus("Electron_mass", true);
    tree->SetBranchStatus("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_DZ", true);
    tree->SetBranchStatus("HLT_Ele32_WPTight_Gsf", true);

    Int_t Electron_charge[100];
    UInt_t nElectron;
    Float_t Electron_phi[100], Electron_eta[100], Electron_pt[100], Electron_pfRelIso03_all[100], Electron_mass[100];
    Int_t Electron_cutBased[100];
    Bool_t HLT_Ele23_Ele12, HLT_Ele32_WPTight;

    tree->SetBranchAddress("Electron_charge", &Electron_charge);
    tree->SetBranchAddress("nElectron", &nElectron);
    tree->SetBranchAddress("Electron_mass", &Electron_mass);
    tree->SetBranchAddress("Electron_pt", &Electron_pt);
    tree->SetBranchAddress("Electron_eta", &Electron_eta);
    tree->SetBranchAddress("Electron_phi", &Electron_phi);
    tree->SetBranchAddress("Electron_cutBased", &Electron_cutBased);
    tree->SetBranchAddress("Electron_pfRelIso03_all", &Electron_pfRelIso03_all);
    tree->SetBranchAddress("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_DZ", &HLT_Ele23_Ele12);
    tree->SetBranchAddress("HLT_Ele32_WPTight_Gsf", &HLT_Ele32_WPTight);

    for (int iEntry = 0; tree->LoadTree(iEntry) >= 0; ++iEntry){    
        
        tree->GetEntry(iEntry);
    
        if (!HLT_Ele32_WPTight) continue;

        vector<int> goodIndices;

        for (UInt_t i = 0; i < nElectron; ++i){
            Bool_t isGoodElectron = (Electron_cutBased[i] >= 3) && (Electron_pfRelIso03_all[i] < isoCut);
            if (isGoodElectron) {
                goodIndices.push_back(i);
            }
        }

        if (goodIndices.size() > 1){

            sort(goodIndices.begin(), goodIndices.end(), [&](int a, int b) {
                return Electron_pt[a] > Electron_pt[b];
            });
            
            int leading = goodIndices[0];
            int subleading = goodIndices[1];

            if (Electron_charge[leading] * Electron_charge[subleading] < 0) {

                if (Electron_pt[leading] > 35 && Electron_pt[subleading] > 35 && 
                    fabs(Electron_eta[leading]) < 2.5 && fabs(Electron_eta[subleading]) < 2.5 && 
                    !(fabs(Electron_eta[leading]) > 1.4442 && fabs(Electron_eta[leading]) < 1.566) &&
                    !(fabs(Electron_eta[subleading]) > 1.4442 && fabs(Electron_eta[subleading]) < 1.566)) {

                    PtEtaPhiMVector electron1(Electron_pt[leading], Electron_eta[leading], Electron_phi[leading], Electron_mass[leading]);
                    PtEtaPhiMVector electron2(Electron_pt[subleading], Electron_eta[subleading], Electron_phi[subleading], Electron_mass[subleading]);

                    auto dilepton = electron1 + electron2;
                     
                    if (dilepton.M() < 100 && dilepton.M() > 80){

                        hDenLeading.Fill(Electron_eta[leading]);
                        hDenSub.Fill(Electron_eta[subleading]);

                        hDen2D.Fill(Electron_eta[leading], Electron_eta[subleading]);
                        if (HLT_Ele23_Ele12){
                            hNumLeading.Fill(Electron_eta[leading]);
                            hNumSub.Fill(Electron_eta[subleading]);

                            hNum2D.Fill(Electron_eta[leading], Electron_eta[subleading]);
                        }
                    }
                }
            }    
        }
    }
    
    TH1D hEffLeading = hNumLeading;
    hEffLeading.Divide(&hDenLeading);
    hEffLeading.SetTitle((sampleName + " Trigger Efficiency; #eta; Efficiency").c_str());
    hEffLeading.GetXaxis()->SetTitle("#eta");
    hEffLeading.GetYaxis()->SetTitle("Efficiency");
    hEffLeading.GetYaxis()->SetTitleOffset(1.3);
    hEffLeading.SetLineColor(kRed);
    hEffLeading.SetStats(0);

    TH1D hEffSub = hNumSub;
    hEffSub.Divide(&hDenSub);
    hEffSub.SetTitle((sampleName + " Subleading Electron Efficiency; #eta; Efficiency").c_str());
    hEffSub.GetXaxis()->SetTitle("#eta");
    hEffSub.GetYaxis()->SetTitle("Efficiency");
    hEffSub.GetYaxis()->SetTitleOffset(1.3);
    hEffSub.SetLineColor(kBlue);
    hEffSub.SetStats(0);

    TH2D hEff2D = hNum2D;
    hEff2D.Divide(&hDen2D);
    hEff2D.SetTitle((sampleName + " Trigger Efficiency 2D;Leading #eta;Subleading #eta").c_str());
    hEff2D.GetXaxis()->SetTitle("Leading electron #eta");
    hEff2D.GetYaxis()->SetTitle("Subleading electron #eta");
    hEff2D.SetStats(0); 

    double totalDen = hDenLeading.Integral()+hDenSub.Integral();
    double totalNum = hNumLeading.Integral()+hNumSub.Integral();
    double eff = (totalDen > 0) ? (totalNum / totalDen) : 0.0;

    cout << "Overall efficiency for " << sampleName << " = " << eff << " (" 
         << totalNum << "/" << totalDen << ")" << endl;

    TCanvas* c2 = new TCanvas("c2", "2D Efficiency", 800, 700);
    c2->SetTicks(1,1);
    hEff2D.Draw("COLZ TEXT"); 
    c2->SaveAs((sampleName + "_efficiency2D.png").c_str());
    delete c2;

    TCanvas* c = new TCanvas("c", "Trigger efficiency", 800, 600);
    c->SetTicks(1, 1);
    hEffLeading.Draw("HIST"); 
    hEffSub.Draw("HIST SAME");
    TLegend *l = new TLegend(0.35, 0.14, 0.65, 0.27);
    l->AddEntry(&hEffLeading, "Leading Electron", "l");
    l->AddEntry(&hEffSub, "Subleading Electron", "l");
    l->SetFillStyle(0);      
    l->Draw("same");
 
    c->SaveAs((sampleName + "_efficiency.png").c_str());
    delete c;

}

TChain* createChainFromFileList(const char* txtFileName, int maxFiles = -1) {
    TChain* chain = new TChain("Events");
    ifstream infile(txtFileName);
    string line;
    int count = 0;
    while (getline(infile, line)) {
        if (maxFiles != -1 && count >= maxFiles) break;
        chain->Add(line.c_str());
        count++;
    }
    return chain;
}

int main(){

    TChain* tree = createChainFromFileList("files/CMS_Run2016H_DoubleEG_NANOAOD_UL2016_MiniAODv2_NanoAODv9-v1_100000_file_index.txt", 4);  
    TChain* mcTree = createChainFromFileList("files/CMS_mc_RunIISummer20UL16NanoAODv9_DYJetsToLL_M-50_TuneCP5_13TeV-amcatnloFXFX-pythia8_NANOAODSIM_106X_mcRun2_asymptotic_v17-v1_30000_file_index.txt", 4);  
    
    fillHistogramFromTree(tree, "Real Data");

    fillHistogramFromTree(mcTree, "MC DY");    

    return 0;
}