TrigGlobEffCorrExample3c.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 3c"

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

Macro Definition Documentation

MSGSOURCE

#define MSGSOURCE   "Example 3c"

Definition at line 94 of file TrigGlobEffCorrExample3c.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 electron(s) as 'Signal'

To emulate PID selection (90% loose-to-medium efficiency)

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

<key in map file>, <PID WP>, <iso WP> Single electron trigger: electrons tagged 'Signal'

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

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

2015 periods D-H, J

Get a random run number, and decorate the event info

electron must be above softest trigger threshold (e12 here)

Add 'Signal' decorations to random electrons also count 'Signal' electrons above e24_xxx threshold

Events must contain enough leptons to trigger

single-electron trigger

dielectron

Finally retrieve the global trigger scale factor

Definition at line 99 of file TrigGlobEffCorrExample3c.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_signal("Signal");
    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 = {
        {"e24_lhmedium_L1EM20VH_OR_e60_lhmedium_OR_e120_lhloose", "Signal", "2015_e24_lhmedium_L1EM20VH_OR_e60_lhmedium_OR_e120_lhloose", "Medium", ""}, 
        {"e12_lhloose_L1EM10VH", "*,Signal", "2015_e12_lhloose_L1EM10VH", "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 global trigger SF tool");
    asg::AnaToolHandle<ITrigGlobalEfficiencyCorrectionTool> myTool("TrigGlobalEfficiencyCorrectionTool/TrigGlobal");
    myTool.setProperty("ElectronEfficiencyTools", electronEffTools).ignore();
    myTool.setProperty("ElectronScaleFactorTools", electronSFTools).ignore();
    const char* triggers2015 = 
        "e24_lhmedium_L1EM20VH_OR_e60_lhmedium_OR_e120_lhloose"
        "|| 2e12_lhloose_L12EM10VH";
    myTool.setProperty("TriggerCombination2015", triggers2015).ignore();
    myTool.setProperty("LeptonTagDecorations", "Signal").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[] = {
        276073, 278727, 279932, 280423, 281130, 282625
    };
    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;
        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 < 13e3f) continue;
 
            myTriggeringElectrons.push_back(electron);
        }
 
        vector<const xAOD::Muon*> myTriggeringMuons;

        unsigned nTrig1L = 0;
        for(auto electron : myTriggeringElectrons)
        {
            bool signal = bernoulliPdf(randomEngine);
            dec_signal(*electron) = signal? 1 : 0;
            if(signal && electron->pt()>25e3f) ++nTrig1L;
        }
        
        if(nTrig1L < 1 
            && myTriggeringElectrons.size() < 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;
}