DAS::JetEnergy Namespace Reference

Classes
struct	AbstractFit
class	Balance
struct	ResolutionFit
struct	ResponseFit

Functions
void	applyJERsmearing (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	applyJEScorrections (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
vector< DAS::Uncertainties::Variation >	GetScaleVariations (const DT::MetaInfo &metainfo, bool applySyst)
void	FitResolution (TDirectory *dir, unique_ptr< TH1 > h, const int steering)
void	loopDirsFromFitResponse (TDirectory *dIn, TDirectory *dOut, const int steering)
void	fitJetResolution (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering)
vector< double >	GetMergedBinning (TH1 *h, const int threshold=30)
pair< float, float >	findDLogExtrema (const unique_ptr< TH1 > &DLog, const double mu, const pair< float, float > Range)
double	round_to (double value, double precision=1.0)
TH1 *	CumulativeResponse (TH1 *res1D)
void	FitResponseByBin (TDirectory *dir, unique_ptr< TH2 > res2D, const int steering)
void	loopDirsFromGetResponse (TDirectory *dIn, TDirectory *dOut, const int steering, const DT::Slice slice)
void	fitJetResponse (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
bool	passDijetSelection (const auto &dijet, const auto &jets, const Variation &v, const double alpha=0.3)
void	getDijetBalance (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	getJetResponse (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	getJMEtable (const fs::path input, const fs::path output, const pt::ptree &config, const int steering)
vector< double >	GetBinning (TAxis *axis)
int	GetNRhoBins (int year)
std::vector< TString >	rhoBins (int year)
std::string	GetAlgo (const Darwin::Tools::UserInfo &metainfo)
std::string	GetShortCampaign (const std::string &campaign)
std::string	ScanCorrections (const auto &corrections, const std::string &key="")
template<typename CorrectionType >
CorrectionType	GetCorrection (const auto &corrections, const std::string &campaign, const std::string &type, const std::string &level, const std::string &suffix)
float	Evaluate (const auto &correction, std::ostream &cout, const RecJet &recJet, const std::optional< float > &rho={}, const std::optional< GenJet > &genJet={}, const std::optional< RecEvent > &recEvt={}, const std::optional< std::string > &systematic={}, const std::optional< float > &JER={}, const std::optional< float > &JERSF={})
int	GetR (const Darwin::Tools::UserInfo &metainfo)

Variables
static const auto	deps = numeric_limits<double>::epsilon()
static const std::vector< double >	pt_JERC_edges = {15, 17, 20, 23, 27, 30, 35, 40, 45, 57, 72, 90, 120, 150, 200, 300, 400, 550, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500}
static const int	nPtJERCbins = pt_JERC_edges.size()-1
static const std::vector< double >	abseta_edges = { 0.000, 0.261, 0.522, 0.783, 1.044, 1.305, 1.566, 1.740, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 2.964, 3.139, 3.489, 3.839, 5.191 }
static const int	nAbsEtaBins = abseta_edges.size()-1
static const std::vector< TString >	absetaBins = MakeTitle(abseta_edges, "|#eta|", false, true, [](double v) { return Form("%.3f", v);})
static const std::map< int, std::vector< double > >	rho_edges
static const std::vector< double >	eta_unc_edges = MakeBinning({-5.4, 5.4}, 0.2, BinningOptions::step)
static const int	nEtaUncBins = eta_unc_edges.size()-1
static const std::vector< std::string >	JES_variations {"AbsoluteStat", "AbsoluteScale", "AbsoluteMPFBias", "Fragmentation", "SinglePionECAL", "SinglePionHCAL", "FlavorQCD", "RelativeJEREC1", "RelativeJEREC2", "RelativeJERHF", "RelativePtBB", "RelativePtEC1", "RelativePtEC2", "RelativePtHF", "RelativeBal", "RelativeSample", "RelativeFSR", "RelativeStatFSR", "RelativeStatEC", "RelativeStatHF", "PileUpDataMC", "PileUpPtRef", "PileUpPtBB", "PileUpPtEC1", "PileUpPtEC2", "PileUpPtHF"}

Function Documentation

applyJERsmearing()

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

Get JER smearing corrections in the form of factors to store in the scales member of RecJet.

Todo:

investigate JER uncertainty breakdown (see Core#35)

Todo:

fix in case of offset (e.g. JES applied on MC)

Parameters

inputs	input ROOT files (n-tuples)
output	name of output ROOT files (n-tuples)
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");
    auto R = GetR(metainfo);
    if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
                                      make_unique<TFile>(inputs.front().c_str() )) );
    auto year = metainfo.Get<int>("flags", "year");
    JMEmatching<vector<RecJet>, vector<GenJet>>::maxDR = (R/10.)/2;
 
    auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto pu     = flow.GetBranchReadOnly<PileUp  >("pileup"  );
    auto genJets = flow.GetBranchReadOnly <vector<GenJet>>("genJets");
    auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
 
    const auto& tables = config.get<fs::path>("corrections.JER.tables");
    metainfo.Set<fs::path>("corrections", "JER", "tables", tables);
    if (!fs::exists(tables))
        BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exist.",
                    tables, make_error_code(errc::no_such_file_or_directory)) );
    const auto campaign = config.get<string>("corrections.JER.campaign"); // correction campaign
    metainfo.Set<string>("corrections", "JER", "campaign", campaign);
    const auto& core_width = config.get<float>("corrections.JER.core_width");
 
    metainfo.Set<float>("corrections", "JER", "core_width", core_width);
    string algo = GetAlgo(metainfo);
 
    ControlPlots::isMC = isMC;
    ControlPlots raw("raw");
    vector<ControlPlots> calib { ControlPlots("nominal") };
 
    const bool applySyst = steering & DT::syst;
    vector<string> vars {"nom"}; // keyword for correctionlib
    if (applySyst) {
        metainfo.Set<string>("variations", RecJet::ScaleVar, "JER" + SysDown);
        metainfo.Set<string>("variations", RecJet::ScaleVar, "JER" + SysUp);
 
        calib.push_back(ControlPlots("JER" + SysDown));
        calib.push_back(ControlPlots("JER" + SysUp));
 
        vars.push_back("up");
        vars.push_back("down");
    }
 
    // preparing correctionlib object
    auto cset = correction::CorrectionSet::from_file(tables.string());
    string suffix = "AK" + to_string(R) + "PF" + algo;
 
    auto SF = GetCorrection<correction::Correction::Ref>(
                            *cset, campaign, "MC", "ScaleFactor", suffix),
         JER = GetCorrection<correction::Correction::Ref>(
                            *cset, campaign, "MC", "PtResolution", suffix);
 
    const auto& smear_tables = config.get<fs::path>("corrections.JER.smear");
    metainfo.Set<fs::path>("corrections", "JER", "smear", smear_tables);
 
    // preparing the JER smearing from JME
    if (!fs::exists(smear_tables))
        BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exist.",
                    smear_tables, make_error_code(errc::no_such_file_or_directory)) );
    correction::Correction::Ref smear;
    bool JME_smearing = smear_tables != "/dev/null";
    if (JME_smearing) {
        auto cset_smear = correction::CorrectionSet::from_file(smear_tables.string());
        smear = cset_smear->at("JERSmear");
    }
 
    // a few lambda functions
    auto getSFs = [year,&SF,&vars](const RecJet& recJet, ostream& cout) {
        vector<float> sfs;
        sfs.reserve(vars.size());
        for (const string& var: vars)
            sfs.push_back(Evaluate(SF, cout, recJet, {}, {}, {}, var));
        return sfs;
    };
    auto applySFs = [applySyst](RecJet& recJet, const vector<float>& sfs, ostream& cout) {
        cout << "applying SFs: " << recJet.CorrPt();
        float nom_scale = recJet.scales.front();
        for (float& scale: recJet.scales)
            scale *= sfs.front();
        float recpt = recJet.CorrPt();
        cout << " -> " << recpt;
 
        if (applySyst) {
            recJet.scales.push_back(nom_scale*sfs[1]);
            recJet.scales.push_back(nom_scale*sfs[2]);
            cout << ' ' << recpt * sfs[1]
                 << ' ' << recpt * sfs[2] << endl;
        }
    };
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
        cout << "===============" << endl;
 
        raw(*genJets, genEvt->weights.front());
        raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
 
        cout << "Running matching" << endl;
        JMEmatching matching(*recJets, *genJets);
        auto fake_its = matching.fake_its; // we extract it because we will modify it
 
        cout << "GenJets:";
        for (const auto& genJet: *genJets)
            cout << '\t' << genJet.p4.Pt();
        cout << "\nRecJets:";
        for (const auto& recJet: *recJets)
            cout << '\t' << recJet.CorrPt();
        cout << "\nCouples:";
        for (const auto& [recJet,genJet]: matching.match_its)
            cout << '\t' << genJet->p4.Pt() << ',' << recJet->CorrPt();
        cout << "\nMisses:";
        for (const auto& genJet: matching.miss_its)
            cout << '\t' << genJet->p4.Pt();
        cout << "\nFakes:";
        for (const auto& recJet: matching.fake_its)
            cout << '\t' << recJet->CorrPt();
        cout << endl;
 
        cout << "Applying scaling smearing for jets matched in the core of the response" << endl;
        for (const auto& [rec_it, gen_it]: matching.match_its) {
            cout << "---------------" << endl;
 
            float genpt = gen_it->p4.Pt(),
                  recpt = rec_it->CorrPt();
 
            float response = recpt / genpt,
                  resolution = Evaluate(JER, cout, *rec_it, pu->rho);;
 
            if (abs(response - 1) > core_width * resolution) {
                cout << "match in the tail" << endl;
                fake_its.push_back(rec_it); // left for later, with stochastic smearing
                continue;
            }
            cout << "match in the core" << endl;
 
            vector<float> sfs = getSFs(*rec_it, cout);
            for (float& sf: sfs) {
                cout << "res sf = " << sf;
                sf = max(0.f, 1+(sf-1)*(recpt-genpt)/recpt);
                cout << " -> pt sf = " << sf << endl;
            }
 
            applySFs(*rec_it, sfs, cout);
        }
 
        cout << "Fakes after identification of matches in tails: ";
        for (const auto& recJet: fake_its)
            cout << '\t' << recJet->CorrPt();
        cout << endl;
 
        cout << "Applying stochastic smearing for all other jets "
                "(either unmatched, or matched in the tails of the response)" << endl;
        for (auto& fake_it: fake_its) {
            cout << "---------------" << endl;
 
            float resolution = Evaluate(JER, cout, *fake_it, pu->rho);;
            vector<float> sfs = getSFs(*fake_it, cout);
 
            for (float& sf: sfs) {
                if (JME_smearing) // recJet, rho, genJet, recEvt, syst, JER, JERSF
                    sf = Evaluate(smear, cout, *fake_it, pu->rho, {}, *recEvt, {}, resolution, sf);
                else {
                    normal_distribution<float> d{0,resolution};
                    static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
                    float response = d(m_random_generator);
                    sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
                    cout << "smearing "
                         << "RecPt=" << fake_it->CorrPt() << ' '
                         << "JER=" << resolution << ' '
                         << "response = " << response << ' '
                         << "-> " << sf << endl;
                }
            }
 
            applySFs(*fake_it, sfs, cout);
        }
 
        cout << "Sorting by descending pt" << endl;
        sort(recJets->begin(), recJets->end(), pt_sort);
 
        cout << "RecJets after smearing:";
        for (const auto& recJet: *recJets)
            cout << '\t' << recJet.CorrPt();
        cout << endl;
 
        for (size_t i = 0; i < calib.size(); ++i) {
            calib.at(i)(*genJets, genEvt->weights.front());
            calib.at(i)(*recJets, genEvt->weights.front() * recEvt->weights.front(), i); 
        }
 
        cout << "Filling" << endl;
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
    raw.Write(fOut);
    for (auto& c: calib)
        c.Write(fOut);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

