DAS Namespace Reference

Namespaces
	BTagging
	Chebyshev
	DoubleCrystalBall
	JetEnergy
	Jets
	JetVeto
	MissingET
	MN
	Muon
	Muons
	Normalisation
	Photon
	Prefiring
	PUprofile
	PUstaub
	Uncertainties
	Unfolding

Classes
class	AbstractEvent
struct	AbstractPhysicsObject
struct	ControlPlots
struct	Di
class	GenericSFApplier
struct	GenEvent
struct	GenJet
struct	GenMuon
struct	GenPhoton
struct	Greta
struct	GretaPolynomial
class	GretaPolynomialReader
struct	Helper
struct	JMEmatching
struct	MET
struct	Parameters
class	PhysicsObject
struct	PileUp
struct	PrimaryVertex
struct	RecEvent
struct	RecJet
struct	RecMuon
struct	RecPhoton
class	Teddy
struct	Thunberg
struct	Trigger
struct	TriggerLumi
struct	Weight

Typedefs
using	GenDijet = Di< GenJet, GenJet >
using	GenDimuon = Di< GenMuon, GenMuon >
using	GenZJet = Di< GenDimuon, GenJet >
using	RecDijet = Di< RecJet, RecJet >
using	RecDimuon = Di< RecMuon, RecMuon >
using	RecZJet = Di< RecDimuon, RecJet >
typedef ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< float > >	FourVector
typedef std::vector< Weight >	Weights

