DAS::JetEnergy Namespace Reference

Classes
struct	AbstractFit
struct	BalanceHists
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})
std::vector< double >	getBinning (int nBins, float first, float last)
void	assertValidBinning (const std::vector< double > &v)
std::vector< float >	GetLogBinning (float minpt, float maxpt, int nbins)
std::string	getAlgo (const Darwin::Tools::MetaInfo &metainfo)
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})
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)
void	ProcessVariations (const vector< TH3 * > &variations, TDirectory &dOut)
void	loopDirsFromGetBalance (TDirectory *dIn, TDirectory *dOut)
void	getPtAveBalance (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering)
template<typename Jet >
bool	DijetSelection (const std::vector< Jet > *jets, int iJES=0, const double alpha=0.3)
vector< TString >	GetScaleVariations (const DT::MetaInfo &metainfo)
void	getPtBalance (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})

Variables
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::map< int, int >	nRhoBins
static const std::map< int, std::vector< TString > >	rhoBins
static const std::vector< double >	eta_unc_edges = {-5.4, -5.0, -4.4, -4.0, -3.5, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.2 , 2.4 , 2.6 , 2.8 , 3.0, 3.5 , 4.0 , 4.4 , 5.0 , 5.4}
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"}
static const float	w = 0.8
static const auto	deps = numeric_limits<double>::epsilon()

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;
 
    unique_ptr<TChain> tIn = DT::GetChain(inputs);
    unique_ptr<TFile> fOut(DT_GetOutput(output));
    auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    auto R = metainfo.Get<int>("flags", "R");
    if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
                                      make_unique<TFile>(inputs.front().c_str() )) );
    JMEmatching<vector<RecJet>, vector<GenJet>>::maxDR = (R/10.)/2;
 
    GenEvent * genEvt = nullptr;
    RecEvent * recEvt = nullptr;
    PileUp * pu = nullptr;
    tIn->SetBranchAddress("genEvent", &genEvt);
    tIn->SetBranchAddress("recEvent", &recEvt);
    tIn->SetBranchAddress("pileup", &pu);
 
    vector<GenJet> * genJets = nullptr;
    vector<RecJet> * recJets = nullptr;
    tIn->SetBranchAddress("genJets", &genJets);
    tIn->SetBranchAddress("recJets", &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 exists.",
                    tables, make_error_code(errc::no_such_file_or_directory)) );
    const auto tag = config.get<string>("corrections.JER.tag"); // correction tag
    metainfo.Set<string>("corrections", "JER", "tag", tag);
    const auto& core_width = config.get<float>("corrections.JER.core_width");
 
    metainfo.Set<float>("corrections", "JER", "core_width", core_width);
 
    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 algo = getAlgo(metainfo);
    string keySF  = tag +  "_ScaleFactor_AK" + to_string(R < 6 ? 4 : 8) + "PF" + algo,
           keyJER = tag + "_PtResolution_AK" + to_string(R < 6 ? 4 : 8) + "PF" + algo;
    correction::Correction::Ref SF, JER;
    try {
        SF = cset->at(keySF );
        JER = cset->at(keyJER);
    }
    catch (out_of_range) {
        BOOST_THROW_EXCEPTION( invalid_argument("`" + keySF  + "` or `" + keyJER +
                                                "` could not be found in " + tables.string()) );
    }
 
    // a few lambda functions
    auto getSFs = [&SF,&vars](const RecJet& recjet) {
        vector<float> sfs;
        for (const string var: vars)
            sfs.push_back(SF->evaluate({recjet.p4.Eta(), var}));
        return sfs;
    };
    auto applySFs = [applySyst](RecJet& recjet, const vector<float>& sfs) {
        float nom_scale = recjet.scales.front();
        for (float& scale: recjet.scales)
            scale *= sfs.front();
 
        if (applySyst) {
            recjet.scales.push_back(nom_scale*sfs[1]);
            recjet.scales.push_back(nom_scale*sfs[2]);
        }
    };
    auto getJER = [pu,&JER](const RecJet& recjet) {
        return JER->evaluate({recjet.p4.Eta(), recjet.CorrPt(), pu->rho});
    };
 
    for (DT::Looper looper(tIn, slice); 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 << "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 = getJER(*rec_it);
 
            cout << "|- genpt = " << genpt << '\n'
                 << "|- recpt = " << recpt << '\n'
                 << "|- response = " << response << '\n'
                 << "|- resolution = " << resolution << endl;
 
            vector<float> sfs = getSFs(*rec_it);
            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;
 
            for (float& sf: sfs) {
                cout << "|- " << sf;
                sf = max(0.f, 1+(sf-1)*(recpt-genpt)/recpt);
                cout << " -> " << sf << endl;
            }
 
            cout << "|- applying SFs" << endl;
            applySFs(*rec_it, sfs);
        }
 
        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) {
 
            float resolution = getJER(*fake_it);
            normal_distribution<float> d{0,resolution};
            static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
            float response = d(m_random_generator);
 
            cout << "|- recpt = " << fake_it->CorrPt() << '\n'
                 << "|- response = " << response << '\n'
                 << "|- resolution = " << resolution << endl;
 
            vector<float> sfs = getSFs(*fake_it);
            for (float& sf: sfs) {
                cout << "|- " << sf;
                sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
                cout << " -> " << sf << endl;
            }
 
            cout << "|- applying SFs" << endl;
 
            applySFs(*fake_it, sfs);
        }
 
        cout << "Sorting by descending pt" << endl;
        sort(recJets->begin(), recJets->end(), pt_sort);
 
        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);
    fOut->cd();
    tOut->Write();
    raw.Write(fOut.get());
    for (auto& c: calib)
        c.Write(fOut.get());
 
    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)

