""" This file provides a parameterized implementation of the probability of a LLP decay to generate an out-of-time trackless jet, described in the CMS analysis "Search for long-lived particles using out-of-time trackless jets in proton-proton collisions at sqrt(s) = 13 TeV". The provided parameterization can reproduce the number of out-of-time trackless jets from full simulation to within 7% for various LLP masses (127-1800 GeV), lifetimes (0.5 - 3 m), and decay modes (a Higgs or a Z boson, a pair of b quarks or a pair of light quarks) considered. A detailed instruction on the usage of the code is provided in the HEPData record: http://dx.doi.org/10.17182/hepdata.135827 """ import numpy as np import math import ROOT as rt def deltaPhi(phi1, phi2): """ Function that computes the delta phi between two particles. """ return (phi1 - phi2 + np.pi) % (2 * np.pi) - np.pi def deltaR(eta1, phi1, eta2, phi2): """ Function that computes the delta R between two particles. """ deta = eta1 - eta2 dphi = deltaPhi(phi1, phi2) return np.sqrt(deta*deta + dphi*dphi) def merged(LLP_r,LLP_z,boson_pt,boson_eta,quark1_pt,quark1_eta,quark1_phi,quark2_pt,quark2_eta,quark2_phi): """ Function that takes as input the LLP r and z in meters, the boson and quarks pt in GeV, and the boson and quarks eta and phi. It calculates the acceptance that the LLP decay has a merged topology. It outputs the LLP r and z satisfying the merged topology requirement. """ dR = deltaR(quark1_eta,quark1_phi,quark2_eta,quark2_phi) merged = (dR<0.8) acc_m = np.logical_and( np.logical_and(boson_pt>30.,abs(boson_eta)<1.) , merged) return LLP_r[acc_m], LLP_z[acc_m] def resolved_q(LLP_r,LLP_z,quark1_pt,quark1_eta,quark1_phi,quark2_pt,quark2_eta,quark2_phi): """ Function that takes as input the LLP r and z in meters, the quarks pt in GeV, and the quarks eta and phi. It calculates the acceptance that the LLP decay has a resolved topology, with one and only one quark in acceptance. It outputs the LLP r and z satisfying the resolved topology requirement, with one and only one quark in acceptance. """ dR = deltaR(quark1_eta,quark1_phi,quark2_eta,quark2_phi) resolved = (dR>=0.8) acc_r_q = np.logical_and( np.logical_xor( np.logical_and(quark1_pt>30.,abs(quark1_eta)<1.,), np.logical_and(quark2_pt>30.,abs(quark2_eta)<1.,) ), resolved) return LLP_r[acc_r_q], LLP_z[acc_r_q] def resolved_qq(LLP_r,LLP_z,quark1_pt,quark1_eta,quark1_phi,quark2_pt,quark2_eta,quark2_phi): """ Function that takes as input the LLP r and z in meters, the quarks pt in GeV, and the quarks eta and phi. It calculates the acceptance that the LLP decay has a resolved topology, with both quarks in acceptance. It outputs the LLP r and z satisfying the resolved topology requirement, with both quarks in acceptance. """ dR = deltaR(quark1_eta,quark1_phi,quark2_eta,quark2_phi) resolved = (dR>=0.8) acc_r_qq = np.logical_and( np.logical_and( np.logical_and(quark1_pt>30.,abs(quark1_eta)<1.), np.logical_and(quark2_pt>30.,abs(quark2_eta)<1.) ), resolved) return LLP_r[acc_r_qq], LLP_z[acc_r_qq] def predict(LLP_r,LLP_z,eff_map): """ Function that predicts the probability that the LLP decay generates one out-of-time trackless jet. It takes as input the LLP decay position z and r in m, and the corresponding efficiency map (a TH2D root histogram). It outputs the probability and its uncertainty as numpy arrays. More details in the HEPData entry. """ weights = [] weights_unc = [] for i in range(len(LLP_r)): b = eff_map.FindBin(LLP_r[i],LLP_z[i]) weights.append(eff_map.GetBinContent(b)) weights_unc.append(eff_map.GetBinError(b)) weights = np.array(weights) weights_unc = np.array(weights_unc) return weights, weights_unc def load_maps(LLP_mass): """ Function that loads the efficiency maps taking as input the particle mass. Path to ROOT file needs to be replaced. """ assert(type(LLP_mass)==int or type(LLP_mass)==float) masses = np.array([127,150,175,200,250,300,400,600,800,1000,1250,1500,1800]) #Select the map derived with the closest LLP mass w.r.t. the input value mass = masses[ np.argmin( np.abs(masses - LLP_mass) ) ] dict_fig_mass_m = {} dict_fig_mass_rq = {} dict_fig_mass_rqq = {} dict_fig_mass_m[127] = "006" dict_fig_mass_rq[127] = "007" dict_fig_mass_rqq[127] = "008" dict_fig_mass_m[150] = "009" dict_fig_mass_rq[150] = "010" dict_fig_mass_rqq[150] = "011" dict_fig_mass_m[175] = "012" dict_fig_mass_rq[175] = "013" dict_fig_mass_rqq[175] = "014" dict_fig_mass_m[200] = "015" dict_fig_mass_rq[200] = "016" dict_fig_mass_rqq[200] = "017" dict_fig_mass_m[250] = "018" dict_fig_mass_rq[250] = "019" dict_fig_mass_rqq[250] = "020" dict_fig_mass_m[300] = "021" dict_fig_mass_rq[300] = "022" dict_fig_mass_rqq[300] = "023" dict_fig_mass_m[400] = "024" dict_fig_mass_rq[400] = "025" dict_fig_mass_rqq[400] = "026" dict_fig_mass_m[600] = "027" dict_fig_mass_rq[600] = "028" dict_fig_mass_rqq[600] = "029" dict_fig_mass_m[800] = "030" dict_fig_mass_rq[800] = "031" dict_fig_mass_rqq[800] = "032" dict_fig_mass_m[1000] = "033" dict_fig_mass_rq[1000] = "034" dict_fig_mass_rqq[1000] = "035" dict_fig_mass_m[1250] = "036" dict_fig_mass_rq[1250] = "037" dict_fig_mass_rqq[1250] = "038" dict_fig_mass_m[1500] = "039" dict_fig_mass_rq[1500] = "040" dict_fig_mass_rqq[1500] = "041" dict_fig_mass_m[1800] = "042" dict_fig_mass_rq[1800] = "043" dict_fig_mass_rqq[1800] = "044" print("Loading efficiency map for LLP mass ", mass) PATH = "/PATH_TO_EFFICIENCY_MAPS/" rootFile_m = PATH+'/Figure-aux_'+ dict_fig_mass_m[mass]+'.root' rootFile_r_q = PATH+'/Figure-aux_'+ dict_fig_mass_rq[mass]+'.root' rootFile_r_qq = PATH+'/Figure-aux_'+ dict_fig_mass_rqq[mass]+'.root' hist_name = 'h2' file_m = rt.TFile(rootFile_m, 'READ') file_r_q = rt.TFile(rootFile_r_q, 'READ') file_r_qq = rt.TFile(rootFile_r_qq, 'READ') eff_map_m = file_m.Get(hist_name) eff_map_r_q = file_r_q.Get(hist_name) eff_map_r_qq = file_r_qq.Get(hist_name) eff_map_m.SetDirectory(0) eff_map_r_q.SetDirectory(0) eff_map_r_qq.SetDirectory(0) file_m.Close() file_r_q.Close() file_r_qq.Close() return eff_map_m, eff_map_r_q, eff_map_r_qq if __name__ == "__main__": """ Dummy IO of input lists and efficiency parameterization (Additional Figure XXX). More details in the HEPData entry. Input lists z, r, boson and quarks pt, phi, eta need to be replaced. LLP r and z must be provided in meters (not cm). """ LLP_mass = 275 LLP_r = np.array([ 2.87450266, 0.1952939, 0.46330371, 0.11667854, 1.35574925, 0.03328017, 0.2940428, 0.28128052, 0.02864034]) LLP_z = np.array([ 1.35661197, 0.22584859, 0.14422756, 0.48841801, 2.43715215, 0.14440206, 0.13775492, 0.53336173, 0.15120745]) boson_pt = np.array([ 606.480896, 530.89263916, 322.73410034, 351.94665527, 133.66444397, 119.462883, 287.01843262, 22.0531044, 88.31693268]) boson_eta = np.array([ 0.57100648, 0.59867728, 0.24541342, 1.70422399, 1.05397666, 1.23075068, 0.06160917, 0.07102091, 1.29030764]) quark1_pt = np.array([ 524.51318359, 439.63412476, 186.05549622, 313.06686401, 148.00857544, 87.64628601, 235.77085876, 50.70792007, 110.85513306]) quark1_eta = np.array([ 0.63173032, 0.69264203, 0.06470954, 1.77668536, 0.67247403, 1.22006834, 0.07688265, 0.76620543, 0.98105001]) quark1_phi = np.array([ 0.82929206, 2.28655791, 0.94022733, 0.52706158, 0.0355283, 1.11143875, 2.58550024, 2.48103976, 1.06563103]) quark2_pt = np.array([ 85.37778473, 91.36347961, 154.35249329, 41.38116074, 23.88053703, 77.5853653, 76.88776398, 43.96770477, 38.72476959]) quark2_eta = np.array([ 0.13542943, 0.08433751, 0.42713669, 0.79854375, 1.6709702, 0.62330586, 0.00585777, 0.83615202, 0.52454126]) quark2_phi = np.array([ 1.13482451, 2.23388171, 1.59046102, 0.89850473, 2.31684518, 0.41743258, 1.6284169, 0.21210209, 2.87484264]) assert(len(LLP_z) == len(LLP_r) == len(boson_pt) == len(boson_eta) == len(quark1_pt) == len(quark1_eta) == len(quark1_phi) == len(quark2_pt) == len(quark2_eta) == len(quark2_phi)) eff_map_m, eff_map_rq, eff_map_rqq = load_maps(LLP_mass) LLP_r_m,LLP_z_m = merged(LLP_r,LLP_z,boson_pt,boson_eta,quark1_pt,quark1_eta,quark1_phi,quark2_pt,quark2_eta,quark2_phi) weights_m, weights_unc_m = predict(LLP_r_m,LLP_z_m,eff_map_m) LLP_r_rq,LLP_z_rq = resolved_q(LLP_r,LLP_z,quark1_pt,quark1_eta,quark1_phi,quark2_pt,quark2_eta,quark2_phi) weights_rq, weights_unc_rq = predict(LLP_r_rq,LLP_z_rq,eff_map_rq) LLP_r_rqq,LLP_z_rqq = resolved_qq(LLP_r,LLP_z,quark1_pt,quark1_eta,quark1_phi,quark2_pt,quark2_eta,quark2_phi) weights_rqq, weights_unc_rqq = predict(LLP_r_rqq,LLP_z_rqq,eff_map_rqq) n_m = np.sum(weights_m) n_m_unc = np.sqrt( np.sum( weights_unc_m*weights_unc_m ) ) n_rq = np.sum(weights_rq) n_rq_unc = np.sqrt( np.sum( weights_unc_rq*weights_unc_rq ) ) n_rqq = np.sum(weights_rqq) n_rqq_unc = np.sqrt( np.sum( weights_unc_rqq*weights_unc_rqq ) ) print("Predicted merged: ", n_m, " +- ", n_m_unc) print("Predicted resolved q: ", n_rq, " +- ", n_rq_unc) print("Predicted resolved qq: ", n_rqq, " +- ", n_rqq_unc)