applyJEScorrections()

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

Apply JES corrections to data and MC *n*tuples, as well as uncertainties in data (all sources are considered separately)

Parameters

inputs	input ROOT files (n-tuples)
output	name of output root file containing the 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");
    auto R = GetR(metainfo);
    auto year = metainfo.Get<int>("flags", "year");
 
    auto genEvt = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto pu     = flow.GetBranchReadOnly<PileUp  >("pileup"  );
    auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
 
    const auto tables = config.get<fs::path>("corrections.JES.tables"); // path to directory containing corrections in JSON format
    metainfo.Set<fs::path>("corrections", "JES", "tables", tables);
    if (!fs::exists(tables))
        BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exists.",
                    tables, make_error_code(errc::no_such_file_or_directory)) );
    auto campaign = config.get<string>("corrections.JES.campaign"); // calibration campaign
    metainfo.Set<string>("corrections", "JES", "campaign", campaign); // e.g. `Summer19UL18_RunA` (everything before "MC" or "DATA")
    string algo = GetAlgo(metainfo);
 
    ControlPlots::isMC = isMC;
    ControlPlots::verbose = steering & DT::verbose;
    ControlPlots raw("raw");
    vector<ControlPlots> calib { ControlPlots("nominal") };
 
    const bool applySyst = steering & DT::syst;
 
    if (applySyst)
    for (string source: JES_variations) {
        metainfo.Set<string>("variations", RecJet::ScaleVar, source + SysDown);
        metainfo.Set<string>("variations", RecJet::ScaleVar, source + SysUp);
 
        calib.push_back(ControlPlots(source + SysDown));
        calib.push_back(ControlPlots(source + SysUp));
    }
 
    // prepare correctionlib
    auto cset = correction::CorrectionSet::from_file(tables.string());
 
    string suffix = "AK" + to_string(R) + "PF" + algo;
    auto nomSF = GetCorrection<correction::CompoundCorrection::Ref>(
            cset->compound(), campaign, isMC ? "MC" : "DATA", "L1L2L3Res", suffix);
    cout << ScanCorrections(cset->compound(), nomSF->name()) << endl;
    // Example:
    // - Summer19UL18_RunA_V6_DATA_L1L2L3Res_AK4PFchs
    // - Summer19UL18_RunB_V6_DATA_L1L2L3Res_AK4PFchs
    // - Summer19UL18_RunC_V6_DATA_L1L2L3Res_AK4PFchs
    // - Summer19UL18_RunD_V6_DATA_L1L2L3Res_AK4PFchs
    // - Summer19UL18_V5_MC_L1L2L3Res_AK4PFchs
 
    if (steering & DT::verbose)
        cout << ScanCorrections(*cset) << endl;
 
    vector<correction::Correction::Ref> varSFs;
    varSFs.reserve(JES_variations.size());
 
    string short_campaign = GetShortCampaign(campaign);
    if (applySyst)
    for (string source: JES_variations)
        varSFs.push_back(GetCorrection<correction::Correction::Ref>(
                        *cset, short_campaign, "MC", source, suffix));
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        double weight = recEvt->weights.front();
        if (isMC) weight *= genEvt->weights.front();
 
        // CP before JES
        raw(*recJets, weight);
 
        for (auto& recJet: *recJets) {
 
            // sanity check
            if (recJet.scales.size() != 1 || recJet.scales.front() != 1)
                BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected JES corrections", tIn) );
 
            // nominal value
            float corr = Evaluate(nomSF, cout, recJet, pu->rho, {}, *recEvt);
            recJet.scales.front() = corr;
 
            if (!applySyst) continue;
 
            // load uncertainties in JECs
            recJet.scales.reserve(JES_variations.size()*2+1);
            for (const correction::Correction::Ref& varSF: varSFs) {
                float var = Evaluate(varSF, cout, recJet);
                recJet.scales.push_back(corr*(1-var));
                recJet.scales.push_back(corr*(1+var));
            }
        }
 
        // CP after corrections
        for (size_t i = 0; i < calib.size(); ++i)
            calib.at(i)(*recJets, weight, i);
 
        sort(recJets->begin(), recJets->end(), pt_sort); // sort again the jets by pt
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
    raw.Write(fOut);
    for (auto& c: calib)
        c.Write(fOut);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