Todo:

CorrPt?

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;
 
    unique_ptr<TChain> tIn = DT::GetChain(inputs);
    unique_ptr<TFile> fOut(DT_GetOutput(output));
    auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    auto R = metainfo.Get<int>("flags", "R");
 
    RecEvent * recEvt = nullptr;
    GenEvent * genEvt = nullptr;
    PileUp * pu = nullptr;
    tIn->SetBranchAddress("recEvent", &recEvt);
    if (isMC)
        tIn->SetBranchAddress("genEvent", &genEvt);
    tIn->SetBranchAddress("pileup", &pu);
 
    vector<RecJet> * recJets = nullptr;
    tIn->SetBranchAddress("recJets", &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)) );
    const auto tag = config.get<string>("corrections.JES.tag"); // correction tag
    metainfo.Set<string>("corrections", "JES", "tag", tag);
 
    ControlPlots::isMC = isMC;
    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 algo = getAlgo(metainfo);
    string key = tag + "_L1L2L3Res_AK" + to_string(R < 6 ? 4 : 8) + "PF" + algo;
    correction::CompoundCorrection::Ref nomSF = cset->compound().at(key);
    vector<correction::Correction::Ref> varSFs;
    varSFs.reserve(JES_variations.size());
    for (string source: JES_variations) {
        key = tag + '_' + source + "_AK" + to_string(R < 6 ? 4 : 8) + "PF" + algo;
        varSFs.push_back(cset->at(key));
    }
 
    for (DT::Looper looper(tIn, slice); 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
            auto corr = nomSF->evaluate({recJet.area, recJet.p4.Eta(), recJet.p4.Pt(), pu->rho});
            cout << recJet.p4.Pt() << ' ' << corr << '\n';
            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) {
                auto var = varSF->evaluate({recJet.p4.Eta(), recJet.p4.Pt()}); 
                recJet.scales.push_back(corr*(1-var));
                recJet.scales.push_back(corr*(1+var));
            }
        }
        cout << flush;
 
        // 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);
    fOut->cd();
    tOut->Write();
    raw.Write(fOut.get());
    for (auto& c: calib)
        c.Write(fOut.get());
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

assertValidBinning()

