DAS::Normalisation Namespace Reference

Classes
struct	Efficiency
struct	Functor
struct	LumiUnc
struct	TriggerEff
struct	TriggerLumi

Enumerations
enum	{ WEIGHT = 0b001, SIGN = 0b010, CUTOFF = 0b100 }

Functions
void	getHLTJetResponse (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
unique_ptr< TH3 >	makeRespHist (TString name, TString study)
void	applyDataNormalisation (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	reset (Weights &wgts, const double v=0)
float	GetNormFactor (const vector< fs::path > &sumWgts, const float xsec)
void	applyMClumi (const vector< fs::path > &inputs, const vector< fs::path > &sumWgts, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	applyNormFactor (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	getSumWeights (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
std::vector< DAS::RecJet >::iterator	phaseSel (std::vector< DAS::RecJet > &recjets)
std::ostream &	operator<< (std::ostream &Stream, const TriggerLumi &tr_lu)
float	getPrescale (Trigger *trigger, size_t indx)
void	getTriggerCurves (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	getTriggerTurnons (const fs::path &input, const fs::path &outputRoot, const fs::path &outputTxt, const pt::ptree &config, const int steering)
DAS::FourVector	match (const DAS::FourVector &jet, const std::vector< DAS::FourVector > *hltJets)
float	GetX (TF1 *f)
template<typename STREAM >
TF1 *	FitTriggerEfficiency (int trigger, TH1 *h, STREAM &Stream, int from, const char *options="NQERS")
template<typename STREAM >
float	GetTriggerEfficiency (int trigger, TH1 *h, STREAM &Stream, float from, float minEff)

Variables
vector< double >	thresholds
int	nThresholds
const float	threshold = 0.995

Enumeration Type Documentation

anonymous enum

anonymous enum

Enumerator
WEIGHT	use the weights
SIGN	use the sign
CUTOFF	apply a selection on the weight

     {
    WEIGHT = 0b001, 
    SIGN   = 0b010, 
    CUTOFF = 0b100  
};

Function Documentation

applyDataNormalisation()

void DAS::Normalisation::applyDataNormalisation	(	const vector< fs::path > &	inputs,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Apply the the kinematic selection, ensuring that events in the dataset have fired a trigger and correcting the weight of the events with the associated trigger prescales, trigger efficiency, and total luminosity. The weight $\frac{ \text{Prescale} * \mathcal{L}^{-1} }{ \epsilon_\text{trigger} } $ is applied by event. The argument 'strategy' defines the selection criterion of the events and the calculation of the weight (see functor in applyDataNormalisation.h).

Todo:

add lumi unc variations

Todo:

revisit the logic (having inputs modified + a return object may be confusing)

Todo:

reimplement the logic to avoid a selection on the weight

Parameters

inputs	input ROOT files (n-tuple)
output	output ROOT file (n-tuple)
config	config handled with `Darwin::Tools::Options`
steering	parameters obtained from explicit options
slice	number and index of slice

                               {1,0} 
        )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    if (isMC)
        BOOST_THROW_EXCEPTION( DE::BadInput("Only real data may be used as input.", metainfo) );
 
    // event information
    auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
    auto trigger = flow.GetBranchReadOnly<Trigger>("jetTrigger");
 
    // physics object information
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto hltJets = flow.GetBranchReadOnly<vector<FourVector>>("hltJets");
 
    int year = metainfo.Get<int>("flags", "year");
 
    // defining functor to apply data prescales
    auto lumi_file    = config.get<fs::path>("corrections.normalisation.luminosities" ),
         unc_file     = config.get<fs::path>("corrections.normalisation.uncertainties"),
         turnon_file  = config.get<fs::path>("corrections.normalisation.turnons"      ),
         efficiencies = config.get<fs::path>("corrections.normalisation.efficiencies" );
    auto strategy = config.get<string>("corrections.normalisation.strategy"),
         method   = config.get<string>("corrections.normalisation.method"  );
    Functor normalisation(lumi_file, turnon_file, efficiencies, strategy, method, year);
    LumiUnc lumi_unc(unc_file, year);
    fOut->cd();
 
    bool applySyst = steering & DT::syst;
    if (applySyst)
        for (string source: lumi_unc.sources)
            metainfo.Set<string>("variations", RecEvent::WeightVar, source);
 
    // a few control plots
    ControlPlots::isMC = false;
    ControlPlots corrNoTrigEff("corrNoTrigEff");
 
    vector<ControlPlots> corr { ControlPlots("nominal") };
    TList * ev_wgts = metainfo.List("variations", RecEvent::WeightVar);
    for (TObject * obj: *ev_wgts) {
        auto name = dynamic_cast<TObjString*>(obj)->GetString();
        corr.push_back(ControlPlots(name));
    }
    //if (applySyst)
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (applySyst)
            lumi_unc(recEvt);
 
        if (recJets->size() == 0) continue;
 
        RecJet leadingJet = normalisation(*recEvt, *recJets, *hltJets, *trigger);
        cout << "leadingJet: " << leadingJet << endl;
        // the normalisation itself is done in `Normalisation::Functor`
 
        if (recEvt->weights.front() <= 0) continue;
        if (steering & DT::fill) tOut->Fill();
 
        float evtWgts = recEvt->weights.front();
        float efficiency = normalisation.eff(leadingJet);
 
        cout << "evtWgts = " << evtWgts << '\n'
             << "efficiency = " << efficiency << endl;
        if (efficiency <= 0) continue;
 
        corrNoTrigEff(*recJets, evtWgts*efficiency); // compensate
        for (size_t i = 0; i < corr.size(); ++i)
            corr.at(i)(*recJets, recEvt->weights.at(i));
 
        // Exclusive curves are filled (i.e. one event can populate only a trigger curve).
        // The ibit value is determined by
        // the leading jet but also the other jets of the event can populate the trigger curve.
        for (auto& jet: *recJets) {
            float y = jet.AbsRap();
            float w = recEvt->weights.front();
            normalisation.eff.contribs.at(normalisation.ibit)->Fill(jet.CorrPt(), y, w);
        }
    }
 
    fOut->cd();
    corrNoTrigEff.Write(fOut);
    for (size_t i = 0; i < corr.size(); ++i)
        corr.at(i).Write(fOut);
    TDirectory * controlplots = fOut->mkdir("controlplots");
    controlplots->cd();
    for (TH2 * h: normalisation.eff.contribs) {
        h->SetDirectory(controlplots);
        h->Write();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

applyMClumi()

void DAS::Normalisation::applyMClumi	(	const vector< fs::path > &	inputs,
		const vector< fs::path > &	sumWgts,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Normalise with sum of weights (to be computed) and cross section (given)

Note

A cut-off for events with hard scale above 5 TeV is applied, since these events are not realistic

Todo:

Events with scale > 5TeV are removed because Pythia used to generate too many of them. Check that this is still needed and that the cut is still appropriate.

Parameters

inputs	input ROOT files (n-tuple)
sumWgts	input ROOT files (histogram)
output	output ROOT file (n-tuple)
config	config handled with `Darwin::Tools::Options`
steering	parameters obtained from explicit options
slice	number and index of slice

                               {1,0} 
        )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
                                      make_unique<TFile>(inputs.front().c_str() )) );
 
    // normalisation factor
    auto xsection = config.get<float>("corrections.MCnormalisation.xsection");
    float factor = GetNormFactor(sumWgts, xsection);
    metainfo.Set<float>("corrections", "MCnormalisation", "xsection", xsection);
 
    // declaring branches
    auto genEvt = flow.GetBranchReadWrite<GenEvent>("genEvent");
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto genJets = flow.GetBranchReadOnly<vector<GenJet>>("genJets", DT::facultative);
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets", DT::facultative);
 
    // control plot
    ControlPlots::isMC = true;
    ControlPlots plots("controlplots");
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (genEvt->hard_scale > 5000 /* GeV */) continue;
 
        // sanity check
        if (genEvt->weights.size() != 1 || recEvt->weights.size() != 1)
            BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected event weights", tIn) );
        if (recJets != nullptr && recJets->size() > 0 && recJets->front().scales.size() > 1)
            BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected jet energy scale variations", tIn) );
 
        // renormalisation
        genEvt->weights.front() *= factor;
        if (steering & DT::fill) tOut->Fill();
 
        // control plot
        if (genJets != nullptr)
            plots(*genJets, genEvt->weights.front()                          );
        if (recJets != nullptr)
            plots(*recJets, genEvt->weights.front() * recEvt->weights.front());
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    plots.Write(fOut);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

applyNormFactor()

void DAS::Normalisation::applyNormFactor	(	const vector< fs::path > &	inputs,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Multiplies the gen weight of every event by a factor provided on the command line.

Parameters

inputs	input ROOT files (n-tuple)
output	output ROOT file (n-tuple)
config	config handled with `Darwin::Tools::Options`
steering	parameters obtained from explicit options
slice	number and index of slice

                               {1,0} 
        )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
                                      make_unique<TFile>(inputs.front().c_str() )) );
 
    const float factor = config.get<float>("corrections.normFactor");
    metainfo.Set<float>("corrections", "normFactor", factor);
 
    // declaring branches
    auto evnt = flow.GetBranchReadOnly<GenEvent>("genEvent");
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        // sanity check
        if (evnt->weights.size() != 1)
            BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected event weights", tIn) );
 
        // normalise
        for (auto& weight: evnt->weights)
            weight *= factor;
 
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

FitTriggerEfficiency()

TF1* DAS::Normalisation::FitTriggerEfficiency	(	int	trigger,
		TH1 *	h,
		STREAM &	Stream,
		int	from,
		const char *	options = "NQERS"
	)

Fit the trigger efficiency with sigmoid function in [trigger, 3* trigger], where trigger is the HLT pt threshold.

The sigmoid function is defined as follows:

\epsilon (p_T) = a + 0.5 \times (1-a) \times \left( 1 + \erf \left( \frac{p_T - \mu}{\sigma} \right) \right)

where a, $\mu$ and $\sigma$ are the three parameters.

NOTE:

• the output is very sensitive to the fit range
• unclear if/how JES uncertainties should be taken into account

Parameters

trigger	trigger HLT threshold
h	efficiency histogram
Stream	`cout` or file
from	last turnon, to know from where the ref is efficient
options	fit options for TF1 function

{
    assert(from > 0);
    int mbin = h->GetXaxis()->FindBin(  from),
        Mbin = h->GetXaxis()->FindBin(3*from);
    float m = h->GetBinLowEdge(mbin+1), M = h->GetBinLowEdge(Mbin);
    for (int bin = 0; bin <= mbin; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
    for (int bin = Mbin+1; bin <= h->GetNbinsX()+1; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
 
    cout << "Fitting trigger " << trigger << " from " << m << " to " << M << endl;
 
    // function
    TF1 * f = new TF1(Form("sigmoid%d", trigger),
            "[0]+0.5*(1-[0])*(1+erf((x-[1])/[2]))",
            m, M);
 
    f->SetParNames("a", "#mu", "#sigma");
 
    //f->FixParameter(0,0);
    f->SetParameter(0,0.1);
    f->SetParLimits(0,0,0.9);
 
    f->SetParameter(1,1.1*trigger);
    f->SetParLimits(1,m,M);
 
    f->SetParameter(2,10);
    f->SetParLimits(2,1,150);
 
    // fit
    h->Fit(f, options, "", m, M);
 
    float turnon = GetX(f);
 
    Stream << trigger << '\t' << turnon << '\n';
 
    return f;
}

getHLTJetResponse()

void DAS::Normalisation::getHLTJetResponse	(	const vector< fs::path > &	inputs,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Get HLT single-jet response in bins of kinematics in real data. The output can be used as input to fitJetResponse to determine and understand the deviations from the Gaussian expectations, which prevent from fitting the trigger efficiency with the error function.

Parameters

inputs	input ROOT files (n-tuples)
output	name of output ROOT file (histograms)
config	config file from `DTOptions`
steering	steering parameters from `DTOptions`
slice	slices for running

                               {1,0} 
        )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    if (isMC)
        BOOST_THROW_EXCEPTION( DE::BadInput("Only DATA may be used as input.", metainfo) );
 
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto recjets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto hltjets = flow.GetBranchReadOnly<vector<FourVector>>("hltJets");
 
    unique_ptr<TH3> inclResp = makeRespHist("inclusive", "HLT"); // Inclusive response
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        auto recWgt = rEv->weights.front();
        // Matching & Fill response
        for (const RecJet& recjet: *recjets) {
            FourVector hltjet = match(recjet.p4, hltjets);
            auto ptrec = recjet.CorrPt();
            auto pthlt = hltjet.Pt();
 
            int pthltbin = inclResp->GetXaxis()->FindBin(pthlt);
            double hlt_lowedge = inclResp->GetXaxis()->GetBinLowEdge(pthltbin);
 
            // Fill response
            auto response = ptrec/hlt_lowedge;
            auto yrec = recjet.AbsRap();
 
            inclResp->Fill(pthlt, yrec, response, recWgt);
 
        } // End of recjets' loop
    } // End of event loop
 
    // Creating directory hierarchy and saving histograms
    inclResp->SetDirectory(fOut);
    inclResp->SetTitle("Response");
    inclResp->Write("Response");
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

GetNormFactor()

float DAS::Normalisation::GetNormFactor	(	const vector< fs::path > &	sumWgts,
		const float	xsec
	)

Retrieve the sum of the weights obtained with getSumWeights and combines it with the input cross section to calculate the normalisation factor.

Parameters

sumWgts	file containing the hist
xsec	cross section value

{
    auto h = DT::GetHist<TH2>(sumWgts, "hSumWgt");
    auto sumw = h->Integral(0, -1, // duplicates in overflow cancel overweighted entries
                            1, 2); // both positive and negative weights
    if (sumw <= 0)
        BOOST_THROW_EXCEPTION( DE::BadInput("Negative sum of weights.", h) );
    if (xsec <= 0)
        BOOST_THROW_EXCEPTION( invalid_argument("Negative cross section.") );
    auto factor = xsec/sumw;
    cout << xsec << '/' << sumw << '=' << factor << endl;
    return factor;
}

getPrescale()

float DAS::Normalisation::getPrescale	(	Trigger *	trigger,
		size_t	indx
	)

{
    float preHLT   = trigger->PreHLT[indx];
    float preL1min = trigger->PreL1min[indx];
    float preL1max = trigger->PreL1max[indx];
    assert(preL1min == preL1max);
    float prescale = preHLT * preL1min;
    return prescale;
}

getSumWeights()

void DAS::Normalisation::getSumWeights	(	const vector< fs::path > &	inputs,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Sum the weights from ROOT n-tuple. Mostly useful for MC to estimate the effective number of events, which is necessary to normalise.

The histogram is filled in a way so that it can be used for two purposes:

• provide a distribution of the weights and choose a max weight,
• extract the normalisation factor by takaing its integral.

Note

The values beyond maxAbsValue are filled a 2nd time in the overflow with the opposite weight. The truncated sum of the weights is then obtained from the integral by including the overflow.

Parameters

inputs	input ROOT files (n-tuple)
output	output ROOT file (histograms)
config	config handled with `Darwin::Tools::Options`
steering	parameters obtained from explicit options
slice	number and index of slice

                               {1,0} 
        )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
                                      make_unique<TFile>(inputs.front().c_str() )) );
 
    const auto maxAbsValue = config.get<float>("corrections.weights.maxAbsValue");
    const auto mode = config.get<string>("corrections.weights.mode");
 
    cout << "maxAbsValue = " << maxAbsValue << endl;
 
    auto evnt = flow.GetBranchReadOnly<GenEvent>("genEvent");
 
    TString title = ";log_{10}(w);sign;#sum_{i=1}^{events} w_{i}";
    if (maxAbsValue > 0) title = Form("w < %d", maxAbsValue);
    auto h = make_unique<TH2D>("hSumWgt", title, 100, 3, 40,
                                                 2, -1, 1);
    const double overflow = h->GetXaxis()->GetXmax()+1.;
 
    uint32_t wm = 0;
    if      (mode == "weight"  ) wm = WEIGHT | SIGN;
    else if (mode == "signOnly") wm = SIGN;
    else if (mode != "count"   )
         BOOST_THROW_EXCEPTION( invalid_argument(mode + " is unknown") );
 
    if (maxAbsValue > 0) wm |= CUTOFF;
 
    cout << "wm = " << bitset<8>(wm) << endl;
 
    if (steering & DT::verbose)
        cout << setw(20) << "w"
             << setw(20) << "absw"
             << setw(20) << "x [log10(abs(w))]"
             << setw(20) << "y [sign]"
             << setw(20) << "z [eff. weight]" << endl;
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        const double w = evnt->weights.front();
        if (w == 0.) continue;
 
        if (!isfinite(w))
            BOOST_THROW_EXCEPTION( runtime_error("Infinite or NaN weight:"s + w) );
 
        // the weight is used twice:
        // - once for the bin to fill, which should always be [log(abs(w)),+/-1]
        // - once for the actual weight in the histogram, which depends on the options
 
        const double absw = std::abs(w),
                     x = log10(absw),
                     y = copysign(0.5,w); // bin center
        double z = wm & WEIGHT ? absw : 1;
        if (wm & SIGN) z = copysign(z,w);
 
        cout << setw(20) << w
             << setw(20) << absw
             << setw(20) << x
             << setw(20) << (2*y) // we want the sign only (bin center is irrelevant)
             << setw(20) << z << endl;
 
        h->Fill(x, y, z);
 
        if (x >= overflow)
            cerr << orange << "Warning: a good weight is being filled in the overflow.\n" << def;
 
        if ((wm & CUTOFF) && absw >= maxAbsValue)
            h->Fill(overflow, y, -z);
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    h->Write();
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

getTriggerCurves()

void DAS::Normalisation::getTriggerCurves	(	const vector< fs::path > &	inputs,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Get trigger curves from the data n-tuples.

Todo:

fix error

Todo:

fix error

Parameters

inputs	input ROOT file (n-tuple)
output	output ROOT file (n-tuple)
config	config handled with `Darwin::Tools::options`
steering	parameters obtained from explicit options
slice	number and index of slice

                                       {1,0} 
                )
{
    //cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
 
    const auto usePrescales = config.get<bool>("corrections.normalisation.use_prescales");
    metainfo.Set<bool>("corrections", "normalisation", "use_prescales", usePrescales);
 
    auto lumi_file = config.get<fs::path>("corrections.normalisation.luminosities");
    pt::ptree triggers_lumi;
    pt::read_info(lumi_file.string(), triggers_lumi);
    for (auto const& trigger_lumi: triggers_lumi)
        thresholds.push_back(stoi(trigger_lumi.first));
    nThresholds = thresholds.size()-1;
    metainfo.Set<fs::path>("corrections", "normalisation", "luminosities", lumi_file);
 
    auto trigger = flow.GetBranchReadOnly<Trigger>("jetTrigger");
    auto recjets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto hltjets = flow.GetBranchReadOnly<vector<FourVector>>("hltJets");
 
    map<TString, size_t> methods { // methods[name] = shift (first trigger for which the method is used)
        {"emulation", 1}, // default method
        {"emulationShift", 2}, // emulation but with a lower-threshold reference trigger (e.g. trig 40 as reference for trig 80 as test)
        {"TnP", 0} // tag-and-probe method (used in any case of the lowest reference trigger(s)
    };
 
    TH2 * pt_correlation = new TH2F("ptCorr","; p_{T}^{PF,corr}; p_{T}^{hlt}",nPtBins,pt_edges.data(),nPtBins,pt_edges.data());
    TH1 * Deta = new TH1F("deta","; #Delta #eta; nevt", 200, -5, 5);
    TH1 * Dphi = new TH1F("dphi","; #Delta #phi; nevt", 200, -M_PI, M_PI);
    TH1 * dR   = new TH1F("dR"  ,"; #Delta R; nevt"   , 200, 0, 10);
 
    map<TString, map<int, Efficiency>> curves;
    for (auto method: methods)
    for (size_t i = method.second; i < thresholds.size(); ++i) {
        int t = thresholds.at(i);
        cout << method.first << '\t' << t << '\n';
        curves[method.first].insert( {t, Efficiency(method.first, t)} );
    }
    cout << flush;
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
 
        // sanity check
        size_t nbits = trigger->Bit.size();
        if (nbits < thresholds.size())
            BOOST_THROW_EXCEPTION( DE::AnomalousEvent( Form("nbits = %ld < thresholds.size() = %ld",
                                                             nbits,         thresholds.size()), tIn) );
 
        // avoid empty events
        if (recjets->size() == 0) continue;
 
        for (const RecJet& pfjet: *recjets){
            FourVector hltjet = match(pfjet.p4, hltjets);
            if ( hltjet == FourVector()) continue;
 
            Deta->Fill(pfjet.p4.Eta() - hltjet.Eta(), pfjet.weights.front());
            Dphi->Fill(DeltaPhi(pfjet.p4,hltjet), pfjet.weights.front());
            dR->Fill(DeltaR(pfjet.p4, hltjet), pfjet.weights.front());
            pt_correlation->Fill(pfjet.CorrPt(), hltjet.Pt(), pfjet.weights.front());
        }
 
        // NOTE: from now on, only the leading jet is considered
 
        auto leadingInTk = phaseSel(*recjets); // iterator
 
        if (leadingInTk == recjets->end()) continue;
 
        FourVector pfjet0 = leadingInTk->p4;
        pfjet0 *= leadingInTk->scales.front();
 
        FourVector hltjet0 = match(pfjet0, hltjets);
        if ( hltjet0 == FourVector()) continue;
 
        /***  using emulation method(s) ***/
        for (auto method: methods) {
 
            if (method.first.Contains("TnP")) continue;
            size_t shift = method.second;
 
            for (size_t i = shift; i < thresholds.size(); ++i) {
 
                int t = thresholds[i];
                auto& eff = curves.at(method.first).at(t);
 
                eff.hAll->Fill(pfjet0.Pt(), hltjet0.Pt());
 
                // test if previous trigger has fired
                if (!trigger->Bit[i-shift]) continue;
                eff.hFired->Fill(pfjet0.Pt(), hltjet0.Pt());
 
                if (hltjet0.Pt() < thresholds.at(i-shift)) continue;
                eff.hFiredThreshold->Fill(pfjet0.Pt(), hltjet0.Pt());
 
                float wgt = usePrescales ? getPrescale(trigger, i-shift) : 1;
                bool fired = hltjet0.Pt() > t;
                eff.Fill(pfjet0, fired, wgt*recjets->front().weights.front());
 
                eff.hHLTmap->Fill( pfjet0.Phi(),pfjet0.Eta(),wgt*recjets->front().weights.front() );
            }
        } // end of loop on methods
 
        /*** using tag & probe method ***/
        for (size_t i = 0; i < thresholds.size(); ++i) {
 
            int t = thresholds[i];
            auto& eff = curves.at("TnP").at(t);
 
            eff.hAll->Fill(pfjet0.Pt(), hltjet0.Pt());
 
            // first, it needs to have fired
            if (!trigger->Bit[i])
                continue;
            eff.hFired->Fill(pfjet0.Pt(), hltjet0.Pt());
 
            // then we need dijet configurations at PF level:
            // -> 1) at least 2 jets
            if (recjets->size() == 1)
                continue;
            // -> 2) back-to-back jets
            if (DeltaPhi(recjets->at(0).p4, recjets->at(1).p4) < 2.4) continue;
            // -> 3) a possible 3rd jet should not be significant
            if (recjets->size() > 2) {
                const double pt0 = recjets->at(0).CorrPt(),
                             pt1 = recjets->at(1).CorrPt(),
                             pt2 = recjets->at(2).CorrPt();
                if (pt2 > 0.15*(pt0 + pt1))
                    continue;
            }
 
            // we take the probe and the tag randomly
            static TRandom3 Rand(/* seed = */ metainfo.Seed<10989678>(slice));
            double r = Rand.Uniform();
            int itag = (r<0.5),
                iprobe = (r>=0.5);
 
            FourVector pftag   = recjets->at(itag  ).p4,
                       pfprobe = recjets->at(iprobe).p4;
 
            float pfprobe_wgt = recjets->at(iprobe).weights.front();
            float pftag_wgt = recjets->at(itag).weights.front();
 
            pftag   *= recjets->at(itag  ).scales.front();
            pfprobe *= recjets->at(iprobe).scales.front();
 
            // -> 4) selection on the rapidity of the system
            if (std::abs(recjets->at(itag).p4.Eta()) > 1.3) continue;
 
            // match PF at HLT
            // 1) test if tag can be matched
            //    - if not, continue
            // 2) then test if the probe can be match
 
            const FourVector& hlttag = match(pftag, hltjets);
            const FourVector& hltprobe = match(pfprobe, hltjets);
            if ( hlttag == FourVector() || hltprobe == FourVector()) continue;
 
            if (hlttag.Pt() > t) {
                eff.hFiredThreshold->Fill(pfjet0.Pt(), hltjet0.Pt());
 
                const FourVector& hltprobe = match(pfprobe, hltjets);
                bool fired = hltprobe.Pt() > t;
 
                float wgt = usePrescales ? getPrescale(trigger, i) : 1;
                eff.Fill(pfprobe, fired, wgt * pftag_wgt * pfprobe_wgt);
                eff.hHLTmap->Fill( pfjet0.Phi(),pfjet0.Eta(),wgt*recjets->front().weights.front() );
            }
        }
    }
 
    cout << "Writing to file" << endl;
    for (auto& method: curves) {
        fOut->cd();
        auto d = fOut->mkdir(method.first);
        for (auto& curve: method.second)
            curve.second.Write(d);
    }
    fOut->cd();
    dR->Write();
    Deta->Write();
    Dphi->Write();
    pt_correlation->Write();
 
    metainfo.Set<bool>("git", "complete", true);
 
    //cout << __func__ << ' ' << slice << " end" << endl;
}