Functions
void	mergeNtuples (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
std::vector< TString >	MakeTitle (const std::vector< double > &edges, const char *v, bool lowEdge, bool highEdge, std::function< const char *(double)> format)
Darwin::Tools::Options	Options (const char *, int=Darwin::Tools::none)
std::unique_ptr< TH1 >	getHistSafe (std::unique_ptr< TFile > &f, const std::string &name)
template<class Object >
Weights &	weightsRef (Object &obj)
double	identity (double x)
template<const size_t d>
TF1 *	GetSmoothFit (TH1 *h, int ifm, int ifM, int nSigmaStop=0, const char *options="Q0SNRE", bool autoStop=false)
template<const size_t d>
TF1 *	GetSmoothFit (TH1 *h, double m, double M, int nSigmaStop=0, const char *options="Q0SNRE", bool autoStop=false)
bool	pt_sort (const PhysicsObject &j1, const PhysicsObject &j2)
template<typename TTreePtr >
bool	branchExists (const TTreePtr &tree, TString brName)
template<typename T >
std::vector< T * >	GetObjects (TDirectory *dir)
TDirectory *	GetDirectory (TDirectory *dir, const std::vector< const char * > &names)
void	EventLoop (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	SimpleExec (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering)
vector< double >	GetTriggerTurnons (const fs::path &f)
void	getMETfraction (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
std::map< int, TriggerLumi >	GetLumiFromFiles (const std::filesystem::path &lumi_file, const std::filesystem::path &turnon_file)
map< int, TriggerLumi >	GetLumiFromFiles (const fs::path &lumi_file, const fs::path &turnon_file)
bool	operator== (const PhysicsObject &l, const PhysicsObject &r)
bool	operator< (const PhysicsObject &l, const PhysicsObject &r)
bool	operator> (const PhysicsObject &l, const PhysicsObject &r)
bool	operator== (const Weight &w, const int v)
bool	operator== (const Weight &w, const float v)
bool	operator== (const Weight &w, const double v)
bool	operator== (const Weight &l, const Weight &r)
double	operator* (const Weight &w, const int v)
double	operator* (const int v, const Weight &w)
double	operator* (const Weight &w, const float v)
double	operator* (const float v, const Weight &w)
double	operator* (const Weight &w, const double v)
double	operator* (const double v, const Weight &w)
Weight &	operator*= (Weight &w, const int v)
Weight &	operator/= (Weight &w, const int v)
Weight &	operator*= (Weight &w, const float v)
Weight &	operator/= (Weight &w, const float v)
Weight &	operator*= (Weight &w, const double v)
Weight &	operator/= (Weight &w, const double v)
double	operator* (const Weight &w1, const Weight &w2)
Weights &	operator*= (Weights &wgts, const int v)
Weights &	operator/= (Weights &wgts, const int v)
Weights &	operator*= (Weights &wgts, const float v)
Weights &	operator/= (Weights &wgts, const float v)
Weights &	operator*= (Weights &wgts, const double v)
Weights &	operator/= (Weights &wgts, const double v)
void	getHighScalePUeventIDs (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	getHighScalePUeventIDsHT (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
pair< int, int >	GetBinRange (TH1 *h)
void	ScaleWidth (TH1 *h, TH2 *cov)
void	getSmoothFits (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
Eigen::VectorXd	identity (const Eigen::VectorXd &x)

Variables
static const std::vector< double >	pthat_edges {30,34,39,44,50,56,63,71,80,89,98,108,120,131,143,156,170,196,226,260,300,336,375,420,470,500,531,564,600,645,693,744,800,846,894,946,1000,1088,1183,1287,1400,1491,1587,1690,1800,1934,2078,2233,2400,2579,2771,2978,3200,3578,4000,4472,5000}
static const std::vector< double >	pt_edges {9, 12, 15,18,21,24,28,32,37,43,49,56,64,74,84,97,114,133,153,174,196,220,245,272,300,330,362,395,430,468,507,548,592,638,686,737,790,846,905,967,1032,1101,1172,1248,1327,1410,1497,1588,1684,1784,1890,2000,2116,2238,2366,2500,2640,2787,2941,3103,3273,3450,3637,3832,4037}
static const std::vector< double >	Mjj_edges {160, 200, 249, 306, 372, 449, 539, 641, 756, 887, 1029, 1187, 1361, 1556, 1769, 2008, 2273, 2572, 2915, 3306, 3754, 4244, 4805, 5374, 6094, 6908, 7861, 8929, 10050 }
static const std::vector< double >	pt_av_edges {147, 175, 207, 243, 284, 329, 380, 437, 499, 569, 646, 732, 827, 931, 1046, 1171, 1307, 1458, 1621, 1806, 2003, 2217, 2453, 2702}
static const std::vector< double >	y_edges {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0}
static const std::vector< double >	n_edges {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}
static const int	nYbins = y_edges.size()-1
static const int	nPtBins = pt_edges.size()-1
static const int	nMjjBins = Mjj_edges.size()-1
static const int	nPtDijetBins = pt_av_edges.size()-1
static const int	maxMult = n_edges.size()-1
static const std::vector< TString >	yBins = MakeTitle(y_edges, "|y|", false, true, [](double v) { return Form("%.1f", v);} )
static const double	ptmin = pt_edges.front()
static const double	ptmax = pt_edges.back()
static const double	ymin = y_edges.front()
static const double	ymax = y_edges.back()
static const int	nPUbins = 100
const std::string	SysUp = "Up"
const std::string	SysDown = "Down"

Typedef Documentation

FourVector

typedef ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< float > > FourVector

GenDijet

using GenDijet = Di<GenJet,GenJet>

GenDimuon

using GenDimuon = Di<GenMuon,GenMuon>

GenZJet

using GenZJet = Di<GenDimuon,GenJet>

RecDijet

using RecDijet = Di<RecJet,RecJet>

RecDimuon

using RecDimuon = Di<RecMuon,RecMuon>

RecZJet

using RecZJet = Di<RecDimuon,RecJet>

Weights

typedef std::vector<Weight> Weights

Function Documentation

branchExists()

bool DAS::branchExists	(	const TTreePtr &	tree,
		TString	brName
	)

Check if branch exists

Loop over branches of tree and check if brName exists

{
    auto brList = tree->GetListOfBranches();
    bool brFound = false;
    for (auto it = brList->begin(); it != brList->end(); ++it) {
        TString name = (*it)->GetName();
        if (name == brName) brFound = true;
    }
    return brFound;
}

EventLoop()

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

Template for function (TODO)

Parameters

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

                                  {1,0} 
            )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    // TODO: fetch info from MetaInfo and/or fill it with new information
 
    auto genJets = isMC ? flow.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
    auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets", DT::facultative);
    // TODO: load any branch you may need
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        // TODO: modify the branches
 
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

GetBinRange()

pair<int,int> DAS::GetBinRange

(

TH1 *

h

)

Find range for non-empty bins in an input histogram.

{
    int N = h->GetNbinsX();
    int im, iM;
    for (im =  1; h->GetBinContent(im) <= 0 && im <= N; ++im);
    for (iM = im; h->GetBinContent(iM) >  0 && iM <= N; ++iM);
    --iM;
    return {im,iM};
}

GetDirectory()

TDirectory* DAS::GetDirectory ( TDirectory *  dir, const std::vector< const char * > &  names  )

inline

Get ((...)sub)subdirectory of dir.

{
    for (const char * name: names) {
        TDirectory * subdir = dir->GetDirectory(name);
        if (!subdir)
            BOOST_THROW_EXCEPTION( std::invalid_argument(
                    Form("`%s` could not be found in `%s`", name, dir->GetName())) );
        dir = subdir;
    }
    return dir;
}

getHighScalePUeventIDs()

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

Make a 2-column list of MB event IDs with the max genpt (if any).

Unlike most commands, this command should be run on the output of the n-tupliser directly, and only makes sense (even though it would technically work) for MB samples without pileup.

Todo:

include test in CI

Parameters

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

                                  {1,0} 
            )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice, "ntupliser/events");
 
    auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto genJets = flow.GetBranchReadOnly<vector<GenJet>>("genJets");
 
    const auto minpt = config.get<float>("corrections.PUstaub.minpt");
 
    unordered_map<unsigned long long, float> maxgenpts;
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (genJets->size() == 0) continue;
        float pthat = genEvt->hard_scale,
              maxgenpt = genJets->front().p4.Pt();
 
        if (pthat > 0 ? (maxgenpt > 3*pthat) // refinement to remove crazy parton shower
                      : (maxgenpt > minpt)) // refinment to remove diffractive events
                maxgenpts[recEvt->evtNo] = maxgenpt;
    }
 
    pt::ptree prop_tree;
    for (const auto& evtID_pt: maxgenpts) {
        auto evtID = Form("%llu", evtID_pt.first);
        auto pt    =              evtID_pt.second;
        prop_tree.put<float>(evtID, pt);
    }
    write_info(output.c_str(), prop_tree);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

getHighScalePUeventIDsHT()

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

Make a 2-column list of MB event IDs with the max genht (if any).

Unlike most commands, this command should be run on the output of the n-tupliser directly, and only makes sense (even though it would technically work) for MB samples without pileup.

Todo:

include test in CI

Parameters

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

                                  {1,0} 
            )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice, "ntupliser/events");
 
    auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto genJets = flow.GetBranchReadOnly<vector<GenJet>>("genJets");
 
    const auto minht = config.get<float>("corrections.PUstaub.minht");
 
    unordered_map<unsigned long long, float> hts;
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (genJets->size() == 0) continue;
 
        float ht = 0;
        for (const auto& jet : *genJets) {
            ht += jet.p4.Pt();
        }
 
        if (ht > minht)
            hts[recEvt->evtNo] = ht;
 
    }
 
    pt::ptree prop_tree;
    for (const auto& evtID_ht: hts) {
        auto evtID = Form("%llu", evtID_ht.first);
        auto ht    =              evtID_ht.second;
        prop_tree.put<float>(evtID, ht);
    }
    write_info(output.c_str(), prop_tree);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