void DAS::JetEnergy::assertValidBinning ( const std::vector< double > &  v )

inline

{
    assert(v.size() > 1);
    for (size_t i = 1; i<v.size(); ++i) {
        if (v[i] > v[i-1]) continue;
        std::cerr << i << ' ' << v[i] << ' ' << v[i-1] << '\n';
    }
}

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

DijetSelection()

bool DAS::JetEnergy::DijetSelection	(	const std::vector< Jet > *	jets,
		int	iJES = 0,
		const double	alpha = 0.3
	)

Apply dijet topology selection independently from the level.

Parameters

jets	full jet vector
iJES	JES variation index (only for rec-level)
alpha	reject significant extra jets

{
    if (jets->size() < 2) return false;
 
    // back-to-back
    using ROOT::Math::VectorUtil::DeltaPhi;
    if (DeltaPhi(jets->at(0).p4, jets->at(1).p4) < 2.7) return false;
 
    // ECAL acceptance
    if (std::abs(jets->at(0).p4.Eta()) >= 3.0 ||
        std::abs(jets->at(1).p4.Eta()) >= 3.0) return false;
 
    // impact of extra jets
    if (jets->size() > 2) {
        if constexpr (is_same<Jet,RecJet>()) {
            const auto pt0 = jets->at(0).CorrPt(iJES),
                       pt1 = jets->at(1).CorrPt(iJES),
                       pt2 = jets->at(2).CorrPt(iJES);
            if (pt2 > alpha*(pt0 + pt1)/2) return false;
        }
        else {
            const auto pt0 = jets->at(0).p4.Pt(),
                       pt1 = jets->at(1).p4.Pt(),
                       pt2 = jets->at(2).p4.Pt();
            if (pt2 > alpha*(pt0 + pt1)/2) return false;
        }
    }
 
    return true;
}

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");
    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();
        int nQuantiles = 1;
        double median;
        double probability[1] = {0.5};
        h->GetQuantiles(nQuantiles, &median, probability);
 
        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(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);
    }
 
    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::MetaInfo &  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>("flag", "year");
    if (year > 2019) return "puppi"; // run 3
    if (year > 2014) return "chs"; // run 2
                     return ""; // run 1
}

getBinning()

std::vector<double> DAS::JetEnergy::getBinning ( int  nBins, float  first, float  last  )

inline