CumulativeResponse()

TH1* DAS::JetEnergy::CumulativeResponse

(

TH1 *

res1D

)

Creating cumulative response histogram. Given a response 1D histogram, it adds to a bin the content of all previous bins.

{
    double cont = 0;
    const int NptBins = res1D->GetNbinsX();
    for (int n = 1; n <= NptBins; ++n) {
        cont+=res1D->GetBinContent(n);
        res1D->SetBinContent(n, cont);
    }
    return res1D;
}

Evaluate()

float DAS::JetEnergy::Evaluate	(	const auto &	correction,
		std::ostream &	cout,
		const RecJet &	recJet,
		const std::optional< float > &	rho = {},
		const std::optional< GenJet > &	genJet = {},
		const std::optional< RecEvent > &	recEvt = {},
		const std::optional< std::string > &	systematic = {},
		const std::optional< float > &	JER = {},
		const std::optional< float > &	JERSF = {}
	)

Wrapper to evaluate scale factors with DAS objects from correctionlib with exceptions.

The implementation is generic, namely for all years, corrections, and versions of the corrections, which may depend on various variables. The variables names (e.g. JetPt) are only applicable to JERC. This brute-force approach is necessary because of the unpredictable changes of conventions.

Exceptions

A	std::invalid_argument is thrown when an unrecognised variable is requested by the tables. It should then be sufficient to edit function.

Todo:

Use boost::core::demangle() to provide a more readable function name?

Parameters

correction	simple or compound correction
cout	output stream, possible redirected
recJet	for all corrections
rho	e.g. for smearing
genJet	for stoch smearing
recEvt	for the run or event number
systematic	e.g. for JER
JER	e.g. for stoch smearing
JERSF	e.g. for stoch smearing

                                                {}, 
                const std::optional<GenJet>& genJet = {}, 
                const std::optional<RecEvent>& recEvt = {}, 
                const std::optional<std::string>& systematic = {}, 
                const std::optional<float>& JER = {}, 
                const std::optional<float>& JERSF = {}  
               )
try {
    using namespace std;
 
    cout << correction->name();
    vector<correction::Variable::Type> inputs;
    inputs.reserve(10);
 
    auto push_back = [&cout,&inputs](const auto& v) {
        inputs.push_back(v);
        cout << '=' << v;
    };
 
    for (const correction::Variable& input: correction->inputs()) {
        const string& n = input.name();
        cout << ' ' << n;
        if      (n == "JetPt"     ) push_back(recJet.CorrPt());
        else if (n == "JetEta"    ) push_back(recJet.p4.Eta());
        else if (n == "JetPhi"    ) push_back(recJet.p4.Phi());
        else if (n == "JetEta"    ) push_back(recJet.p4.Phi());
        else if (n == "JetA"      ) push_back(recJet.area);
        else if (n == "Rho"       ) push_back(*rho);
        else if (n == "systematic") push_back(*systematic);
        else if (n == "GenPt"     ) push_back(genJet ? genJet->CorrPt() : -1);
        else if (n == "EventID"   ) push_back((int) recEvt->evtNo);
        else if (n == "run"       ) push_back((float) recEvt->runNo);
        else if (n == "JER"       ) push_back(*JER);
        else if (n == "JERSF"     ) push_back(*JERSF);
        else BOOST_THROW_EXCEPTION( invalid_argument("`"s + n + "` is needed by "s
                                    + correction->name() + " but not recognized."s) );
    }
    auto corr = correction->evaluate(inputs);
    cout << " -> " << corr << endl;
    return corr;
}
catch (std::runtime_error& e) {
    const char * name = typeid(*correction).name();
    auto location = boost::source_location(__FILE__, __LINE__, name);
    boost::throw_exception(e, location);
}

findDLogExtrema()

pair<float,float> DAS::JetEnergy::findDLogExtrema	(	const unique_ptr< TH1 > &	DLog,
		const double	mu,
		const pair< float, float >	Range
	)

Find the positions of the extrema, starting from the mean and within a given range (to avoid regions with too large fluctuations).

{
    const int m = DLog->FindBin(mu);
 
    // find maximum on the LHS
    int imax = m;
    for (int i = m; i >= DLog->FindBin(Range.first); --i)
        if (DLog->GetBinContent(i) > DLog->GetBinContent(imax)) imax = i;
    float x1 = DLog->GetBinLowEdge(imax+1);
 
    // find minimum on the RHS
    int imin = m;
    for (int i = m; i <= DLog->FindBin(Range.second); ++i)
        if (DLog->GetBinContent(i) < DLog->GetBinContent(imin)) imin = i;
    float x2 = DLog->GetBinLowEdge(imin);
 
    return {x1,x2};
}

fitJetResolution()

void DAS::JetEnergy::fitJetResolution	(	const fs::path &	input,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering
	)

Fit resolution curves from fitJetResponse

Parameters

input	input ROOT file (histograms)
output	name of output ROOT (histograms)
config	config file from `DTOptions`
steering	steering parameters from `DTOptions`

{
    cout << __func__ << " start" << endl;
 
    auto fIn  = make_unique<TFile>(input .c_str(), "READ"),
         fOut = make_unique<TFile>(output.c_str(), "RECREATE"); // TODO
    JetEnergy::loopDirsFromFitResponse(fIn.get(), fOut.get(), steering);
 
    cout << __func__ << " stop" << endl;
}

fitJetResponse()

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

Fit response distributions from getJetResponse

Parameters

input	input ROOT file (histograms)
output	name of output ROOT (histograms)
config	config file from `DTOptions`
steering	steering parameters from `DTOptions`
slice	number and index of slice

                                   {1,0} 
            )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    auto fIn  = make_unique<TFile>(input .c_str(), "READ"),
         fOut = make_unique<TFile>(output.c_str(), "RECREATE");
    JetEnergy::loopDirsFromGetResponse(fIn.get(), fOut.get(), steering, slice);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

FitResolution()

void DAS::JetEnergy::FitResolution	(	TDirectory *	dir,
		unique_ptr< TH1 >	h,
		const int	steering
	)

Fit the resolution obtained from fitJetResponse in a given eta/rho bin.

Todo:

parameters in bins of eta (and rho/mu?)

--> one should probably give the directory with all eta bins as input

