2018 CMS Muon HATS at the LPC

Start with creating an empty analyzer:

cd $CMSSW_BASE/src/HATSExercises
mkedanlzr MuonExercise3
cd MuonExercise3

Then copy the template code MuonExercise3.cc, the BuildFile, the configuration file MuonExercise3.py and some prepared tools:

xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise3/plugins/MuonExercise3.cc plugins/MuonExercise3.cc
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise3/MuonExercise3.py MuonExercise3.py
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise3/plugins/BuildFile.xml plugins/BuildFile.xml
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise3/test/overlayHistos.cc test/overlayHistos.cc

We will consider muons satisfying some very basic criteria:

• p_T > 20 GeV
• |η| < 2.4, i.e. within the tracker acceptance
• the muon must have an inner-tracker track

if (!vertex->isFake() && vertex->ndof()>4 && vertex->position().Rho()position().Z())<24) { ... }

To better understand the selection criteria, it is useful to classify each muon according to how it was produced:

• prompt, i.e. from the decay of a W or Z boson or from a τ lepton produced in the hard proton-proton interaction
• heavy flavor decay, i.e. from the decay of a b-quark or c-quark hadron
• light flavor decay, i.e. from the decay of a light-quark hadron, such as a pion or a kaon
• fake muons, such as punch-through or matching of random tracks and muon segments.

These are some possible variables that we will investigate:

1. whether the muon is a global muon, i.e. it was reconstructed with a combined fit of tracker and muon chambers measurements
   • muon->isGlobalMuon()
2. whether the muon is a tracker muon, i.e. it was identified by geometrically matching an inner track with segments in the muon chambers
   • muon->isTrackerMuon()
3. Normalized χ² of the global track fit
   • if(muon->isGlobalMuon()) muon->globalTrack()->normalizedChi2()
4. Number of muon chamber hits included in the global-muon track fit
   • if(muon->isGlobalMuon()) muon->globalTrack()->hitPattern().numberOfValidMuonHits()
5. Number of muon stations with matched segments
   • muon->numberOfMatchedStations()
6. Number of hits in the pixel detector
   • muon->innerTrack()->hitPattern().numberOfValidPixelHits()
7. Number of hits in the tracker layers
   • muon->innerTrack()->hitPattern().trackerLayersWithMeasurement()
8. Transverse impact parameter of the track with respect to the vertex from which the muon originated
   • muon->muonBestTrack()->dxy(firstGoodVertex->position())
9. Isolation based on the sum of p_T of charged-hadron PFCandidates from the leading primary vertex in the event, in a cone of ΔR < 0.4 around the muon
   • muon->pfIsolationR04().sumChargedHadronPt
10. Isolation calculated with neutral-hadron PFCandidates in a cone of ΔR < 0.4 around the muon
    • muon->pfIsolationR04().sumNeutralHadronEt
11. Isolation calculated with photon PFCandidates in a cone of ΔR < 0.4 around the muon
    • muon->pfIsolationR04().sumPhotonEt
12. Isolation calculated with all charged particles in a cone of ΔR < 0.4 around the muon, but not from the leading primary vertex (i.e. pileup contribution to the isolation sum)
    • muon->pfIsolationR04().sumPUPt
13. PF-based combined relative isolation, Δβ-corrected for pileup
    • (muon->pfIsolationR04().sumChargedHadronPt + max(0., mu->pfIsolationR04().sumNeutralHadronEt + mu->pfIsolationR04().sumPhotonEt - 0.5*mu->pfIsolationR04().sumPUPt)) / muon->pt()
14. Tracker-based relative isolation
    • muon->isolationR03().sumPt / muon->pt()

The methods to access these informations you can see in the list above and on the Muon ID page (look in the table of tight id requirements). For an even more in depth look you can see the documentation of pat::Muon. More information about the types of isolations can be found here.

Now we need to plot all these variables in histograms. Given the large number of variables and muon classes, a convenient way to go is to use the TFileService utility — we already tried it before — and create an array of histograms, one for each category: prompt muons, heavy flavor decays, light flavor decays, and fake muons. E.g., see the example below for the case of muon p_T.
In the class declaration:

enum MuonParentage {PROMPT, HF, LF, OTHER, N_MUPAR_TYPES};
const char* enumNames[MuonParentage::N_MUPAR_TYPES] = {"prompt", "hf", "lf", "other"};
TH1F* h_pt[MuonParentage::N_MUPAR_TYPES];

for (size_t idx=0; idx < MuonParentage::N_MUPAR_TYPES; idx++)
    h_pt[idx] = fs->make(Form("pt_%s", enumNames[idx]), Form("pt_%s", enumNames[idx]), 20,  0., 100.);

// After determining the parentage of the muon
h_pt[parentage]->Fill(muon->pt());

Complete the analyzer adding histograms for all the identification variables, then run it separately on the three samples. Then merge the three output files using the following command:

mergeTFileServiceHistograms -i histos1.root -i histos2.root -i histos3.root -o merged_histos.root

Finally, for each variable, overlay the normalized distributions from each source onto the same plot, for comparison. You can use the provided tool overlayHistos.cc, as follows:

root -l -b -q 'overlayHistos.cc("merged_histos.root")'

MuonSelectors

• Loose Muons

bool pat::Muon::isLooseMuon() 
bool muon::isLooseMuon(const reco::Muon& muon) {
    return muon.isPFMuon() && (muon.isGlobalMuon() || muon.isTrackerMuon());
}

• Medium Muons

bool pat::Muon::isMediumMuon() 
bool muon::isMediumMuon(const reco::Muon& muon) {
    if( !( isLooseMuon(muon) && muon.innerTrack()->validFraction()>0.8 )) return false; 

    bool goodGlb = muon.isGlobalMuon() && 
        muon.globalTrack()->normalizedChi2()<3. && 
        muon.combinedQuality().chi2LocalPosition<12. && 
        muon.combinedQuality().trkKink<20.; 

    return (segmentCompatibility(muon) > (goodGlb ? 0.303 : 0.451)); 
}

• Tight Muons

bool pat::Muon::isTightMuon(const reco::Vertex& vtx) 
bool muon::isTightMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
    if(!muon.isPFMuon() || !muon.isGlobalMuon()) return false;
    bool muID = muon.isGlobalMuon() && muon.globalTrack()->normalizedChi2()<10. && muon.globalTrack()->hitPattern().numberOfValidMuonHits()>0 && muon.numberOfMatchedStations()>1;  
    bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement()>5 && muon.innerTrack()->hitPattern().numberOfValidPixelHits()>0; 
    bool ip = fabs(muon.muonBestTrack()->dxy(vtx.position()))<0.2 && fabs(muon.muonBestTrack()->dz(vtx.position()))<0.5;
  
    return muID && hits && ip;
} 

• Soft Muons

bool pat::Muon::isSoftMuon(const reco::Vertex& vtx) 
bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
    bool muID = muon::isGoodMuon(muon, TMOneStationTight);
    if(!muID) return false;
  
    bool layers = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement()>5 &&
        muon.innerTrack()->hitPattern().pixelLayersWithMeasurement()>0;
    bool ishighq = muon.innerTrack()->quality(reco::Track::highPurity); 
    bool ip = fabs(muon.innerTrack()->dxy(vtx.position()))<0.3 && fabs(muon.innerTrack()->dz(vtx.position()))<20.;
  
    return layers && ip && ishighq;
}

• HighPt Muons

bool pat::Muon::isHighPtMuon(const reco::Vertex& vtx) 
bool muon::isHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx){
    bool muID = muon.isGlobalMuon() && muon.globalTrack()->hitPattern().numberOfValidMuonHits()>0 && (muon.numberOfMatchedStations()>1);
    if(!muID) return false;

    bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement()>5 &&
        muon.innerTrack()->hitPattern().numberOfValidPixelHits()>0; 
    bool momQuality = muon.tunePMuonBestTrack()->ptError()/muon.tunePMuonBestTrack()->pt()<0.3;
    bool ip = fabs(muon.innerTrack()->dxy(vtx.position()))<0.2 && fabs(muon.innerTrack()->dz(vtx.position()))<0.5;
  
  return muID && hits && momQuality && ip;
}

mu->pfIsolationR04().sumChargedHadronPt   // pT sum of charged hadrons from the main primary vertex of the event
mu->pfIsolationR04().sumNeutralHadronEt   // pT sum of neutral hadrons
mu->pfIsolationR04().sumPhotonEt          // pT sum of photons