{
    std::vector<double> bins(nBins+1);
    for(int i = 0; i <= nBins; ++i)
        bins[i] = first + ((last-first)/nBins) * i;
    return bins;
}

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;
 
    unique_ptr<TChain> tIn = DT::GetChain(inputs);
    unique_ptr<TFile> fOut(DT_GetOutput(output));
    auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
 
    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;
 
    RecEvent * rEv = nullptr;
    tIn->SetBranchAddress("recEvent", &rEv);
    const auto hiPU = config.get<bool>("skims.pileup");
    PileUp * pileUp = nullptr;
    if (hiPU)
        tIn->SetBranchAddress("pileup", &pileUp);
 
    vector<RecJet> * recjets = nullptr;
    vector<GenJet> * genjets = nullptr;
    tIn->SetBranchAddress("recJets", &recjets);
    tIn->SetBranchAddress("genJets", &genjets);
 
    unique_ptr<TH3> inclResp = makeRespHist("inclusive", "JER"); // 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");
    if (hiPU) {
        rhobinsResp.reserve(nRhoBins.at(year));
         mubinsResp.reserve(nRhoBins.at(year));
        for (int rhobin = 1; rhobin <= nRhoBins.at(year); ++rhobin) {
            rhobinsResp.push_back( makeRespHist( Form("rhobin%d", rhobin) , "JER") );
             mubinsResp.push_back( makeRespHist( Form( "mubin%d", rhobin) , "JER") );
        }
    }
 
    for (DT::Looper looper(tIn, slice); 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) {
            static const size_t nrho = static_cast<size_t>(nRhoBins.at(year));
            for (irho = 0; pileUp->rho > rho_edges.at(year).at(*irho+1) && *irho < nrho; ++*irho);
            if (*irho >= nrho)
                BOOST_THROW_EXCEPTION( DE::AnomalousEvent(Form("irho = %lu > nRhoBins.at(%d) = %d",
                                                            *irho, year, nRhoBins.at(year)), 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
    inclResp->SetDirectory(fOut.get());
    inclResp->SetTitle("Response (inclusive)");
    inclResp->Write("Response");
    if (hiPU)
    for (int irho = 0; irho < nRhoBins.at(year); ++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.get());
        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.get());
        mubinsResp.at(irho)->SetTitle("Response (" + title + ")");
        mubinsResp.at(irho)->Write("Response");
    }
 
    fOut->cd();
    metainfo.Set<bool>("git", "complete", true);
    tOut->Write();
 
    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:

identify algo

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;
 
    unique_ptr<TFile> fOut(DT_GetOutput(output));
 
    auto R = config.get<int>("flags.R");
    auto year = config.get<int>("flags.year");
    string algo = "chs"; 
    cout << R << ' ' << year << ' ' << algo << endl;
 
    const auto tag = config.get<string>("corrections.JES.tag"), // correction tag
               level = config.get<string>("corrections.JES.level"); // correction level
    string key = tag + "_" + level + "_AK" + to_string(R) + "PF" + algo;
    cout << key << endl;
 
    auto cset = correction::CorrectionSet::from_file(input.string());
    correction::Correction::Ref sf;
    try {
        sf = cset->at(key);
    }
    catch (std::out_of_range& e) {
        BOOST_THROW_EXCEPTION( std::invalid_argument("Nothing corresponding to that key") );
    }
 
    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
    }
    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 == "ScaleFactor") {
        h = make_unique<TH1D>(level.c_str(), level.c_str(),
                                   nAbsEtaBins, abseta_edges .data()); 
        named_inputs["systematic"] = "nom";
    }
    else // assuming 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) {
        // 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")
            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;
}

GetLogBinning()

std::vector<float> DAS::JetEnergy::GetLogBinning ( float  minpt, float  maxpt, int  nbins  )

inline

Make equidistant binning in log pt.

Todo:

this accumulates imprecision --> improve

{
    assert(maxpt > minpt);
    assert(minpt > 0);
    assert(nbins > 1);
 
    std::vector<float> edges;
 
    float R = pow(maxpt/minpt,1./nbins);
    for (float pt = minpt; pt <= maxpt; pt *= R) edges.push_back(pt);
    std::cout << edges.size() << std::endl;
 
    return edges;
}

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

getPtAveBalance()

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

This function iterate over the histograms produced by getPtBalance and produces the the momentum balance histograms in bins of ptave and eta. It calls the loopDir() function defined above.

Todo:

use config and steering (e.g. for verbosity)

Parameters

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

{
    auto fIn  = make_unique<TFile>(input .c_str(), "READ"),
         fOut = make_unique<TFile>(output.c_str(), "RECREATE");
    JetEnergy::loopDirsFromGetBalance(fIn.get(), fOut.get());
}

getPtBalance()

