Exercise-4-Muon-Efficiencies.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tag-and-Probe is a data-driven technique used to calculate lepton reconstruction, identification, and trigger efficiencies. This technique uses narrow dilepton resonances, such as Z (for muons with relatively high pT) or J/ψ (for muons with lower pT). Almost-unbiased estimates of the efficiencies can be obtained at the different stages of muon trigger and offline reconstruction. Events are selected with strict requirements on one muon (the tag), and with a more relaxed selection on the other muon (the probe), such that the probe muon can be used to measure the efficiency in question without large biases. The fraction of probe muons that pass the selection under study gives an estimate of its efficiency. The invariant mass of the tag-probe pair is used to select Z→μμ or J/ψ→μμ events.\n",
    "The Tag-and-Probe technique is generally used to measure and compare efficiencies in data and in MC simulation, and thus to compute a correction scale factor that can be applied to MC events to match the efficiency observed in data. These scale factors are typically determined as functions of pT and η. If necessary, their dependence on other kinematic variables can be investigated too — e.g. vs the number of vertices, in case of strong pileup dependence. In some cases, customized scale factors are necessary for some analyses, depending on their specific trigger and offline thresholds.\n",
    "\n",
    "Despite the tight selection on the tag muon and the invariant mass constraints, the selected Z→μμ or J/ψ→μμ sample generally contains background events, which appear as a nonresonant continuum underneath the resonance peak. Therefore the background must be subtracted, to ensure that the efficiency is measured with signal muons only. This can be achieved by fitting the invariant mass spectrum to signal + background shapes (e.g. analytical functions or MC templates). Finding proper functions or templates for signal and background is often the most challenging part of the process.\n",
    "\n",
    "The total lepton efficiency is generally factorized in multiple steps as follows:\n",
    "total lepton efficiency = (tracking) × (reconstruction/tracking) × (ID/reconstruction) × (isolation/ID) × (trigger/isolation)\n",
    "In each efficiency step, the denominator determines the selection of the probe. In the last steps (in part. isolation and trigger), the probe selection is tighter and, therefore, the background level is quite low and the background subtraction is generally easier or not even needed, in some cases. In this exercise, you will measure efficiencies using simulated Z→μμ events:\n",
    "\n",
    "DYJets.root\n",
    "\n",
    "Since you are using a pure Z sample, you won't need background subtraction nor fitting. The efficiency will simply be computed by counting probes before and after the selection under study.\n",
    "\n",
    "You will start by computing \"true\" efficiencies using the generator-level information. This will be your reference efficiencies. Next, you will implement a simple tag-and-probe algorithm to measure the efficiencies with a data-driven approach, and you will compare your results to the \"true\" efficiencies. Finally, you will try to use the same algorithm on a real single-muon data sample:\n",
    "\n",
    "MuonData.root<br>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "import ROOT\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import awkward as ak\n",
    "from coffea.nanoevents import NanoEventsFactory, NanoAODSchema\n",
    "from math import pi, sqrt, cosh, cos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dR=0.1 #set delta R requirement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: \"MC-true\" efficiency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First,  select the prompt final state gen level muons with pT > 20 GeV and |η| < 2.4, and fill histograms with the main kinematic variables. Include a distribution of the number of reconstructed vertices\n",
    "\n",
    "Next, select reco muons from Muon collections and check that they are matched to one of the GenParticle muons used . Then you can start filling histograms (using the generator-level information from the associated GenParticle) at different selection stages:\n",
    "* all reconstructed muons with a valid inner track (i.e. skip muons that are reconstructed exclusively as standalone muons)\n",
    "* reconstructed muons passing some standard ID: e.g. loose, medium, tight\n",
    "* tight muons that pass a requirement on the Δβ-corrected combined relative isolation (e.g. < 0.12)\n",
    "* isolated muons that pass a single-muon trigger or a combination of single-muon triggers, e.g. HLT_IsoMu24, HLT_IsoTkMu24 This step requires matching the muon candidate to a trigger object.\n",
    "\n",
    "Finally, create efficiency graphs by computing the ratio of the histograms above. You can use the TGraphAsymmErrors class, which gives you different options to compute the errors on the efficiency.\n",
    "\n",
    "Compute the following factorized efficiencies:\n",
    "* reconstructed muons / generated muons\n",
    "* muon ID / reconstructed muons\n",
    "* isolation / muon ID\n",
    "* single-muon trigger / isolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events = NanoEventsFactory.from_root(\n",
    "    '/eos/uscms/store/user/cmsdas/2025/short_exercises/muons/DYJets.root',\n",
    "     treepath = \"Events\", \n",
    "    schemaclass=NanoAODSchema,\n",
    "    metadata={\"dataset\": \"DYJets\"},\n",
    "    entry_stop = 100000,\n",
    ").events()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "selelct gen particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events[\"Muon\",\"PV_npvsGood\"] = events.PV.npvsGood\n",
    "events[\"GenPart\",\"PV_npvsGood\"] = events.PV.npvsGood\n",
    "muon = events.Muon    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gen_muon= events.GenPart[((abs(events.GenPart.pdgId)==13) & (events.GenPart.status==1) & (abs(events.GenPart.eta)< 2.4))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trigger_muon= events.TrigObj[abs(events.TrigObj.id)==13]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deltaR(probe, trigger):\n",
    "    return np.sqrt((probe.eta-trigger.eta)**2+((probe.phi-trigger.phi+pi)%(2*pi)-pi)**2)\n",
    "\n",
    "def trigger_match(leptons, trigobjs, pdgid, pt, filterbit):\n",
    "        pass_pt = trigobjs.pt > pt\n",
    "        pass_id = abs(trigobjs.id) == pdgid\n",
    "        pass_filterbit = (trigobjs.filterBits & (0x1 << filterbit)) != 0\n",
    "        trigger_cands = trigobjs[pass_pt & pass_id & pass_filterbit]\n",
    "        delta_r = leptons.metric_table(trigger_cands, metric=deltaR)\n",
    "        pass_delta_r = delta_r < 0.1\n",
    "        n_of_trigger_matches = ak.sum(pass_delta_r, axis=2)\n",
    "        trig_matched_locs = n_of_trigger_matches >= 1\n",
    "\n",
    "        return trig_matched_locs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c={}\n",
    "def TH1_fill_and_plot(th1, array, name): #fill and plot TH1\n",
    "    \n",
    "    for a in ak.flatten(array):\n",
    "        th1.Fill(a)\n",
    "        \n",
    "    c[f'{name}']= ROOT.TCanvas(f'{name}',\"\",900,600)\n",
    "    th1.Draw()\n",
    "    c[f'{name}'].Draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c1={}\n",
    "def eff_fill_and_plot(th1, name): #plot the efficiency graph\n",
    "     \n",
    "    c1[f'{name}']= ROOT.TCanvas(f'{name}',\"\",900,600)\n",
    "    th1.Draw()\n",
    "    c1[f'{name}'].Draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets set a few cut on the event\n",
    "\n",
    "* Event cut: must has at least a good vertex\n",
    "* Particle cut: \n",
    "\n",
    "* 1. correspond genpart must be good muon \n",
    "* 2. fired a trigger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "good_gen_mask=(((muon.isGlobal==1)|(muon.isTracker==1))& (muon.genPartIdx>-1) & (muon.matched_gen.pt>20)& (abs(muon.matched_gen.eta)<2.4)& (abs(muon.matched_gen.pdgId)==13)& (muon.matched_gen.status==1)& ((muon.matched_gen.statusFlags%2)==1))\n",
    "good_gen_mask = ak.fill_none(good_gen_mask, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_mask_goodpv= events.PV.npvsGood>0 #at least a good vertex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_mask_goodpv_HLT_isomu24= (events.PV.npvsGood>0) & (events.HLT.IsoMu24==1) #fired the iso muon trigger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trigger_match_mask=trigger_match(muon, trigger_muon, 13, 0, 8 )\n",
    "print (trigger_match_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "muon_mask={}\n",
    "muon_mask['rec']= good_gen_mask\n",
    "muon_mask['rec_loose']= good_gen_mask& (muon.looseId==1) \n",
    "muon_mask['rec_medium']= good_gen_mask& (muon.mediumId==1) \n",
    "muon_mask['rec_tight']= good_gen_mask& (muon.tightId==1)\n",
    "muon_mask['rec_tight_iso']= good_gen_mask& (muon.tightId==1)& (muon.pfRelIso04_all< 0.12)\n",
    "muon_mask['rec_tight_iso_hlt']= good_gen_mask& (muon.tightId==1)& (muon.pfRelIso04_all< 0.12)& trigger_match_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "th1dic={} #dictionary of TH1D\n",
    "efficiencydic={} #dictionary of TGraphAsymmErrors\n",
    "for name in ['gen', 'rec', 'rec_loose', 'rec_medium', 'rec_tight', 'rec_tight_iso', 'rec_tight_iso_hlt']:\n",
    "    if name== 'gen':           #Cannot plot efficiency for gen itself, so if \n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0) #constrcut\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], gen_muon.pt, f'h_{name}_pt') #fill and plot\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], gen_muon.eta, f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], gen_muon.PV_npvsGood, f'h_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "        continue\n",
    "        \n",
    "    if name== 'rec_tight_iso_hlt':\n",
    "        print (name, name2)\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.matched_gen.pt[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.matched_gen.eta[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_nvtx')\n",
    "\n",
    "        efficiencydic[f'gae_{name}_pt']   = ROOT.TGraphAsymmErrors(20) \n",
    "        efficiencydic[f'gae_{name}_pt'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetXaxis().SetTitle(f'h_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_pt'].Divide(th1dic[f'h_{name}_pt']  , th1dic[f'h_{name2}_pt']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_pt'], f'gae_{name}_pt')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_eta']   = ROOT.TGraphAsymmErrors(24 ) \n",
    "        efficiencydic[f'gae_{name}_eta'].SetName(f'gae_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetXaxis().SetTitle(f'h_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_eta'].Divide(th1dic[f'h_{name}_eta']  , th1dic[f'h_{name2}_eta']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_eta'], f'gae_{name}_eta')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_nvtx']   = ROOT.TGraphAsymmErrors(90 ) \n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetName(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetXaxis().SetTitle(f'h_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_nvtx'].Divide(th1dic[f'h_{name}_nvtx']  , th1dic[f'h_{name2}_nvtx']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_nvtx'], f'gae_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "    else:\n",
    "        print (name, name2)\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.matched_gen.pt[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.matched_gen.eta[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_nvtx')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_pt']   = ROOT.TGraphAsymmErrors(20 ) \n",
    "        efficiencydic[f'gae_{name}_pt'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetXaxis().SetTitle(f'h_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_pt'].Divide(th1dic[f'h_{name}_pt']  , th1dic[f'h_{name2}_pt']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_pt'], f'gae_{name}_pt')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_eta']   = ROOT.TGraphAsymmErrors(24 ) \n",
    "        efficiencydic[f'gae_{name}_eta'].SetName(f'gae_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetXaxis().SetTitle(f'h_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_eta'].Divide(th1dic[f'h_{name}_eta']  , th1dic[f'h_{name2}_eta']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_eta'], f'gae_{name}_eta')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_nvtx']   = ROOT.TGraphAsymmErrors(90 ) \n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetName(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetXaxis().SetTitle(f'h_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_nvtx'].Divide(th1dic[f'h_{name}_nvtx']  , th1dic[f'h_{name2}_nvtx']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_nvtx'], f'gae_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2 \"Tag-and-Probe\" efficiency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, in this case, you will not be accessing any generator-level information. This analyzer could be used on real data\n",
    "\n",
    "First, select the tag with the following cuts:\n",
    "* Tight muon ID\n",
    "* pT > 25 GeV and |η| < 2.4\n",
    "* Combined relative PF isolation < 0.12 (in a cone ΔR < 0.4)\n",
    "* |dZ| < 0.1\n",
    "* Match a trigger object passing a single-muon trigger (e.g. HLT_IsoMu24)\n",
    "\n",
    "Once you have found a tag, you can search for a probe. To be consistent with the previous part of the exercise, the probe should be any muon having pT > 20 GeV and |η| < 2.4 with a valid inner track. Additionally, to ensure that the tag-probe pair is compatible with a Z→μμ decay, require that tag and probe have opposite electric charges, and that their invariant mass is within a narrow window around the nominal Z mass. Since you are not performing any background subtraction here, you should choose a very narrow window, to reduce the level of background (in real data). E.g. you can choose (91 ± 5) GeV.\n",
    "\n",
    "Once you have found both a tag and a probe, you can proceed to fill the plots using the probe information, at different stages of probe selection. Make the same plots as in Part 1, starting from the reconstructed muons with a valid track, obviously you don't have a collection of \"all the generated muons\" in this case! Create the efficiency TGraphAsymmErrors as in Part 1, starting from the ID efficiencies.\n",
    "\n",
    "Finally, run the code and compare the tag-and-probe efficiencies with the MC-truth efficiencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tag_muon_mask=(muon.pt>25)& (abs(muon.eta)<2.4)& (muon.tightId==1)& (abs(muon.dz)<0.1)& (muon.pfRelIso04_all< 0.12)& trigger_match_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tag_muon=muon[tag_muon_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_muon_mask=[]\n",
    "# make sure the z mass is within 91+-5 window\n",
    "for tag4s, probe4s, in zip(tag_muon, muon):\n",
    "    list=[]\n",
    "    if len(probe4s)==0:\n",
    "        probe_muon_mask.append([])\n",
    "        continue\n",
    "    for probe4 in probe4s:\n",
    "        bestmass=0\n",
    "        if len(tag4s)==0:\n",
    "            list.append(False)\n",
    "            continue\n",
    "        for tag4 in tag4s:\n",
    "            if tag4.charge*probe4.charge>0:\n",
    "                continue\n",
    "            mass=sqrt(2*tag4.pt*probe4.pt*(cosh(tag4.eta - probe4.eta) - cos(tag4.phi - probe4.phi)))\n",
    "            if abs(mass-91.)>5.:\n",
    "                continue\n",
    "            else:\n",
    "                bestmass= mass\n",
    "                break\n",
    "        if bestmass!=0:\n",
    "            list.append(True)\n",
    "        else:\n",
    "            list.append(False)\n",
    "    probe_muon_mask.append(list)\n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_muon_masking=ak.from_iter(probe_muon_mask)& (muon.pt>20)& (abs(muon.eta)<2.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_muon= muon[probe_muon_masking]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "muon_mask={}\n",
    "muon_mask['rec']= probe_muon_masking\n",
    "muon_mask['rec_loose']= probe_muon_masking& (muon.looseId==1) \n",
    "muon_mask['rec_medium']= probe_muon_masking& (muon.mediumId==1) \n",
    "muon_mask['rec_tight']=  probe_muon_masking& (muon.tightId==1)\n",
    "muon_mask['rec_tight_iso']= probe_muon_masking& (muon.tightId==1)& (muon.pfRelIso04_all< 0.12)\n",
    "muon_mask['rec_tight_iso_hlt']= probe_muon_masking& (muon.tightId==1)& (muon.pfRelIso04_all< 0.12)& trigger_match_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "th1dic={}\n",
    "efficiencydic={}\n",
    "for name in ['rec', 'rec_loose', 'rec_medium', 'rec_tight', 'rec_tight_iso', 'rec_tight_iso_hlt']:\n",
    "    if name== 'rec':\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.pt[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.eta[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "        continue\n",
    "        \n",
    "    if name== 'rec_tight_iso_hlt':\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.pt[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.eta[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_nvtx')\n",
    "\n",
    "        nbinpt   = th1dic['h_rec_pt'].GetNbinsX() \n",
    "        nbineta  = th1dic['h_rec_eta'].GetNbinsX() \n",
    "        nbinnvtx = th1dic['h_rec_nvtx'].GetNbinsX()\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_pt']   = ROOT.TGraphAsymmErrors(nbinpt  ) \n",
    "        efficiencydic[f'gae_{name}_pt'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetXaxis().SetTitle(f'h_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_pt'].Divide(th1dic[f'h_{name}_pt']  , th1dic[f'h_{name2}_pt']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_pt'], f'gae_{name}_pt')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_eta']   = ROOT.TGraphAsymmErrors(nbineta ) \n",
    "        efficiencydic[f'gae_{name}_eta'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_eta'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetXaxis().SetTitle(f'h_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_eta'].Divide(th1dic[f'h_{name}_eta']  , th1dic[f'h_{name2}_eta']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_eta'], f'gae_{name}_eta')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_nvtx']   = ROOT.TGraphAsymmErrors(nbinnvtx ) \n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetName(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetTitle(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetXaxis().SetTitle(f'h_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_nvtx'].Divide(th1dic[f'h_{name}_nvtx']  , th1dic[f'h_{name2}_nvtx']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_nvtx'], f'gae_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "    else:\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.pt[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.eta[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_nvtx')\n",
    "        \n",
    "        nbinpt   = th1dic['h_rec_pt'].GetNbinsX() \n",
    "        nbineta  = th1dic['h_rec_eta'].GetNbinsX() \n",
    "        nbinnvtx = th1dic['h_rec_nvtx'].GetNbinsX()\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_pt']   = ROOT.TGraphAsymmErrors(nbinpt  ) \n",
    "        efficiencydic[f'gae_{name}_pt'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetXaxis().SetTitle(f'h_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_pt'].Divide(th1dic[f'h_{name}_pt']  , th1dic[f'h_{name2}_pt']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_pt'], f'gae_{name}_pt')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_eta']   = ROOT.TGraphAsymmErrors(nbineta ) \n",
    "        efficiencydic[f'gae_{name}_eta'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_eta'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetXaxis().SetTitle(f'h_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_eta'].Divide(th1dic[f'h_{name}_eta']  , th1dic[f'h_{name2}_eta']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_eta'], f'gae_{name}_eta')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_nvtx']   = ROOT.TGraphAsymmErrors(nbinnvtx ) \n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetName(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetTitle(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetXaxis().SetTitle(f'h_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_nvtx'].Divide(th1dic[f'h_{name}_nvtx']  , th1dic[f'h_{name2}_nvtx']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_nvtx'], f'gae_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Questions 1:\n",
    "* In the tag selection, you requested that the tag muon be matched to a trigger object. Why is it important to make sure that the tag has fired a trigger?\n",
    "* In the MC-truth approach, you can use both muons in the event to measure the efficiency. In the tag-and-probe, you only use the probe. Does this mean that you can use only one muon per event?\n",
    "* Do you observe significant differences between the \"true\" efficiencies and the tag-and-probe ones? In what plots and in what regions? What is the maximum size of the differences?\n",
    "What factors might impact or bias the tag-and-probe efficiency estimate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Run the \"Tag-and-Probe\" analyzer on real data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As already mentioned, the tag-and-probe approach is data driven. The analyzer does not access any generator-level information, thus you can safely run it on real data. The following sample--Run2016G_SingleMuon_NANOAOD.root is taken from the 2016G CMS data set, SingleMuon stream, requiring at least a good tag candidate and a good probe candidate per event."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events = NanoEventsFactory.from_root(\n",
    "    \"/eos/uscms/store/user/cmsdas/2025/short_exercises/muons/MuonData.root\",\n",
    "    treepath = \"Events\",\n",
    "    schemaclass=NanoAODSchema,\n",
    "    metadata={\"dataset\": \"Data\"},\n",
    "    entry_stop = 100000,\n",
    ").events()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events[\"Muon\",\"PV_npvsGood\"] = events.PV.npvsGood\n",
    "muon = events.Muon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trigger_muon= events.TrigObj[abs(events.TrigObj.id)==13]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deltaR(probe, trigger):\n",
    "    return np.sqrt((probe.eta-trigger.eta)**2+((probe.phi-trigger.phi+pi)%(2*pi)-pi)**2)\n",
    "\n",
    "def trigger_match(leptons, trigobjs, pdgid, pt, filterbit):\n",
    "        pass_pt = trigobjs.pt > pt\n",
    "        pass_id = abs(trigobjs.id) == pdgid\n",
    "        pass_filterbit = (trigobjs.filterBits & (0x1 << filterbit)) != 0\n",
    "        trigger_cands = trigobjs[pass_pt & pass_id & pass_filterbit]\n",
    "        delta_r = leptons.metric_table(trigger_cands, metric=deltaR)\n",
    "        pass_delta_r = delta_r < 0.1\n",
    "        n_of_trigger_matches = ak.sum(pass_delta_r, axis=2)\n",
    "        trig_matched_locs = n_of_trigger_matches >= 1\n",
    "\n",
    "        return trig_matched_locs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trigger_match_mask = trigger_match(muon, trigger_muon, 13, 0, 8 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tag_muon_mask=(muon.pt>25)& (abs(muon.eta)<2.4)& (muon.tightId==1)& (muon.dz<0.1)& (muon.pfRelIso04_all< 0.12)& trigger_match_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tag_muon= muon [tag_muon_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_muon_mask=[]\n",
    "# make sure the z mass is within 91+-5 window\n",
    "for tag4s, probe4s, in zip(tag_muon, muon):\n",
    "    list=[]\n",
    "    if len(probe4s)==0:\n",
    "        probe_muon_mask.append([])\n",
    "        continue\n",
    "    for probe4 in probe4s:\n",
    "        bestmass=0\n",
    "        if len(tag4s)==0:\n",
    "            list.append(False)\n",
    "            continue\n",
    "        for tag4 in tag4s:\n",
    "            if tag4.charge*probe4.charge>0:\n",
    "                continue\n",
    "            mass=sqrt(2*tag4.pt*probe4.pt*(cosh(tag4.eta - probe4.eta) - cos(tag4.phi - probe4.phi)))\n",
    "            if abs(mass-91.)>5.:\n",
    "                continue\n",
    "            else:\n",
    "                bestmass= mass\n",
    "                break\n",
    "        if bestmass!=0:\n",
    "            list.append(True)\n",
    "        else:\n",
    "            list.append(False)\n",
    "    probe_muon_mask.append(list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_mask_goodpv= events.PV.npvsGood>0 #at least a good vertex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_mask_goodpv_HLT_isomu24= (events.PV.npvsGood>0)& (events.HLT.IsoMu24==1)#fired the iso muon trigger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_muon_masking= (ak.from_iter(probe_muon_mask)& (muon.pt>20)& (abs(muon.eta)<2.4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe_muon= muon[probe_muon_masking]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "muon_mask={}\n",
    "muon_mask['rec']= probe_muon_masking\n",
    "muon_mask['rec_loose']= probe_muon_masking& (muon.looseId==1) \n",
    "muon_mask['rec_medium']= probe_muon_masking& (muon.mediumId==1) \n",
    "muon_mask['rec_tight']=  probe_muon_masking& (muon.tightId==1)& (muon.mediumId==1)\n",
    "muon_mask['rec_tight_iso']= probe_muon_masking& (muon.tightId==1)& (muon.pfRelIso04_all< 0.12)& (muon.mediumId==1)\n",
    "muon_mask['rec_tight_iso_hlt']= probe_muon_masking& (muon.tightId==1)& (muon.pfRelIso04_all< 0.12)& (muon.mediumId==1) & trigger_match_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "th1dic={}\n",
    "efficiencydic={}\n",
    "for name in ['rec', 'rec_loose', 'rec_medium', 'rec_tight', 'rec_tight_iso', 'rec_tight_iso_hlt']:\n",
    "    if name== 'rec':\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.pt[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.eta[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "        continue\n",
    "        \n",
    "    if name== 'rec_tight_iso_hlt':\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.pt[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.eta[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv_HLT_isomu24], f'h_{name}_nvtx')\n",
    "\n",
    "        nbinpt   = th1dic['h_rec_pt'].GetNbinsX() \n",
    "        nbineta  = th1dic['h_rec_eta'].GetNbinsX() \n",
    "        nbinnvtx = th1dic['h_rec_nvtx'].GetNbinsX()\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_pt']   = ROOT.TGraphAsymmErrors(nbinpt  ) \n",
    "        efficiencydic[f'gae_{name}_pt'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetXaxis().SetTitle(f'h_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_pt'].Divide(th1dic[f'h_{name}_pt']  , th1dic[f'h_{name2}_pt']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_pt'], f'gae_{name}_pt')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_eta']   = ROOT.TGraphAsymmErrors(nbineta ) \n",
    "        efficiencydic[f'gae_{name}_eta'].SetName(f'gae_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetXaxis().SetTitle(f'h_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_eta'].Divide(th1dic[f'h_{name}_eta']  , th1dic[f'h_{name2}_eta']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_eta'], f'gae_{name}_eta')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_nvtx']   = ROOT.TGraphAsymmErrors(nbinnvtx ) \n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetName(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetTitle(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetXaxis().SetTitle(f'h_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_nvtx'].Divide(th1dic[f'h_{name}_nvtx']  , th1dic[f'h_{name2}_nvtx']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_nvtx'], f'gae_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "    else:\n",
    "        th1dic[f'h_{name}_pt']=ROOT.TH1D(\"pt\"  , f'h_{name}_pt;muon p_T [GeV];Events',  20,  20.0, 120.0)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_pt'], muon.pt[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_pt')\n",
    "        th1dic[f'h_{name}_eta']=ROOT.TH1D(\"eta\" , f'h_{name}_eta;muon #eta;Events',  24,  -2.4,   2.4)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_eta'], muon.eta[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_eta')\n",
    "        th1dic[f'h_{name}_nvtx'] =ROOT.TH1D(\"nvtx\", f'h_{name}_nvtx;N. vertices;Events',  90,   0.5,  90.5)\n",
    "        TH1_fill_and_plot(th1dic[f'h_{name}_nvtx'], muon.PV_npvsGood[muon_mask[f'{name}']][event_mask_goodpv], f'h_{name}_nvtx')\n",
    "        \n",
    "        nbinpt   = th1dic['h_rec_pt'].GetNbinsX() \n",
    "        nbineta  = th1dic['h_rec_eta'].GetNbinsX() \n",
    "        nbinnvtx = th1dic['h_rec_nvtx'].GetNbinsX()\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_pt']   = ROOT.TGraphAsymmErrors(nbinpt  ) \n",
    "        efficiencydic[f'gae_{name}_pt'].SetName(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].SetTitle(f'gae_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetXaxis().SetTitle(f'h_{name}_pt')\n",
    "        efficiencydic[f'gae_{name}_pt'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_pt'].Divide(th1dic[f'h_{name}_pt']  , th1dic[f'h_{name2}_pt']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_pt'], f'gae_{name}_pt')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_eta']   = ROOT.TGraphAsymmErrors(nbineta ) \n",
    "        efficiencydic[f'gae_{name}_eta'].SetName(f'gae_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].SetTitle(f'gae_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetXaxis().SetTitle(f'h_{name}_eta')\n",
    "        efficiencydic[f'gae_{name}_eta'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_eta'].Divide(th1dic[f'h_{name}_eta']  , th1dic[f'h_{name2}_eta']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_eta'], f'gae_{name}_eta')\n",
    "        \n",
    "        efficiencydic[f'gae_{name}_nvtx']   = ROOT.TGraphAsymmErrors(nbinnvtx ) \n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetName(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].SetTitle(f'gae_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetXaxis().SetTitle(f'h_{name}_nvtx')\n",
    "        efficiencydic[f'gae_{name}_nvtx'].GetYaxis().SetTitle('Efficiency')\n",
    "        try:\n",
    "            efficiencydic[f'gae_{name}_nvtx'].Divide(th1dic[f'h_{name}_nvtx']  , th1dic[f'h_{name2}_nvtx']  , \"cl=0.683 b(1,1) mode\")\n",
    "        except:\n",
    "            print(\"An exception occurred\")\n",
    "        eff_fill_and_plot(efficiencydic[f'gae_{name}_nvtx'], f'gae_{name}_nvtx')\n",
    "        \n",
    "        name2=name\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Questions 2:\n",
    "* Are there very significant differences between the data and MC efficiencies? In what plots and regions?\n",
    "* What factors may cause differences between data and MC efficiencies? Is it possible to eliminate these differences?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Muon efficiency corrections\n",
    "\n",
    "Muon efficiencies (see Exercise four for how these are measured) be different for data and simulated events. It is therefore necessary to correct the efficiencies in MC with so-called scale factors (SF). \n",
    "For standard ID, isolation, and trigger selections, these efficiencies and SFs are provided centrally by MUO POG and can be found among the recommendations provided on this [Twiki](https://twiki.cern.ch/twiki/bin/view/CMS/TWikiPAGsMUO). The corrections are provided in a JSON format that has been standardized accross CMS by the XPOG group. These files can be read and the corrections be applied using the [correctionlib](https://github.com/cms-nanoAOD/correctionlib) package. The JSON files themselves are collected centrally and availabe in the [sonpog-integration](https://gitlab.cern.ch/cms-nanoAOD/jsonpog-integration/-/tree/master/POG/MUO?ref_type=heads) repository. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import correctionlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us apply the efficiency SFs for the tight muon ID efficiency with respect to reconstructed muons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ceval = correctionlib.CorrectionSet.from_file(\"muon_Z.json.gz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The correction is parametrized as a function of muon pT and eta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SFs = ceval[\"NUM_TightID_DEN_TrackerMuons\"].evaluate(ak.flatten(muon.eta[muon_mask['rec_tight']]), ak.flatten(muon.pt[muon_mask['rec_tight']]),'nominal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SFs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the correction to be applied to the muon efficiency in MC is rather small and in the oder of a few percent, depending on muon kinematics. These values now have to be applied as event weights in the analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}