DAS::Muon Namespace Reference

Classes
class	IDApplier
class	TriggerApplier

Enumerations
enum	Selector {   None = 0u, CutBasedIdLoose = 1u<< 0, CutBasedIdMedium = 1u<< 1, CutBasedIdMediumPrompt = 1u<< 2,   CutBasedIdTight = 1u<< 3, CutBasedIdGlobalHighPt = 1u<< 4, CutBasedIdTrkHighPt = 1u<< 5, PFIsoVeryLoose = 1u<< 6,   PFIsoLoose = 1u<< 7, PFIsoMedium = 1u<< 8, PFIsoTight = 1u<< 9, PFIsoVeryTight = 1u<<10,   TkIsoLoose = 1u<<11, TkIsoTight = 1u<<12, SoftCutBasedId = 1u<<13, SoftMvaId = 1u<<14,   MvaLoose = 1u<<15, MvaMedium = 1u<<16, MvaTight = 1u<<17, MiniIsoLoose = 1u<<18,   MiniIsoMedium = 1u<<19, MiniIsoTight = 1u<<20, MiniIsoVeryTight = 1u<<21 }

Functions
bool	PassTriggerSelection (const Trigger &trg, const vector< RecMuon > &recMuons, int year)
void	applyDimuonTriggerStrategy (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	applyMuonSelection (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	applyRochesterCorr (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
vector< float >	GetMassEdges (int nbins=100, double M=3000, double m=30)
void	getDimuonSpectrum (const vector< fs::path > &inputs, const fs::path &output, const int steering, const DT::Slice slice={1, 0})
void	getZmmgControlPlots (const vector< fs::path > &inputs, const fs::path &output, const int steering, const DT::Slice slice={1, 0})
std::tuple< std::filesystem::path, std::string >	getLocation (const std::string &location)

Enumeration Type Documentation

Selector

enum Selector

Adapted from CMSSW

Enumerator
None
CutBasedIdLoose
CutBasedIdMedium
CutBasedIdMediumPrompt
CutBasedIdTight
CutBasedIdGlobalHighPt
CutBasedIdTrkHighPt
PFIsoVeryLoose
PFIsoLoose
PFIsoMedium
PFIsoTight
PFIsoVeryTight
TkIsoLoose
TkIsoTight
SoftCutBasedId
SoftMvaId
MvaLoose
MvaMedium
MvaTight
MiniIsoLoose
MiniIsoMedium
MiniIsoTight
MiniIsoVeryTight

              {
    None                   = 0u, // added
    CutBasedIdLoose        = 1u<< 0,
    CutBasedIdMedium       = 1u<< 1,
    CutBasedIdMediumPrompt = 1u<< 2,  // medium with IP cuts
    CutBasedIdTight        = 1u<< 3,
    CutBasedIdGlobalHighPt = 1u<< 4,  // high pt muon for Z',W' (better momentum resolution)
    CutBasedIdTrkHighPt    = 1u<< 5,  // high pt muon for boosted Z (better efficiency)
    PFIsoVeryLoose         = 1u<< 6,  // reliso<0.40
    PFIsoLoose             = 1u<< 7,  // reliso<0.25
    PFIsoMedium            = 1u<< 8,  // reliso<0.20
    PFIsoTight             = 1u<< 9,  // reliso<0.15
    PFIsoVeryTight         = 1u<<10,  // reliso<0.10
    TkIsoLoose             = 1u<<11,  // reliso<0.10
    TkIsoTight             = 1u<<12,  // reliso<0.05
    SoftCutBasedId         = 1u<<13,
    SoftMvaId              = 1u<<14,
    MvaLoose               = 1u<<15,
    MvaMedium              = 1u<<16,
    MvaTight               = 1u<<17,
    MiniIsoLoose           = 1u<<18,  // reliso<0.40
    MiniIsoMedium          = 1u<<19,  // reliso<0.20
    MiniIsoTight           = 1u<<20,  // reliso<0.10
    MiniIsoVeryTight       = 1u<<21   // reliso<0.05
};

Function Documentation

applyDimuonTriggerStrategy()

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

Code to apply the Double Muon trigger SFs documented at the following link: https://indico.cern.ch/event/1106050/contributions/4653418/ Expecting one the following triggers:

• HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass8_v <– expected from n-tupliser
• HLT_Mu19_TrkIsoVVL_Mu9_TrkIsoVVL_DZ_Mass8_v
• HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass3p8_v <– expected from n-tupliser
• HLT_Mu19_TrkIsoVVL_Mu9_TrkIsoVVL_DZ_Mass3p8_v

The code was written assuming the calibration files found in /eos/cms/store/group/phys_muon/jmijusko/DoubleMuonTrigger_SF_UL/ which contains the following TH2Ds with axes corresponding to the respective muon pseudorapidities and bin content to the efficiency.

• ScaleFactor[WP]_UL[year]
• ScaleFactor[WP]_UL[year]_stat
• ScaleFactor[WP]_UL[year]_syst
• ScaleFactor[WP]_UL[year]_alt (except for 2018)
• ScaleFactor[WP]_UL[year]_massBin (except for 2018)
• ScaleFactor[WP]_UL[year]_massRange (except for 2018) where [WP] = Tight or Medium (for both ID and ISO) [year] = 2016_HIPM, 2016, 2017, or 2018

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 tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    auto year = metainfo.Get<int>("flags", "year");
 
    auto efficiency = config.get<string>("corrections.muons.trigger_eff");
    auto [filename, histname] = getLocation(efficiency);
 
    TriggerApplier applier(filename, histname, year,
                           isMC && filename != "/dev/null", steering & DT::syst);
    metainfo.Set<string>("corrections", "muons", "trigger_eff", efficiency);
    for (const auto& name: applier.weightNames())
        metainfo.Set<string>("variations", RecEvent::WeightVar, name);
 
    auto trg = flow.GetBranchReadOnly<Trigger>("muonTrigger");
    auto recEvent = flow.GetBranchReadWrite<RecEvent>("recEvent");
    auto recMuons = flow.GetBranchReadOnly<vector<RecMuon>>("recMuons");
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        applier(*recEvent, *trg, *recMuons);
 
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

applyMuonSelection()

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

Generic executable for all types of muon efficiency corrections. The command should be called once for each table.

TWiki Tables

==== RECO NUM = TrackerMuons DEN = genTracks

• binning: abseta_pt
• uncertainties: – (?)

==== Muon ID NUM = HighPtID LooseID MediumID MediumPromptID SoftID TightID TrkHighPtID DEN = TrackerMuons

• binning: abseta_pt
• uncertainties: stat, syst

==== Muon ISO DEN \ NUM LooseRelIso LooseRelTkIso TightRelIso TightRelTkIso LooseID x MediumID x x MediumPromptID x x TightID & IPCut x x HighPtID & IPCut x x TrkHighPtID & IPCut x x

• binning: abseta_pt
• uncertainties: stat, syst

==== trigger DEN \ NUM IsoMu24 IsoMu24_or_Mu50 Mu50_or_OldMu100_or_TkMu100 IdMedium & PFIsoMedium x x IdTight & PFIsoTight x x IdGlobalHighPt & TkIsoLoose x

• binning: abseta_pt, charge_abseta_pt, charge_eta_pt, eta_pt
• uncertainties: stat, syst

==== tracking Quoting the TWiki:

Official tracking POG SFs, where the probe is a standalone muon and it is required to be matched to a generalTrack in its vicinity, are independent of momentum. The SFs are found to be above 0.99, and are considered equivalent to unity. You can find the detailed plots per eta and pT bins in this TWiki. The recommendation are to not apply them, unless you target a precison analysis.

Todo:

ControlPlots

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 tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
 
    auto recMuons = flow.GetBranchReadWrite<vector<RecMuon>>("recMuons");
 
    auto efficiency = config.get<string>("corrections.muons.efficiency");
    if (efficiency == "") {
        auto& efficiencies = config.get_child("corrections.muons.efficiency");
        if (efficiencies.size() == 0)
            BOOST_THROW_EXCEPTION( invalid_argument("Missing efficiency") );
        auto it = efficiencies.begin();
        while (metainfo.Find("corrections", "muons", "efficiency", it->first))
            ++it;
        efficiency = it->first;
    }
    auto [filename, histname] = getLocation(efficiency);
    cout << filename << ' ' << histname << endl;
 
    IDApplier applier(filename, histname, filename != "/dev/null", steering & DT::syst);
    if (metainfo.Find("corrections", "muons", "efficiency")) {
        // 0 = previous, former
        auto efficiency0 = metainfo.Get<string>("corrections", "muons", "efficiency");
        auto [filename0, histname0] = getLocation(efficiency0);
        if (filename == filename0 && filename != "/dev/null")
            cerr << orange << "The same file has been used twice in a row, "
                              "which is unexpected: `" << filename << "`. "
                              "Proceed at your own risks.\n" << def;
        IDApplier applier0(filename0, histname0, false, false);
        if (applier.DEN != applier0.NUM)
            cerr << orange << "The former and present muon efficiency corrections are "
                              "unexpected combinations: `" << applier0.name << "` and `"
                           << applier.name <<"`. Proceed at your own risks.\n" << def;
    }
    metainfo.Set<string>("corrections", "muons", "efficiency", efficiency);
 
    for (const auto& name: applier.weightNames())
        metainfo.Set<string>("variations", RecMuon::WeightVar, name);
 
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        applier(*recMuons);
 
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

applyRochesterCorr()

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

TWiki

 set   nmembers  comment

Default 0 1 default, reference based on madgraph sample, with adhoc ewk (sw2eff and Z width) and Z pt (to match data) weights. Stat 1 100 pre-generated stat. replicas; Zpt 2 1 derived without reweighting reference pt to data. Ewk 3 1 derived without applying ad-hoc ewk weights deltaM 4 1 one representative set for alternative profile deltaM mass window

Ewk2 5 1 reweight the reference from constant to s-dependent Z width

For statistical replicas, std. dev. gives uncertainty. For the rest, difference wrt the cental is assigned as syst. error. Total uncertainty is calculated as a quadrature sum of all components.

Todo:

Implement alternative nominal correction

Todo:

Zpt: rather provide as alternative? (issue #76)

Todo:

Ewk2: understand 3rd remark from MUON POG

Todo:

clear scales?

Todo:

should we force up (down) in the first (second) place?

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 tOut = flow.GetOutputTree(output);
 
    DT::MetaInfo metainfo(tOut);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
 
    auto recMuons = flow.GetBranchReadWrite<vector<RecMuon>>("recMuons");
    auto genMuons = isMC ? flow.GetBranchReadWrite<vector<GenMuon>>("genMuons") : nullptr;
 
    auto tables = config.get<fs::path>("corrections.muons.rochester.tables");
    if (!fs::exists(tables))
        BOOST_THROW_EXCEPTION( fs::filesystem_error("Tables can't be found.",
                        tables, make_error_code(errc::no_such_file_or_directory)));
    metainfo.Set<fs::path>("corrections", "muons", "rochester", "tables", tables);
    RoccoR rc(tables.c_str());
 
    auto alternative = config.get<bool>("corrections.muons.rochester.alternative");
    if (alternative)
        BOOST_THROW_EXCEPTION( logic_error("Alternative nominal correction not yet implemented") );
    metainfo.Set<bool>("corrections", "muons", "rochester", "alternative", alternative);
 
    const bool applySyst = steering & DT::syst;
    if (applySyst) {
        for (const string& syst: {"stat", "Zpt", "Ewk", "deltaM", "Ewk2"}) {
            metainfo.Set<string>("variations", RecMuon::WeightVar, syst + SysUp);
            metainfo.Set<string>("variations", RecMuon::WeightVar, syst + SysDown);
        }
    }
 
    vector<GenMuon> mcands; // matching candidates
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (isMC) mcands = *genMuons; // copy to avoid modifying the event
 
        for (auto& recMuon: *recMuons) {
 
            float pt = recMuon.p4.Pt(),
                  eta = recMuon.p4.Eta(),
                  phi = recMuon.p4.Phi();
            int Q = recMuon.Q;
            cout << Q << '\t' << pt << '\t' << eta << '\t' << phi << endl;
 
            auto& nl = recMuon.nTkHits;
 
            // 1) get scale factor (lambda)
            function<double(int /* unc set */,int /* member */)> SF;
            if (isMC) {
                // find the closest gen muon within DR <= 0.1
                auto DR = [&recMuon](const GenMuon& m) { return DeltaR(m.p4, recMuon.p4) < 0.1; };
                auto getClosest = [&DR](const GenMuon& m1, const GenMuon& m2) { return DR(m1) < DR(m2); };
                auto closestGenMuon_it = min_element(genMuons->begin(), genMuons->end(), getClosest);
 
                if (closestGenMuon_it != mcands.end() && DR(*closestGenMuon_it))
                    SF = [&](int s, int m) {
                        float genpt = closestGenMuon_it->p4.Pt();
                        return rc.kSpreadMC(Q, pt, eta, phi, genpt, s, m);
                    };
                else {
                    static TRandom3 r(metainfo.Seed<897532>(slice));
                    const auto u = r.Rndm();
                    SF = [&](int s, int m) { return rc.kSmearMC(Q, pt, eta, phi, nl, u, s, m); };
                }
            }
            else
                SF = [&](int s, int m) { return rc.kScaleDT(Q, pt, eta, phi, s, m); };
 
            // 2) store it in scale
            recMuon.scales.clear(); 
            auto nom = SF(0,0);
            recMuon.scales.push_back(nom);
 
            // 3) variations
 
            if (!applySyst) continue;
 
            static const int nvars = 5;
            recMuon.scales.reserve(1+2*nvars);
            auto push_back = [&recMuon,&nom](float unc) {
                recMuon.scales.push_back(unc);
                unc -= nom; unc *= -1; unc += nom;
                recMuon.scales.push_back(unc);
            };
 
            // 3a) first stat unc (from replicas)
            float sum1 = 0, sum2 = 0;
            for (int m = 0; m < 100; ++m) {
                float sf = SF(1, m);
                sum1 += sf;
                sum2 += sf*sf;
            }
            double unc = sqrt(sum1*sum1 - sum2) / 100;
 
            push_back(unc);
 
            // 3b) then the syst unc
            for (int s = 2; s <= nvars; ++s) {
                unc = SF(s,0);
                push_back(unc);
            }
        }
        if (steering & DT::fill) tOut->Fill();
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

getDimuonSpectrum()

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

Obtain dimuon spectrum.

Parameters

inputs	input ROOT files (n-tuples)
output	output ROOT file (n-tuple)
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);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
 
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto genEvt = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
 
    auto recMuons = flow.GetBranchReadOnly<vector<RecMuon>>("recMuons");
    auto recPhotons = flow.GetBranchReadOnly<vector<RecPhoton>>("recPhotons", DT::facultative);
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets", DT::facultative);
 
    vector<float> mass_edges = GetMassEdges();
    unique_ptr<TH1F> h_dimuon = make_unique<TH1F>("dimuon", "", mass_edges.size() - 1, mass_edges.data());
    unique_ptr<TH1F> h_mmg = make_unique<TH1F>("mumugamma", "", mass_edges.size() - 1, mass_edges.data());
 
    unique_ptr<TH1F> h_dimuon_peak = make_unique<TH1F>("dimuon_peak", "", 60, 71, 121);
    unique_ptr<TH1F> h_mmg_peak = make_unique<TH1F>("mumugamma_peak", "", 60, 71, 121);
 
    unique_ptr<TH2F> h_dimuon_mmg = make_unique<TH2F>("dimuon_mmg", "", mass_edges.size() - 1, mass_edges.data(),
                                                                        mass_edges.size() - 1, mass_edges.data());
    unique_ptr<TH2F> h_dimuon_mmg_peak = make_unique<TH2F>("dimuon_mmg_peak", "", 60, 71, 121,
                                                                                  60, 71, 121);
 
    vector<double> pt_edges = {25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 110., 130., 150., 170., 190., 220., 250., 400., 1000.};
    unique_ptr<TH3D> h_Zjet = make_unique<TH3D>("Zjet", "", (int)pt_edges.size() - 1, pt_edges.data(),
                                                            (int)y_edges.size() - 1, y_edges.data(),
                                                            (int)y_edges.size() - 1, y_edges.data());
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (recMuons->size() < 2) continue;
        if (recMuons->at(0).CorrPt() < 20) continue;
        if (recMuons->at(1).CorrPt() < 15) continue;
 
        FourVector m0 = recMuons->at(0).CorrP4(),
                   m1 = recMuons->at(1).CorrP4();
        FourVector diMuon = m0 + m1;
        float w_Z = recEvt->weights.front()
                * recMuons->at(0).weights.front()
                * recMuons->at(1).weights.front();
        if (isMC)
            w_Z *= genEvt->weights.front();
        h_dimuon->Fill(diMuon.M(), w_Z);
        h_dimuon_peak->Fill(diMuon.M(), w_Z);
 
        if (recJets != nullptr && !recJets->empty() && // at least one jet
                abs(diMuon.M() - 91.2) < 20 && diMuon.Pt() > 25) {
            const auto& jet = recJets->front().CorrP4();
            if (jet.Pt() > 20 && abs(jet.Rapidity()) < 2.4 && // jet selection
                    min(m0.Pt(), m1.Pt()) > 29 && max(abs(m0.Eta()), abs(m1.Eta())) < 2.4 && // muon selection
                    DeltaR(jet, m0) > 0.4 && DeltaR(jet, m1) > 0.4) { // opening angle
                float w_j = recJets->front().weights.front(),
                      y_b = 0.5 * (diMuon.Rapidity() + recJets->front().Rapidity()),
                      y_s = 0.5 * (diMuon.Rapidity() - recJets->front().Rapidity());
                h_Zjet->Fill(diMuon.Pt(), y_b, y_s, w_Z * w_j);
            }
        }
 
        if (diMuon.M() >= 76) continue;
        if (recPhotons != nullptr && !recPhotons->empty()) {
            FourVector ph = recPhotons->front().CorrP4();
            if (ph.Pt() < 20) continue;
            FourVector mumugamma = diMuon + ph;
            float w_ph = recPhotons->front().weights.front();
 
            h_mmg->Fill(mumugamma.M(), w_Z * w_ph);
            h_mmg_peak->Fill(mumugamma.M(), w_Z * w_ph);
 
            h_dimuon_mmg->Fill(diMuon.M(), mumugamma.M(), w_Z * w_ph);
            h_dimuon_mmg_peak->Fill(diMuon.M(), mumugamma.M(), w_Z * w_ph);
        }
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    h_dimuon->SetDirectory(fOut);
    h_dimuon_peak->SetDirectory(fOut);
    h_mmg->SetDirectory(fOut);
    h_mmg_peak->SetDirectory(fOut);
    h_Zjet->SetDirectory(fOut);
 
    h_dimuon->Write();
    h_dimuon_peak->Write();
    h_mmg->Write();
    h_mmg_peak->Write();
    h_Zjet->Write();
 
    h_dimuon_mmg->SetDirectory(fOut);
    h_dimuon_mmg_peak->SetDirectory(fOut);
    h_dimuon_mmg->Write();
    h_dimuon_mmg_peak->Write();
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

getLocation()

std::tuple<std::filesystem::path, std::string> DAS::Muon::getLocation

(

const std::string &

location

)

{
    using namespace std;
    namespace al = boost::algorithm;
    namespace fs = filesystem;
 
    vector<string> words;
    al::split(words, location, al::is_any_of(":"), al::token_compress_on);
    if (words.size() != 2)
        BOOST_THROW_EXCEPTION( invalid_argument("Expecting `/path/to/file.ext:histname`") );
    fs::path filename = words.front();
    string histname = words.back();
    return {filename, histname};
}

GetMassEdges()

vector<float> DAS::Muon::GetMassEdges	(	int	nbins = 100,
		double	M = 3000,
		double	m = 30
	)

                                                                            {
    float R = pow(M / m, 1. / nbins);
    vector<float> edges;
    edges.reserve(nbins + 1);
    for (float edge = m; edge <= M; edge *= R)
        edges.push_back(edge);
    return edges;
}

getZmmgControlPlots()

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

Obtain angular separation plots and fake/miss rates of photons/muons.

Parameters

inputs	input ROOT files (n-tuples)
output	output ROOT file (n-tuple)
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);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
 
    auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto genEvt = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
 
    auto recMuons = flow.GetBranchReadOnly<vector<RecMuon>>("recMuons");
    auto genMuons = isMC ? flow.GetBranchReadOnly<vector<GenMuon>>("genMuons") : nullptr;
    auto recPhotons = flow.GetBranchReadOnly<vector<RecPhoton>>("recPhotons");
    auto genPhotons = isMC ? flow.GetBranchReadOnly<vector<GenPhoton>>("genPhotons") : nullptr;
 
    array h_deltaR {
        make_unique<TH1F>("deltaR0", "", 100, 0, 5.0),
            make_unique<TH1F>("deltaR1", "", 100, 0, 5.0)
    };
    auto h_deltaRmin = make_unique<TH1F>("deltaRmin", "", 100, 0, 5.0);
 
    auto h_fake_mu = make_unique<TH1F>("fake_mu", "", 2, -0.5, 1.5);
    auto h_miss_mu = make_unique<TH1F>("miss_mu", "", 2, -0.5, 1.5);
    auto h_fake_ph = make_unique<TH1F>("fake_ph", "", 2, -0.5, 1.5);
    auto h_miss_ph = make_unique<TH1F>("miss_ph", "", 2, -0.5, 1.5);
    auto h_fake_mmg = make_unique<TH1F>("diff_fake_mmg", "", 50, 0, 100);
    auto h_miss_mmg = make_unique<TH1F>("diff_miss_mmg", "", 100, 0, 200);
 
    array h_gen_reco_mu_dR {
        make_unique<TH1F>("gen_reco_mu0_dR", "", 50, 0, 0.5),
            make_unique<TH1F>("gen_reco_mu1_dR", "", 50, 0, 0.5)
    };
    auto h_gen_reco_ph_dR = make_unique<TH1F>("gen_reco_ph_dR", "", 50, 0, 0.5);
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[maybe_unused]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        // muon selection
        if (recMuons->size() < 2) continue;
        if (recMuons->at(0).CorrPt() < 20) continue;
        if (recMuons->at(1).CorrPt() < 15) continue;
 
        // photon selection
        if (recPhotons->empty()) continue;
    auto rph = recPhotons->front().CorrP4();
        if (rph.Pt() < 20) continue;
 
        // \todo Use Di<RecMuon,RecMuon>
        // \todo Move phase space selection here
        auto recdimuon = recMuons->at(0) + recMuons->at(1);
        if (recdimuon.CorrP4().M() >= 76) continue;
 
        double w_Z = recEvt->weights.front() * recdimuon.Weight();
        if (isMC)
            w_Z *= genEvt->weights.front();
 
        bool match_mmg = true; // will use logical ands to check the matching to individual objects
        if (isMC) {
            const float maxDr = 0.1;
 
            if (genMuons->size() > 1) {
                for (size_t irec = 0; irec < 2; irec++) {
                // the leading and subleading muons may swap due to resolution effects
                float dR0 = DeltaR(recMuons->at(irec).CorrP4(), genMuons->at(0).CorrP4()),
                      dR1 = DeltaR(recMuons->at(irec).CorrP4(), genMuons->at(1).CorrP4());
                size_t igen = dR0 < dR1 ? 0 : 1;
                float dR = min(dR0, dR1);
 
                bool match_mu = dR < maxDr;
                double rW = recMuons->at(irec).weights.front(),
                       gW = genMuons->at(igen).weights.front();
                h_fake_mu->Fill(match_mu, rW);
                h_miss_mu->Fill(match_mu, gW);
                match_mmg &= match_mu;
 
                if (match_mu)
                    h_gen_reco_mu_dR[igen]->Fill(dR, gW*rW);
                }
            }
            else match_mmg = false;
 
            // \todo Is the leading photon the right photon? (or should we look for the photon closest to one of the two leading muons?)
            if (recPhotons->size() > 0) {
                bool match_ph = DeltaR(recPhotons->front().CorrP4(), genPhotons->front().CorrP4()) < maxDr;
                h_fake_ph->Fill(match_ph, recPhotons->front().weights.front());
                h_miss_ph->Fill(match_ph, genPhotons->front().weights.front());
                match_mmg &= match_ph;
            }
            else match_mmg = false;
        }
 
        // \todo Use Di<Di<RecMuon,RecMuon>, RecPhoton>
        double w_ph = recPhotons->front().weights.front();
 
        float dR0 = DeltaR(rph, recMuons->at(0).CorrP4()),
              dR1 = DeltaR(rph, recMuons->at(1).CorrP4());
        h_deltaR[0]->Fill(dR0, w_Z * w_ph);
        h_deltaR[1]->Fill(dR1, w_Z * w_ph);
        h_deltaRmin->Fill(min(dR0, dR1), w_Z * w_ph);
 
        if (isMC) {
 
            if (match_mmg) {
                FourVector gph = genPhotons->front().CorrP4();
                h_gen_reco_ph_dR->Fill(DeltaR(gph, rph), w_Z * w_ph);
            }
            else {
                auto mumugamma = recdimuon.CorrP4() + rph;
                h_fake_mmg->Fill(mumugamma.M(), w_Z * w_ph);
                h_miss_mmg->Fill(mumugamma.M(), w_Z * w_ph);
            }
        }
    }
 
    metainfo.Set<bool>("git", "complete", true);
 
    h_deltaRmin->SetDirectory(fOut);
    h_deltaRmin->Write();
    for (auto& h: h_deltaR) {
        h->SetDirectory(fOut);
        h->Write();
    }
 
    h_fake_mu->SetDirectory(fOut);
    h_miss_mu->SetDirectory(fOut);
    h_fake_ph->SetDirectory(fOut);
    h_miss_ph->SetDirectory(fOut);
    h_fake_mmg->SetDirectory(fOut);
    h_miss_mmg->SetDirectory(fOut);
    h_gen_reco_ph_dR->SetDirectory(fOut);
    h_fake_mu->Write();
    h_miss_mu->Write();
    h_fake_ph->Write();
    h_miss_ph->Write();
    h_fake_mmg->Write();
    h_miss_mmg->Write();
    h_gen_reco_ph_dR->Write();
    for (auto& h: h_gen_reco_mu_dR) {
        h->SetDirectory(fOut);
        h->Write();
    }
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

PassTriggerSelection()

bool DAS::Muon::PassTriggerSelection	(	const Trigger &	trg,
		const vector< RecMuon > &	recMuons,
		int	year
	)

{
    if (recMuons.size() < 2) return false;
 
    if (recMuons.at(0).CorrPt() < 19) return false;
    if (recMuons.at(1).CorrPt() <  9) return false;
 
    auto dimuon = recMuons.at(0).p4 + recMuons.at(1).p4;
    float m = dimuon.M();
    bool pass = false;
    switch (year) {
        case 2017:
            pass |= (m > 3.8 && trg.Bit.at(1));
            [[fallthrough]];
        case 2018:
            pass |= (m >   8 && trg.Bit.at(0));
            break;
        default:
            BOOST_THROW_EXCEPTION( std::runtime_error(to_string(year) + " has not been implemented yet"s) );
    }
 
    return pass;
}