getHistSafe()

std::unique_ptr<TH1> DAS::getHistSafe ( std::unique_ptr< TFile > &  f, const std::string &  name  )

inline

Gets an histogram from a TFile.

Exceptions

Darwin::Exceptions::BadInput

if the histogram doesn't exist

Todo:

Move to Darwin?

{
    auto h = std::unique_ptr<TH1>(f->Get<TH1>(name.c_str()));
    if (!h) {
        std::cerr << name << std::endl;
        BOOST_THROW_EXCEPTION( Darwin::Exceptions::BadInput("Can't find required histogram", f) );
    }
    h->SetDirectory(nullptr);
    return h;
}

GetLumiFromFiles() [1/2]

map<int, TriggerLumi> DAS::GetLumiFromFiles	(	const fs::path &	lumi_file,
		const fs::path &	turnon_file
	)

GetLumiFromFiles returns a map containing key value pairs corresponding to the trigger threshold at HLT (e.g. 40, 60, etc.) and the object containing the turn-on point and its effective luminosity.

{
    map<int, TriggerLumi> trigger_luminosities;
 
    ifstream infile; // will be used twice consecutively
 
    auto checkFile = [&infile]() {
        if (!infile.fail()) return;
        cerr << "Failure at opening the file\n";
        exit(EXIT_FAILURE);
    };
 
    // 1) get effective luminosity
    assert(fs::exists(lumi_file));
    infile.open(lumi_file);
    cout << "Opening lumi file: " << lumi_file << endl;
    checkFile();
    while (infile.good()) { 
        int trigger;
        float lumi;
        infile >> trigger >> lumi;
        if (infile.eof()) break;
        TriggerLumi trigger_lumi(trigger);
        trigger_lumi.weight = 1./lumi;
        cout << trigger << ' ' << lumi << ' ' << trigger_lumi.weight << '\n';
        trigger_luminosities.insert({trigger, trigger_lumi});
    } 
    infile.close(); 
 
    // 2) get turn-on points
    cout << "Opening turn-on file: " << turnon_file << endl;
    assert(fs::exists(turnon_file));
    infile.open(turnon_file);
    checkFile();
    while (infile.good()) { 
        int trigger;
        double turnon;
        infile >> trigger >> turnon;
        if (infile.eof()) break;
        cout << trigger << ' ' << turnon << '\n';
        TriggerLumi& trigger_lumi = trigger_luminosities.at(trigger);
        trigger_lumi.turnon = ceil(turnon); //goes to int 
    } 
    infile.close(); 
 
    // then remove pairs with dummy entry
    for (auto trigger_lumi = trigger_luminosities.begin();
              trigger_lumi != trigger_luminosities.end(); /* nothing */) {
 
        // remove or increment
        if (   trigger_lumi->second.turnon == DUMMY
            || trigger_lumi->second.weight == DUMMY) 
            trigger_lumi = trigger_luminosities.erase(trigger_lumi);
        else ++trigger_lumi;
 
    }
 
    // display
    for (auto trigger_lumi: trigger_luminosities) 
        cout << trigger_lumi.first << '\t'
             << trigger_lumi.second.turnon << '\t'
             << trigger_lumi.second.weight << '\n';
 
    cout << flush;
 
    return trigger_luminosities;
}