GetTriggerEfficiency()

float DAS::Normalisation::GetTriggerEfficiency	(	int	trigger,
		TH1 *	h,
		STREAM &	Stream,
		float	from,
		float	minEff
	)

DON'T fit the trigger efficiency but just find the first bin above threshold

Todo:

Where does the factor 6 come from?

Parameters

trigger	trigger HLT threshold
h	efficiency histogram
Stream	`cout` or file
from	last turnon, to know from where the ref is efficient
minEff	min efficiency (e.g. different for barrel and forward region)

{
    int mbin = h->GetXaxis()->FindBin(  from),
        Mbin = h->GetXaxis()->FindBin(6*from); 
    for (int bin = 0; bin < mbin; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
    for (int bin = Mbin+1; bin <= h->GetNbinsX()+1; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
    for (int bin = mbin; bin <= Mbin; ++bin) {
        float content = h->GetBinContent(bin);
        if (content < minEff) continue;
        float turnon = h->GetBinLowEdge(bin);
        Stream << trigger << '\t' << turnon << '\n';
        return turnon;
    }
    return -1;
}

getTriggerTurnons()

void DAS::Normalisation::getTriggerTurnons	(	const fs::path &	input,
		const fs::path &	outputRoot,
		const fs::path &	outputTxt,
		const pt::ptree &	config,
		const int	steering
	)

Get the trigger turn-on points from the trigger curves.

Todo:

Use metainfo

Todo:

Get year from metainfo?

Todo:

Get lumis from metainfo?

Todo:

implement a more general logic to use only rapidity <3 turnons to determine the overall turnon

Parameters

input	input ROOT file (n-tuple)
outputRoot	output ROOT file with efficiency
outputTxt	output text file with turn-ons
config	config handled with `Darwin::Tools::options`
steering	parameters obtained from explicit options

{
    cout << __func__ << " start" << endl;
 
    auto fOut = DT::GetOutputFile(outputRoot);
 
    //DT::MetaInfo metainfo(tOut);
    //metainfo.Check(config);
    int year = config.get<int>("flags.year"); 
 
    cout << setw(10) << "trigger";
    for (auto yBin: yBins) cout << setw(18) << yBin;
    cout << endl;
 
    pt::ptree triggers_lumi;
    auto lumi_file = config.get<fs::path>("corrections.normalisation.luminosities"); 
    pt::read_info(lumi_file.c_str(), triggers_lumi);
    vector<int> triggerThresholds;
    for (auto trigger_lumi: triggers_lumi)
        triggerThresholds.push_back(stoi(trigger_lumi.first));
 
    auto fIn = make_unique<TFile>(input.c_str(), "READ");
    map<int, int> turnonsFinal;
    TIter Next(fIn->GetListOfKeys()) ;
    TKey* key = nullptr;
    TString defMethod = "TnP"; // only for the first trigger(s)
    while ( (key = dynamic_cast<TKey*>(Next())) ) { // looping over methods
        if (TString(key->ReadObj()->ClassName()) != TString("TDirectoryFile"))
            continue;
 
        fIn->cd();
        auto dIn = dynamic_cast<TDirectory*>(key->ReadObj());
 
        int lastturnon = triggerThresholds.front();
        cout << dIn << endl;
 
        TString method = dIn->GetName();
        cout << method << endl;
 
        fOut->cd();
        auto dOut = fOut->mkdir(method);
 
        TIter Next2(dIn->GetListOfKeys()) ;
        while ( (key = dynamic_cast<TKey*>(Next2())) ) { // looping of triggers
            if (TString(key->ReadObj()->ClassName()) != TString("TDirectoryFile"))
                continue;
 
            dIn->cd();
            auto ddIn = dynamic_cast<TDirectory*>(key->ReadObj());
 
            TString trigName = ddIn->GetName();
            assert(trigName.Contains("HLT"));
 
            dOut->cd();
            auto ddOut = dOut->mkdir(trigName);
 
            trigName.ReplaceAll("HLT","");
            int trigger = trigName.Atoi();
            if (method == "TnP"            && trigger < triggerThresholds.at(0)) continue;
            if (method == "emulation"      && trigger < triggerThresholds.at(1)) continue;
            if (method == "emulationShift" && trigger < triggerThresholds.at(2)) continue;
            cout << setw(10) << trigger;
 
            auto h2  = unique_ptr<TH2>(ddIn->Get<TH2>("test"));
            {
                auto den = unique_ptr<TH2>(ddIn->Get<TH2>("ref"));
                h2->Divide(h2.get(), den.get(), 1, 1, "B");
            }
 
            ddOut->cd();
 
            h2->SetDirectory(ddOut);
            h2->SetName("efficiency");
 
            vector<int> turnonsNow;
            for (int y = 1; y <= h2->GetNbinsY(); ++y) {
 
                // efficiency curve
                auto h1 = unique_ptr<TH1>(h2->ProjectionX(Form("efficiency_ybin%d", y), y, y));
                //cout << y << '\t';
                static ostream bitBucket(0); // printing to nowhere
                double threshold = 0.995;
                switch (year) {
                    case 2017:
                        if (y > 5) threshold = 0.99;
                        break;
                    case 2018:
                        if (y > 3) threshold = 0.99;
                        break;
                    case 2016:
                    default:
                        if (y > 5) threshold = 0.99; //2016 Pre
                        //if (y > 5) threshold = 0.99; 2016 Post
                }
                float from = h1->GetBinLowEdge(h1->FindBin(lastturnon)+1);
                int turnon = Normalisation::GetTriggerEfficiency(trigger, h1.get(), bitBucket, from, threshold);
                h1->SetTitle(Form("%d (%.3f)", turnon, threshold));
                h1->SetDirectory(dOut);
                h1->Write();
                turnonsNow.push_back(turnon);
            }
 
            auto maxturnon = max_element(turnonsNow.begin(), turnonsNow.begin()+6);
            h2->SetTitle(Form("%d", *maxturnon));
            h2->Write();
            for (auto turnon: turnonsNow) {
                if (turnon == *maxturnon) cout << "\x1B[32m\e[1m";
                if (turnon < 0) cout << "\x1B[31m\e[1m";
                cout << setw(18) << turnon;
                cout << "\x1B[30m\e[0m";
            }
            cout << endl;
 
            // final choice
            if (method == defMethod)
                turnonsFinal[trigger] = *maxturnon;
 
            defMethod = "emulation";
            lastturnon = *maxturnon;
 
            if (lastturnon == -1) {
                    cerr << "Didn't find any turnon for trigger " << trigger << '\n';
                    lastturnon = trigger;
            }
        } // end of loop over triggers
    } // end of loop over method
 
    pt::ptree turnon_file;
    for (auto& turnon: turnonsFinal)
        turnon_file.put<int>(to_string(turnon.first), turnon.second);
    pt::write_info(outputTxt.string(), turnon_file);
 
    cout << __func__ << " end" << endl;
}

GetX()

float DAS::Normalisation::GetX

(

TF1 *

f

)

Replaces TF1::GetX() which does not return the same value when run after compilation or in interactive mode...

Principle: just scan the pt axis with step = 1 until eff > 0.995

{
    double m = 0, M = 7000;
    f->GetRange(m,M);
 
    for (float pt = m; pt < M; ++pt)
        if (f->Eval(pt) > threshold)
            return pt;
 
    cerr << "No threshold was found from " << m << " to " << M << ".\nAborting.\n";
    exit(EXIT_FAILURE);
}

makeRespHist()

unique_ptr<TH3> DAS::Normalisation::makeRespHist	(	TString	name,
		TString	study
	)

Create TH3 histograms for the Jet response and HLT response.

Parameters

name	abitrary name
study	expecting JER or HLT

{
    using namespace DAS::JetEnergy;
 
    static int nResBins = 200;
    static auto resBins = getBinning(nResBins, 0, 2);
    unique_ptr<TH3D> h;
    if      (study == "HLT"){
        // Filling in bins of ptHLT and eta
        static const char * axistitles = ";p_{T}^{HLT};|#eta^{rec}|;#frac{p_{T}^{rec}}{p_{T}^{HLT bin low edge}}";
        static const vector<double> HLT_binning = {40., 60., 80., 140., 200., 260., 320., 400., 450., 500.};
        static int nHLTbins = HLT_binning.size()-1;
        h = make_unique<TH3D>(name, axistitles,
                nHLTbins, HLT_binning.data(),
                nYbins, y_edges.data(),
                nResBins, resBins.data());
    }
    else if (study == "JER"){
        // Filling in bins of ptgen results in better conditioned response distributions
        // Pure rec level weight is used
        static const char * axistitles = ";p_{T}^{gen};|#eta^{rec}|;#frac{p_{T}^{rec}}{p_{T}^{gen}}";
        h = make_unique<TH3D>(name, axistitles,
                nPtJERCbins, pt_JERC_edges.data(),
                nAbsEtaBins, abseta_edges.data(),
                nResBins, resBins.data());
    }
 
    h->SetDirectory(nullptr);
    return h;
}

match()

DAS::FourVector DAS::Normalisation::match	(	const DAS::FourVector &	jet,
		const std::vector< DAS::FourVector > *	hltJets
	)

                                                                                         {
    for (const auto& hltjet: *hltJets){
        using ROOT::Math::VectorUtil::DeltaR;
        if (DeltaR(hltjet, jet) < 0.3)
            return hltjet;
    }
    return DAS::FourVector();
}

operator<<()

std::ostream& DAS::Normalisation::operator<<	(	std::ostream &	Stream,
		const TriggerLumi &	tr_lu
	)

{
    Stream << tr_lu.pt << '\t' << tr_lu.turnon << '\t' << tr_lu.weight << '\t' << 1./tr_lu.weight;
    return Stream;
}

phaseSel()

std::vector<DAS::RecJet>::iterator DAS::Normalisation::phaseSel

(

std::vector< DAS::RecJet > &

recjets

)

                                                                      {
    auto leadingInTk = recjets.begin();
    while (leadingInTk != recjets.end()
            && std::abs(leadingInTk->Rapidity()) >= 3.0)
        ++leadingInTk;
 
    return leadingInTk;
}

reset()

void DAS::Normalisation::reset ( Weights &  wgts, const double  v = 0  )

inline

{ std::fill(wgts.begin(), wgts.end(), Weight{v,0}); }

Variable Documentation

nThresholds

int nThresholds

threshold

const float threshold = 0.995

thresholds

vector<double> thresholds