TrigGlobEffCorrExample3e.cxx File Reference

#include <TFile.h>
#include <TError.h>
#include "AthAnalysisBaseComps/AthAnalysisHelper.h"
#include "POOLRootAccess/TEvent.h"
#include "AsgMessaging/MessageCheck.h"
#include "AsgTools/AnaToolHandle.h"
#include "EgammaAnalysisInterfaces/IAsgElectronEfficiencyCorrectionTool.h"
#include "MuonAnalysisInterfaces/IMuonTriggerScaleFactors.h"
#include "TriggerAnalysisInterfaces/ITrigGlobalEfficiencyCorrectionTool.h"
#include "xAODEventInfo/EventInfo.h"
#include "xAODEgamma/ElectronContainer.h"
#include "xAODMuon/MuonContainer.h"
#include "PATCore/PATCoreEnums.h"
#include "AthContainers/AuxElement.h"
#include "AthContainers/ConstAccessor.h"
#include "AthContainers/Decorator.h"
#include <sstream>
#include <random>
#include <vector>
#include <array>

Go to the source code of this file.

Macros
#define	MSGSOURCE   "Example 3e"

Functions
int	main (int argc, char *argv[])

Macro Definition Documentation

MSGSOURCE

#define MSGSOURCE   "Example 3e"

Definition at line 84 of file TrigGlobEffCorrExample3e.cxx.

Function Documentation

main()

int main	(	int	argc,
		char *	argv[] )

For property 'ElectronEfficiencyTools':

For property 'ElectronScaleFactorTools':

For property 'ListOfLegsPerTool':

For property 'ListOfTagsPerTool':

For property 'ElectronLegsPerTag':

To tag electrons according to the PID criteria they fulfil

To emulate PID selection (90% loose-to-medium/medium-to-tight eff.)

RAII on-the-fly creation of electron CP tools:

<key in map file>, <PID WP>, <iso WP> Single-electron trigger: only electrons tagged 'PID2' (TightLH+iso)

Electron-muon trigger: electrons tagged 'PID2' or 'PID1' (MediumLH)

Dielectron trigger: all electrons (tagged or not)

one instance per trigger leg x working point

two instances: 0 -> MC efficiencies, 1 -> SFs

Safer to retrieve the name from the final ToolHandle, it might be prefixed (by the parent tool name) when the handle is copied

For property 'MuonTools':

Special character '?' indicates that the decorated value is to be suffixed to the name (=> 'PID1' for medium, 'PID2' for tight)

Uniform random run number generation spanning the target dataset. In real life, use the PileupReweightingTool instead!

2016 periods A-L

Get a random run number, and decorate the event info

