| 1 | #!/usr/bin/env python
|
|---|
| 2 | """
|
|---|
| 3 | This script computes and prints the signal efficiency when reinterpreting the CMS analysis searching for LLPs that decay in the endcap muon detectors (https://arxiv.org/abs/2107.04838)
|
|---|
| 4 | The event-level and cluster-level selections follow the exact selections applied in the paper and the recasting instructions provided in the HEPData entry (https://www.hepdata.net/record/104408)
|
|---|
| 5 | The user would need to normalize to the correct signal cross section and luminosity to get the expected signal yield.
|
|---|
| 6 | """
|
|---|
| 7 |
|
|---|
| 8 | import sys
|
|---|
| 9 | import ROOT
|
|---|
| 10 | import math
|
|---|
| 11 | def deltaR(eta1, phi1, eta2, phi2):
|
|---|
| 12 | return (dPhi(phi1,phi2)**2+(eta1-eta2)**2)**0.5
|
|---|
| 13 |
|
|---|
| 14 | def dPhi(phi1, phi2):
|
|---|
| 15 | delta = phi1-phi2
|
|---|
| 16 | while delta > math.pi: delta -= 2* math.pi
|
|---|
| 17 | while delta < math.pi: delta += 2* math.pi
|
|---|
| 18 | return delta
|
|---|
| 19 |
|
|---|
| 20 | if __name__ == '__main__':
|
|---|
| 21 | try:
|
|---|
| 22 | input = raw_input
|
|---|
| 23 | except:
|
|---|
| 24 | pass
|
|---|
| 25 |
|
|---|
| 26 | if len(sys.argv) < 2:
|
|---|
| 27 | print(" Usage: Example1.py input_file")
|
|---|
| 28 | sys.exit(1)
|
|---|
| 29 |
|
|---|
| 30 | ROOT.gSystem.Load("libDelphes")
|
|---|
| 31 |
|
|---|
| 32 | try:
|
|---|
| 33 | ROOT.gInterpreter.Declare('#include "classes/DelphesClasses.h"')
|
|---|
| 34 | ROOT.gInterpreter.Declare('#include "external/ExRootAnalysis/ExRootTreeReader.h"')
|
|---|
| 35 | except:
|
|---|
| 36 | pass
|
|---|
| 37 |
|
|---|
| 38 | inputFile = sys.argv[1]
|
|---|
| 39 |
|
|---|
| 40 | # Create chain of root trees
|
|---|
| 41 | chain = ROOT.TChain("Delphes")
|
|---|
| 42 | chain.Add(inputFile)
|
|---|
| 43 |
|
|---|
| 44 | # Create object of class ExRootTreeReader
|
|---|
| 45 | treeReader = ROOT.ExRootTreeReader(chain)
|
|---|
| 46 | numberOfEntries = treeReader.GetEntries()
|
|---|
| 47 |
|
|---|
| 48 | # Get pointers to branches used in this analysis
|
|---|
| 49 | branchCluster = treeReader.UseBranch("CscCluster")
|
|---|
| 50 | branchMET = treeReader.UseBranch("MissingET")
|
|---|
| 51 | branchElectron = treeReader.UseBranch("Electron")
|
|---|
| 52 | branchMuon = treeReader.UseBranch("Muon")
|
|---|
| 53 | branchJet = treeReader.UseBranch("Jet")
|
|---|
| 54 | branchWeight = treeReader.UseBranch("Weight")
|
|---|
| 55 | branchEvent = treeReader.UseBranch("Event")
|
|---|
| 56 |
|
|---|
| 57 |
|
|---|
| 58 | signal_yield = 0
|
|---|
| 59 | total_weight = 0
|
|---|
| 60 |
|
|---|
| 61 | # Loop over all events
|
|---|
| 62 | for entry in range(0, numberOfEntries):
|
|---|
| 63 | # Load selected branches with data from specified event
|
|---|
| 64 | treeReader.ReadEntry(entry)
|
|---|
| 65 |
|
|---|
| 66 | ## main MC event weight
|
|---|
| 67 | w = branchWeight[0].Weight
|
|---|
| 68 | total_weight += w
|
|---|
| 69 |
|
|---|
| 70 | ################################
|
|---|
| 71 | # Event-level selections
|
|---|
| 72 | ################################
|
|---|
| 73 | # Require MET > 200 GeV
|
|---|
| 74 | if branchMET.At(0).MET < 200: continue
|
|---|
| 75 |
|
|---|
| 76 | # Require at least 1 jet with pT > 50 GeV and abs(eta) < 2.4
|
|---|
| 77 | nJet = 0
|
|---|
| 78 | for i in range(branchJet.GetEntries()):
|
|---|
| 79 | jet = branchJet.At(i)
|
|---|
| 80 | if jet.PT > 50 and abs(jet.Eta)< 2.4: nJet+=1
|
|---|
| 81 | if nJet == 0: continue
|
|---|
| 82 |
|
|---|
| 83 | # Require 0 lepton
|
|---|
| 84 | nLeptons = 0
|
|---|
| 85 | for i in range(branchElectron.GetEntries()):
|
|---|
| 86 | if branchElectron.At(i).PT > 35 and abs(branchElectron.At(i).Eta)< 2.5: nLeptons+=1
|
|---|
| 87 | for i in range(branchElectron.GetEntries()):
|
|---|
| 88 | if branchElectron.At(i).PT > 25 and abs(branchElectron.At(i).Eta)< 2.4: nLeptons+=1
|
|---|
| 89 | if nLeptons > 0:continue
|
|---|
| 90 | ################################
|
|---|
| 91 | # Cluster-level selections
|
|---|
| 92 | ################################
|
|---|
| 93 | nCscCluster = 0
|
|---|
| 94 | for i in range(branchCluster.GetEntries()):
|
|---|
| 95 | cluster = branchCluster.At(i)
|
|---|
| 96 |
|
|---|
| 97 | # check for jet veto
|
|---|
| 98 | maxJetVetoPt = 0
|
|---|
| 99 | for j in range(branchJet.GetEntries()):
|
|---|
| 100 | jet = branchJet.At(j)
|
|---|
| 101 | if deltaR(cluster.Eta, cluster.Phi, jet.Eta, jet.Phi) < 0.4:
|
|---|
| 102 | maxJetVetoPt = max(maxJetVetoPt, jet.PT)
|
|---|
| 103 | nCscCluster+=1
|
|---|
| 104 | if (maxJetVetoPt<10 and abs(dPhi(cluster.Phi, branchMET.At(0).Phi)) < 0.75 and cluster.T < 12.5 and cluster.T > -5): nCscCluster+=1
|
|---|
| 105 | if nCscCluster == 0:continue
|
|---|
| 106 |
|
|---|
| 107 | signal_yield+= w
|
|---|
| 108 |
|
|---|
| 109 | print("final signal efficiency is:" + str(signal_yield/total_weight))
|
|---|