{
    dir->cd();
    TString stdPrint = h->GetTitle();
    stdPrint.ReplaceAll(TString( h->GetName() ), "");
    cout << h->GetName() << '\t' << stdPrint << endl;
    auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
 
    ResolutionFit res(h, cout);
 
    res.NSC();
    res.NSCd();
 
}

FitResponseByBin()

void DAS::JetEnergy::FitResponseByBin	(	TDirectory *	dir,
		unique_ptr< TH2 >	res2D,
		const int	steering
	)

Projection 2D histograms onto 1D histograms and fit response

Fit is currently performed in 3 stages:

1. Fit the Gaussian core to get a better estimation of the mean and sigma.
2. Extend the fit right range, transposing to a Single-Sided Crystal Ball fit, using the mean and sigma better define the fit parameters.
3. Extend the fit left range, transposing to a Double-Sided Crystal Ball fit.

Fit iteratively on each stage by adjusting the fit range so that you minimize the amount of failed fits.

{
    dir->cd();
 
    auto h_ptBins = unique_ptr<TH1>( dynamic_cast<TH1*>(res2D->ProjectionX("h_ptBins", 1, -1)) ); // Only for pt binnning extraction
    h_ptBins->GetXaxis()->SetTitle("p_{T}^{gen}");
 
    // Declaring histograms (control plots)
    map<TString, unique_ptr<TH1>> hists;
    for (TString name: {"meanCore", "sigmaCore", // Obtained from DCB fit
                        "kLtail", "kRtail", // DCB tails transitions
                        "nLtail", "nRtail", // DCB tails slope
                        "Lcore", "Rcore", // estimates from derivative
                        "chi2ndf", "Median"}) {
        // for successful fits
        hists.insert( {name, unique_ptr<TH1>(dynamic_cast<TH1*>( h_ptBins->Clone(name) )) } );
        hists[name]->Reset();
        hists[name]->SetTitle(name);
        hists[name]->SetDirectory(nullptr);
 
        // for failed fits
        name += "Fail";
        hists.insert( {name, unique_ptr<TH1>(dynamic_cast<TH1*>( h_ptBins->Clone(name) )) } );
        hists[name]->Reset();
        hists[name]->SetTitle(name);
        hists[name]->SetDirectory(nullptr);
    }
 
    int Npt = res2D->GetXaxis()->GetNbins();
    for (int ptbin = 1; ptbin <= Npt; ++ptbin) { // Repeat fit procedure for all response distributions
        double lowEdge = res2D->GetXaxis()->GetBinLowEdge(ptbin),
               upEdge  = res2D->GetXaxis()->GetBinUpEdge (ptbin);
        const char * level     = TString( res2D->GetXaxis()->GetTitle() ).Contains("gen") ? "gen" : "rec";
        const char * hlt_level = TString( res2D->GetXaxis()->GetTitle() ).Contains("HLT") ? "HLT" : level;
        TString title    = Form("%.0f < p_{T}^{%s} < %.0f", lowEdge, hlt_level, upEdge),
                stdPrint = TString( res2D->GetTitle() );
        stdPrint.ReplaceAll("zx projection", "");
        stdPrint.ReplaceAll(")", "");
        stdPrint += ", " + title + ")";
        cout << stdPrint << def << endl;
        static auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
 
        auto original = unique_ptr<TH1>( dynamic_cast<TH1*>(res2D->ProjectionY("original", ptbin, ptbin)) ); // Response distribution
        if (original->GetEntries() < 100) {
            cout << orange << "Skipped due to low number of events!" << def << endl;
            continue;
        }
 
        vector<double> edges = GetMergedBinning(original.get());
        auto h = unique_ptr<TH1>( original->Rebin(edges.size()-1, "h", edges.data()) );
        const char * h_title = TString( res2D->GetXaxis()->GetTitle() ).Contains("HLT") ? "HLT bin low edge" : "gen";
        h->GetXaxis()->SetTitle(Form("#frac{p_{T}^{rec}}{p_{T}^{%s}}", h_title));
        //h->Rebin(2);
 
        TDirectory * d_pt = dir->mkdir(Form("ptbin%d", ptbin), title);
        d_pt->cd();
        h->SetTitle(title);
        h->SetDirectory(d_pt);
        float normalisation = abs( h->Integral() );
        h->Scale(1.f/normalisation, "width");
        //cout << "normalisation = " << h->Integral() << endl;
        h->Write();
 
        // Calculating median reponse and Interquartile range (IQR) error
        // While the mean has standard error estimates, the median uses the
        // interquartile range to estimate the robustness of the median against outliers or skewed data.
        // https://en.wikipedia.org/wiki/Interquartile_range
        vector<double> probabilities{0.25, 0.5, 0.75};
        const int nQuantiles = probabilities.size();
        vector<double> quantiles(nQuantiles);
        h->GetQuantiles(nQuantiles, quantiles.data(), probabilities.data());
        double& median = quantiles[1];
        double IQRError = M_PI / 2.0 * (quantiles[2] - quantiles[0]) / sqrt(h->GetEntries());
 
        TString hlt_title = Form("%.0f < p_{T}^{%s} < %.0f", lowEdge, hlt_level, upEdge);
 
        // Calculating Cumulative response histogram and storing it to directory
        auto hc_tmp = dynamic_cast<TH1*>( h->Clone("hc_tmp") );
        auto hc = CumulativeResponse(std::move(hc_tmp)); // Cumulative response
        hc->SetNameTitle("h_cumulative", hlt_title);
        hc->SetDirectory(d_pt);
        hc->Write();
 
        // the complex fit is split over several steps using a ResponseFit object
        ResponseFit resp(h, cout);
 
        hists.at("Lcore")->SetBinContent(ptbin, resp.p[ResponseFit::KL]);
        hists.at("Rcore")->SetBinContent(ptbin, resp.p[ResponseFit::KR]);
 
        resp.gausCore();
        resp.SSCB();
        resp.DSCB();
 
        TString suffix = resp.good() ? "" : "Fail";
 
        // Save fit parameters into control plots
        hists.at( "meanCore" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::MU]);
        hists.at( "meanCore" + suffix)->SetBinError  (ptbin, resp.e[ResponseFit::MU]);
        hists.at("sigmaCore" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::SIGMA]);
        hists.at("sigmaCore" + suffix)->SetBinError  (ptbin, resp.e[ResponseFit::SIGMA]);
 
        hists.at("kLtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::KL]);
        hists.at("kLtail" + suffix)->SetBinError  (ptbin, resp.e[ResponseFit::KL]);
        hists.at("nLtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::NL]);
        hists.at("nLtail" + suffix)->SetBinError  (ptbin, resp.e[ResponseFit::NL]);
        hists.at("kRtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::KR]);
        hists.at("kRtail" + suffix)->SetBinError  (ptbin, resp.e[ResponseFit::KR]);
        hists.at("nRtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::NR]);
        hists.at("nRtail" + suffix)->SetBinError  (ptbin, resp.e[ResponseFit::NR]);
 
        hists.at("chi2ndf" + suffix)->SetBinContent(ptbin, resp.chi2ndf   .value_or(0));
        hists.at("chi2ndf" + suffix)->SetBinError  (ptbin, resp.chi2ndfErr.value_or(0));
 
        hists.at("Median" + suffix)->SetBinContent(ptbin, median);
        hists.at("Median" + suffix)->SetBinError(ptbin, IQRError);
    }
 
    dir->cd();
    TString title = res2D->GetTitle();
    title.ReplaceAll("Response distribution", "");
    title.ReplaceAll("zx projection", "");
    for (auto&& h: hists) {
        h.second->SetDirectory(dir);
        h.second->SetTitle(h.second->GetTitle() + title);
        h.second->Write();
    }
}

GetAlgo()

std::string DAS::JetEnergy::GetAlgo ( const Darwin::Tools::UserInfo &  metainfo )

inline

Determine from metainfo if CHS or PUPPI has been used to reconstruct the jets.

{
    if (metainfo.Find("flags", "labels", "CHS"  )) return "chs";
    if (metainfo.Find("flags", "labels", "PUPPI")) return "Puppi";
 
    std::cerr << orange << "Couldn't identify CHS or PUPPI. Running default "
                           "CHS for Run 2 (PUPPI foir Run 3)\n" << def;
 
    auto year = metainfo.Get<int>("flags", "year");
    if (year > 2019) return "Puppi"; // run 3
    if (year > 2014) return "chs"; // run 2
                     return ""; // run 1
}

GetBinning()

vector<double> DAS::JetEnergy::GetBinning

(

TAxis *

axis

)

{
    int n = axis->GetNbins();
    vector<double> bins(n+1);
    for (int i = 0; i <= n; ++i)
        bins[i] = axis->GetBinLowEdge(i+1);
    bins[n] = axis->GetBinUpEdge(n);
    return bins;
}

GetCorrection()

CorrectionType DAS::JetEnergy::GetCorrection	(	const auto &	corrections,
		const std::string &	campaign,
		const std::string &	type,
		const std::string &	level,
		const std::string &	suffix
	)

Returns the correction corresponding to a key.

< expecting correction::[Compound]Correction::Ref

Parameters

corrections	from `correctionlib`
campaign	e.g. Summer19UL18_RunA
type	MC or DATA
level	e.g. L1L2L3Res, TimePtEta
suffix	e.g. AK4chs

{
    namespace al = boost::algorithm;
    for (const auto& correction: corrections) {
        if (!al::starts_with(correction.first, campaign) ||
            !al::contains   (correction.first, type) ||
            !al::contains   (correction.first, level) ||
            !al::ends_with  (correction.first, suffix))
            continue;
 
        return correction.second;
    }
    using namespace std;
    BOOST_THROW_EXCEPTION(
            invalid_argument(
                "No `"s + campaign + "*"s + type + "*" + level + "*"s + suffix
                + "` correction can be found in the given tables. "s
                + ScanCorrections(corrections)
            )
        );
}

getDijetBalance()

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

Get Pt balance from dijet events.

Parameters

inputs	input ROOT files (n-tuples)
output	name of output root file containing the 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");
    const auto alpha = config.get<float>("skim.dijet.alpha");
    metainfo.Set<float>("skim", "dijet", "alpha", alpha);
 
    auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent"): nullptr;
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    
    auto genJets = isMC ? flow.GetBranchReadOnly<vector<GenJet>>("genJets") : nullptr;
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
 
    auto genDijet = isMC ? flow.GetBranchReadOnly<GenDijet>("genDijet") : nullptr;
    auto recDijet = flow.GetBranchReadOnly<RecDijet>("recDijet");
 
    const bool applySyst = (steering & DT::syst) && !isMC;
    const auto variations = GetScaleVariations(metainfo, applySyst);
 
    // detector level
    vector<Balance> recBalances;
    TDirectory * rec = fOut->mkdir("rec");
    for (const Variation& v: variations)
        recBalances.emplace_back(rec, v, pt_JERC_edges, abseta_edges);
 
    // particle level
    vector<Balance> genBalances;
    TDirectory * gen = nullptr;
    if (isMC) {
        gen = fOut->mkdir("gen");
        genBalances.emplace_back(gen, variations.front(), pt_JERC_edges, abseta_edges);
    }
 
    TRandom3 r3(metainfo.Seed<39856>(slice));
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        // we take the probe and the tag randomly
        auto r = r3.Uniform();
        int itag =   (r <  0.5),
            iprobe = (r >= 0.5);
 
        for (Balance& bal: recBalances) {
            auto gw = isMC ? gEv->Weight(bal.v) : Weight{1.,0},
                 rw =        rEv->Weight(bal.v);
 
            if (passDijetSelection(*recDijet, *recJets, bal.v, alpha) && 
                (abs((itag ? recDijet->first : recDijet->second)->p4.Eta()) < 1.3))
                bal((itag   ? *recDijet->first : *recDijet->second),
                    (iprobe ? *recDijet->first : *recDijet->second), gw*rw);
        }
 
        if (!isMC) continue;
 
        for (Balance& bal: genBalances) {
            auto gw = gEv->Weight(bal.v);
            if (passDijetSelection(*genDijet, *genJets, bal.v, alpha) && 
                (abs((itag ? genDijet->first : genDijet->second)->p4.Eta()) < 1.3))
                bal((itag   ? *genDijet->first : *genDijet->second),
                    (iprobe ? *genDijet->first : *genDijet->second), gw);
        }
 
    } // end of event loop
 
    for (auto& p: recBalances) p.Write();
    for (auto& p: genBalances) p.Write();
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

getJetResponse()

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

Get JER differential resolution and balance.

Todo:

with gen level weight too?

Todo:

remove upper boundaries from titles of last bins

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 MC may be used as input.", metainfo) );
 
    const auto R = metainfo.Get<int>("flags", "R");
 
    JMEmatching<>::maxDR = R/10./2.;
    cout << "Radius for matching: " << JMEmatching<>::maxDR << endl;
 
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    const auto hiPU = config.get<bool>("skims.pileup");
    auto pileUp = hiPU ? flow.GetBranchReadOnly<PileUp>("pileup") : nullptr;
 
    auto recjets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto genjets = flow.GetBranchReadOnly<vector<GenJet>>("genJets");
 
    unique_ptr<TH3> inclResp = makeRespHist("inclusive"); // Response inclusive in rho
    vector<unique_ptr<TH3>> rhobinsResp, // Response in bins of rho
                            mubinsResp;  // Response in bins of mu (true pileup)
 
    const auto year = metainfo.Get<int>("flags", "year");
    const size_t nRhoBins = GetNRhoBins(year);
    if (hiPU) {
        rhobinsResp.reserve(nRhoBins);
         mubinsResp.reserve(nRhoBins);
        for (int rhobin = 1; rhobin <= nRhoBins; ++rhobin) {
            rhobinsResp.push_back( makeRespHist( Form("rhobin%d", rhobin) ) );
             mubinsResp.push_back( makeRespHist( Form( "mubin%d", rhobin) ) );
        }
    }
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        auto recWgt = rEv->weights.front();
        optional<size_t> irho, imu;
        if (hiPU) {
            for (irho = 0; pileUp->rho > rho_edges.at(year).at(*irho+1)
                    && *irho < nRhoBins; ++*irho);
            if (*irho >= nRhoBins)
                BOOST_THROW_EXCEPTION(
                        DE::AnomalousEvent(Form("irho = %lu > nRhoBins = %zu",
                                                *irho, nRhoBins), tIn) );
            static const size_t imuMax = mubinsResp.size()-1;
            imu = static_cast<size_t>(pileUp->GetTrPU())/10;
            *imu = min(*imu,imuMax);
        }
 
        // Matching
        if (genjets->size() > 3) genjets->resize(3); // Keep only three leading gen jets
        JMEmatching matching(*recjets, *genjets);
        matching.match_its.resize(min(matching.match_its.size(),3lu));
 
        for (auto& [rec_it,gen_it]: matching.match_its) {
            auto ptrec = rec_it->CorrPt();
            auto ptgen = gen_it->p4.Pt();
 
            // Fill resolution from smearing
            auto response = ptrec/ptgen;
            auto etarec   = abs( rec_it->p4.Eta() );
 
            inclResp->Fill(ptgen, etarec, response, recWgt);
            if (irho) rhobinsResp.at(*irho)->Fill(ptgen, etarec, response, recWgt);
            if (imu )  mubinsResp.at(*imu )->Fill(ptgen, etarec, response, recWgt);
 
        } // End of for (pairs) loop
    } // End of while (event) loop
 
    // Creating directory hierarchy and saving histograms
    fOut->cd();
    inclResp->SetDirectory(fOut);
    inclResp->SetTitle("Response (inclusive)");
    inclResp->Write("Response");
    if (hiPU)
    for (int irho = 0; irho < nRhoBins; ++irho) {
 
 
        TString title = Form("%.2f < #rho < %.2f", rho_edges.at(year).at(irho), rho_edges.at(year).at(irho+1));
        fOut->cd();
        TDirectory * d_rho = fOut->mkdir(Form("rhobin%d", irho+1), title);
        d_rho->cd();
        rhobinsResp.at(irho)->SetDirectory(fOut);
        rhobinsResp.at(irho)->SetTitle("Response (" + title + ")");
        rhobinsResp.at(irho)->Write("Response");
 
        title = Form("%d < #mu < %d", irho*10, (irho+1)*10);
        fOut->cd();
        TDirectory * d_mu = fOut->mkdir(Form("mubin%d", irho+1), title);
        d_mu->cd();
        mubinsResp.at(irho)->SetDirectory(fOut);
        mubinsResp.at(irho)->SetTitle("Response (" + title + ")");
        mubinsResp.at(irho)->Write("Response");
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

getJMEtable()

void DAS::JetEnergy::getJMEtable	(	const fs::path	input,
		const fs::path	output,
		const pt::ptree &	config,
		const int	steering
	)

Obtain corrections and uncertainties in ROOT format for plotting and comparison.

Todo:

Compound corrections?

Todo:

Use THn to get one more dimension

Todo:

check binning

Todo:

check binning

Todo:

check binning

Parameters

input	JetMET JES tables in JSON format
output	name of output root file
config	config file from `DTOptions`
steering	steering parameters from `DTOptions`

{
    cout << __func__ << " start" << endl;
 
    auto fOut = DT::GetOutputFile(output);
 
    auto isMC = config.get<bool>("flags.isMC");
    auto R = config.get<int>("flags.R");
    auto year = config.get<int>("flags.year");
    DT::UserInfo userinfo(config);
    string algo = GetAlgo(userinfo);
    cout << R << ' ' << year << ' ' << algo << endl;
 
    auto campaign = config.get<string>("corrections.JME.campaign"); // calibration campaign
    auto level = config.get<string>("corrections.JME.level"); // calibration level
    string suffix = "AK" + to_string(R) + "PF" + algo;
 
    auto cset = correction::CorrectionSet::from_file(input.string());
    auto sf = GetCorrection<correction::Correction::Ref>(
                    *cset, campaign, isMC ? "MC" : "DATA", level, suffix);
 
    cout << "Declaring histogram" << endl;
    unique_ptr<TH1> h;
    map<string, correction::Variable::Type> named_inputs;
    if (level == "L1FastJet" || level == "PtResolution") {
        h = make_unique<TH3D>(level.c_str(), level.c_str(),
                                   nAbsEtaBins, abseta_edges .data(),
                                   nPtJERCbins, pt_JERC_edges.data(),
                                   rho_edges.at(year).size()-1,
                                   rho_edges.at(year).data());
        named_inputs["JetA"] = M_PI * pow(R*0.1, 2); // Jet area estimation
        named_inputs["JetPhi"] = 0; 
    }
    else if (level == "L2Relative" || level == "L2L3Residual")
        h = make_unique<TH2D>(level.c_str(), level.c_str(),
                                   nAbsEtaBins, abseta_edges .data(),
                                   nPtJERCbins, pt_JERC_edges.data());
    else if (level == "PtResolution") {
    // ├── 📈 Summer23BPixPrompt23_RunD_JRV1_MC_PtResolution_AK4PFPuppi (v1)
    // │   PtResolution for AK4PFPuppi jets, created from Summer23BPixPrompt23_RunD_JRV1_MC by using
    // │   https://gitlab.cern.ch/cms-jetmet/jerc2json
    // │   Node counts: Binning: 39, Formula: 305
    // │   ╭────────────────────────────────────────── ▶ input ───────────────────────────────────────────╮
    // │   │ JetEta (real)                                                                                │
    // │   │ pseudorapidity of the jet                                                                    │
    // │   │ Range: [-5.191, 5.191), overflow ok                                                          │
    // │   ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
    // │   ╭────────────────────────────────────────── ▶ input ───────────────────────────────────────────╮
    // │   │ JetPt (real)                                                                                 │
    // │   │ pT of the jet before specific correction (for JER and uncertainties: after all corrections   │
    // │   │ applied)                                                                                     │
    // │   │ Range: [-inf, inf), overflow ok                                                              │
    // │   ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
    // │   ╭────────────────────────────────────────── ▶ input ───────────────────────────────────────────╮
    // │   │ Rho (real)                                                                                   │
    // │   │ energy density rho (as measure of PU)                                                        │
    // │   │ Range: [0.0, 52.05), overflow ok                                                             │
    // │   ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
    // │   ╭──── ◀ output ─────╮
    // │   │ correction (real) │
    // │   │ No description    │
    // │   ╰───────────────────╯
        h = make_unique<TH3D>(level.c_str(), level.c_str(),
                                   nAbsEtaBins, abseta_edges .data(), 
                                   nPtJERCbins, pt_JERC_edges.data(),
                                   rho_edges.at(year).size()-1,
                                   rho_edges.at(year).data());
    }
    else if (level == "ScaleFactor") {
    // ├── 📈 Summer23BPixPrompt23_RunD_JRV1_MC_ScaleFactor_AK4PFPuppi (v1)
    // │   ScaleFactor for AK4PFPuppi jets, created from Summer23BPixPrompt23_RunD_JRV1_MC by using
    // │   https://gitlab.cern.ch/cms-jetmet/jerc2json
    // │   Node counts: Binning: 43, Category: 1470
    // │   ╭────────────────────────────────────────── ▶ input ───────────────────────────────────────────╮
    // │   │ JetEta (real)                                                                                │
    // │   │ pseudorapidity of the jet                                                                    │
    // │   │ Range: [-5.191, 5.191), overflow ok                                                          │
    // │   ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
    // │   ╭────────────────────────────────────────── ▶ input ───────────────────────────────────────────╮
    // │   │ JetPt (real)                                                                                 │
    // │   │ pT of the jet before specific correction (for JER and uncertainties: after all corrections   │
    // │   │ applied)                                                                                     │
    // │   │ Range: [10.0, 5373.0), overflow ok                                                           │
    // │   ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
    // │   ╭────────────────────────────────────────── ▶ input ───────────────────────────────────────────╮
    // │   │ systematic (string)                                                                          │
    // │   │ systematics: nom, up, down                                                                   │
    // │   │ Values: down, nom, up                                                                        │
    // │   ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
    // │   ╭──── ◀ output ─────╮
    // │   │ correction (real) │
    // │   │ No description    │
    // │   ╰───────────────────╯
    //
        if (year < 2020) // run 2
            h = make_unique<TH1D>(level.c_str(), level.c_str(),
                                   nAbsEtaBins, abseta_edges.data()); 
        else // run 3
            h = make_unique<TH2D>(level.c_str(), level.c_str(),
                                   nAbsEtaBins, abseta_edges.data(), 
                                   nPtJERCbins, pt_JERC_edges.data());
        named_inputs["systematic"] = "nom";
    }
    else // assuming JES uncertainties
        h = make_unique<TH2D>(level.c_str(), level.c_str(),
                                   nEtaUncBins, eta_unc_edges.data(),
                                   nPtJERCbins, pt_JERC_edges.data());
 
    cout << "dim = " << h->GetDimension() << endl;
    for (int etabin = 1; etabin <= h->GetNbinsX(); ++etabin)
    for (int  ptbin = 1;  ptbin <= h->GetNbinsY(); ++ ptbin)
    for (int rhobin = 1; rhobin <= h->GetNbinsZ(); ++rhobin) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        // set values by name
        named_inputs["JetEta"] = h->GetXaxis()->GetBinCenter(etabin);
        named_inputs["JetPt" ] = h->GetYaxis()->GetBinCenter( ptbin);
        named_inputs["Rho"   ] = h->GetZaxis()->GetBinCenter(rhobin);
 
        // order the input value in the necessary order
        vector<correction::Variable::Type> inputs;
        for (const correction::Variable& input: sf->inputs())
            inputs.push_back(named_inputs.at(input.name()));
 
        // get and fill the value
 
        double value = sf->evaluate(inputs);
 
        cout << setw(5) << ptbin << setw(5) << etabin << setw(5) << rhobin << setw(15) << value << endl;
 
        if (level == "L1FastJet" || level == "PtResolution")
            h->SetBinContent(etabin, ptbin, rhobin, value);
        else if (level == "ScaleFactor")
            h->SetBinContent(etabin, value);
        else
            h->SetBinContent(etabin, ptbin, value);
    }
    h->Write();
 
    cout << __func__ << " end" << endl;
}