GetLumiFromFiles() [2/2]

std::map<int, TriggerLumi> DAS::GetLumiFromFiles	(	const std::filesystem::path &	lumi_file,
		const std::filesystem::path &	turnon_file
	)

path to text file with turn-on points

Parameters

lumi_file

path to text file with effective luminosities

getMETfraction()

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

Get MET fraction for each jet trigger separately.

Parameters

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

                                     {1,0} 
               )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice);
    auto tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
 
    auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto met = flow.GetBranchReadOnly<MET>("met");
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
 
    vector<double> MET_edges;
    for (float edge = 0; edge <= 1; edge += 0.01)
        MET_edges.push_back(edge);
    int nMETbins = MET_edges.size()-1;
 
    auto turnon_file = config.get<fs::path>("corrections.normalisation.turnons");
    vector<double> trig_edges = GetTriggerTurnons(turnon_file);
    int nTrigBins = trig_edges.size() - 1;
 
    auto METfraction = new TH2F("METfraction", ";MET fraction;p^{leading}_{T}   (GeV)",
                                    nMETbins, MET_edges.data(),
                                    nTrigBins, trig_edges.data());
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (recJets->empty()) continue;
        static auto rapAcc = y_edges.back();
        if (recJets->front().AbsRap() >= rapAcc) continue;
        if (std::abs(met->SumEt) < feps) continue;
 
        auto pt = recJets->front().CorrPt();
        auto frac = met->Et/met->SumEt;
        double w = (isMC ? gEv->weights.front().v : 1) * rEv->weights.front().v;
 
        cout << setw(15) << pt << setw(15) << frac << setw(15) << w << endl;
 
        // we do not use the jet weight,
        // because the pt is only determined
        // to find the right trigger window
        // -> only the event weight matters
        METfraction->Fill(frac, pt, w);
    }
 
    METfraction->Write();
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

GetObjects()

std::vector<T*> DAS::GetObjects

(

TDirectory *

dir

)

Find all objects of type T directly in a TDirectory (i.e. not recursive)

{
    if (!dir)
        BOOST_THROW_EXCEPTION( std::runtime_error("Trying to access inexistent `TDirectory`") );
    std::vector<T*> objects;
    for (const auto&& element: *(dir->GetListOfKeys())) {
        auto key = dynamic_cast<TKey*>(element);
        if (auto obj = dynamic_cast<T*>(key->ReadObj()); obj)
            objects.push_back(obj);
    }
    return objects;
}

GetSmoothFit() [1/2]

TF1* DAS::GetSmoothFit	(	TH1 *	h,
		double	m,
		double	M,
		int	nSigmaStop = 0,
		const char *	options = "Q0SNRE",
		bool	autoStop = false
	)

Parameters

h	histogram to fit (NB: assumed to be normalised to bin width)
m	min
M	max
nSigmaStop	stop automatically before maxdegree if good chi2/ndf +/- nSigma ~ 1
options	fit options

{
    assert(m < M);
 
    int ifm = h->FindBin(m+feps),
        ifM = h->FindBin(M-feps);
    cout << "Axis interval (" << m << ',' << M << ") to bin interval (" << ifm << ',' << ifM << ')' << endl;
 
    assert(ifm < ifM);
 
    return GetSmoothFit<d>(h, ifm, ifM, nSigmaStop, options, autoStop);
}

GetSmoothFit() [2/2]

TF1* DAS::GetSmoothFit	(	TH1 *	h,
		int	ifm,
		int	ifM,
		int	nSigmaStop = 0,
		const char *	options = "Q0SNRE",
		bool	autoStop = false
	)

Parameters

h	histogram to fit (NB: assumed to be normalised to bin width)
ifm	min index
ifM	max index
nSigmaStop	stop automatically before maxdegree if good chi2/ndf +/- nSigma ~ 1
options	fit options