const reco::MuonPFIsolation &pfR04 = mu->pfIsolationR04();
double combRelIso = (pfR04.sumChargedHadronPt + pfR04.sumNeutralHadronEt + pfR04.sumPhotonE) / mu->pt();   // combined relative isolation

1. the Σp_T of charged particles from pileup vertices (mu->pfIsolationR04().sumPUPt), and
2. the ratio of neutral-to-charged particles expected in LHC proton-proton collisions. From simulation studies, this ratio results to be about 0.5.

const reco::MuonPFIsolation &pfR04 = mu->pfIsolationR04();
double corrCombRelIso = (pfR04.sumChargedHadronPt + std::max(0.0, pfR04.sumNeutralHadronEt + pfR04.sumPhotonEt - 0.5*pfR04.sumPUPt)) / mu->pt();

Solutions

Answers to the questions raised during the exercise can be found in this slides.
If you are running out of time, you may copy a solution from below, look at the implementation, run it, and look at the histograms. You can also find the complete analyzer for the exercise here:
/uscms_data/d1/hats/Muon2018/Exercises/MuonExercise3/plugins/MuonExercise3_solutions.cc

cd $CMSSW_BASE/src/HATS
mkedanlzr MuonExercise4
cd MuonExercise4

Then copy the template code MuonExercise4.cc, the BuildFile and the configuration file MuonExercise4.py:

xrdcp -f root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/plugins/MuonExercise4McTruth.cc plugins/MuonExercise4McTruth.cc
xrdcp -f root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/MuonExercise4.py MuonExercise4.py
xrdcp -f root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/plugins/BuildFile.xml plugins/BuildFile.xml
xrdcp -f root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/test/muonExercise4mcTruth_cfg.py test/muonExercise4mcTruth_cfg.py
xrdcp -f root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/test/muonExercise4TagProbe_cfg.py test/muonExercise4TagProbe_cfg.py
xrdcp -f  root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/test/muonExercise4TagProbe_data_cfg.py test/muonExercise4TagProbe_data_cfg.py

Edit the analyzer MuonExercise4McTruth.cc .

First, loop over the reco::GenParticle collection, select the prompt muons with p_T > 20 GeV and |η| < 2.4, and fill histograms with the main kinematic variables. Include a distribution of the number of reconstructed vertices, to check the pileup dependence of the efficiencies. E.g.

for(auto&& genPart : *(genColl.product())) {
  // Check if it's a muon from Drell-Yan process
  if(genPart.isPromptFinalState() && std::abs(genPart.pdgId()) == 13) {
    // Only muons within acceptance and pt>20
    if(genPart.pt() > 20. && std::abs(genPart.eta())<2.4) {
      h_pt->Fill(genPart.pt());
      h_nvtx->Fill(nGoodVtx);
    }
  }
}

• all reconstructed muons with a valid inner track (i.e. skip muons that are reconstructed exclusively as standalone muons)
• reconstructed muons passing some standard ID: e.g. loose, medium, tight
• tight muons that pass a requirement on the Δβ-corrected combined relative isolation (e.g. < 0.15)
• 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 (see example in the template file).

// In the class definition
TGraphAsymmErrors *gae_rec_pt;

// In the constructor
gae_rec_pt = fs->make<TGraphAsymmErrors>();

// In the endJob or destructor
gae_rec_pt->Divide(h_rec_pt, h_pt, "cl=0.683 b(1,1) mode"); // 1-sigma errors computed with a Bayesian approach

try{gae_rec_pt->Divide(h_rec_pt, h_pt, "cl=0.683 b(1,1) mode");} catch(cms::Exception& ex) {}  

• reconstructed muons / generated muons
• muon ID / reconstructed muons
• isolation / muon ID
• single-muon trigger / isolation

cd $CMSSW_BASE/src/CMSDASExercises/MuonExercise4
cmsenv
xrdcp -f root://cmseos.fnal.gov///store/user/cmsdas/2018/short_exercises/Muons/Exercises/MuonExercise4/plugins/MuonExercise4TagProbe.cc plugins/MuonExercise4TagProbe.cc

