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