{
    assert(ifm > 0);
    assert(ifM <= h->GetNbinsX());
    assert(h->GetDimension() == 1);
 
    double contentfm = h->GetBinContent(ifm),
           contentfM = h->GetBinContent(ifM);
 
    double m = h->GetBinLowEdge(ifm), M = h->GetBinLowEdge(ifM+1);
 
    cout << "From bin " << ifm << " (m = " << m << ", content = " << contentfm << ")"
         << " to bin "  << ifM << " (M = " << M << ", content = " << contentfM << ")" << endl;
 
    if (abs(contentfm) <= deps || abs(contentfM) <= deps) {
        cerr << "Fit cannot start with empty bins\n";
        //h->Print("all");
        return nullptr;
    }
 
    double errorfm = h->GetBinError(ifm),
           errorfM = h->GetBinError(ifM);
 
    if (abs(errorfm) <= deps || abs(errorfM) <= deps) {
        cerr << "There must be uncertainties for the fit\n";
        //h->Print("all");
        return nullptr;
    }
 
    double fm = log(contentfm),
           fM = log(contentfM);
 
    double x2n = dmax, error = dmax;
 
    static int iname = 0; ++iname;
    //Greta greta(maxdegree, m, M);
    Thunberg<d, log, exp> greta(m, M);
    const size_t maxdegree = greta.maxdegree;
    TF1 * f = new TF1(Form("greta%d", iname), greta, m, M, greta.npars);
 
    // setting initial parameters
    f->SetParameter(0,(fM+fm)/2);
    f->SetParameter(1,(fM-fm)/2);
    for (size_t degree = 2; degree <= maxdegree; ++degree)
        f->FixParameter(degree,0);
 
    auto dump = [f](int degree) {
        cout << "parameters:\t";
        for (int i = 0; i <= degree; ++i) cout << f->GetParameter(i) << ' ';
        cout << endl;
    };
 
    //running
    for (size_t degree = 2; degree <= maxdegree; ++degree) {
 
        dump(degree-1);
        f->ReleaseParameter(degree);
 
        TFitResultPtr r = h->Fit(f, options);
 
        double chi2 = r->Chi2();
        if (!isnormal(chi2)) {
            h->Print("all");
            cerr << "Fit aborted (chi2 = " << chi2 << ")\n";
            continue;
        }
        double ndf = r->Ndf();
 
        int status = r->Status();
        if (status != 0) {
            h->Print("all");
            cerr << "Fit aborted (status = " << status << ")\n";
            continue;
        }
 
        auto lastx2n = x2n;
        x2n = chi2/ndf;
        error = sqrt(2./ndf);
        cout << "\e[1m" << degree << ' ' << status << '\t' << chi2 << ' ' << ndf << '\t' << x2n << "\u00B1" << error << "\e[0m" << endl;
 
        if (ndf == 2) break;
        if (autoStop && x2n > lastx2n) {
            cerr << "Chi2/ndf has increased -- stepping back and stopping!\n";
            f->FixParameter(degree,0);
            --degree;
            h->Fit(f, options);
            break;
        }
        if (abs(x2n-1) <= nSigmaStop*error) break;
    }
    dump(maxdegree);
 
    return f;
}

getSmoothFits()

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

Smooth relative uncertainties on the output of getUnfPhysDist.

