TrigGlobEffCorrExample3a.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 3a"

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

Macro Definition Documentation

MSGSOURCE

#define MSGSOURCE   "Example 3a"

Definition at line 99 of file TrigGlobEffCorrExample3a.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 the leading electron as 'Signal'

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

<key in map file>, <PID WP>, <iso WP> For leading electron:

For subleading electron(s):

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)

Events must contain at least two lepton above trigger threshold

Let's pretend that the leading electron passes TightLH+isolation, thus the event passes the selection; now we tag the electrons:

Leading electron tagged as 'Signal' -> decorated value set to 1 Subleading electron(s) not tagged -> decorated value set to 0

Finally retrieve the global trigger scale factor

Definition at line 104 of file TrigGlobEffCorrExample3a.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");

    vector<asg::AnaToolHandle<IAsgElectronEfficiencyCorrectionTool>> factory;
    enum{ cLEGS, cTAG, cKEY, cPID, cISO };
    std::vector<std::array<std::string,5> > toolConfigs = {
        {"e12_lhloose_L1EM10VH", "Signal", "2015_e12_lhloose_L1EM10VH", "Tight", "FixedCutTightTrackOnly"},
        {"e12_lhloose_L1EM10VH", "*", "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 = "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;

        if(myTriggeringElectrons.size() < 2) continue;

        auto compareByPt = [](const xAOD::Electron*e1, const xAOD::Electron*e2)
            { return e1->pt() < e2->pt(); };
        auto leadingElectron = *std::max_element(myTriggeringElectrons.begin(),
            myTriggeringElectrons.end(), compareByPt);
        for(auto electron : myTriggeringElectrons)
        {
            dec_signal(*electron) = (electron==leadingElectron)? 1 : 0;
        }
 

        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;
}