GetMergedBinning()

vector<double> DAS::JetEnergy::GetMergedBinning	(	TH1 *	h,
		const int	threshold = 30
	)

Provide a new binning ensuring a sufficient number of entries in each bin to assume Gaussian uncertainties.

For a gaussian distribution, as it is assumed for the number of events in each bin, the error is given as sqrt(n) where n the number of events in the bin. Assuming both the bin content and error have been normalised, by taking the division of the (BinContent/BinError)^2 ~ n, one can estimate the number of events in that bin.

There is no special treatment for negative weights, as the samples used to develop the code have no such weight. A warning is however displayed in case negative bins are found.

{
    const int N = h->GetNbinsX();
    vector<double> edges(1,0);
    edges.reserve(N);
    int sum = 0;
    for (int i = 1; i <= N; ++i) {
        auto content = h->GetBinContent(i);
        if (content == 0) continue;
        if (content < 0)
            cerr << orange << "Found negative bins in "
                           << h->GetName() << def << '\n';
        auto error = h->GetBinError(i);
        sum += pow(content/error,2);
        if (sum < threshold) continue;
        sum = 0;
        edges.push_back(h->GetBinLowEdge(i+1));
    }
    if (edges.back() != h->GetBinLowEdge(N+1))
        edges.push_back(h->GetBinLowEdge(N+1));
 
    return edges;
}