void DAS::JetEnergy::getPtBalance	(	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.

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;
 
    unique_ptr<TChain> tIn = DT::GetChain(inputs);
    unique_ptr<TFile> fOut(DT_GetOutput(output));
    auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
 
    RecEvent * rEv = nullptr;
    GenEvent * gEv = nullptr;
    tIn->SetBranchAddress("recEvent", &rEv);
    if (isMC)
        tIn->SetBranchAddress("genEvent", &gEv);
    vector<RecJet> * recjets = nullptr;
    vector<GenJet> * genjets = nullptr;
    tIn->SetBranchAddress("recJets", &recjets);
    if (isMC)
        tIn->SetBranchAddress("genJets", &genjets);
 
    const bool applySyst = (steering & DT::syst) == DT::syst;
    const auto variations = applySyst && !isMC ? GetScaleVariations(metainfo)
                                               : vector<TString>{"nominal"};
 
    // detector level
    vector<BalanceHists<RecJet>> recBalances;
    recBalances.emplace_back("rec_barrel_vs_any", variations);
 
    // hadron level
    vector<BalanceHists<GenJet>> genBalances;
    if (isMC)
        genBalances.emplace_back("gen_barrel_vs_any", vector<TString>{"nominal"});
 
    TRandom3 r3(metainfo.Seed<39856>(slice));
    for (DT::Looper looper(tIn, slice); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
 
        auto gw = isMC ? gEv->weights.front() : Weight{1.,0},
             rw = rEv->weights.front();
 
        // we take the probe and the tag randomly
        auto r = r3.Uniform();
        int itag = (r<0.5),
            iprobe = (r>=0.5);
 
        if (DijetSelection<RecJet>(recjets))
        for (auto bal: recBalances)
            bal(recjets->at(itag), recjets->at(iprobe), gw*rw);
 
        if (!isMC) continue;
 
        if (DijetSelection<GenJet>(genjets))
        for (auto bal: genBalances)
            bal(genjets->at(itag), genjets->at(iprobe), gw);
 
    } // end of event loop
 
    for (auto p: recBalances) p.Write(fOut.get());
    for (auto p: genBalances) p.Write(fOut.get());
 
    fOut->cd();
    metainfo.Set<bool>("git", "complete", true);
    tOut->Write();
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

GetScaleVariations()

vector<TString> DAS::JetEnergy::GetScaleVariations

(

const DT::MetaInfo &

metainfo

)

{
    vector<TString> variations{"nominal"};
    const TList * groupContents = metainfo.List("variations", RecJet::ScaleVar);
    if (!groupContents)
        BOOST_THROW_EXCEPTION( invalid_argument(groupContents->GetName()) );
 
    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.push_back(os->GetString());
    }
    return variations;
}

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

loopDirsFromGetBalance()

void DAS::JetEnergy::loopDirsFromGetBalance	(	TDirectory *	dIn,
		TDirectory *	dOut
	)

This function is called by getPtAveBalance(). Recusively iterates over the internal structure of a file/directory and each it fins a serie of TH3 in the file, it interpret it as the many variations of the jet energy balance and us it to calculate the mean of the balance and also the balance curves in bins of pt and eta.

NOTE: additional remark on the implementation: Once the TH3 representing the binned in pt_balance, eta, y (presumably in the order (x, y, z), the pt_balance is projected into a 1-D histogram for each bin of rapidity and pt_ave (this done for each of the variations). For each bin of rapidity, the mean of the pt_balance is computed (see formula in th comments of ProcessVariations() the standard deviations from the balance are calculated from the statistical uncertainities and from the variations (presumably the variations of the JES).

{
    // Recursively explore the directory present in the file.
    // At each level where there are TH3, they are stored in a vector.
    vector<TH3*> variations;
    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) {
                cerr << red << "ddIn = " << ddIn << '\n' << def;
                continue;
            }
            TDirectory * ddOut = dOut->mkdir(ddIn->GetName(), ddIn->GetTitle());
            ddOut->cd();
            cout << "creating directory " << bold << ddIn->GetPath() << normal << endl;
 
            loopDirsFromGetBalance(ddIn, ddOut);
        }
        else if ( ( TString( key->ReadObj()->ClassName() ).Contains("TH3") ) ) {
            auto bal3D = dynamic_cast<TH3*>( key->ReadObj() );
            bal3D->SetDirectory(0);
            variations.push_back(bal3D);
        }
        else cout << "Found " << orange << key->ReadObj()->GetName() << normal << ".. Is of "
                 << underline << key->ReadObj()->ClassName() << normal << " type, ignoring it" << endl;
    } // End of for (listof keys) loop
 
    if (variations.size() > 0)
        ProcessVariations(variations, *dOut);
}

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, move(hist2D), steering);
            }
        }
        else cout << orange << "Ignoring " << key->ReadObj()->GetName() << " of `"
                  << key->ReadObj()->ClassName() << "` type" << def << endl;
    } // End of for (list of keys) loop
}

