applyDataNormalisation.h

Go to the documentation of this file.

#include <cstdlib>
#include <cassert>
#include <iostream>
#include <vector>
#include <map>
#include <filesystem>
 
#include "Core/Objects/interface/Jet.h"
#include "Core/Objects/interface/Event.h"
 
#include "Core/CommonTools/interface/toolbox.h"
 
#include <TFile.h>
#include <TChain.h>
#include <TH2.h>
#include <TString.h>
 
#include "Math/VectorUtil.h"
 
#include "Core/Normalisation/interface/TriggerEff.h"
#include "Core/Normalisation/interface/TriggerLumi.h"
#include "Core/Normalisation/interface/PhaseSelection.h"
 
namespace fs = std::filesystem;
 
namespace DAS::Normalisation {
 
inline void reset (Weights& wgts, const double v = 0) { std::fill(wgts.begin(), wgts.end(), Weight{v,0}); }
 
struct Functor {
    const std::map<int, TriggerLumi> triggers_lumi;
    TriggerEff eff;
    const size_t ntriggers;
    const double inv_total_lumi;
    const float minpt;
    size_t ibit;
    std::map<int, TriggerLumi>::const_reverse_iterator it;
 
private:
    const RecJet& (Functor::* NormalisationPtr)(RecEvent&, std::vector<RecJet>&, std::vector<FourVector>&, Trigger&);
    const RecJet& (Functor::* NormalisationCommonPtr)(RecEvent&, const RecJet&, const Trigger&);
 
public:
    Functor (fs::path lumi_file,        // !< path to the a two column file containg trigger lumis
             fs::path turnon_file,      // !< file containing turnons output of `getTriggerTurnons`
             fs::path trigger_curves,   // !< .root trigger curves, output of `getTriggerTurnons`
             const std::string& strategy,    // !< strategy 'pt' or 'eta'
             const std::string& method,      // !< method 'presc' or 'lumi'
             const int year) :          // !< year of the dataset (201?)
        triggers_lumi(GetLumiFromFiles(lumi_file, turnon_file)),
        eff(trigger_curves, triggers_lumi, year),
        ntriggers(triggers_lumi.size()),
        inv_total_lumi((prev(triggers_lumi.end())->second).weight),
        minpt(triggers_lumi.begin()->second.turnon)
    {
        using namespace std;
 
        cout << "ntriggers = " << ntriggers << '\n'
             << "tot lumi = " << 1./inv_total_lumi << '\n'
             << "minpt = " << minpt << endl;
 
        if (method == "prescales")
            NormalisationCommonPtr = &Functor::NormalisationCommonPresc;
        else if (method == "lumi")
            NormalisationCommonPtr = &Functor::NormalisationCommonLumi;
        else {
            cerr << "method ('" << method << "') is not defined, aborting\n";
            exit(EXIT_FAILURE);
        }
 
        if (strategy == "pt")
            NormalisationPtr = &Functor::NormalisationStd;
        else if (strategy == "eta")
            NormalisationPtr = &Functor::NormalisationFwd;
        else {
            cerr << "strategy of selection and normalisation ('" << strategy << "') is not defined, aborting\n";
            exit(EXIT_FAILURE);
        }
    }
 
    const RecJet& operator() (RecEvent & evnt,
                              std::vector<RecJet>& recJets,
                              std::vector<FourVector>& hltJets,
                              Trigger& trigger)
    {
        return (this->*NormalisationPtr)(evnt, recJets, hltJets, trigger);
    }
 
private:
    const RecJet& NormalisationCommonPresc (RecEvent & evnt,
                                            const RecJet& leadingJet,
                                            const Trigger& trigger)
    {
        using namespace std;
 
        // get prescale
        double preHLT   = trigger.PreHLT[ibit];
        double preL1min = trigger.PreL1min[ibit];
        double preL1max = trigger.PreL1max[ibit];
        if (preL1min != preL1max)
            cerr << "\x1B[31m\e[1m" << preL1min << ' ' <<  preL1max << "\x1B[30m\e[0m\n";
 
        double prescale = preHLT * preL1min; 
        // sanity checks:
        // 1) we assume that the same prescale is indeed constant for a given LS
        {
            static vector<map<pair<int,int>,int>> prescales(ntriggers);
            pair<int,int> runlum = {evnt.runNo, evnt.lumi};
            if (prescales.at(ibit).count(runlum))
                assert(prescales.at(ibit).at(runlum) == prescale);
            else
                prescales.at(ibit)[runlum] = prescale;
        }
 
        // setting the inverse of the eff lumi to the event including correction of for the trigger efficiency
        double efficiency = eff(leadingJet);
        //cout << __LINE__ << ' ' << efficiency << endl; /// \todo verbose option
        if (efficiency <= 0)
            reset(evnt.weights);
        else {
            double norm = prescale * inv_total_lumi / efficiency;
            //cout << __LINE__ << ' ' << norm << endl; /// \todo verbose option
            evnt.weights *= norm;
        }
        return leadingJet;
    }
 
    const RecJet& NormalisationCommonLumi (RecEvent& evnt,
                                           const RecJet& leadingJet,
                                           const Trigger& trigger)
    {
        // retrieve the efficiency and the effective luminosity
        double efficiency = eff(leadingJet);
        double inv_eff_lumi = (it->second).weight;
 
 
        if (efficiency <= 0)
            reset(evnt.weights);
        else {
            double norm = inv_eff_lumi / efficiency;
            evnt.weights *= norm;
        }
        return leadingJet;
    }
 
