TrigGlobEffCorrExample5a.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/

/*
 *    This example shows how to configure the tool for photon triggers. 
 *    Here, the minimal configuration is shown (one symmetric diphoton trigger)
 */
 
// ROOT include(s):
#include <TFile.h>
#include <TError.h>
 
// Infrastructure include(s):
#ifdef XAOD_STANDALONE
    #include "xAODRootAccess/Init.h"
    #include "xAODRootAccess/TEvent.h"
    #include "xAODRootAccess/TStore.h"
#else
    #include "AthAnalysisBaseComps/AthAnalysisHelper.h"
    #include "POOLRootAccess/TEvent.h"
#endif
 
#include "AsgMessaging/MessageCheck.h"
// EDM include(s):
#include "AsgTools/AnaToolHandle.h"
#include "EgammaAnalysisInterfaces/IAsgPhotonEfficiencyCorrectionTool.h"
#include "TriggerAnalysisInterfaces/ITrigGlobalEfficiencyCorrectionTool.h"
#include "xAODEventInfo/EventInfo.h"
#include "xAODEgamma/PhotonContainer.h"
#include "PATCore/PATCoreEnums.h"
#include "AthContainers/ConstAccessor.h"
#include "AthContainers/Decorator.h"
 
// stdlib include(s):
#include <random>
#include <vector>
#include <array>
using std::vector;
using std::string;
 
#define MSGSOURCE "Example 5a"
 
ANA_MSG_HEADER(Test)
ANA_MSG_SOURCE(Test, MSGSOURCE)
using namespace Test;
int main(int argc, char* argv[])
{
    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 photon CP tools");
    ToolHandleArray<IAsgPhotonEfficiencyCorrectionTool> photonEffTools;
    ToolHandleArray<IAsgPhotonEfficiencyCorrectionTool> photonSFTools; 

    vector<asg::AnaToolHandle<IAsgPhotonEfficiencyCorrectionTool>> factory;
 
    const char* mapPath = "PhotonEfficiencyCorrection/2015_2018/"
            "rel21.2/Summer2018_Rec_v1/map1.txt";
    for(int j=0;j<2;++j) 
    {
        string name = "AsgPhotonEfficiencyCorrectionTool/" + std::string(j? "PhTrigEff" : "PhTrigSF");
        auto t = factory.emplace(factory.end(), name);
        t->setProperty("MapFilePath", mapPath).ignore();
        t->setProperty("TriggerKey", string(j?"":"Eff_") + "DI_PH_2015_2016_g20_tight_2016_g22_tight_2017_2018_g20_tight_icalovloose_L1EM15VHI").ignore();
        t->setProperty("IsoKey", "Loose").ignore();
        t->setProperty("ForceDataType", (int)PATCore::ParticleDataType::Full).ignore();
        if(t->initialize() != StatusCode::SUCCESS)
        {
            Error(MSGSOURCE, "Unable to initialize the photon CP tool <%s>!",
                    t->name().c_str());
            return 3;
        }
        auto& handles = (j? photonSFTools : photonEffTools);
        handles.push_back(t->getHandle());
    }
    
    /* ********************************************************************** */
    
    Info(MSGSOURCE, "Configuring the global trigger SF tool");
    asg::AnaToolHandle<ITrigGlobalEfficiencyCorrectionTool> myTool("TrigGlobalEfficiencyCorrectionTool/TrigGlobal");
    myTool.setProperty("PhotonEfficiencyTools", photonEffTools).ignore();
    myTool.setProperty("PhotonScaleFactorTools", photonSFTools).ignore();
    std::map<std::string, std::string> triggers;
    triggers["266904-302872"] = "2g20_tight"; 
    triggers["302919-311481"] = "2g22_tight"; 
    triggers["2017"] = "2g22_tight_L12EM15VHI";
    triggers["2018"] = "2g22_tight_L12EM15VHI";
    myTool.setProperty("TriggerCombination", triggers).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[] = {
        276262, 278727, 279932, 280423, 281130, 282625,
        297730, 300345, 301912, 302737, 303638, 303943, 305291, 307124, 
        305359, 309311, 310015,
        325713, 329385, 330857, 332720, 334842, 336497, 336832, 338183,
        348885, 349534, 350310, 352274, 354107, 354826, 355261, 355331,
        355529, 357050, 359191
    };
    std::uniform_int_distribution<unsigned> uniformPdf(0,
            sizeof(periodRuns)/sizeof(*periodRuns) - 1);
    std::default_random_engine randomEngine;
    
    /* ********************************************************************** */
    
    Info(MSGSOURCE, "Starting the event loop");
    unsigned errors = 0;
    double nSuitableEvents = 0., sumW = 0.;
    static const SG::Decorator<unsigned> RandomRunNumberDec("RandomRunNumber");
    static const SG::ConstAccessor<int> truthTypeAcc("truthType");
    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::Photon*> myTriggeringPhotons;
        const xAOD::PhotonContainer* photons = nullptr;
        event.retrieve(photons,"Photons").ignore();
        for(auto photon : *photons)
        {
            if(!photon->caloCluster()) continue;
            float eta = fabs(photon->caloCluster()->etaBE(2));
            float pt = photon->pt();
            if(pt<10e3f || eta>=2.37 || (eta>1.37 && eta<1.52)) continue;
            int t = truthTypeAcc(*photon);
            if(t!=14) continue;
            if(pt < (runNumber>=302919? 23e3f : 21e3f)) continue;
            myTriggeringPhotons.push_back(photon);
        }

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

        double sf = 1.;
        auto cc = myTool->getEfficiencyScaleFactor(myTriggeringPhotons, 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;
}