GetNRhoBins()

int DAS::JetEnergy::GetNRhoBins ( int  year )

inline

{ return rho_edges.at(year).size()-1; }

GetR()

int DAS::JetEnergy::GetR

(

const Darwin::Tools::UserInfo &

metainfo

)

Determine the R values for which JetMET provides corrections.

{
    auto year = metainfo.Get<int>("flags", "year"),
         R = metainfo.Get<int>("flags", "R");
 
    std::string warning = "Not a standard jet size.\n";
    if (year > 2014) { // run 2 and run 3
        if (R != 4 && R != 8)
            std::cerr << orange << warning << def;
        return R < 6 ? 4 : 8;
    }
    else { // run 1
        if (R != 5 && R != 7)
            std::cerr << orange << warning << def;
        return R < 6 ? 5 : 7;
    }
}

GetScaleVariations()

vector<DAS::Uncertainties::Variation> DAS::JetEnergy::GetScaleVariations	(	const DT::MetaInfo &	metainfo,
		bool	applySyst
	)

Get scale variations

Note

one may want to get the Rochester variations as well for Z+jet balance

this function may be redundant with DAS::Unfolding::GetVariations()

{
    vector<DAS::Uncertainties::Variation> variations;
    variations.emplace_back("nominal","nominal");
    if (!applySyst) return variations;
 
    const TList * groupContents = metainfo.List("variations", DAS::RecJet::ScaleVar);
    if (!groupContents)
        BOOST_THROW_EXCEPTION( invalid_argument(groupContents->GetName()) );
 
    size_t i = 0;
    for (const TObject * obj: *groupContents) {
        const auto os = dynamic_cast<const TObjString*>(obj);
        if (!os) BOOST_THROW_EXCEPTION( invalid_argument(obj->GetName()) );
        TString s = os->GetString();
        s.ReplaceAll(SysUp, "");
        s.ReplaceAll(SysDown, "");
        if (find(begin(JES_variations), end(JES_variations), s) != end(JES_variations))
            variations.emplace_back(DAS::RecJet::ScaleVar, os->GetString(), ++i);
    }
    return variations;
}