analyzer MuonExercise4TagProbe.cc .
Note that, in this case, you will not be accessing any generator-level information. This analyzer could be used on real data.

You need two nested loops on the pat::Muon collection, one to select the tag and one to select the probe. First, select the tag with the following cuts:

• Tight muon ID
• p_T > 25 GeV and |η| < 2.4
• Combined relative PF isolation < 0.15 (in a cone ΔR < 0.4)
• Match a trigger object passing a single-muon trigger (e.g. HLT_IsoMu24 || HLT_IsoTkMu24)

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 p_T > 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.

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! In the endJob or in the destructor, create the efficiency TGraphAsymmErrors as in Part 1, starting from the ID efficiencies.

Finally, run the code and compare the tag-and-probe efficiencies with the MC-truth efficiencies.

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 is taken from the 2016G CMS data set, SingleMuon stream, requiring at least a good tag candidate and a good probe candidate per event (with tags firing " HLT_IsoMu24 || HLT_IsoTkMu24 "):
root://cmseos.fnal.gov///store/user/hats/2018/Muon/Samples/SingleMuon_Run2016G_18Apr2017_3.root
root://cmseos.fnal.gov///store/user/hats/2018/Muon/Samples/SingleMuon_Run2016G_18Apr2017_4.root
root://cmseos.fnal.gov///store/user/hats/2018/Muon/Samples/SingleMuon_Run2016G_18Apr2017_9.root
root://cmseos.fnal.gov///store/user/hats/2018/Muon/Samples/SingleMuon_Run2016G_18Apr2017_10.root
You will need to check out a new CMSSW release, the same used to process the data:

cd $CMSSW_BASE/src/CMSDASExercises/MuonExercise4
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise4/plugins/MuonExercise4McTruth.cc plugins/MuonExercise4McTruth.cc

Remove the following lines from the BuildFile.xml:

<library   file="MuonExercise4McTruth.cc" name="MuonExercise4McTruth">
  <flags   EDM_PLUGIN="1"/>
</library>

Finally, run your analyzer on the data skim above, check the efficiency plots, and compare them to the efficiencies in the simulated sample. Check also the invariant mass distributions at the different stages of the probe selection.

xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/CRAFT15.root .
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/Startup100.root .
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/Startup500.root .
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/Interfill.root .
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/Asym100.root .
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/Asym500.root .
xrdcp -f root://cmseos.fnal.gov///store/user/hats/2018/Muon/Exercises/MuonExercise5/muonExercise5a.cc .

TFile *f = TFile::Open(TString("$CMSSW_BASE/src/WSUDiLeptons/CosmicEndpoint/test/CRAFT15.root"));
TFile *g = TFile::Open(TString("$CMSSW_BASE/src/WSUDiLeptons/CosmicEndpoint/test/Startup100.root"));
TFile *h = TFile::Open(TString("$CMSSW_BASE/src/WSUDiLeptons/CosmicEndpoint/test/Startup500.root"))

First use CRAFT and the two MC_Startup files.

root -x -b -q muonExercise5a.cc++()

This will output a .root file containing all of the histograms that we are interested in.

The name of the file should be: ROOTHistograms.root

Look in the chi^2 plot to see if you notice anything interesting.

Note what value of bias the minimum of this distribution is.

std::string const& outFile = "ROOTHistograms";

When you have completed this step, you should have two .ROOT files, one for CRAFT vs Startup MC and one for Interfill vs Asym MC.

To start, copy the muonExercise5b.cc file to your working area:

xrdcp -f root://cmseos.fnal.gov///store/user/hats/2017/muon/Exercises/MuonExercise5/muonExercise5b.cc .

TFile *file = new TFile(TString("$CMSSW_BASE/src/WSUDiLeptons/CosmicEndpoint/test/"+rootFile+".root"), "READ");

Also change the input file name to match the file you would like to analyze:

<pre>std::string const& rootFile = "ROOTHistograms";

Start with CRAFT and MC Startup.

With the proper files specified, run the code using the following command:

root -x -b -q muonExercise5b.cc++()

This produces an output .root file called fittedHistograms.root. Open this with a TBrowser and look for a TCanvas called "Canvas" and check to see if the bias is consistent with zero bias.