electron must be above softest trigger threshold (e7 here

muon must be above softest trigger threshold (mu24 here)

Add 'PID' decorations to random electrons also count 'Tight' electrons above e26_xxx threshold and 'Medium' electrons above e7_xxx threshold

Set decorated value to 2 for TightLH+iso or 1 for MediumLH

Events must contain enough leptons to trigger

single-electron trigger

electron-muon

dielectron

Finally retrieve the global trigger scale factor

Definition at line 89 of file TrigGlobEffCorrExample3e.cxx.

{
    ANA_CHECK_SET_TYPE(bool)
    const char* filename = nullptr;
    bool debug = false, cmdline_error = false, toys = false;
    for(int i=1;i<argc;++i)
    {
        if(string(argv[i]) == "--debug") debug = true;
        else if(string(argv[i]) == "--toys") toys = true;
        else if(!filename && *argv[i]!='-') filename = argv[i];
        else cmdline_error = true;
    }
    if(!filename || cmdline_error)
    {
        Error(MSGSOURCE, "No file name received!");
        Error(MSGSOURCE, "  Usage: %s [--debug] [--toys] [DxAOD file name]", argv[0]);
        return 1;
    }
    #ifdef XAOD_STANDALONE
        xAOD::Init(MSGSOURCE).ignore();
        TFile* file = TFile::Open(filename, "READ");
        if(!file)
        {
            Error(MSGSOURCE, "Unable to open file!");
            return 2;
        }
        xAOD::TEvent event(xAOD::TEvent::kClassAccess);
        xAOD::TStore store;
        StatusCode::enableFailure();
    #else
       IAppMgrUI* app = POOL::Init();
       POOL::TEvent event(POOL::TEvent::kClassAccess);
       TString file(filename);
    #endif
    event.readFrom(file).ignore();
    Long64_t entries = event.getEntries();
    Info(MSGSOURCE, "Number of events in the file: %lli", entries);
 
    /* ********************************************************************** */
    
    Info(MSGSOURCE, "Configuring the electron CP tools");
    ToolHandleArray<IAsgElectronEfficiencyCorrectionTool> electronEffTools;
    ToolHandleArray<IAsgElectronEfficiencyCorrectionTool> electronSFTools; 
    std::map<std::string,std::string> legsPerTool;
    std::map<std::string,std::string> tagsPerTool;
    std::map<std::string,std::string> legsPerTag;
    static const SG::Decorator<char> dec_pid("PID");
    std::bernoulli_distribution bernoulliPdf(0.9);

    vector<asg::AnaToolHandle<IAsgElectronEfficiencyCorrectionTool>> factory;
    enum{ cLEGS, cTAG, cKEY, cPID, cISO };
    std::vector<std::array<std::string,5> > toolConfigs = {
        {"e26_lhtight_nod0_ivarloose_OR_e60_lhmedium_nod0_OR_e140_lhloose_nod0", "PID2", 
            "2016_e26_lhtight_nod0_ivarloose_OR_e60_lhmedium_nod0_OR_e140_lhloose_nod0", "Tight", "GradientLoose"},
        {"e7_lhmedium_nod0", "PID1,PID2", "2016_e7_lhmedium_nod0", "Medium", ""},
        {"e17_lhvloose_nod0", "*,PID1,PID2", "2016_e17_lhvloose_nod0", "LooseBLayer", ""},
     };
 
    const char* mapPath = "ElectronEfficiencyCorrection/2015_2017/"
            "rel21.2/Moriond_February2018_v2/map6.txt";
    for(auto& cfg : toolConfigs) 
    for(int j=0;j<2;++j) 
    {
        string name = "AsgElectronEfficiencyCorrectionTool/"
                + ((j? "ElTrigEff_" : "ElTrigSF_")
                + std::to_string(factory.size()/2));
        auto t = factory.emplace(factory.end(), name);
        t->setProperty("MapFilePath", mapPath).ignore();
        t->setProperty("TriggerKey", string(j?"":"Eff_") + cfg[cKEY]).ignore();
        t->setProperty("IdKey", cfg[cPID]).ignore();
        t->setProperty("IsoKey", cfg[cISO]).ignore();
 
        t->setProperty("CorrelationModel", "TOTAL").ignore();
        t->setProperty("ForceDataType", (int)PATCore::ParticleDataType::Full).ignore();
        if(t->initialize() != StatusCode::SUCCESS)
        {
            Error(MSGSOURCE, "Unable to initialize the electron CP tool <%s>!",
                    t->name().c_str());
            return 3;
        }
        auto& handles = (j? electronSFTools : electronEffTools);
        handles.push_back(t->getHandle());
        name = handles[handles.size()-1].name();
        legsPerTool[name] = cfg[cLEGS];
        tagsPerTool[name] = cfg[cTAG];
        if(!j)
        {
            for(auto& tag : ::split_comma_delimited(cfg[cTAG]))
            {
                if(legsPerTag[tag]=="") legsPerTag[tag] = cfg[cLEGS];
                else legsPerTag[tag] += "," + cfg[cLEGS];
            }
        }
 
    }
 
    /* ********************************************************************** */
    
    Info(MSGSOURCE, "Configuring the muon CP tools");
    ToolHandleArray<CP::IMuonTriggerScaleFactors> muonTools;
    asg::AnaToolHandle<CP::IMuonTriggerScaleFactors> muonTool("CP::MuonTriggerScaleFactors/MuonTrigEff");
    muonTool.setProperty("MuonQuality", "Tight").ignore();
    muonTool.setProperty("useRel207", false).ignore();
    if(muonTool.initialize() != StatusCode::SUCCESS)
    {
        Error(MSGSOURCE, "Unable to initialize the muon CP tool!");
        return 3;
    }
    muonTools.push_back(muonTool.getHandle());
    
    /* ********************************************************************** */
    
    Info(MSGSOURCE, "Configuring the global trigger SF tool");
    asg::AnaToolHandle<ITrigGlobalEfficiencyCorrectionTool> myTool("TrigGlobalEfficiencyCorrectionTool/TrigGlobal");
    myTool.setProperty("ElectronEfficiencyTools", electronEffTools).ignore();
    myTool.setProperty("ElectronScaleFactorTools", electronSFTools).ignore();
    myTool.setProperty("MuonTools", muonTools).ignore();
    const char* triggers2016 = 
        "e26_lhtight_nod0_ivarloose_OR_e60_lhmedium_nod0_OR_e140_lhloose_nod0"
        "|| e7_lhmedium_nod0_mu24"
        "|| 2e17_lhvloose_nod0";
    myTool.setProperty("TriggerCombination2016", triggers2016).ignore();
    myTool.setProperty("LeptonTagDecorations", "PID?").ignore();
    myTool.setProperty("ListOfLegsPerTool", legsPerTool).ignore();
    myTool.setProperty("ListOfTagsPerTool", tagsPerTool).ignore();
    myTool.setProperty("ListOfLegsPerTag", legsPerTag).ignore();
 
    if(debug) myTool.setProperty("OutputLevel", MSG::DEBUG).ignore();
    if(toys) myTool.setProperty("NumberOfToys", 1000).ignore();
    if(myTool.initialize() != StatusCode::SUCCESS)
    {
        Error(MSGSOURCE, "Unable to initialize the TrigGlob tool!");
        return 3;
    }
    
    const unsigned periodRuns[] = {
        296939, 300345, 301912, 302737, 303638, 303943, 305291, 307124, 
        305359, 309311, 310015
    };
    std::uniform_int_distribution<unsigned> uniformPdf(0,
            sizeof(periodRuns)/sizeof(*periodRuns) - 1);
    std::default_random_engine randomEngine;
    
    static const SG::ConstAccessor<int> truthType("truthType");
    static const SG::ConstAccessor<int> truthOrigin("truthOrigin");
    
    /* ********************************************************************** */
    
    Info(MSGSOURCE, "Starting the event loop");
    unsigned errors = 0;
    double nSuitableEvents = 0., sumW = 0.;
    static const SG::Decorator<unsigned> RandomRunNumberDec("RandomRunNumber");
    for(Long64_t entry = 0; entry < entries; ++entry)
    {
        event.getEntry(entry);
        
        const xAOD::EventInfo* eventInfo = nullptr;
        event.retrieve(eventInfo,"EventInfo").ignore();
        unsigned runNumber = periodRuns[uniformPdf(randomEngine)];
        RandomRunNumberDec(*eventInfo) = runNumber;
 
        unsigned nTrig_e26 = 0, nTrig_e7 = 0, nTrig_e17 = 0;
        
        vector<const xAOD::Electron*> myTriggeringElectrons;
        const xAOD::ElectronContainer* electrons = nullptr;
        event.retrieve(electrons,"Electrons").ignore();
        for(auto electron : *electrons)
        {
            if(!electron->caloCluster()) continue;
            float eta = fabs(electron->caloCluster()->etaBE(2));
            float pt = electron->pt();
            if(pt<10e3f || eta>=2.47) continue;
            if(!truthType.isAvailable(*electron)) continue;
            if(!truthOrigin.isAvailable(*electron)) continue;
            int t = truthType(*electron), o = truthOrigin(*electron);
            if(t!=2 || !(o==10 || (o>=12 && o<=22) || o==43)) continue;
            if(pt < 7e3f) continue;
            if(pt >= 18e3f) ++nTrig_e17;
 
            myTriggeringElectrons.push_back(electron);
        }
 
        vector<const xAOD::Muon*> myTriggeringMuons;
        const xAOD::MuonContainer* muons = nullptr;
        event.retrieve(muons,"Muons").ignore();
        for(auto muon : *muons)
        {
            if(runNumber >= 324320) break; // delete line once all SFs available for 2017
            float pt = muon->pt();
            if(pt<10e3f || fabs(muon->eta())>=2.5) continue;
            auto mt = muon->muonType();
            if(mt!=xAOD::Muon::Combined && mt!=xAOD::Muon::MuonStandAlone) continue;
            auto& mtp = *(muon->primaryTrackParticle());
            if(!truthType.isAvailable(mtp)) continue;
            if(!truthOrigin.isAvailable(mtp)) continue;
            int t = truthType(mtp), o = truthOrigin(mtp);
            if(t!=6 || !(o==10 || (o>=12 && o<=22) || o==43)) continue;
            if(pt < 25.2e3f) continue;
 
            myTriggeringMuons.push_back(muon);
        }

        for(auto electron : myTriggeringElectrons)
        {
            bool medium = bernoulliPdf(randomEngine);
            bool tight = medium && bernoulliPdf(randomEngine);
            dec_pid(*electron) = tight? 2 : medium? 1 : 0;
            if(medium && electron->pt()>8e3f) ++nTrig_e7;
            if(tight && electron->pt()>27e3f) ++nTrig_e26;
        }
        
        if(nTrig_e26 < 1 
            && (nTrig_e7==0 || myTriggeringMuons.size()==0) 
            && nTrig_e17 < 2) 
        {
            continue;
        }
 

        double sf = 1.;
        auto cc = myTool->getEfficiencyScaleFactor(myTriggeringElectrons,
            myTriggeringMuons, sf);
        if(cc==CP::CorrectionCode::Ok)
        {
            nSuitableEvents += 1;
            sumW += sf;
        }
        else
        {
            Warning(MSGSOURCE, "Scale factor evaluation failed");
            ++errors;
        }
        if(errors>10)
        {
            Error(MSGSOURCE, "Too many errors reported!");
            break;
        }
    }
    Info(MSGSOURCE, "Average scale factor: %f (over %ld events)",
            sumW / nSuitableEvents, long(nSuitableEvents));
    #ifndef XAOD_STANDALONE
        ANA_CHECK(app->finalize());
    #endif
    return errors? 4 : 0;
}