DAS::Normalisation Namespace Reference

Classes
struct	Efficiency
struct	Functor
struct	LumiUnc
class	Strategy
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})
template<typename T = TH3D>
std::unique_ptr< T >	makeRespHist (TString name, std::vector< double > pt_edges=DAS::JetEnergy::pt_JERC_edges, TString title=";p_{T}^{gen};|#eta^{rec}|;#frac{p_{T}^{rec}}{p_{T}^{gen}}")
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 int steering, 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)
void	loopOverDirs (TDirectory *dIn, TDirectory *dOut, const int steering, const DT::Slice slice)
void	getHistoEfficiencyCurves (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
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	getTriggerEfficiencyEmulation (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)
vector< double >	extractBinning (const unique_ptr< TH1 > &h3D, const TString &axis)
unique_ptr< TH1F >	makeTriggerEffHist (const unique_ptr< TEfficiency > &t_eff, int iHLT, int iEta)
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, inputs);
    auto tIn = flow.GetInputTree(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

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, inputs);
    auto tIn = flow.GetInputTree(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(steering, 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;
 
        // sanity check
        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, inputs);
    auto tIn = flow.GetInputTree(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;
}

extractBinning()

vector<double> DAS::Normalisation::extractBinning	(	const unique_ptr< TH1 > &	h3D,
		const TString &	axis
	)

Extract the axis binning from a TH1 obgect and return it in a vector.

{
    const TAxis *naxis = nullptr;
    if      (axis == "X") naxis = h3D->GetXaxis();
    else if (axis == "Y") naxis = h3D->GetYaxis();
    else if (axis == "Z") naxis = h3D->GetZaxis();
 
    int ndiv = naxis->GetNbins();
    vector<double> binning(ndiv+1);
    naxis->GetLowEdge(&binning[0]);
    binning[ndiv] = naxis->GetBinUpEdge(ndiv);
 
    return binning;
}

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

getHistoEfficiencyCurves()

void DAS::Normalisation::getHistoEfficiencyCurves	(	const fs::path &	input,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Extract trigger efficiency distributions from getTriggerEfficiency*Method*

Parameters

input	input ROOT file (output of getTriggerEfficiency*Method*)
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__ << " start" << endl;
    
    auto fIn  = make_unique<TFile>(input .c_str(), "READ"),
         fOut = make_unique<TFile>(output.c_str(), "RECREATE");
    Normalisation::loopOverDirs(fIn.get(), fOut.get(), steering, slice);
    
    cout << __func__ << ' ' << slice << " stop" << endl;
} // end of void

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, inputs);
    auto tIn = flow.GetInputTree(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"); // 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 etarec = abs(recjet.p4.Eta());
 
            inclResp->Fill(pthlt, etarec, 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 int	steering,
		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::Flow(steering, sumWgts).GetInputHist<TH2>("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;
    if (steering & DT::verbose)
        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, inputs);
    auto tIn = flow.GetInputTree(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 < %f", 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, inputs);
    auto tIn = flow.GetInputTree(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;
}

getTriggerEfficiencyEmulation()

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

Get the efficiency in 3D TEfficiency objects from data n-tuples using different selection criteria for the matching.

Todo:

: what if there are HLT jets but no rec jets? -> global inefficiency?

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__ << " start" << endl;
 
    DT::Flow flow(steering, inputs);
    auto tIn = flow.GetInputTree(slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    pt::ptree triggers_lumi;
    auto lumi_file = config.get<fs::path>("corrections.normalisation.luminosities");
    pt::read_info(lumi_file.c_str(), triggers_lumi);
    vector<double> thresholds;
    for (auto [trigger, _]: triggers_lumi)
        thresholds.push_back(stoi(trigger));
    if (thresholds.size() == 0)
        BOOST_THROW_EXCEPTION(invalid_argument("Empty list of triggers"));
    thresholds.push_back(3*thresholds.back()); // Added one more boundary to have an upper HLT limit
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    const auto R = metainfo.Get<int>("flags", "R");
 
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto trigger = flow.GetBranchReadOnly<Trigger>("jetTrigger");
    auto recjets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto hltjets = flow.GetBranchReadOnly<vector<FourVector>>("hltJets");
 
    using Matching = JMEmatching<vector<FourVector>, vector<RecJet>>;
    Matching::maxDR = R/10./2.;
    cout << "Radius for matching: " << Matching::maxDR << endl;
 
    Strategy leading_pt   ("leading_pt"   , thresholds),
             best_match   ("best_match"   , thresholds),
             most_central ("most_central" , thresholds),
             isolated_jets("isolated_jets", thresholds);
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        // sanity check
        size_t nbits = trigger->Bit.size();
        if (nbits < thresholds.size()-1)  // manually added an extra HLT boundary, -1 to satisfy the condition
            BOOST_THROW_EXCEPTION( DE::AnomalousEvent( Form("nbits = %ld < thresholds.size() = %ld",
                                                             nbits,        thresholds.size()), tIn) );
 
        if (recjets->empty()) continue;
 
        cout << "Running matching" << endl;
        Matching matching(*hltjets, *recjets);
 
        auto ev_w = rEv->Weight();
 
        auto get_hlt_match = [&matching,&hltjets](const RecJet& recjet) {
            // match_it = { hlt jet iterator, rec jet iterator }
            auto criterion = [&recjet](const auto& match_it) {
                return *(match_it.second) == recjet;
            };
            auto it = ranges::find_if(matching.match_its, criterion);
            return (it == matching.match_its.end()) ? hltjets->end() : it->first;
        };
 
        // leading pt
        {
            cout << "Leading pt" << endl;
            const auto& recjet = recjets->front();
            const auto hltjet_it = get_hlt_match(recjet);
            if (hltjet_it != hltjets->end())
                leading_pt.Fill(recjet, *hltjet_it, trigger->Bit, ev_w);
        }
 
        // most central
        {
            cout << "Most central" << endl;
            auto criterion = [](const RecJet& jet1, const RecJet& jet2) {
                return jet1.AbsRap() < jet2.AbsRap();
            };
            auto recjet_it = ranges::min_element(*recjets, criterion);
            const auto hltjet_it = get_hlt_match(*recjet_it);
            if (hltjet_it != hltjets->end())
                most_central.Fill(*recjet_it, *hltjet_it, trigger->Bit, ev_w);
        }
 
        // best match
        if (matching.match_its.size() > 0) {
            cout << "Best match" << endl;
            const auto& [hltjet_it, recjet_it] = matching.match_its.front();
            best_match.Fill(*recjet_it, *hltjet_it, trigger->Bit, ev_w);
        }
 
        // all isolated jets
        {
            cout << "All isolated jets" << endl;
            for (auto ref_it = recjets->begin(); ref_it != recjets->end(); ++ref_it) {
                const auto hltjet_it = get_hlt_match(*ref_it);
 
                auto is_neighbor = [ref_it, R](const RecJet& test_jet) {
                    return &test_jet != &(*ref_it) && DeltaR(ref_it->p4, test_jet.p4) < R / 10.; // Avoid comparing jet to itself and check proximity.
                };
                // Look for the first neighboring jet. If one is found, continue (not isolated).
                if (any_of(recjets->begin(), recjets->end(), is_neighbor)) continue;
                isolated_jets.Fill(*ref_it, *hltjet_it, trigger->Bit, ev_w);
            }
        }
    }
 
    leading_pt   .Write(fOut);
    most_central .Write(fOut);
    best_match   .Write(fOut);
    isolated_jets.Write(fOut);
 
    cout << __func__ << " stop" << endl;
}

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

loopOverDirs()

void DAS::Normalisation::loopOverDirs	(	TDirectory *	dIn,
		TDirectory *	dOut,
		const int	steering,
		const DT::Slice	slice
	)

Considering the possible strategies to determine the trigger efficiency of jets systems stored in 3D ROOT::TEfficiency objects, they need to be extarcted to other more convenient ROOT objects used in plotting and fitting:

• TH1
• TGraphAsymmErrors Based on whether the errors are symmetrical or not, the appropriate class should be used. In this case, the errors are symmetrical, so the TH1 class will be used.

This class includes the extraction of the corresponding histograms that are filled with the efficiency, for all HLT and |eta| bins.

{
    for (const auto&& obj: *(dIn->GetListOfKeys())) {
if (!obj->InheritsFrom("TKey")) continue;
        auto const key = dynamic_cast<TKey*>(obj);
        //dIn->cd();    // redundant?
        auto ddIn = dynamic_cast<TDirectory*>( key->ReadObj() );
        if (ddIn == nullptr) continue;
        ddIn->cd();
        cout << ddIn << endl;
 
        TString match = ddIn->GetName();
        cout << "Matching method: " << match << endl;
        cout << endl; 
 
        dOut->cd();
        auto ddOut = dOut->mkdir(match);  // Create a folder for each matching method
        ddOut->cd();
 
        auto t_eff = unique_ptr<TEfficiency>( dynamic_cast<TEfficiency*>( ddIn->Get("efficiency") ) );  // ROOT::TEfficiency from input
        if (t_eff == nullptr) continue;
 
        auto h3D = t_eff->GetTotalHistogram();  // used only to extract HLT thresholds
 
        // Extracting HLT bins (thresholds)
        auto h_HLTBins = unique_ptr<TH1>( dynamic_cast<TH1*>(h3D->Clone("h_HLTBins")) );
        vector<double> HLTBins = extractBinning(h_HLTBins, "X");
        const int nHLTBins = HLTBins.size()-1;
 
        // Extracting Eta bins
        auto h_EtaBins = unique_ptr<TH1>( dynamic_cast<TH1*>(h3D->Clone("h_EtaBins")) );
        vector<double> EtaBins = extractBinning(h_EtaBins, "Y");
        const int nEtaBins = EtaBins.size()-1;
 
        for (int iHLT = 1; iHLT <= nHLTBins; ++iHLT){
            if (t_eff->GetDimension() != 3) continue;  // works only for 3D TEfficiency plots
            
            auto dddOut = ddOut->mkdir(Form("HLT%.0f", HLTBins[iHLT-1]));  // Create a sub-folder for each HLT threshold
            dddOut->cd();
 
            vector<unique_ptr<TH1>> h1D_eff;
            for (int iEta = 1; iEta <= nEtaBins; ++iEta){
                auto ddddOut = dddOut->mkdir(Form("etabin%d", iEta));  // Create a sub-sub-folder for each Eta bin
                ddddOut->cd();
 
                cout << "HLT" << HLTBins[iHLT-1] << std::setw(10) << EtaBins[iEta-1] << "<|#eta|<" << EtaBins[iEta] << endl;
 
                // Fill 1D histogram with efficiency for each HLT and Eta bin
                h1D_eff.push_back(makeTriggerEffHist(t_eff, iHLT, iEta));
                h1D_eff.at(iEta-1)->SetDirectory(0);
 
                h1D_eff.at(iEta-1)->SetDirectory(ddddOut);
                h1D_eff.at(iEta-1)->Write();
            }  // end of eta loop
        }  // end of HLT loop
    }  // end of list-of-keys loop
}

makeRespHist()

std::unique_ptr<T> DAS::Normalisation::makeRespHist	(	TString	name,
		std::vector< double >	pt_edges = DAS::JetEnergy::pt_JERC_edges,
		TString	title = ";p_{T}^{gen};|#eta^{rec}|;#frac{p_{T}^{rec}}{p_{T}^{gen}}"
	)

Create a TH3 histogram for the Jet response and HLT jet response or trigger efficiency curve.

                         {T}^{gen};|#eta^{rec}|;#frac{p_{T}^{rec}}{p_{T}^{gen}}")
{
    using namespace std;
    using namespace DAS::JetEnergy;
 
    static int nResBins = 200;
    static auto resBins = getBinning(nResBins, 0, 2);
    int nPtBins = pt_edges.size()-1;
    unique_ptr<T> h = make_unique<T>(name, title,
                                     nPtBins, pt_edges.data(),
                                     nAbsEtaBins, abseta_edges.data(),
                                     nResBins, resBins.data());
 
    h->SetDirectory(nullptr);
    return h;
}

makeTriggerEffHist()

unique_ptr<TH1F> DAS::Normalisation::makeTriggerEffHist	(	const unique_ptr< TEfficiency > &	t_eff,
		int	iHLT,
		int	iEta
	)

Extract the HLT efficiency curves in 1D histograms from TEfficiency objects.

{
    auto  h3D = t_eff->GetTotalHistogram();  // used only to extract binning
 
    // Extracting HLT bins (thresholds)
    auto h3D_HLT = unique_ptr<TH1>( dynamic_cast<TH1*>(h3D->Clone("h3D_HLT")) );  // create one clone for HLT bins
    vector<double> HLTBins = extractBinning(h3D_HLT, "X");
 
    // Extracting resolution bins
    auto h3D_res = unique_ptr<TH1>( dynamic_cast<TH1*>(h3D->Clone("h3D_res")) );  // a separate clone for Res bins
    vector<double> ResBins = extractBinning(h3D_res, "Z");
    const int nResBins = ResBins.size()-1;
 
    TString name = "efficiency";
    TString title = Form("HLT%.0f_eta%d", HLTBins[iHLT-1], iEta+1);
    unique_ptr<TH1F> h = make_unique<TH1F>(name, title, nResBins, 0, 2*HLTBins[iHLT-1]);
    for (int ipt = 1; ipt <= nResBins; ++ipt){        // loop over pt bins (z axis of TEfficiency)
        h->SetBinContent( ipt, t_eff->GetEfficiency       (t_eff->GetGlobalBin(iHLT, iEta, ipt)) );
        h->SetBinError  ( ipt, t_eff->GetEfficiencyErrorUp(t_eff->GetGlobalBin(iHLT, iEta, ipt)) );
    } // end of pt loop
 
    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