Parameters

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

                                  {1,0} 
            )
{
    cout << __func__ << " start" << endl;
 
    auto fIn = make_unique<TFile>(input.c_str(), "READ");
 
    auto fOut = DT::GetOutputFile(output);
 
    const TString nDist = config.get<string>("uncertainties.smoothness.distribution"),
                  nCov  = config.get<string>("uncertainties.smoothness.covariance");
    cout << nDist << ' ' << nCov << endl;
 
    Step::verbose = steering & DT::verbose;
    for (const auto& var: GetObjects<TDirectory>(fIn.get())) {
        // parallelisation (option `-j`)
        {
            static int ivar = 0; // used for parallelisation
            ++ivar;
            if (slice.second != ivar % slice.first) continue;
        }
 
        TString varname = var->GetName();
        if (varname.Contains("fake")) continue; // TODO: calculate cov matrix
        if (varname.Contains("miss")) continue; // TODO: calculate cov matrix
 
        // looping over all available observables
        for (const auto& obs: GetObjects<TDirectory>(var->GetDirectory(nDist))) {
            TString obsname = obs->GetName();
            cout << varname << ' ' << obsname << endl;
 
            [[ maybe_unused ]] // trick to only activate `cout` in the loop if option `-v` has been given
            static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
            var->cd();
            TDirectory * dInDist = GetDirectory(var, {nDist, obsname}),
                       * dInCov  = GetDirectory(var, {nCov, obsname, obsname}),
                       * dOutDist = fOut->mkdir(var->GetName(), "", true)
                                        ->mkdir(nDist, "", true)
                                        ->mkdir(obsname);
            for (const auto& h: GetObjects<TH1>(dInDist)) {
                if (h->Integral() <= 0) {
                    delete h;
                    continue;
                }
 
                auto [im, iM] = GetBinRange(h);
                if (iM - im <= 2) {
                    delete h;
                    continue;
                }
 
                TString n = h->GetName(),
                        covName = n + "__" + n;
 
                auto cov = dInCov->Get<TH2>(covName);
                if (!cov) {
                    h->Print("all");
                    BOOST_THROW_EXCEPTION( DE::BadInput("Covariance matrix "s + covName
                                    + " could not be found for "s + varname, *dInCov) );
                }
 
                ScaleWidth(h, cov);
 
                using namespace Step;
                auto f = GetSmoothFit<Basis::Chebyshev,log,exp>(h, cov, im, iM, 10,
                                                EarlyStopping::Chi2ndf, 2, 0, cout);
                // NOTE: degree 10 and 2 sigmas are inspired from personal experience with inclusive jet
                //       --> we may want to make options from the command line in the future
                dOutDist->cd();
                for (int i = im; i <= iM; ++i) {
                    double err2 = cov->GetBinContent(i,i);
                    h->SetBinError(i, sqrt(err2));
                }
                h->SetDirectory(dOutDist);
                h->Write();
                f->Write(Form("fit_%s", h->GetName()));
 
                delete h;
                delete cov;
                delete f;
            }
        }
    }
 
    cout << __func__ << " stop" << endl;
}

GetTriggerTurnons()

vector<double> DAS::GetTriggerTurnons

(

const fs::path &

f

)

Load the two-column file with the HLT and corresponding offline PF thresholds.

{
    if (!fs::exists(f))
        BOOST_THROW_EXCEPTION(fs::filesystem_error("Bad input", f,
                    make_error_code(errc::no_such_file_or_directory)));
 
    pt::ptree thresholds;
    read_info(f.string(), thresholds);
 
    vector<double> edges(1, 30);
    for (const auto& threshold: thresholds)
        edges.push_back(threshold.second.get_value<double>());
    edges.push_back(6500);
 
    return edges;
}

identity() [1/2]

Eigen::VectorXd DAS::identity

(

const Eigen::VectorXd &

x

)

identity (Eigen)

Trivial function returning argument, used as default function in the constructor of Teddy

identity() [2/2]

double DAS::identity ( double  x )

inline

{ return x; };

MakeTitle()

std::vector<TString> DAS::MakeTitle ( const std::vector< double > &  edges, const char *  v, bool  lowEdge, bool  highEdge, std::function< const char *(double)>  format  )

inline

Make a vector of properly formated titles from bin edges.

{
    assert(edges.size() > 1);
    std::vector<TString> titles;
    for (size_t i = 1; i < edges.size(); ++i) {
        TString title;
        if (i >              1 ||  lowEdge) title += Form("%s < ", format(edges.at(i-1)));
        title += v;
        if (i < edges.size()-1 || highEdge) title += Form(" < %s", format(edges.at(i  )));
        titles.push_back(title);
    }
    return titles;
}

mergeNtuples()

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

Get fractioned n-tuples after the n-tuplisation and merge them into a single file. Complementary information about the pile-up may be included. A minimal selection is applied, e.g. on the primary vertex (PV).

Parameters

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

                                       {1,0} 
                )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering);
    auto tIn = flow.GetInputTree(inputs, slice, "ntupliser/events");
    auto tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    bool isMC = metainfo.Get<bool>("flags", "isMC");
    int year = metainfo.Get<int>("flags", "year");
 
    // event information
    auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
    auto vtx = flow.GetBranchReadWrite<PrimaryVertex>("primaryvertex");
 
    // physics object information
    auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets", DT::facultative);
    auto recMuons = flow.GetBranchReadWrite<vector<RecMuon>>("recMuons", DT::facultative);
    auto recPhotons = flow.GetBranchReadWrite<vector<RecPhoton>>("recPhotons", DT::facultative);
 
    // jet veto business
    const auto& jetveto_file = config.get<fs::path>("corrections.jetvetomap.filename");
    optional<JetVeto::Conservative> jetveto;
    if (jetveto_file != "/dev/null")
        jetveto = make_optional<JetVeto::Conservative>(jetveto_file);
    else
        cout << orange << "No jet veto file was provided." << def << endl;
    metainfo.Set<fs::path>("corrections", "jetvetomap", "filename", jetveto_file);
 
    // MET filter business
    auto met = flow.GetBranchReadOnly<MET>("met");
    optional<MissingET::Filters> metfilters;
    if (config.get<bool>("corrections.METfilters"))
        metfilters = make_optional<MissingET::Filters>(year);
    else
        cout << orange << "No MET filters will be applied." << def << endl;
    metainfo.Set<bool>("corrections", "METfilters", (bool) metfilters);
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        const bool badPV = abs(vtx->z) > 24 /* cm */ || abs(vtx->Rho) > 2 /* cm */ || vtx->fake;
        if (badPV) {
            recEvt->weights *= 0;
            if (recJets) recJets->clear();
            if (recMuons) recMuons->clear();
            if (recPhotons) recPhotons->clear();
        }
 
        if (recJets && jetveto)
            for (auto& recJet: *recJets)
                (*jetveto)(recJet);
 
        if (metfilters) (*metfilters)(met, recEvt);
 
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

