Exercise-3-Muon-identification-and-isolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous exercises we applied cuts on various quantities embedded in the muon object without understanding what they were or why we were imposing the suggested requirements. In this exercise, we will build on the discussion in the introduction and the practice with the nanoAOD in the previous exercise to analyze the main properties and quality variables of muon tracks and how they can be used to identify muons from different sources. We will use the following NANOAODSIM files, containing simulated Drell–Yan, top-pair, and QCD multijet events:\n",
    "* DYJetsToMuMu_M-50_TuneCP5_14TeV-madgraphMLM-pythia8.root\n",
    "* TTJets_DiLept_TuneCP5_13TeV-madgraphMLM-pythia8.root\n",
    "* QCD_Pt-20toInf_MuEnrichedPt15_TuneCP5_13TeV_pythia8.root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import awkward as ak\n",
    "from coffea.nanoevents import NanoEventsFactory, NanoAODSchema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: Retrive and plot some variables\n",
    "We will consider muons satisfying some very basic criteria:\n",
    "\n",
    "* pT > 20 GeV\n",
    "* |η| < 2.4, i.e. within the tracker acceptance \n",
    "* the muon must have an inner-tracker track\n",
    "<br>\n",
    "\n",
    "We will also retrieve and select the primary proton-proton vertices with some very basic requirements:\n",
    "\n",
    "To begin, we can creat a few lists of file names, variable names and muon flavor to loop over those variables and plot histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = '/eos/uscms/store/user/cmsdas/2025/short_exercises/muons/'\n",
    "filenames= ['DYJets.root', 'TTToLL.root', 'QCDMuEnriched.root']\n",
    "variables= ['pt', 'eta', 'isGlobal', 'isTracker', 'nStations', 'nTrackerLayers', 'dxy', 'pfRelIso04_all', 'tkRelIso']\n",
    "muonid= ['looseId', 'mediumId', 'tightId', 'softId', 'highPtId']\n",
    "muonflavs= ['prompt', 'tau', 'lf', 'fake']\n",
    "bins=[20, 25, 2, 2, 6, 25, 40, 30, 30]\n",
    "binranges=[[0,100], [-2.5, 2.5], [-0.5, 1.5], [-0.5, 1.5], [-0.5, 5.5], [-0.5, 24.5], [-0.10, 0.11], [0, 0.3], [0, 0.3]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index=[0,1,2]\n",
    "n=len(filenames)\n",
    "file=[None]*n\n",
    "muon=[None]*n\n",
    "muon_genpartflav=[None]*n\n",
    "npvs=[None]*n\n",
    "npvsgood=[None]*n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eventsDY = NanoEventsFactory.from_root(\n",
    "    filepath + filenames[0],\n",
    "    treepath = \"Events\",\n",
    "    schemaclass=NanoAODSchema,\n",
    "    metadata={\"dataset\": \"DYJets\"},\n",
    "    entry_stop = 100000,\n",
    ").events()\n",
    "\n",
    "eventsTT = NanoEventsFactory.from_root(\n",
    "    filepath + filenames[1],\n",
    "    treepath = \"Events\",\n",
    "    schemaclass=NanoAODSchema,\n",
    "    metadata={\"dataset\": \"TT\"},\n",
    "    entry_stop = 100000,    \n",
    ").events()\n",
    "\n",
    "\n",
    "eventsQCD = NanoEventsFactory.from_root(\n",
    "    filepath + filenames[2],\n",
    "    treepath = \"Events\",\n",
    "    schemaclass=NanoAODSchema,\n",
    "    metadata={\"dataset\": \"QCD\"},\n",
    "    entry_stop = 100000,    \n",
    ").events()\n",
    "\n",
    "events = [eventsDY, eventsTT, eventsQCD]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#getting the primary vertices\n",
    "for idx in index:\n",
    "    npvsgood[idx]= events[idx].PV.npvsGood\n",
    "    events[idx][\"Muon\",\"PV_npvs\"] = events[idx].PV.npvs\n",
    "    print (events[idx].Muon.PV_npvs)\n",
    "    muon[idx] = events[idx].Muon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to proceed, each event must have at least one good vertex(selection:!isFake && ndof > 4 && abs(z) <= 24 && position.Rho <= 2), this information is stored in PV_npvsGood "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So every event have at lest one good vertex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To better understand the selection criteria, it is useful to classify each muon according to how it was produced:\n",
    "\n",
    "* prompt, i.e. from the decay of a W or Z boson or from a τ lepton produced in the hard proton-proton interaction\n",
    "* heavy flavor decay, i.e. from the decay of a b-quark or c-quark hadron\n",
    "* light flavor decay, i.e. from the decay of a light-quark hadron, such as a pion or a kaon\n",
    "* fake muons, such as punch-through or matching of random tracks and muon segments.\n",
    "\n",
    "In nanoAOD file, the above information is stored in Muon_genPartFlav, here is its description:<br>\n",
    "Flavour of genParticle for MC matching to status==1 muons: 1 = prompt muon (including gamma*->mu mu), 15 = muon from prompt tau, 5 = muon from b, 4 = muon from c, 3 = muon from light or unknown, 0 = unmatched <br>\n",
    "\n",
    "These are some possible variables that we will investigate:<br>\n",
    "whether the muon is a global muon, i.e. it was reconstructed with a combined fit of tracker and muon chambers measurements\n",
    "* Muon_isGlobal\n",
    "\n",
    "whether the muon is a tracker muon, i.e. it was identified by geometrically matching an inner track with segments in the muon chambers\n",
    "* Muon_isTracker\n",
    "\n",
    "Number of muon stations with matched segments\n",
    "* Muon_nStations\n",
    "\n",
    "Number of hits in the tracker layers\n",
    "* Muon_nTrackerLayers\n",
    "\n",
    "Transverse impact parameter of the track with respect to the vertex from which the muon originated\n",
    "* Muon_dxy\n",
    "\n",
    "PF-based combined relative isolation, Δβ-corrected for pileup\n",
    "* Muon_pfRelIso04_all\n",
    "\n",
    "Tracker-based relative isolation\n",
    "* Muon_tkRelIso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot them "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diclist=[]\n",
    "for idx in index:\n",
    "    muon_mask=(muon[idx].pt>20)& (abs(muon[idx].eta)<2.4)\n",
    "    dic={}\n",
    "    for flav, type in zip([1, 15 , 3, 0], muonflavs):\n",
    "        muon_mask_all= muon_mask& (muon[idx].genPartFlav==flav)\n",
    "        event_mask= npvsgood[idx]>0 #at least one good primary vertex\n",
    "        for var in variables:\n",
    "            dic[f'{var}_{type}']=(muon[idx])[f'{var}'][muon_mask_all][event_mask]\n",
    "        for id in muonid:        \n",
    "            id_mask_combined=muon_mask_all& (((muon[idx])[f'{id}'])!=0)\n",
    "            dic[f'{id}_{type}_pt']=(muon[idx])['pt'][id_mask_combined][event_mask]\n",
    "    \n",
    "    muon_mask_tightid=muon_mask& (muon[idx].tightId!=0)\n",
    "    muon_mask_tightid_pfreliso04_all=muon_mask_tightid& (muon[idx].pfRelIso04_all<0.15)\n",
    "    dic['npvs_muon_mask']= (muon[idx].PV_npvs)[event_mask]\n",
    "    dic['npvs_muon_mask_tightid']=(muon[idx].PV_npvs)[muon_mask_tightid][event_mask]\n",
    "    dic['npvs_muon_mask_tightid_pfreliso04_all']=(muon[idx].PV_npvs)[muon_mask_tightid_pfreliso04_all][event_mask]\n",
    "    diclist.append(dic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_variables={}\n",
    "for type in muonflavs:\n",
    "        for var in variables:\n",
    "            combined_variables[f'{var}_{type}']=np.concatenate((ak.flatten(diclist[0][f'{var}_{type}']), ak.flatten(diclist[1][f'{var}_{type}']), ak.flatten(diclist[2][f'{var}_{type}'])))\n",
    "        for id in muonid:\n",
    "            combined_variables[f'{id}_{type}_pt']=np.concatenate((ak.flatten(diclist[0][f'{id}_{type}_pt']), ak.flatten(diclist[1][f'{id}_{type}_pt']), ak.flatten(diclist[2][f'{id}_{type}_pt'])))\n",
    "for name in ['npvs_muon_mask','npvs_muon_mask_tightid','npvs_muon_mask_tightid_pfreliso04_all']:\n",
    "    combined_variables[name]=np.concatenate((ak.flatten(diclist[0][name]), ak.flatten(diclist[1][name]), ak.flatten(diclist[2][name])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for var, bin, binrange in zip(variables, bins, binranges):\n",
    "    plt.figure(figsize=(9,6))\n",
    "    for type, color in zip(muonflavs, ['red', 'blue', 'green', 'black']):\n",
    "        plt.hist((combined_variables[f'{var}_{type}']), bins=bin, range=binrange, density=True, histtype='step', linewidth=2, color=f'{color}', label=f'{var}_{type}')\n",
    "    plt.legend()\n",
    "    plt.xlabel(f'{var}')\n",
    "    plt.ylabel('Events')\n",
    "    #plt.yscale('log')       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions 1**:\n",
    "* What types of muons would you expect from each sample?\n",
    "* For each variable, overlay the histograms from different sources. Can you explain the differences you observe?\n",
    "* Explain the differences in the number of matched muon stations, number of hit in the tracker layers and isolation variables.\n",
    "* By comparing the various distributions, could you come up with a possible set of cuts to isolate each source of muons?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: Standard muon definitions in CMS\n",
    "The CMS Muon Physics Object Group (MUO POG) takes care of everything that concerns muon reconstruction, identification, high-level triggering, performance evaluation and monitoring, corrections, use in physics analysis, etc. Among other tasks, it develops and maintains a number of standard identification and isolation criteria, which are broadly used in analysis across all CMS. The full list of official definitions can be found on this [TWiki page](https://twiki.cern.ch/twiki/bin/view/CMS/SWGuideMuonIdRun2). Here you can find a summary of the most common criteria. These definitions are conveniently implemented as selectors for reco::Muon objects in the [MuonSelectors](https://github.com/cms-sw/cmssw/blob/master/DataFormats/MuonReco/interface/MuonSelectors.h) class, as well as functions in the pat::Muon class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* *Loose Muons*\n",
    "```c++\n",
    "bool pat::Muon::isLooseMuon() \n",
    "bool muon::isLooseMuon(const reco::Muon& muon) {\n",
    "    return muon.isPFMuon() && (muon.isGlobalMuon() || muon.isTrackerMuon());\n",
    "}\n",
    "```\n",
    "\n",
    "* *Medium Muons*\n",
    "```c++\n",
    "bool pat::Muon::isMediumMuon() \n",
    "bool muon::isMediumMuon(const reco::Muon& muon) {\n",
    "    if( !( isLooseMuon(muon) && muon.innerTrack()->validFraction()>0.8 )) return false; \n",
    "\n",
    "    bool goodGlb = muon.isGlobalMuon() && \n",
    "        muon.globalTrack()->normalizedChi2()<3. && \n",
    "        muon.combinedQuality().chi2LocalPosition<12. && \n",
    "        muon.combinedQuality().trkKink<20.; \n",
    "\n",
    "    return (segmentCompatibility(muon) > (goodGlb ? 0.303 : 0.451)); \n",
    "}\n",
    "```\n",
    " \n",
    "* *Tight Muons*\n",
    "```c++\n",
    "bool pat::Muon::isTightMuon(const reco::Vertex& vtx) \n",
    "bool muon::isTightMuon(const reco::Muon& muon, const reco::Vertex& vtx) {\n",
    "    if(!muon.isPFMuon() || !muon.isGlobalMuon()) return false;\n",
    "    bool muID = muon.isGlobalMuon() && muon.globalTrack()->normalizedChi2()<10. && muon.globalTrack()->hitPattern().numberOfValidMuonHits()>0 && muon.numberOfMatchedStations()>1;  \n",
    "    bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement()>5 && muon.innerTrack()->hitPattern().numberOfValidPixelHits()>0; \n",
    "    bool ip = fabs(muon.muonBestTrack()->dxy(vtx.position()))<0.2 && fabs(muon.muonBestTrack()->dz(vtx.position()))<0.5;\n",
    "  \n",
    "    return muID && hits && ip;\n",
    "} \n",
    "```\n",
    "\n",
    "* *Soft Muons* \n",
    "```c++ \n",
    "bool pat::Muon::isSoftMuon(const reco::Vertex& vtx) \n",
    "bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx) {\n",
    "    bool muID = muon::isGoodMuon(muon, TMOneStationTight);\n",
    "    if(!muID) return false;\n",
    "  \n",
    "    bool layers = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement()>5 &&\n",
    "        muon.innerTrack()->hitPattern().pixelLayersWithMeasurement()>0;\n",
    "    bool ishighq = muon.innerTrack()->quality(reco::Track::highPurity); \n",
    "    bool ip = fabs(muon.innerTrack()->dxy(vtx.position()))<0.3 && fabs(muon.innerTrack()->dz(vtx.position()))<20.;\n",
    "  \n",
    "    return layers && ip && ishighq;\n",
    "}\n",
    "```\n",
    "\n",
    "* *HighPt Muons* \n",
    "\n",
    "```c++ \n",
    "bool pat::Muon::isHighPtMuon(const reco::Vertex& vtx) \n",
    "bool muon::isHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx){\n",
    "    bool muID = muon.isGlobalMuon() && muon.globalTrack()->hitPattern().numberOfValidMuonHits()>0 &(muon.numberOfMatchedStations()>1);\n",
    "    if(!muID) return false;\n",
    "\n",
    "    bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement()>5 &&\n",
    "        muon.innerTrack()->hitPattern().numberOfValidPixelHits()>0; \n",
    "    bool momQuality = muon.tunePMuonBestTrack()->ptError()/muon.tunePMuonBestTrack()->pt()<0.3;\n",
    "    bool ip = fabs(muon.innerTrack()->dxy(vtx.position()))<0.2 && fabs(muon.innerTrack()->dz(vtx.position()))<0.5;\n",
    "  \n",
    "  return muID && hits && momQuality && ip;\n",
    "}\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions 2**:\n",
    "* Given the implementation of the functions above, can you guess which types of muons each selection is meant to address, and why?\n",
    "* The HighPt muon definition is very similar to the Tight muon one, but it adds a requirement on the transverse momentum resolution, σ(pT)/pT < 0.3. Can you find an explanation?\n",
    "* Verify your guesses by adding these selections in your analyzer. Produce distributions for some variables (e.g. pT) for each muon category, after applying each of the selections above, one at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for type in muonflavs:\n",
    "    plt.figure(figsize=(9,6))\n",
    "    for id, color in zip(muonid, ['red', 'purple', 'blue', 'green', 'black']):\n",
    "        plt.hist((combined_variables[f'{id}_{type}_pt']), bins=20, range=[0, 100], density=False, histtype='step',linewidth=2, color=f'{color}', label=f'{id}_{type}_pt')\n",
    "    plt.legend()\n",
    "    plt.xlabel(f'{type}_pt [GeV]')\n",
    "    plt.ylabel('Events')\n",
    "    #plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Muon isolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now take a detailed look at the isolation variables mentioned at the beginning of Exercise 3. The most common muon isolation algorithm in CMS makes use of the PF candidates found in a region of &Delta;R < 0.4 around the muon track: \n",
    "```c++ \n",
    "mu->pfIsolationR04().sumChargedHadronPt   // pT sum of charged hadrons from the main primary vertex of the event \n",
    "mu->pfIsolationR04().sumNeutralHadronEt   // pT sum of neutral hadrons \n",
    "mu->pfIsolationR04().sumPhotonEt          // pT sum of photons \n",
    "```\n",
    "In order to exploit the full-detector information, these variables can be combined in a single isolation variable: \n",
    "```c++ \n",
    "const reco::MuonPFIsolation &pfR04 = mu->pfIsolationR04();\n",
    "double combRelIso = (pfR04.sumChargedHadronPt + pfR04.sumNeutralHadronEt + pfR04.sumPhotonE) / mu->pt();   // combined relative isolation \n",
    "```\n",
    "The combined isolation turns out to perform better than the individual components separately in terms of efficiency _vs_ background rejection.<BR> \n",
    "Note that for neutral particles (photons and neutral hadrons) it is impossible to determine the vertex they originated from, since they don't have a track. Therefore neutral particles from pileup vertices contribute to the $P_T$ sum, and the performance of the combined isolation results to be strongly dependent on the pileup level. Corrections are available to mitigate such effect. The most common in CMS is called \"&Delta;&beta; correction\": it estimates the &Sigma;$P_T$ of neutral particles coming from pileup vertices using \n",
    "   1. the &Sigma;$P_T$ of charged particles from pileup vertices (mu->pfIsolationR04().sumPUPt), and \n",
    "   1. the ratio of neutral-to-charged particles expected in LHC proton-proton collisions. From simulation studies, this ratio results to be about 0.5. \n",
    "\n",
    "We can now define a &Delta;&beta;-corrected combined relative isolation_, less sensitive to the number of pileup vertices: \n",
    "\n",
    "```c++ \n",
    "const reco::MuonPFIsolation &pfR04 = mu->pfIsolationR04();\n",
    "double corrCombRelIso = (pfR04.sumChargedHadronPt + std::max(0.0, pfR04.sumNeutralHadronEt + pfR04.sumPhotonEt - 0.5*pfR04.sumPUPt)) / mu->pt();\n",
    "``` \n",
    "\n",
    "All the variables described above can be found in a minAOD file. In nanoAOD, the _&Delta;&beta;-corrected combined relative isolation is already calculated for you. It is stored as Muon_pfRelIso04_all. In addition, it is the only isolation variable in a region of &Delta;R < 0.4 around the muon track stored in nanoAOD currently.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Questions 3**:\n",
    "* In your analyzer, add distributions of the number of reconstructed primary vertices at different stages of the selection:\n",
    " 1. for all muons,\n",
    " 1. for muons passing the Tight ID,\n",
    " 1. for muons passing the Tight ID plus the same cut on Muon_pfRelIso04_all, e.g. Muon_pfRelIso04_all < 0.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9,6))\n",
    "for name, color in zip(['npvs_muon_mask','npvs_muon_mask_tightid','npvs_muon_mask_tightid_pfreliso04_all'],['red', 'blue', 'green']):\n",
    "    plt.hist((combined_variables[name]), bins=100, range=[0, 100], density=False, histtype='step',linewidth=2, color=f'{color}', label=name)\n",
    "plt.legend()\n",
    "plt.xlabel('number of pvs')\n",
    "plt.ylabel('Events')\n",
    "    "
   ]
  },
  {
   "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
}