    const RecJet& NormalisationStd (RecEvent & evnt,
                                    std::vector<RecJet>& recJets,
                                    std::vector<FourVector>& hltJets,
                                    Trigger& trigger)
    {
        using namespace std;
 
        // we find the leading jet in tracker acceptance
        auto leadingInTk = phaseSel(recJets);
 
        // remove empty events
        if (leadingInTk == recJets.end()) {
            reset(evnt.weights);
            return recJets.front();
        }
        // we assume that JES corrections have been stored in DAS::RecJet::scales
        auto leading_pt = leadingInTk->CorrPt();
        // remove events if it does not contain any jets with pt > first threshold
        if (leading_pt < minpt) {
            reset(evnt.weights);
            return recJets.front();
        }
        // note: we want to keep jets with pt < minpt if at least the leading jets has >= 64
        // (but other jets are allowed to have a arbitrarily low pt)
 
        it = triggers_lumi.rbegin();
        ibit = ntriggers-1;
        // find the range by looping over triggers (from top to bottom)
        while (it != triggers_lumi.rend()) {
            if (leading_pt >= it->second.turnon) break;
            ++it; --ibit;
        }
 
        //cout << __LINE__ << ' ' << ibit << endl; /// \todo verbose option
 
        // if the leading pt was "too low"
        // (i.e. smaller than the lowest-pt turn-on),
        // then the event is discarded
 
        if (it == triggers_lumi.rend()) {
            cerr << "\x1B[31m\e[1mProblem with `trigger_lumi` iterator (leading_pt = " << leading_pt << ")\x1B[30m\e[0m\n";
            reset(evnt.weights);
            return recJets.front();
        }
        if (ibit >= ntriggers+1) {
            cerr << "\x1B[31m\e[1mProblem with `ibit` (leading_pt = " << leading_pt << ")\x1B[30m\e[0m\n";
            reset(evnt.weights);
            return recJets.front();
        }
 
        // check if the trigger has fired
        if (!trigger.Bit.at(ibit)) {
            reset(evnt.weights);
            return recJets.front();
        }
        return (this->*NormalisationCommonPtr)(evnt, *leadingInTk, trigger);
    }
 
    const RecJet& NormalisationFwd (RecEvent & evnt,
                                    std::vector<RecJet>& recJets,
                                    std::vector<FourVector>& hltJets,
                                    Trigger & trigger)
    {
        using namespace std;
 
        // sorte the jets in decreasing order of rapidity
        auto mnorder = [] (const RecJet& recJet1, const RecJet& recJet2) {
            return (recJet1.p4.Eta() > recJet2.p4.Eta());
        };
        sort(recJets.begin(), recJets.end(), mnorder);
 
        // take the most forward and most backward jets
        auto leadingJet = recJets.begin();
        auto subleadingJet = prev(recJets.end());
        FourVector leadingHltTk = match(leadingJet->p4, &hltJets);
        FourVector subleadingHltTk = match(subleadingJet->p4, &hltJets);
 
        if (leadingHltTk == FourVector() && subleadingHltTk == FourVector()) {
            reset(evnt.weights);
            return recJets.front();
        }
 
        // we assume that JES corrections have been stored in DAS::RecJet::scales
        float leading_pt = leadingJet->CorrPt();
        float subleading_pt = subleadingJet->CorrPt();
        // remove events if it does not contained any jet with pt > pt of the first threshold
        if (leading_pt < minpt && subleading_pt < minpt) {
            reset(evnt.weights);
            return recJets.front();
        }
 
        // note: we want to keep jets with pt < minpt if at least the leading jets has > 64
        // (but other jets are allowed to have a arbitrarily low pt)
 
        float leadingHlt_pt = leadingHltTk.Pt();
        float subleadingHlt_pt = subleadingHltTk.Pt();
 
        it = triggers_lumi.rbegin();
        ibit = ntriggers-1;
 
        // Find the range by looping over triggers (from top to bottom)
        while (it != triggers_lumi.rend()) {
            int recturnon = it->second.turnon;
            float hltthreshold = it->first;
            if ( (leading_pt >= recturnon) && (leadingHlt_pt >= hltthreshold)
                    && (subleading_pt >= recturnon) && (subleadingHlt_pt >= hltthreshold) ) {
                if ( std::abs(subleadingJet->p4.Eta()) > std::abs(leadingJet->p4.Eta()) ) {
                    swap(leadingJet, subleadingJet);
                    swap(leading_pt, subleading_pt);
                    swap(leadingHlt_pt, subleadingHlt_pt);
                }
                break;
            }
            if (leading_pt >= recturnon && leadingHlt_pt >= hltthreshold) break;
            if (subleading_pt >= recturnon && subleadingHlt_pt >= hltthreshold) {
                swap(leadingJet, subleadingJet);
                swap(leading_pt, subleading_pt);
                swap(leadingHlt_pt, subleadingHlt_pt);
                break;
            }
            ++it; --ibit;
        }
 
        // if the leading pt was "too low"
        // (i.e. smaller than the lowest-pt turn-on),
        // then the event is discarded
        if (it == triggers_lumi.rend()) {
            reset(evnt.weights);
            return recJets.front();
        }
 
        if (ibit >= ntriggers+1) {
            cerr << "\x1B[31m\e[1mProblem with `ibit` (leading_pt = " << leading_pt << ")\x1B[30m\e[0m\n";
            reset(evnt.weights);
            return recJets.front();
        }
 
        //check if the trigger has fired
        if (!trigger.Bit.at(ibit)) {
            reset(evnt.weights);
            return recJets.front();
        }
        return (this->*NormalisationCommonPtr)(evnt, *leadingJet, trigger);
    }
};
 
} // end of DAS::Normalisation namespace