operator*() [1/7]

double DAS::operator* ( const double  v, const Weight &  w  )

inline

{ return w.v * v; }

operator*() [2/7]

double DAS::operator* ( const float  v, const Weight &  w  )

inline

{ return w.v * v; }

operator*() [3/7]

double DAS::operator* ( const int  v, const Weight &  w  )

inline

{ return w.v * v; }

operator*() [4/7]

double DAS::operator* ( const Weight &  w, const double  v  )

inline

{ return w.v * v; }

operator*() [5/7]

double DAS::operator* ( const Weight &  w, const float  v  )

inline

{ return w.v * v; }

operator*() [6/7]

double DAS::operator* ( const Weight &  w, const int  v  )

inline

{ return w.v * v; }

operator*() [7/7]

double DAS::operator* ( const Weight &  w1, const Weight &  w2  )

inline

{ return w1.v * w2.v; }

operator*=() [1/6]

Weight& DAS::operator*= ( Weight &  w, const double  v  )

inline

{ w.v *= v; return w; }

operator*=() [2/6]

Weight& DAS::operator*= ( Weight &  w, const float  v  )

inline

{ w.v *= v; return w; }

operator*=() [3/6]

Weight& DAS::operator*= ( Weight &  w, const int  v  )

inline

{ w.v *= v; return w; }

operator*=() [4/6]

Weights& DAS::operator*= ( Weights &  wgts, const double  v  )

inline

{ for (auto& w: wgts) w *= v; return wgts; }

operator*=() [5/6]

Weights& DAS::operator*= ( Weights &  wgts, const float  v  )

inline

{ for (auto& w: wgts) w *= v; return wgts; }

operator*=() [6/6]

Weights& DAS::operator*= ( Weights &  wgts, const int  v  )

inline

{ for (auto& w: wgts) w *= v; return wgts; }

operator/=() [1/6]

Weight& DAS::operator/= ( Weight &  w, const double  v  )

inline

{ w.v /= v; return w; }

operator/=() [2/6]

Weight& DAS::operator/= ( Weight &  w, const float  v  )

inline

{ w.v /= v; return w; }

operator/=() [3/6]

Weight& DAS::operator/= ( Weight &  w, const int  v  )

inline

{ w.v /= v; return w; }

operator/=() [4/6]

Weights& DAS::operator/= ( Weights &  wgts, const double  v  )

inline

{ for (auto& w: wgts) w /= v; return wgts; }

operator/=() [5/6]

Weights& DAS::operator/= ( Weights &  wgts, const float  v  )

inline

{ for (auto& w: wgts) w /= v; return wgts; }

operator/=() [6/6]

Weights& DAS::operator/= ( Weights &  wgts, const int  v  )

inline

{ for (auto& w: wgts) w /= v; return wgts; }

operator<()

bool DAS::operator< ( const PhysicsObject &  l, const PhysicsObject &  r  )

inline

{
    return l.CorrPt() < r.CorrPt();
}

operator==() [1/5]

bool DAS::operator== ( const PhysicsObject &  l, const PhysicsObject &  r  )

inline

{
    return l.p4 == r.p4 && l.scales == r.scales && l.weights == r.weights;
}

operator==() [2/5]

bool DAS::operator== ( const Weight &  l, const Weight &  r  )

inline

{ return l.v == r.v && l.i == r.i; }

operator==() [3/5]

bool DAS::operator== ( const Weight &  w, const double  v  )

inline

{ return w.v == v; }

operator==() [4/5]

bool DAS::operator== ( const Weight &  w, const float  v  )

inline

{ return w.v == v; }

operator==() [5/5]

bool DAS::operator== ( const Weight &  w, const int  v  )