GetShortCampaign()

std::string DAS::JetEnergy::GetShortCampaign

(

const std::string &

campaign

)

Extracts for isntance Summer19UL18 from Summer19UL18_RunA

{
    using namespace std;
    regex r("^(Summer|Fall|Autumn|Winter|Spring)[0-9]{2}[A-Za-z0-9]*");
    smatch campaign_short;
    if (!regex_search(campaign, campaign_short, r))
        BOOST_THROW_EXCEPTION( invalid_argument("The campaign could not be identified"));
    return campaign_short.str();
}

loopDirsFromFitResponse()

void DAS::JetEnergy::loopDirsFromFitResponse	(	TDirectory *	dIn,
		TDirectory *	dOut,
		const int	steering
	)

Copy directory structure. Given an input directory dIn recursively loop over all directory hierarchy/contents and reproduce it in the directory dOut of the output file. For each TH1 instance found containing the resolution, perform a NSC fit by calling FitResolution() function.

{
    auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
    for (const auto&& obj: *(dIn->GetListOfKeys())) {
        auto const key = dynamic_cast<TKey*>(obj);
        if ( key->IsFolder() ) {
            auto ddIn = dynamic_cast<TDirectory*>( key->ReadObj() );
            if ( TString( ddIn->GetName() ).Contains("ptbin") ) continue;
            TDirectory * ddOut = dOut->mkdir(ddIn->GetName(), ddIn->GetTitle());
            ddOut->cd();
            cout << bold << ddIn->GetPath() << def << endl;
 
            loopDirsFromFitResponse(ddIn, ddOut, steering);
        }
        else if ( TString( key->ReadObj()->GetName() ).Contains("sigmaCore") ) {
            auto res = unique_ptr<TH1>( dynamic_cast<TH1*>(key->ReadObj()) );
            res->SetDirectory(dOut);
            dOut->cd();
            res->Write();
 
            if ( TString(res->GetName()).Contains("Fail") ) continue;
 
            // Fit resolution
            FitResolution(dOut, std::move(res), steering);
        }
        else cout << orange << "Ignoring " << key->ReadObj()->GetName() << " of `"
                  << key->ReadObj()->ClassName() << "` type" << def << endl;
    } // End of for (list of keys) loop
    cout << def << flush;
}

