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