inline

{ return w.v == v; }

operator>()

bool DAS::operator> ( const PhysicsObject &  l, const PhysicsObject &  r  )

inline

{
    return l.CorrPt() > r.CorrPt();
}

Options()

DT::Options Options	(	const char *	description,
		int	options = Darwin::Tools::none
	)

Constructs Darwin options with the correct commit information.

{
#ifdef DAS_IN_CMSSW
    // Getting the commit hash is not supported by the scram build system. Pass
    // a null hash, which means using the Darwin commit hash. This is the best
    // we can do.
    return DT::Options(description, options, nullptr);
#else
    return DT::Options(description, options, DAS_COMMIT_HASH);
#endif
}

pt_sort()

bool DAS::pt_sort ( const PhysicsObject &  j1, const PhysicsObject &  j2  )

inline

{
    return j1.CorrPt() > j2.CorrPt();
}

ScaleWidth()

void DAS::ScaleWidth	(	TH1 *	h,
		TH2 *	cov
	)

Divide the histogram by bin width and adapt the covariance matrix accordingly.

{
    for (int i = 1; i <= h->GetNbinsX(); ++i) {
        double width = h->GetBinWidth(i);
        double content = h->GetBinContent(i);
        h->SetBinContent(i, content/width);
        for (int j = 1; j <= h->GetNbinsX(); ++j) {
            double width2 = h->GetBinWidth(j);
            double content = cov->GetBinContent(i,j);
            cov->SetBinContent(i, j, content/width/width2);
        }
    }
}

SimpleExec()

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

Template for function (TODO)

Parameters

input	input ROOT files (histograms)
output	output ROOT file (histograms)
config	config handled with `Darwin::Tools::options`
steering	parameters obtained from explicit options

{
    cout << __func__ << " start" << endl;
 
    // TODO
 
    cout << __func__ << " stop" << endl;
}

weightsRef()

Weights& DAS::weightsRef

(

Object &

obj

)

Customization point for how weights are retrieved from an object.

The default implementation tries to use the weights field.

{
    return obj.weights;
}

Variable Documentation

maxMult

const int maxMult = n_edges.size()-1

static

Mjj_edges

const std::vector<double> Mjj_edges {160, 200, 249, 306, 372, 449, 539, 641, 756, 887, 1029, 1187, 1361, 1556, 1769, 2008, 2273, 2572, 2915, 3306, 3754, 4244, 4805, 5374, 6094, 6908, 7861, 8929, 10050 }

static

n_edges

const std::vector<double> n_edges {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}

static

nMjjBins

const int nMjjBins = Mjj_edges.size()-1

static

nPtBins

const int nPtBins = pt_edges.size()-1

static

nPtDijetBins

const int nPtDijetBins = pt_av_edges.size()-1

static

nPUbins

const int nPUbins = 100

static

nYbins

const int nYbins = y_edges.size()-1

static

pt_av_edges

const std::vector<double> pt_av_edges {147, 175, 207, 243, 284, 329, 380, 437, 499, 569, 646, 732, 827, 931, 1046, 1171, 1307, 1458, 1621, 1806, 2003, 2217, 2453, 2702}

static

pt_edges

const std::vector<double> pt_edges {9, 12, 15,18,21,24,28,32,37,43,49,56,64,74,84,97,114,133,153,174,196,220,245,272,300,330,362,395,430,468,507,548,592,638,686,737,790,846,905,967,1032,1101,1172,1248,1327,1410,1497,1588,1684,1784,1890,2000,2116,2238,2366,2500,2640,2787,2941,3103,3273,3450,3637,3832,4037}

static

pthat_edges

const std::vector<double> pthat_edges {30,34,39,44,50,56,63,71,80,89,98,108,120,131,143,156,170,196,226,260,300,336,375,420,470,500,531,564,600,645,693,744,800,846,894,946,1000,1088,1183,1287,1400,1491,1587,1690,1800,1934,2078,2233,2400,2579,2771,2978,3200,3578,4000,4472,5000}

static

ptmax

const double ptmax = pt_edges.back()

static

ptmin

const double ptmin = pt_edges.front()

static

SysDown

const std::string SysDown = "Down"

inline

Suffix used for "down" uncertainties. Follows the Combine convention.

SysUp

const std::string SysUp = "Up"

inline

Suffix used for "up" uncertainties. Follows the Combine convention.

y_edges

const std::vector<double> y_edges {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0}

static

yBins

const std::vector<TString> yBins = MakeTitle(y_edges, "|y|", false, true, [](double v) { return Form("%.1f", v);} )

static

ymax

const double ymax = y_edges.back()

static

ymin

const double ymin = y_edges.front()

static