ProcessVariations()

void DAS::JetEnergy::ProcessVariations	(	const vector< TH3 * > &	variations,
		TDirectory &	dOut
	)

This function is called by loopDirsFromGetBalance(). It takes as argument a vector of TH3 which presumabilly represents the transverse momentum balance ($p_{T}^{bal}$) histograms in bins of ($p_{T}^{bal}$, $p_{T}^{ave}$, $y_{probe}$) for different variations of the JES. This function project the pt_balance histogram in 1-D for each bin of $p_{T}^{ave}$ and $y_{probe}$. This function calculates the histogram of the mean value of the balance in bins of $p_{T}^{ave}$ for each bin, it calculates the incertainty bands deduced from the statistical variations of the mean and the systematics of the JES variations.

NOTE: for an explicit definition of the variables, see the additional docstring in the core of the code.

The second arguments of this function is the directory in which the histograms are written. This function creates the subdirectory structure necessary to store each series of variations corresponding to each $y_{probe}$ and $p_{T}^{ave}$ bins in a given directory. The structure of the output is:

file/ > dOut/ > eta1/ > nominal (TH1, mean value of the $p_{T}^{bal}$ with statistical unc.)
                      > upper (upper limit of the variations statistical+systematics)
                      > lower (lower limit of the variations staistical + systematics)
                      > ptbin1 > rec_any_vs_any_balance_0
                               > ...
                               > rec_any_vs_any_balance_${nVariations}
                      > ptbin2
                      > ...
                      > ptbin${nPtbins}
              > ...
              > eta${nYbins}

The structure at file/dOut/ level is identical to the structure of the input file.

We use the terminology used in AN2019_167_v13:

• The balance A is : \(A = 2 \frac{p_{T,1}-p_{T,2}}{p_{T,1} + p_{T,2}} \)
• The mean balance is the quantity: \( R = \frac{1 - \fac{ \left\langle A \right\rangle }{2}}{1 + \frac{ \left\langle A \right\rangle }{2}}\)