loopDirsFromGetResponse()

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

Copy directory structure. Given an input file dIn recursively loop over all directory hierarchy/contents and reproduce it in output file dOut. For each TH3 instance found containing the response perform a DCB fit by calling FitResponseByBin() function.

{
    for (const auto&& obj: *(dIn->GetListOfKeys())) {
        auto const key = dynamic_cast<TKey*>(obj);
        if ( key->IsFolder() ) {
            auto ddIn = dynamic_cast<TDirectory*>( key->ReadObj() );
            if (ddIn == nullptr) continue;
            TDirectory * ddOut = dOut->mkdir(ddIn->GetName(), ddIn->GetTitle());
            ddOut->cd();
            cout << bold << ddIn->GetPath() << def << endl;
 
            loopDirsFromGetResponse(ddIn, ddOut, steering, slice);
        }
        else if ( TString( key->ReadObj()->ClassName() ).Contains("TH3") ) {
            // Assuming for axis, x: pt bins, y: eta bins, z: resolution bins
            auto hist3D = unique_ptr<TH3>( dynamic_cast<TH3*>(key->ReadObj()) );
            hist3D->SetDirectory(0);
            bool isAbsEta = hist3D->GetYaxis()->GetBinLowEdge(1) >= 0;
            TDirectory * ddOut = dOut->mkdir(key->GetName(), key->ReadObj()->GetTitle());
            for (int etabin = 1; etabin <= hist3D->GetNbinsY(); ++etabin) {
                static int ieta = 0; // Used for parallelisation
                ++ieta;
                if (slice.second != ieta % slice.first) continue;
 
                float lowedge = hist3D->GetYaxis()->GetBinLowEdge(etabin),
                       upedge = hist3D->GetYaxis()->GetBinUpEdge (etabin);
                TString title = isAbsEta ? Form("%.3f < |#eta| < %.3f", lowedge, upedge)
                                         : Form("%.3f < #eta < %.3f"  , lowedge, upedge);
                hist3D->GetYaxis()->SetRange(etabin, etabin);
                auto hist2D = unique_ptr<TH2>( dynamic_cast<TH2*>(hist3D->Project3D("ZX")) ); // Z vs X, where Z is vertical and X is horizontal axis
                hist2D->SetTitle( TString( hist2D->GetTitle() ).ReplaceAll(") ", ", " + title + ")") );
                ddOut->cd();
                TDirectory * dddOut = ddOut->mkdir(Form("etabin%d", etabin), title);
                cout << bold << dddOut->GetPath() << def << endl;
 
                FitResponseByBin(dddOut, std::move(hist2D), steering);
            }
        }
        else cout << orange << "Ignoring " << key->ReadObj()->GetName() << " of `"
                  << key->ReadObj()->ClassName() << "` type" << def << endl;
    } // End of for (list of keys) loop
}

passDijetSelection()

bool DAS::JetEnergy::passDijetSelection	(	const auto &	dijet,
		const auto &	jets,
		const Variation &	v,
		const double	alpha = 0.3
	)

Apply dijet topology selection independently from the gen or rec level. The balance is part of the standard jet calibration procedure at CMS.

Parameters

dijet	dijet system
jets	full jet vector
v	systematic variation
alpha	reject significant extra jets

{
    if (!dijet) return false;
    if (jets.size() < 2) return false;
 
    // back-to-back
    if (dijet.DeltaPhi() < 2.7 ) return false;
 
    // at least one jet in barrel acceptance
    if (std::abs(dijet.first ->p4.Eta()) >= 1.3 &&
        std::abs(dijet.second->p4.Eta()) >= 1.3) return false;
 
    // impact of extra jets
    if (jets.size() > 2) {
        const auto pt2 = jets.at(2).CorrPt(v);
        if (pt2 > alpha * dijet.HT()) return false;
    }
 
    return true;
}

rhoBins()

std::vector<TString> DAS::JetEnergy::rhoBins ( int  year )

inline

{
    return MakeTitle(rho_edges.at(year), "#rho", false, false,
                        [](double v) { return Form("%.2f", v); });
}

round_to()

double DAS::JetEnergy::round_to ( double  value, double  precision = 1.0  )

inline

Round number value to precision. Recover "lost" precision.

{
    return round(value / precision) * precision;
}

ScanCorrections()

std::string DAS::JetEnergy::ScanCorrections	(	const auto &	corrections,
		const std::string &	key = ""
	)

Prints the available corrections to the standard output. If a key is given as 2nd argument, it is highlighted in the list.

Parameters

corrections	corrections from `correctionlib`
key	key to highlight in the list

{
    std::stringstream s;
    s << "Available corrections:";
    for (const auto& correction: corrections) {
        bool found = correction.first == key;
        if (found) s << bold << green;
        s << ' ' << correction.first;
        if (found) s << def;
    }
    return s.str();
}

Variable Documentation

abseta_edges

const std::vector<double> abseta_edges = { 0.000, 0.261, 0.522, 0.783, 1.044, 1.305, 1.566, 1.740, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 2.964, 3.139, 3.489, 3.839, 5.191 }

inlinestatic

JERC binning (taken from JERCProtoLab repo, /macros/common_info/common_binning.hpp)

absetaBins

const std::vector<TString> absetaBins = MakeTitle(abseta_edges, "|#eta|", false, true, [](double v) { return Form("%.3f", v);})

inlinestatic

deps

const auto deps = numeric_limits<double>::epsilon()

static

eta_unc_edges

const std::vector<double> eta_unc_edges = MakeBinning({-5.4, 5.4}, 0.2, BinningOptions::step)

inlinestatic

binning for JES uncertainties (taken from JES files)

JES_variations

const std::vector<std::string> JES_variations {"AbsoluteStat", "AbsoluteScale", "AbsoluteMPFBias", "Fragmentation", "SinglePionECAL", "SinglePionHCAL", "FlavorQCD", "RelativeJEREC1", "RelativeJEREC2", "RelativeJERHF", "RelativePtBB", "RelativePtEC1", "RelativePtEC2", "RelativePtHF", "RelativeBal", "RelativeSample", "RelativeFSR", "RelativeStatFSR", "RelativeStatEC", "RelativeStatHF", "PileUpDataMC", "PileUpPtRef", "PileUpPtBB", "PileUpPtEC1", "PileUpPtEC2", "PileUpPtHF"}

inlinestatic

nAbsEtaBins

const int nAbsEtaBins = abseta_edges.size()-1

inlinestatic

nEtaUncBins

const int nEtaUncBins = eta_unc_edges.size()-1

inlinestatic

nPtJERCbins

const int nPtJERCbins = pt_JERC_edges.size()-1

inlinestatic

pt_JERC_edges

const std::vector<double> pt_JERC_edges = {15, 17, 20, 23, 27, 30, 35, 40, 45, 57, 72, 90, 120, 150, 200, 300, 400, 550, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500}

inlinestatic

rho_edges

const std::map<int,std::vector<double> > rho_edges

inlinestatic

Initial value:

= {
    { 2016, {0, 6.69, 12.39, 18.09, 23.79, 29.49, 35.19, 40.9, 999} }, 
    { 2017, {0, 7.47, 13.49, 19.52, 25.54, 31.57, 37.59, 999} }, 
    { 2018, {0, 7.35, 13.26, 19.17, 25.08, 30.99, 36.9, 999} }, 
 
    { 2024, {0.0, 8.23, 14.49, 20.75, 27.01, 33.27, 39.53, 45.79, 52.05, 999} }, 
}