Rivet analyses referenceATLAS_2018_I1676551Recursive jigsaw chargino-neutralino search with 2 or 3 charged leptons in 36/fb of 13 TeV ppExperiment: ATLAS (LHC) Inspire ID: 1676551 Status: VALIDATED NOTREENTRY SINGLEWEIGHT Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
A search for electroweak production of supersymmetric particles in two-lepton and three-lepton final states using recursive jigsaw reconstruction, a technique that assigns reconstructed objects to the most probable hemispheres of the decay trees, allowing one to construct tailored kinematic variables to separate the signal and background. The search uses data collected in 2015 and 2016 by the ATLAS experiment in $\sqrt{s}=13$ TeV proton-proton collisions at the CERN Large Hadron Collider corresponding to an integrated luminosity of 36.1/fb. Chargino-neutralino pair production, with decays via W/Z bosons, is studied in final states involving leptons and jets and missing transverse momentum for scenarios with large and intermediate mass splittings between the parent particle and lightest supersymmetric particle, as well as for the scenario where this mass splitting is close to the mass of the $Z$ boson. The latter case is challenging since the vector bosons are produced with kinematic properties that are similar to those in Standard Model processes. Results are found to be compatible with the Standard Model expectations in the signal regions targeting large and intermediate mass splittings, and chargino-neutralino masses up to 600 GeV are excluded at 95% confidence level for a massless lightest supersymmetric particle. Excesses of data above the expected background are found in the signal regions targeting low mass splittings, and the largest local excess amounts to 3.0 standard deviations. Source code: ATLAS_2018_I1676551.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/DirectFinalState.hh"
5#include "Rivet/Projections/IndirectFinalState.hh"
6#include "Rivet/Projections/TauFinder.hh"
7#include "Rivet/Projections/FastJets.hh"
8#include "Rivet/Projections/MissingMomentum.hh"
9#include "Rivet/Projections/Smearing.hh"
10#include "Rivet/Tools/Cutflow.hh"
11
12namespace Rivet {
13
14
15 /// Recursive jigsaw chargino-neutralino search with 2 or 3 charged leptons in 36/fb of 13 TeV pp
16 ///
17 /// @author Derek Yeung, Andy Buckley
18 class ATLAS_2018_I1676551 : public Analysis {
19 public:
20
21 /// Constructor
22 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1676551);
23
24
25 /// Analysis initialization
26 void init() {
27
28 PromptFinalState electrons(Cuts::abspid == PID::ELECTRON);
29 SmearedParticles recoelectrons(electrons, ELECTRON_EFF_ATLAS_RUN2_MEDIUM, ELECTRON_SMEAR_ATLAS_RUN2);
30 declare(recoelectrons, "Electrons");
31
32 PromptFinalState muons(Cuts::abspid == PID::MUON);
33 SmearedParticles recomuons(muons, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2);
34 declare(recomuons, "Muons");
35
36 FastJets jets4(IndirectFinalState(Cuts::open()), JetAlg::ANTIKT, 0.4);
37 SmearedJets recojets4(jets4, JET_SMEAR_CMS_RUN2, JET_BTAG_EFFS(0.77, 1/6., 1/134.));
38 declare(recojets4, "Jets");
39
40 MissingMomentum met(FinalState(Cuts::abseta < 4.9));
41 SmearedMET recomet(met, MET_SMEAR_ATLAS_RUN2);
42 declare(recomet, "MET");
43
44 // Cutflow Setup for 2l-High
45 const strings edges = {"2l-high", "2l-int", "2l-low", "2l-ISR",
46 "3l-high", "3l-int", "3l-low", "3l-ISR" };
47 book(_cutflows, edges);
48 const strings cfnames1 = {"Trigger matching & 2 signal leptons", "Preselection",
49 "Fraction 1 > 0.8", "Fraction 2 < 0.05",
50 "Delta Phi in [0.3, 2.9]", "H_PP_4,1 > 800 GeV"};
51 const strings cfnames2 = {"Trigger matching & 2 signal leptons", "Preselection",
52 "Fraction 1 > 0.8", "Fraction 2 < 0.05",
53 "Delta Phi in [0.3, 2.6]", "H_PP_4,1 > 600 GeV"};
54 const strings cfnames3 = {"Trigger matching & 2 signal leptons", "Preselection",
55 "Fraction 1 in [0.35, 0.6]", "Fraction 2 < 0.05",
56 "Min Delta Phi > 2.4", "H_PP_4,1 > 400 GeV"};
57 const strings cfnames4 = {"Trigger matching & 2 signal leptons", "Preselection",
58 "m_Z in [80, 100] GeV", "m_J in [50, 110] GeV",
59 "Delta Phi > 2.8", "R_ISR in [0.4, 0.75]", "p_CM_T_ISR > 180 GeV",
60 "p_CM_T_I > 100 GeV", "p_CM_T < 30 GeV"};
61 const strings cfnames5 = {"Trigger matching & 3 signal leptons", "Preselection",
62 "m_ll in [75,105] GeV", "m_W_T > 150 GeV",
63 "Fraction 1 > 0.75", "Fraction 2 < 0.8",
64 "H_PP_31 > 500 GeV", "Fraction 3 < 0.2"};
65 const strings cfnames6 = {"Trigger matching & 3 signal leptons", "Preselection",
66 "m_ll in [75,105] GeV", "m_W_T > 130 GeV",
67 "Fraction 1 > 0.8", "Fraction 2 < 0.75",
68 "H_PP_31 > 450 GeV", "Fraction 3 < 0.15"};
69 const strings cfnames7 = {"Trigger matching & 3 signal leptons", "Preselection",
70 "m_ll in [75,105] GeV", "m_W_T > 100 GeV",
71 "Fraction 1 > 0.9", "H_PP_31 > 250 GeV", "Fraction 2 < 0.05"};
72 const strings cfnames8 = {"Trigger matching & 3 signal leptons", "Preselection",
73 "m_ll in [75, 105] GeV", "m_W_T > 100 GeV",
74 "Delta Phi > 2.0", "R_ISR in [0.55, 1.0]", "p_CM_T_ISR > 100 GeV",
75 "p_CM_T_I > 80 GeV", "p_CM_T < 25 GeV"};
76 book(_cutflows->bin(1), edges[0], cfnames1); // Cutflow Setup for 2l-High
77 book(_cutflows->bin(2), edges[1], cfnames2); // Cutflow Setup for 2l-Int
78 book(_cutflows->bin(3), edges[2], cfnames3); // Cutflow Setup for 2l-Low
79 book(_cutflows->bin(4), edges[3], cfnames4); // Cutflow Setup for 2l-ISR
80 book(_cutflows->bin(5), edges[4], cfnames5); // Cutflow Setup for 3l-High
81 book(_cutflows->bin(6), edges[5], cfnames6); // Cutflow Setup for 3l-Int
82 book(_cutflows->bin(7), edges[6], cfnames7); // Cutflow Setup for 3l-Low
83 book(_cutflows->bin(8), edges[7], cfnames8); // Cutflow Setup for 3l-ISR
84
85 }
86
87
88 // Per-event analysis
89 void analyze(const Event& event) {
90
91 _cutflows->groupfillinit();
92
93 // Obtain Electrons, Muons and Jets
94 Particles elecs = apply<ParticleFinder>(event, "Electrons").particlesByPt(Cuts::pT > 10*GeV && Cuts::abseta < 2.47);
95 Particles muons = apply<ParticleFinder>(event, "Muons").particlesByPt(Cuts::pT > 10*GeV && Cuts::abseta < 2.4);
96 Particles leptons = sortByPt(elecs + muons);
97 Jets jets = apply<SmearedJets>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.4);
98
99 // Discard jets within DR = 0.4 of prompt leptons
100 idiscardIfAnyDeltaRLess(jets, leptons, 0.4);
101
102 // 2-lepton High (n=0), Int (n=1) and Low (n=2) Selection
103 for (size_t n=0; n<_types.size(); ++n) {
104
105 while (true) {
106
107 // Require the leading electron and the leading muon to have pT > 25GeV
108 if (elecs.size() != 0 && elecs[0].pT() < 25*GeV) break;
109 if (muons.size() != 0 && muons[0].pT() < 25*GeV) break;
110
111 // Obtain the missing transverse momentum vector
112 Vector3 EtMissX = apply<SmearedMET>(event,"MET").vectorPt();
113 Vector3 EtMiss = -EtMissX;
114
115 // Require only 2 opposite charged, same flavoured leptons
116 double n_leptons = leptons.size();
117 if (n_leptons != 2) break;
118 if (leptons[0].abspid() != leptons[1].abspid()) break;
119 if (leptons[0].charge() * leptons[1].charge() >= 0) break;
120
121 // Obtain the 4-momenta of leptons and implement requirements on them
122 vector<FourMomentum> lepton;
123 for (int l = 0; l < 2; ++l) lepton.push_back(leptons[l].mom());
124
125 double p_l1_T = lepton[0].pT();
126 if (p_l1_T < 25*GeV) break;
127
128 double p_l2_T = lepton[1].pT();
129 if (p_l2_T < 25*GeV) break;
130
131 _cutflows->fillnext("2l"s + _types[n]);
132
133 // Requirement on m_ll
134 double m_ll = (lepton[0]+lepton[1]).mass();
135 if (!inRange(m_ll, 80*GeV, 100*GeV)) break;
136
137 // Obtain the 4-momenta of jets and implement requirements on them
138 double n_jets = jets.size();
139 if (n != 2) {
140 if (n_jets < 2) break;
141 if (any(jets,hasBTag())) break;
142 } else{
143 if (n_jets != 2) break;
144 if (any(jets,hasBTag())) break;
145 }
146
147 vector<FourMomentum> jet;
148 for (int l = 0; l < 2; ++l) jet.push_back(jets[l].momentum());
149
150 double p_j1_T = jet[0].pT();
151 if (p_j1_T < 30*GeV) break;
152
153 double p_j2_T = jet[1].pT();
154 if (p_j2_T < 30*GeV) break;
155
156 double m_jj = (jet[0]+jet[1]).mass();
157 if (n != 2) {
158 if (m_jj < 60*GeV || m_jj > 100*GeV) break;
159 } else{
160 if (m_jj < 70*GeV || m_jj > 90*GeV) break;
161 }
162
163 // Pre-selection requirements cut
164 _cutflows->fillnext("2l"s+_types[n]);
165
166 // Invisible mass JR and Invisible Rapidity JR to obtain Invisible System 4-momentum
167 FourMomentum P_V = lepton[0] + lepton[1] + jet[0] + jet[1];
168 double M_I = sqrt(P_V.mass2() - 4*m_ll*m_jj);
169 double M_V = P_V.mass();
170 double p_I_z = P_V.pz() * sqrt(EtMiss.dot(EtMiss) + sqr(M_I)) / sqrt(P_V.pT2() + sqr(M_V));
171
172 FourMomentum P_I;
173 P_I.setPM(EtMiss.x(), EtMiss.y(), p_I_z, M_I);
174
175 // Lorentz boost to CM frame
176 LorentzTransform LT = LorentzTransform::mkFrameTransform(P_I+P_V);
177
178 FourMomentum P_F_I; //4-momentum of invisible system in CM frame
179 FourMomentum P_F_Va; //4-momentum of visible system a in CM frame
180 FourMomentum P_F_Vb; //4-momentum of visible system b in CM frame
181 FourMomentum P_F_V; //4-momentum of visible system in CM frame
182 FourMomentum P_F_Va1; //4-momentum of visible system a1 in CM frame
183 FourMomentum P_F_Va2; //4-momentum of visible system a2 in CM frame
184 FourMomentum P_F_Vb1; //4-momentum of visible system b1 in CM frame
185 FourMomentum P_F_Vb2; //4-momentum of visible system b2 in CM frame
186 if ((jet[0]+jet[1]).mass() > (lepton[0]+lepton[1]).mass()) {
187 P_F_I = LT.transform(P_I);
188 P_F_Va = LT.transform(jet[0]+jet[1]);
189 P_F_Vb = LT.transform(lepton[0]+lepton[1]);
190 P_F_V = LT.transform(P_V);
191 //
192 P_F_Va1 = LT.transform(jet[0]);
193 P_F_Va2 = LT.transform(jet[1]);
194 P_F_Vb1 = LT.transform(lepton[0]);
195 P_F_Vb2 = LT.transform(lepton[1]);
196 } else {
197 P_F_I = LT.transform(P_I);
198 P_F_Va = LT.transform(lepton[0]+lepton[1]);
199 P_F_Vb = LT.transform(jet[0]+jet[1]);
200 P_F_V = LT.transform(P_V);
201 //
202 P_F_Va1 = LT.transform(lepton[0]);
203 P_F_Va2 = LT.transform(lepton[1]);
204 P_F_Vb1 = LT.transform(jet[0]);
205 P_F_Vb2 = LT.transform(jet[1]);
206 }
207
208 // Obtain variables defined in the Contra-boost Invariant JR
209 double M2c = 2*((P_F_Va.E())*(P_F_Vb.E())+(P_F_Va.p3()).dot(P_F_Vb.p3()));
210 double m2Va = sqr(P_F_Va.mass());
211 double m2Vb = sqr(P_F_Vb.mass());
212 double ka = m2Va-m2Vb+M2c-2*sqrt(m2Va)*sqrt(m2Vb);
213 double kb = m2Vb-m2Va+M2c-2*sqrt(m2Va)*sqrt(m2Vb);
214 double kn = ka*m2Va - kb*m2Vb + M2c*(kb-ka)/2 + (0.5)*sqrt(pow(ka+kb,2)*(pow(M2c,2)-4*m2Va*m2Vb));
215 double k2d = sqr(ka)*m2Va+sqr(kb)*m2Vb+ka*kb*M2c;
216 double k = kn/k2d;
217 double ca = (1+k*ka)/2;
218 double cb = (1+k*kb)/2;
219 double c = 0.5*(P_F_V.E()+sqrt(pow(P_F_V.E(),2)+pow(M_I,2)-pow(P_V.mass(),2)))/(ca*(P_F_Va.E())+cb*(P_F_Vb.E()));
220
221 // Apply Contra-boost Invariant JR to obtain invisible particles' 4-momenta
222 double p_F_Iax = (P_F_Va.px())*(c*ca-1)-(P_F_Vb.px())*c*cb;
223 double p_F_Iay = (P_F_Va.py())*(c*ca-1)-(P_F_Vb.py())*c*cb;
224 double p_F_Iaz = (P_F_Va.pz())*(c*ca-1)-(P_F_Vb.pz())*c*cb;
225 double p_F_Ibx = (P_F_Vb.px())*(c*cb-1)-(P_F_Va.px())*c*ca;
226 double p_F_Iby = (P_F_Vb.py())*(c*cb-1)-(P_F_Va.py())*c*ca;
227 double p_F_Ibz = (P_F_Vb.pz())*(c*cb-1)-(P_F_Va.pz())*c*ca;
228 double E_F_Ia = (c*ca-1)*(P_F_Va.E())+c*cb*(P_F_Vb.E());
229 double E_F_Ib = (c*cb-1)*(P_F_Vb.E())+c*ca*(P_F_Va.E());
230
231 FourMomentum P_F_Ia;
232 FourMomentum P_F_Ib;
233 P_F_Ia.setPE(p_F_Iax,p_F_Iay,p_F_Iaz,E_F_Ia);
234 P_F_Ib.setPE(p_F_Ibx,p_F_Iby,p_F_Ibz,E_F_Ib);
235
236 // Lorentz boost from the CM frame to the P_a and P_b frame
237 LorentzTransform LTPa=LorentzTransform::mkFrameTransform(P_F_Va+P_F_Ia);
238 LorentzTransform LTPb=LorentzTransform::mkFrameTransform(P_F_Vb+P_F_Ib);
239
240 // min(H_Pa_11,H_Pb_11)/min(H_Pa_21,H_Pb,21) requirement
241 if (n != 2) {
242 double H_Pa_11 = (LTPa.transform(P_F_Va1+P_F_Va2)).p()+(LTPa.transform(P_F_Ia)).p();
243 double H_Pb_11 = (LTPb.transform(P_F_Vb1+P_F_Vb2)).p()+(LTPb.transform(P_F_Ib)).p();
244 double H_Pa_21 = (LTPa.transform(P_F_Va1)).p()+(LTPa.transform(P_F_Va2)).p()+(LTPa.transform(P_F_Ia)).p();
245 double H_Pb_21 = (LTPa.transform(P_F_Vb1)).p()+(LTPb.transform(P_F_Vb2)).p()+(LTPb.transform(P_F_Ib)).p();
246 vector<double> V1 = {H_Pa_11,H_Pb_11};
247 vector<double> V2 = {H_Pa_21,H_Pb_21};
248 double fraction1 = min(V1)/min(V2);
249 if (fraction1 < 0.8) break;
250 } else {
251 double H_PP_41 = P_F_Va1.p()+P_F_Va2.p()+P_F_Vb1.p()+P_F_Vb2.p()+(P_F_Ia+P_F_Ib).p();
252 double H_PP_11 = (P_F_Va1+P_F_Va2+P_F_Vb1+P_F_Vb2).p()+(P_F_Ia+P_F_Ib).p();
253 double fraction1 = H_PP_11/H_PP_41;
254 if (fraction1 < 0.35 || fraction1 > 0.6) break;
255 }
256 _cutflows->fillnext("2l"s+_types[n]);
257
258 // Lorentz boost from the CM frame to the lab frame
259 LorentzTransform LTR = LorentzTransform::mkObjTransform(P_I+P_V);
260
261 // p_lab_T_PP/(p_lab_T_PP+p_lab_T_41) requirement
262 double p_lab_T_PP = (LTR.transform(P_F_Va1+P_F_Va2+P_F_Vb1+P_F_Vb2+P_F_Ia+P_F_Ib)).pT();
263 double H_PP_T_41 = P_F_Va1.pT()+P_F_Va2.pT()+P_F_Vb1.pT()+P_F_Vb2.pT()+(P_F_Ia+P_F_Ib).pT();
264 double fraction2 = p_lab_T_PP/(p_lab_T_PP+H_PP_T_41);
265 if (fraction2 > 0.05) break;
266 _cutflows->fillnext("2l"s+_types[n]);
267
268 // Delta-Phi requirement
269 if (n == 0) {
270 Vector3 vector1 = (P_F_Va+P_F_Ia).betaVec();
271 Vector3 vector2 = (LTPa.transform(P_F_Va)).betaVec();
272 if (deltaPhi(vector1,vector2) < 0.3 || deltaPhi(vector1,vector2) > 2.8) break;
273 Vector3 vector3 = (P_F_Vb+P_F_Ib).betaVec();
274 Vector3 vector4 = (LTPb.transform(P_F_Vb)).betaVec();
275 if (deltaPhi(vector3,vector4) < 0.3 || deltaPhi(vector3,vector4) > 2.8) break;
276 } else if (n == 1) {
277 Vector3 vector1 = (P_F_Va+P_F_Ia).betaVec();
278 Vector3 vector2 = (LTPa.transform(P_F_Va)).betaVec();
279 if (deltaPhi(vector1,vector2) < 0.6 || deltaPhi(vector1,vector2) > 2.6) break;
280 Vector3 vector3 = (P_F_Vb+P_F_Ib).betaVec();
281 Vector3 vector4 = (LTPb.transform(P_F_Vb)).betaVec();
282 if (deltaPhi(vector3,vector4) < 0.6 || deltaPhi(vector3,vector4) > 2.6) break;
283 } else {
284 Vector3 j1 = jet[0].p3();
285 Vector3 j2 = jet[1].p3();
286 double delta1 = deltaPhi(j1,EtMiss);
287 double delta2 = deltaPhi(j2,EtMiss);
288 double delta = min(delta1,delta2);
289 if (delta < 2.4) break;
290 }
291 _cutflows->fillnext("2l"s+_types[n]);
292
293 // H_PP_41 requirement
294 double H_PP_41 = P_F_Va1.p()+P_F_Va2.p()+P_F_Vb1.p()+P_F_Vb2.p()+(P_F_Ia+P_F_Ib).p();
295 if (n==0) {
296 if (H_PP_41 < 800*GeV) break;
297 } else if (n==1) {
298 if (H_PP_41 < 600*GeV) break;
299 } else {
300 if (H_PP_41 < 400*GeV) break;
301 }
302 _cutflows->fillnext("2l"s+_types[n]);
303
304 break;
305 }
306 }
307
308
309 // 2-lepton ISR Selection
310 while (true) {
311
312 // Require the leading electron and the leading muon to have pT > 25 GeV
313 if (elecs.size() != 0 && elecs[0].pT() < 25*GeV) break;
314 if (muons.size() != 0 && muons[0].pT() < 25*GeV) break;
315
316 // Require only 2 opposite charged, same flavoured leptons
317 double n_leptons = leptons.size();
318 if (n_leptons != 2) break;
319 if (leptons[0].abspid() != leptons[1].abspid()) break;
320 if (leptons[0].charge()*leptons[1].charge() >= 0) break;
321
322 // Obtain the 4-momenta of leptons and implement requirements on them
323 vector<FourMomentum> lepton;
324 double M;
325 for (int n = 0; n < 2; ++n) {
326 lepton.push_back(leptons[n].mom());
327 M = lepton[n].mass();
328 lepton[n].setPz(0);
329 lepton[n].setE(sqrt(pow(lepton[n].px(),2)+pow(lepton[n].py(),2)+pow(M,2)));
330 }
331 if (lepton[0].pT() < 25*GeV || lepton[1].pT() < 25*GeV) break;
332
333 _cutflows->fillnext("2l-ISR");
334
335 // Obtain the 4-momenta of leptons and implement requirements on them
336 double n_jets = jets.size();
337 if (n_jets < 3 || n_jets > 4) break;
338 if (any(jets, hasBTag())) break;
339
340 vector<FourMomentum> jet;
341 for (int n = 0; n < n_jets; ++n) {
342 jet.push_back(jets[n].momentum());
343 M = jet[n].mass();
344 jet[n].setPz(0);
345 jet[n].setE(sqrt(pow(jet[n].px(),2)+pow(jet[n].py(),2)+pow(M,2)));
346 }
347 if (jet[0].pT() < 30*GeV || jet[1].pT() < 30*GeV) break;
348
349 Vector3 EtMissX = apply<SmearedMET>(event,"MET").vectorPt();
350 Vector3 EtMiss = -EtMissX;
351
352 // Obtain the missing transverse momentum vector
353 FourMomentum EtMiss1;
354 EtMiss1.setPM(EtMiss.x(),EtMiss.y(),EtMiss.z(),sqrt(EtMiss.dot(EtMiss)));
355
356 // Combinatoric Minimization JR
357 vector<double> f;
358 int indexToReturn = 0;
359 int indexValue = 0;
360 int newValue = 0;
361 if (n_jets == 3) {
362 f.push_back((EtMiss1+jet[1]+jet[2]).p());
363 f.push_back((EtMiss1+jet[0]+jet[2]).p());
364 f.push_back((EtMiss1+jet[0]+jet[1]).p());
365 f.push_back((EtMiss1+jet[0]).p());
366 f.push_back((EtMiss1+jet[1]).p());
367 f.push_back((EtMiss1+jet[2]).p());
368
369 for (int i = 0; i < 6; i++) {
370 newValue = f[i];
371 if (newValue >= indexValue) {
372 indexToReturn = i;
373 indexValue = newValue;
374 }
375 }
376
377 if (indexToReturn == 3 || indexToReturn == 4 || indexToReturn == 5) break;
378 }
379 else {
380 f.push_back((EtMiss1+jet[0]+jet[1]).p());
381 f.push_back((EtMiss1+jet[0]+jet[2]).p());
382 f.push_back((EtMiss1+jet[0]+jet[3]).p());
383 f.push_back((EtMiss1+jet[1]+jet[2]).p());
384 f.push_back((EtMiss1+jet[1]+jet[3]).p());
385 f.push_back((EtMiss1+jet[2]+jet[3]).p());
386 f.push_back((EtMiss1+jet[0]).p());
387 f.push_back((EtMiss1+jet[1]).p());
388 f.push_back((EtMiss1+jet[2]).p());
389 f.push_back((EtMiss1+jet[3]).p());
390 f.push_back((EtMiss1+jet[1]+jet[2]+jet[3]).p());
391 f.push_back((EtMiss1+jet[0]+jet[2]+jet[3]).p());
392 f.push_back((EtMiss1+jet[0]+jet[1]+jet[3]).p());
393 f.push_back((EtMiss1+jet[0]+jet[1]+jet[2]).p());
394
395 for (int i = 0; i < 14; i++) {
396 newValue = f[i];
397 if (newValue >= indexValue) {
398 indexToReturn = i;
399 indexValue = newValue;
400 }
401 }
402
403 if (indexToReturn == 6 || indexToReturn == 7 || indexToReturn == 8 || indexToReturn == 9 ||
404 indexToReturn == 10 || indexToReturn == 11 || indexToReturn == 12 || indexToReturn == 13) break;
405 }
406
407 _cutflows->fillnext("2l-ISR");
408
409 // g1 and g2 are the set of jets belonging to the ISR and signal system respectively
410 vector<vector<double>> g1, g2;
411 if (n_jets == 3) {
412 g2 = {{1,2},{0,2},{0,1}};
413 } else {
414 g1 = {{2,3},{1,3},{1,2},{0,3},{0,2},{0,1}};
415 g2 = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
416 }
417
418 // m_Z requirement
419 double m_Z = ((leptons[0]).mom()+(leptons[1]).mom()).mass();
420 if (m_Z < 80*GeV || m_Z > 100*GeV) break;
421 _cutflows->fillnext("2l-ISR");
422
423 // m_J requirement
424 double m_J = ((jets[g2[indexToReturn][0]]).momentum()+(jets[g2[indexToReturn][1]]).momentum()).mass();
425 if (m_J < 50*GeV || m_J > 110*GeV) break;
426 _cutflows->fillnext("2l-ISR");
427
428 // Compute variables delta_phi, R_ISR, P_T_ISR, p_T_I, p_T and implement their requirements
429 double p_T_ISR;
430 double p_T_I;
431 double p_T;
432 double R_ISR;
433 double delta_phi;
434 if (n_jets == 3) {
435 FourMomentum CM = EtMiss1+jet[indexToReturn]+jet[g2[indexToReturn][0]]+jet[g2[indexToReturn][1]]+lepton[0]+lepton[1];
436 LorentzTransform LT=LorentzTransform::mkFrameTransform(CM);
437 p_T_ISR = (LT.transform(jet[indexToReturn])).pT();
438 p_T_I = (LT.transform(EtMiss1)).pT();
439 p_T = (CM).pT();
440
441 Vector3 p_I = (LT.transform(EtMiss1)).p3();
442 Vector3 p_S = (LT.transform(EtMiss1+jet[g2[indexToReturn][0]]+jet[g2[indexToReturn][1]]+lepton[0]+lepton[1])).p3();
443 R_ISR = (p_S.dot(p_I))/(p_S.dot(p_S));
444 delta_phi = deltaPhi((LT.transform(EtMiss1)).p3(),(LT.transform(jet[indexToReturn])).p3());
445
446 } else {
447
448 FourMomentum CM = EtMiss1+jet[g1[indexToReturn][0]]+jet[g1[indexToReturn][1]]+jet[g2[indexToReturn][0]]+jet[g2[indexToReturn][1]]+lepton[0]+lepton[1];
449 LorentzTransform LT=LorentzTransform::mkFrameTransform(CM);
450 p_T_ISR = (LT.transform(jet[g1[indexToReturn][0]]+jet[g1[indexToReturn][1]])).pT();
451 p_T_I = (LT.transform(EtMiss1)).pT();
452 p_T = (CM).pT();
453
454 Vector3 p_I = (LT.transform(EtMiss1)).p3();
455 Vector3 p_S = (LT.transform(EtMiss1+jet[g2[indexToReturn][0]]+jet[g2[indexToReturn][1]]+lepton[0]+lepton[1])).p3();
456 R_ISR = (p_S.dot(p_I))/(p_S.dot(p_S));
457 delta_phi = deltaPhi((LT.transform(EtMiss1)).p3(),(LT.transform(jet[g1[indexToReturn][0]]+jet[g1[indexToReturn][1]])).p3());
458 }
459
460 if (delta_phi < 2.8) break;
461 _cutflows->fillnext("2l-ISR");
462
463 if (R_ISR < 0.4 || R_ISR > 0.75) break;
464 _cutflows->fillnext("2l-ISR");
465
466 if (p_T_ISR < 180*GeV) break;
467 _cutflows->fillnext("2l-ISR");
468
469 if (p_T_I < 100*GeV) break;
470 _cutflows->fillnext("2l-ISR");
471
472 if (p_T > 20*GeV) break;
473 _cutflows->fillnext("2l-ISR");
474
475 break;
476 }
477
478
479 // 3-lepton High (n=0), Int (n=1) and Low (n=2) Selection
480 for (int n = 0; n < 3; ++n) {
481
482 while (true) {
483
484 // Require the leading electron and the leading muon to have pT > 25 GeV
485 if (elecs.size() != 0 && elecs[0].pT() < 25*GeV) break;
486 if (muons.size() != 0 && muons[0].pT() < 25*GeV) break;
487
488 // Require 3 leptons and a pair of leptons with opposite charge and same flavour
489 double n_leptons = leptons.size();
490 if (n_leptons != 3) break;
491
492 vector<vector<double>> g = {{1,2},{0,2},{0,1}};
493 int indexToReturn = 3;
494 int indexValue = 100000;
495 int newValue = 0;
496 for (int i = 0; i < 3; i++) {
497 if (!(leptons[g[i][0]].abspid() == leptons[g[i][1]].abspid() && leptons[g[i][0]].charge() != leptons[g[i][1]].charge())) continue;
498 newValue = abs((leptons[g[i][0]].mom()+leptons[g[i][1]].mom()).mass() - 90*GeV);
499 if (newValue <= indexValue) {
500 indexToReturn = i;
501 indexValue = newValue;
502 }
503 }
504
505 if (indexToReturn == 3) break;
506
507 _cutflows->fillnext("3l"s+_types[n]);
508
509 // Requirements on n_jets
510 double n_jets = jets.size();
511 if (n != 2) {
512 if (n_jets >= 3) break;
513 if (any(jets, hasBTag())) break;
514 } else {
515 if (n_jets != 0) break;
516 }
517
518 // Obtain 4-momentum of the leptons and implement the requirements on their transverse momentum
519 vector<FourMomentum> lepton;
520 for (int l = 0; l < 3; ++l) {
521 lepton.push_back(leptons[l].mom());
522 }
523
524 double p_l1_T = lepton[0].pT();
525 if (p_l1_T < 60) break;
526
527 if (n==0) {
528 double p_l2_T = lepton[1].pT();
529 if (p_l2_T < 60*GeV) break;
530 } else if (n==1) {
531 double p_l2_T = lepton[1].pT();
532 if (p_l2_T < 50*GeV) break;
533 } else {
534 double p_l2_T = lepton[1].pT();
535 if (p_l2_T < 40*GeV) break;
536 }
537
538 if (n==0) {
539 double p_l3_T = (lepton[2]).pT();
540 if (p_l3_T < 40*GeV) break;
541 } else {
542 double p_l3_T = (lepton[2]).pT();
543 if (p_l3_T < 30*GeV) break;
544 }
545
546 // Pre-selection Cut
547 _cutflows->fillnext("3l"s+_types[n]);
548
549 // Requirement on m_ll
550 double m_ll = (lepton[g[indexToReturn][0]]+lepton[g[indexToReturn][1]]).mass();
551 if (m_ll < 75*GeV || m_ll > 105*GeV) break;
552
553 _cutflows->fillnext("3l"s+_types[n]);
554
555 // Obtain the missing transverse momentum vector
556 Vector3 EtMissX = apply<SmearedMET>(event,"MET").vectorPt();
557 Vector3 EtMiss = -EtMissX;
558
559 // Requirement on m_W_T
560 double deltaphi = deltaPhi(EtMiss,(lepton[indexToReturn]).p3());
561 double m_W_T = sqrt(2*((lepton[indexToReturn]).pT())*sqrt(EtMiss.dot(EtMiss))*(1-cos(deltaphi)));
562
563 if (n==0) {
564 if (m_W_T < 150*GeV) break;
565 } else if (n==1) {
566 if (m_W_T < 130*GeV) break;
567 } else {
568 if (m_W_T < 100*GeV) break;
569 }
570
571 _cutflows->fillnext("3l"s+_types[n]);
572
573 // Invisible mass JR and Invisible Rapidity JR to obtain Invisible System 4-momentum
574 FourMomentum P_V = lepton[0]+lepton[1]+lepton[2];
575 double M_I = sqrt(P_V.mass2() - 4*m_ll*((lepton[indexToReturn]).mass()));
576 double M_V = P_V.mass();
577 double p_I_z = (P_V.pz())*sqrt(EtMiss.dot(EtMiss)+sqr(M_I)) / sqrt(P_V.pT2() + sqr(M_V));
578
579 FourMomentum P_I;
580 P_I.setPM(EtMiss.x(), EtMiss.y(), p_I_z, M_I);
581
582 // Lorentz boost to CM frame
583 LorentzTransform LT = LorentzTransform::mkFrameTransform(P_I+P_V);
584 FourMomentum P_F_I; //4-momentum of invisible system in CM frame
585 FourMomentum P_F_Va; //4-momentum of visible system a in CM frame
586 FourMomentum P_F_Vb; //4-momentum of visible system b in CM frame
587 FourMomentum P_F_V; //4-momentum of visible system in CM frame
588 FourMomentum P_F_Va1; //4-momentum of visible system a1 in CM frame
589 FourMomentum P_F_Va2; //4-momentum of visible system a2 in CM frame
590 FourMomentum P_F_Vb1; //4-momentum of visible system b1 in CM frame
591 FourMomentum P_F_Vb2; //4-momentum of visible system b2 in CM frame
592 if ((lepton[indexToReturn]).mass() > (lepton[g[indexToReturn][0]]+lepton[g[indexToReturn][1]]).mass()) {
593 P_F_I = LT.transform(P_I);
594 P_F_Va = LT.transform(lepton[indexToReturn]);
595 P_F_Vb = LT.transform(lepton[g[indexToReturn][0]]+lepton[g[indexToReturn][1]]);
596 P_F_V = LT.transform(P_V);
597 //
598 P_F_Va1 = LT.transform(lepton[indexToReturn]);
599 P_F_Vb1 = LT.transform(lepton[g[indexToReturn][0]]);
600 P_F_Vb2 = LT.transform(lepton[g[indexToReturn][1]]);
601
602 } else {
603
604 P_F_I = LT.transform(P_I);
605 P_F_Va = LT.transform(lepton[g[indexToReturn][0]]+lepton[g[indexToReturn][1]]);
606 P_F_Vb = LT.transform(lepton[indexToReturn]);
607 P_F_V = LT.transform(P_V);
608 //
609 P_F_Va1 = LT.transform(lepton[g[indexToReturn][0]]);
610 P_F_Va2 = LT.transform(lepton[g[indexToReturn][1]]);
611 P_F_Vb1 = LT.transform(lepton[indexToReturn]);
612 }
613
614 // Obtain variables defined in the Contra-boost Invariant JR
615 double M2c = 2*((P_F_Va.E())*(P_F_Vb.E())+(P_F_Va.p3()).dot(P_F_Vb.p3()));
616 double m2Va = sqr(P_F_Va.mass());
617 double m2Vb = sqr(P_F_Vb.mass());
618 double ka = m2Va-m2Vb+M2c-2*sqrt(m2Va)*sqrt(m2Vb);
619 double kb = m2Vb-m2Va+M2c-2*sqrt(m2Va)*sqrt(m2Vb);
620 double kn = ka*m2Va - kb*m2Vb + M2c*(kb-ka)/2 + (0.5)*sqrt(pow(ka+kb,2)*(pow(M2c,2)-4*m2Va*m2Vb));
621 double k2d = sqr(ka)*m2Va+sqr(kb)*m2Vb+ka*kb*M2c;
622 double k = kn/k2d;
623 double ca = (1+k*ka)/2;
624 double cb = (1+k*kb)/2;
625 double c = 0.5*(P_F_V.E()+sqrt(pow(P_F_V.E(),2)+pow(M_I,2)-pow(P_V.mass(),2)))/(ca*(P_F_Va.E())+cb*(P_F_Vb.E()));
626
627 // Apply Contra-boost Invariant JR to obtain invisible particles' 4-momenta
628 double p_F_Iax = (P_F_Va.px())*(c*ca-1)-(P_F_Vb.px())*c*cb;
629 double p_F_Iay = (P_F_Va.py())*(c*ca-1)-(P_F_Vb.py())*c*cb;
630 double p_F_Iaz = (P_F_Va.pz())*(c*ca-1)-(P_F_Vb.pz())*c*cb;
631 double p_F_Ibx = (P_F_Vb.px())*(c*cb-1)-(P_F_Va.px())*c*ca;
632 double p_F_Iby = (P_F_Vb.py())*(c*cb-1)-(P_F_Va.py())*c*ca;
633 double p_F_Ibz = (P_F_Vb.pz())*(c*cb-1)-(P_F_Va.pz())*c*ca;
634 double E_F_Ia = (c*ca-1)*(P_F_Va.E())+c*cb*(P_F_Vb.E());
635 double E_F_Ib = (c*cb-1)*(P_F_Vb.E())+c*ca*(P_F_Va.E());
636
637 FourMomentum P_F_Ia;
638 FourMomentum P_F_Ib;
639 P_F_Ia.setPE(p_F_Iax,p_F_Iay,p_F_Iaz,E_F_Ia);
640 P_F_Ib.setPE(p_F_Ibx,p_F_Iby,p_F_Ibz,E_F_Ib);
641
642 // Lorentz Transform to the Pb frame from the CM frame
643 LorentzTransform LTP;
644 if ((lepton[indexToReturn]).mass() > (lepton[g[indexToReturn][0]]+lepton[g[indexToReturn][1]]).mass()) {
645 LTP = LorentzTransform::mkFrameTransform(P_F_Vb+P_F_Ib);
646 } else {
647 LTP = LorentzTransform::mkFrameTransform(P_F_Va+P_F_Ia);
648 }
649
650 // p_lab_T_PP/(p_lab_T_PP+p_lab_T_31) requirement
651 LorentzTransform LTR = LorentzTransform::mkObjTransform(P_I+P_V);
652 double p_lab_T_PP = (LTR.transform(P_F_Va+P_F_Vb+P_F_Ia+P_F_Ib)).pT();
653 double H_PP_T_31;
654 double H_PP_31;
655 double H_Pb_11;
656 double H_Pb_21;
657 if ((lepton[indexToReturn]).mass() > (lepton[g[indexToReturn][0]]+lepton[g[indexToReturn][1]]).mass()) {
658 H_PP_T_31 = P_F_Va1.pT()+P_F_Vb1.pT()+P_F_Vb2.pT()+(P_F_Ia+P_F_Ib).pT();
659 H_PP_31 = P_F_Va1.p()+P_F_Vb1.p()+P_F_Vb2.p()+(P_F_Ia+P_F_Ib).p();
660 H_Pb_11 = (LTP.transform(P_F_Vb1+P_F_Vb2)).p()+(LTP.transform(P_F_Ib)).p();
661 H_Pb_21 = (LTP.transform(P_F_Vb1)).p()+(LTP.transform(P_F_Vb2)).p()+(LTP.transform(P_F_Ib)).p();
662 } else {
663 H_PP_T_31 = P_F_Va1.pT()+P_F_Va2.pT()+P_F_Vb1.pT()+(P_F_Ia+P_F_Ib).pT();
664 H_PP_31 = P_F_Va1.p()+P_F_Va2.p()+P_F_Vb1.p()+(P_F_Ia+P_F_Ib).p();
665 H_Pb_11 = (LTP.transform(P_F_Va1+P_F_Va2)).p()+(LTP.transform(P_F_Ia)).p();
666 H_Pb_21 = (LTP.transform(P_F_Va1)).p()+(LTP.transform(P_F_Va2)).p()+(LTP.transform(P_F_Ia)).p();
667 }
668
669 double fraction1 = H_PP_T_31/H_PP_31;
670 if (n==0) {
671 if (fraction1 < 0.75) break;
672 } else if (n==1) {
673 if (fraction1 < 0.8) break;
674 } else {
675 if (fraction1 < 0.9) break;
676 }
677 _cutflows->fillnext("3l"s+_types[n]);
678
679 double fraction2 = H_Pb_11/H_Pb_21;
680 if (n==0) {
681 if (fraction2 < 0.8) break;
682 _cutflows->fillnext("3l"s+_types[n]);
683 } else if (n==1) {
684 if (fraction2 < 0.75) break;
685 _cutflows->fillnext("3l"s+_types[n]);
686 } else { /* ??? */ }
687
688 if (n==0) {
689 if (H_PP_31 < 550*GeV) break;
690 _cutflows->fillnext("3l"s+_types[n]);
691 } else if (n==1) {
692 if (H_PP_31 < 450*GeV) break;
693 _cutflows->fillnext("3l"s+_types[n]);
694 } else {
695 if (H_PP_31 < 250*GeV) break;
696 _cutflows->fillnext("3l"s+_types[n]);
697 }
698
699 double fraction3 = p_lab_T_PP/(p_lab_T_PP+H_PP_T_31);
700 if (n==0) {
701 if (fraction3 > 0.2) break;
702 _cutflows->fillnext("3l"s+_types[n]);
703 } else if (n==1) {
704 if (fraction3 > 0.15) break;
705 _cutflows->fillnext("3l"s+_types[n]);
706 } else {
707 if (fraction3 > 0.05) break;
708 _cutflows->fillnext("3l"s+_types[n]);
709 }
710
711 break;
712 }
713 }
714
715
716 // 3-lepton ISR Selection
717 while (true) {
718
719 // Require the leading electron and the leading muon to have pT > 25GeV
720 if (elecs.size() != 0 && elecs[0].pT() < 25*GeV) break;
721 if (muons.size() != 0 && muons[0].pT() < 25*GeV) break;
722
723 // Require only 3 leptons
724 double n_leptons = leptons.size();
725 if (n_leptons != 3) break;
726
727 // Require and identify the pair of opposite-charged and same-flavoured leptons
728 vector<vector<double>> g = {{1,2},{0,2},{0,1}};
729 int indexToReturn = 3;
730 int indexValue = 100000;
731 int newValue = 0;
732 for (int i = 0; i < 3; i++) {
733 if (!(leptons[g[i][0]].abspid() == leptons[g[i][1]].abspid() && leptons[g[i][0]].charge() != leptons[g[i][1]].charge())) continue;
734 newValue = abs(((leptons[g[i][0]]).mom()+(leptons[g[i][1]]).mom()).mass()-90);
735 if (newValue <= indexValue) {
736 indexToReturn = i;
737 indexValue = newValue;
738 }
739 }
740
741 if (indexToReturn == 3) break;
742 _cutflows->fillnext("3l-ISR");
743
744 // Requirements on jet number and forbid B-tags
745 double n_jets=jets.size();
746 if (n_jets < 1 || n_jets > 3) break;
747 if (any(jets, hasBTag())) break;
748
749 // Obtain the 4-momenta of leptons and implement requirements on them
750 vector<FourMomentum> lepton;
751 double M;
752 for (int n = 0; n < 3; ++n) {
753 lepton.push_back(leptons[n].mom());
754 M = lepton[n].mass();
755 lepton[n].setPz(0);
756 lepton[n].setE(sqrt(pow(lepton[n].px(),2)+pow(lepton[n].py(),2)+pow(M,2)));
757 }
758 if (lepton[0].pT() < 25*GeV || lepton[1].pT() < 25*GeV || lepton[2].pT() < 20*GeV) break;
759
760 // Pre-selection Cut
761 _cutflows->fillnext("3l-ISR");
762
763 // Obtain the missing momentum 3-vector
764 Vector3 EtMiss = apply<SmearedMET>(event,"MET").vectorMissingPt();
765
766 // Requirement on m_ll
767 double m_ll = (leptons[g[indexToReturn][0]].mom()+leptons[g[indexToReturn][1]].mom()).mass();
768 if (m_ll < 75*GeV || m_ll > 105*GeV) break;
769 _cutflows->fillnext("3l-ISR");
770
771 // Requirement on m_W_T
772 double deltaphi = deltaPhi(EtMiss,(lepton[indexToReturn]).p3());
773 double m_W_T = sqrt(2*((lepton[indexToReturn]).pT())*sqrt(EtMiss.dot(EtMiss))*(1-cos(deltaphi)));
774 if (m_W_T < 100*GeV) break;
775 _cutflows->fillnext("3l-ISR");
776
777 // Obtain the missing transverse 4-momentum
778 FourMomentum EtMiss1;
779 EtMiss1.setPM(EtMiss.x(),EtMiss.y(),EtMiss.z(),sqrt(EtMiss.dot(EtMiss)));
780 double p_T_ISR;
781 double p_T_I;
782 double p_T;
783 double R_ISR;
784 double delta_phi;
785 FourMomentum ISR;
786
787 // Obtain the 4-momentum of the ISR system
788 vector<FourMomentum> jet;
789 for (int n = 0; n < n_jets; ++n) {
790 jet.push_back(jets[n].momentum());
791 M = jet[n].mass();
792 jet[n].setPz(0);
793 jet[n].setE(sqrt(pow(jet[n].px(),2)+pow(jet[n].py(),2)+pow(M,2)));
794 }
795
796 if (n_jets==1) {
797 ISR = jet[0];
798 } else if (n_jets==2) {
799 ISR = jet[0]+jet[1];
800 } else {
801 ISR = jet[0]+jet[1]+jet[2];
802 }
803
804 // Four Momentum of the CM System
805 FourMomentum CM = EtMiss1+ISR+lepton[0]+lepton[1]+lepton[2];
806 LorentzTransform LT = LorentzTransform::mkFrameTransform(CM);
807 p_T_ISR = (LT.transform(ISR)).pT();
808 p_T_I = (LT.transform(EtMiss1)).pT();
809 p_T = (CM).pT();
810
811 Vector3 p_I = (LT.transform(EtMiss1)).p3();
812 Vector3 p_S = (LT.transform(EtMiss1+lepton[0]+lepton[1]+lepton[2])).p3();
813 R_ISR = (p_S.dot(p_I))/(p_S.dot(p_S));
814 delta_phi = angle((LT.transform(EtMiss1)).p3(),(LT.transform(ISR)).p3());
815
816 if (delta_phi < 2.0) break;
817 _cutflows->fillnext("3l-ISR");
818
819 if (R_ISR < 0.55 || R_ISR > 1.0) break;
820 _cutflows->fillnext("3l-ISR");
821
822 if (p_T_ISR < 100*GeV) break;
823 _cutflows->fillnext("3l-ISR");
824
825 if (p_T_I < 80*GeV) break;
826 _cutflows->fillnext("3l-ISR");
827
828 if (p_T > 25*GeV) break;
829 _cutflows->fillnext("3l-ISR");
830
831 break;
832 }
833
834 }
835
836
837 /// Finalise cutflow scaling etc.
838 void finalize() {
839
840 vector<double> scales{ 1673., 4369., 65247 };
841 for (size_t n=0; n < _types.size(); ++n) {
842 _cutflows->binAt("2l"s+_types[n])->normalizeFirst(scales[n]);
843 _cutflows->binAt("3l"s+_types[n])->normalizeFirst(scales[n]);
844 }
845 _cutflows->binAt("2l-ISR")->normalizeFirst(scales[2]);
846 _cutflows->binAt("3l-ISR")->normalizeFirst(scales[2]);
847
848 MSG_INFO("CUTFLOWS:\n\n" << _cutflows);
849 }
850
851 CutflowsPtr _cutflows;
852 const vector<string> _types{"-high", "-int", "-low"};
853
854 };
855
856
857 RIVET_DECLARE_PLUGIN(ATLAS_2018_I1676551);
858
859}
|