{
    const int N = variations.size();
    Int_t nBins_local = variations.front()->GetNbinsY();
 
    // Filling an array with the lower edges of each bins in transverse momentum + the upper edge of the last bin.
    double LowEdges[nBins_local] = {0.};
    variations.front()->GetYaxis()->GetLowEdge(LowEdges);
    double genBins_local[nBins_local + 1];
    for (int i = 0 ; i<nBins_local ; ++i) genBins_local[i]=LowEdges[i];
    genBins_local[nBins_local] = variations.front()->GetYaxis()->GetBinUpEdge(nBins_local);
 
    // Iterating over the rapidity bins
    for (int y = 1 ; y <= variations.front()->GetNbinsZ() ; ++y) {
 
        // Create the directory "eta_bin"
        TDirectory *direta = dOut.mkdir( Form("eta%d", y) );
        cout <<"Creating directory " << bold << direta->GetPath() << normal << endl;
        direta->cd();
 
        // `nominal`, `upper` and `lower` contains average value of $p_{T,1}/p_{T,2}$ in bins of $p_{T,ave}$
        TH1 *nominal = new TH1F("nominal", "", nBins_local, genBins_local),
            *upper = new TH1F("upper", "", nBins_local, genBins_local),
            *lower = new TH1F("lower", "", nBins_local, genBins_local);
 
        // Setting the directory in which to write the histograms
        nominal->SetDirectory(direta);
        upper->SetDirectory(direta);
        lower->SetDirectory(direta);
 
        // Setting the title of the axis
        nominal->SetXTitle("p_{T,ave}");
        nominal->SetYTitle("#frac{1 - #LT #frac{#Delta p_{T}}{p_{T,ave}} #GT}{1 + #LT #frac{#Delta p_{T}}{p_{T,ave}} #GT }");
        upper->SetXTitle("p_{T,ave}");
        upper->SetYTitle("#frac{1 - #LT #frac{#Delta p_{T}}{p_{T,ave}} #GT}{1 + #LT #frac{#Delta p_{T}}{p_{T,ave}} #GT }");
        lower->SetXTitle("p_{T,ave}");
        lower->SetYTitle("#frac{1 - #LT #frac{#Delta p_{T}}{p_{T,ave}} #GT}{1 + #LT #frac{#Delta p_{T}}{p_{T,ave}} #GT }");
 
        // Iterating over the pT_Ave bins.
        for (int ipt = 1 ; ipt <= nBins_local ; ++ipt) {
 
            // `errUp` contains the error coming from the 'above' variations.
            // `errDn` contains the error coming from the 'under' variations.
            double errUp = variations.front()->GetBinError(ipt),
                   errDn = errUp, content=0;
 
            TDirectory *dirpt = direta->mkdir( Form("ptbin%d",ipt) );
            dirpt->cd();
            cout <<"Creating directory " << bold << dirpt->GetPath() << normal << endl;
 
            // Iterating over the variations.
            for (int sv=0 ; sv<N ; ++sv) {
                auto balance = variations[sv]->ProjectionX("py", ipt, ipt, y, y, "e");
                balance->SetName(variations[sv]->GetName());
                balance->SetTitle(variations[sv]->GetTitle());
                balance->SetDirectory(dirpt);
                double I = balance->Integral();
 
                if (std::abs(I)<deps) {
                    cout << "integral is null in bin (pt, eta): (" << ipt << "," << y << ") in variation: " << sv << endl;
                    continue;
                }
                else balance->Scale(1./I);
 
                double mean = balance->GetMean();
                double meanErr = balance->GetMeanError();
 
                if (sv == 0) {
                    double R = (1.-mean/2.)/(1.+mean/2.),
                           Rup = (1.-(mean+meanErr)/2.)/(1.+(mean+meanErr)/2.),
                           Rdn = (1.-(mean-meanErr)/2.)/(1.+(mean-meanErr)/2.);
                    // The nominal value is set once for each bins of \f$p_{T,ave}\f$
                    // while iterating over the variations.
                    nominal->SetBinContent(ipt, R);
                    nominal->SetBinError(ipt, std::abs(Rup-Rdn)/2.);
 
                    // the nominal content is set once while iterating over the variations
                    content=mean;
                }
                else {
                    // Iteration over the 'non-nominal' variations (sv>1).
                    // For each variation, if the variation is above the nominal content, var-content
                    // is added to errUp. Else var-content is added to errDn.
                    // errUp or errDn are reach via the address `err` which contain
                    // errDn or errUp depending on the need.
                    double& err = mean > content ? errUp : errDn;
                    err = hypot(err, mean - content);
                }
                balance->Write();
            }
            // See comment l.136 (close to the filling of `nominal`.
            // Content contains the nominal mean value of balance in each bins of pt.
            upper->SetBinContent(ipt, (1.-(content+errUp)/2)/(1.+(content+errUp)/2) );
            lower->SetBinContent(ipt, (1.-(content-errDn)/2)/(1.+(content-errDn)/2) );
        }
        direta->cd();
        nominal->Write();
        upper->Write();
        lower->Write();
    }
}

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

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 = {-5.4, -5.0, -4.4, -4.0, -3.5, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.2 , 2.4 , 2.6 , 2.8 , 3.0, 3.5 , 4.0 , 4.4 , 5.0 , 5.4}

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

nRhoBins

const std::map<int,int> nRhoBins

inlinestatic

Initial value:

= {
    {2016, rho_edges.at(2016).size()-1 },
    {2017, rho_edges.at(2017).size()-1 },
    {2018, rho_edges.at(2018).size()-1 }
}

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

rhoBins

const std::map<int,std::vector<TString> > rhoBins

inlinestatic

Initial value:

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

w

const float w = 0.8

inlinestatic

Todo:

check exact list of variations with JME group