Rivet analyses referenceATLAS_2012_I1204447Inclusive multi-lepton searchExperiment: ATLAS (LHC) Inspire ID: 1204447 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
A generic search for anomalous production of events with at least three charged leptons is presented. The search uses a pp-collision data sample at a center-of-mass energy of $\sqrt{s}$ = 7 TeV corresponding to 4.6/fb of integrated luminosity collected in 2011 by the ATLAS detector at the CERN Large Hadron Collider. Events are required to contain at least two electrons or muons, while the third lepton may either be an additional electron or muon, or a hadronically decaying tau lepton. Events are categorized by the presence or absence of a reconstructed tau-lepton or Z-boson candidate decaying to leptons. No significant excess above backgrounds expected from Standard Model processes is observed. Results are presented as upper limits on event yields from non-Standard-Model processes producing at least three prompt, isolated leptons, given as functions of lower bounds on several kinematic variables. Fiducial efficiencies for model testing are also provided. This Rivet module implements the event selection and the fiducial efficiencies to test various models for their exclusion based on observed/excluded limits. Source code: ATLAS_2012_I1204447.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/VisibleFinalState.hh"
6#include "Rivet/Projections/VetoedFinalState.hh"
7#include "Rivet/Projections/IdentifiedFinalState.hh"
8#include "Rivet/Projections/UnstableParticles.hh"
9#include "Rivet/Projections/FastJets.hh"
10
11namespace Rivet {
12
13
14 class ATLAS_2012_I1204447 : public Analysis {
15 public:
16
17 /// Constructor
18 ATLAS_2012_I1204447()
19 : Analysis("ATLAS_2012_I1204447") { }
20
21
22 /// Book histograms and initialise projections before the run
23 void init() {
24
25 // To calculate the acceptance without having the fiducial lepton efficiencies included, this part can be turned off
26 _use_fiducial_lepton_efficiency = true;
27
28 // Random numbers for simulation of ATLAS detector reconstruction efficiency
29 srand(160385);
30
31 // Read in all signal regions
32 _signal_regions = getSignalRegions();
33
34 // Set number of events per signal region to 0
35 for (size_t i = 0; i < _signal_regions.size(); i++)
36 book(_eventCountsPerSR[_signal_regions[i]], "_eventCountsPerSR_" + _signal_regions[i]);
37
38 // Final state including all charged and neutral particles
39 const FinalState fs((Cuts::etaIn(-5.0, 5.0) && Cuts::pT >= 1*GeV));
40 declare(fs, "FS");
41
42 // Final state including all charged particles
43 declare(ChargedFinalState(Cuts::abseta < 2.5 && Cuts::pT > 1*GeV), "CFS");
44
45 // Final state including all visible particles (to calculate MET, Jets etc.)
46 declare(VisibleFinalState(Cuts::abseta < 5.0), "VFS");
47
48 // Final state including all AntiKt 04 Jets
49 VetoedFinalState vfs;
50 vfs.addVetoPairId(PID::MUON);
51 declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
52
53 // Final state including all unstable particles (including taus)
54 declare(UnstableParticles(Cuts::abseta < 5.0 && Cuts::pT > 5*GeV), "UFS");
55
56 // Final state including all electrons
57 IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*GeV);
58 elecs.acceptIdPair(PID::ELECTRON);
59 declare(elecs, "elecs");
60
61 // Final state including all muons
62 IdentifiedFinalState muons(Cuts::abseta < 2.5 && Cuts::pT > 10*GeV);
63 muons.acceptIdPair(PID::MUON);
64 declare(muons, "muons");
65
66 // Book histograms
67 book(_h_HTlep_all ,"HTlep_all" , 30, 0, 1500);
68 book(_h_HTjets_all ,"HTjets_all", 30, 0, 1500);
69 book(_h_MET_all ,"MET_all" , 20, 0, 1000);
70 book(_h_Meff_all ,"Meff_all" , 30, 0, 3000);
71
72 book(_h_e_n ,"e_n" , 10, -0.5, 9.5);
73 book(_h_mu_n ,"mu_n" , 10, -0.5, 9.5);
74 book(_h_tau_n ,"tau_n", 10, -0.5, 9.5);
75
76 book(_h_pt_1_3l ,"pt_1_3l", 100, 0, 2000);
77 book(_h_pt_2_3l ,"pt_2_3l", 100, 0, 2000);
78 book(_h_pt_3_3l ,"pt_3_3l", 100, 0, 2000);
79 book(_h_pt_1_2ltau ,"pt_1_2ltau", 100, 0, 2000);
80 book(_h_pt_2_2ltau ,"pt_2_2ltau", 100, 0, 2000);
81 book(_h_pt_3_2ltau ,"pt_3_2ltau", 100, 0, 2000);
82
83 book(_h_excluded ,"excluded", 2, -0.5, 1.5);
84 }
85
86
87 /// Perform the per-event analysis
88 void analyze(const Event& event) {
89
90 // Muons
91 Particles muon_candidates;
92 const Particles charged_tracks = apply<ChargedFinalState>(event, "CFS").particles();
93 const Particles visible_particles = apply<VisibleFinalState>(event, "VFS").particles();
94 for (const Particle& mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt()) {
95 // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of muon itself)
96 double pTinCone = -mu.pT();
97 for (const Particle& track : charged_tracks) {
98 if (deltaR(mu.momentum(), track.momentum()) < 0.3)
99 pTinCone += track.pT();
100 }
101
102 // Calculate eTCone30 variable (pT of all visible particles within dR<0.3)
103 double eTinCone = 0.;
104 for (const Particle& visible_particle : visible_particles) {
105 if (visible_particle.abspid() != PID::MUON && inRange(deltaR(mu.momentum(), visible_particle.momentum()), 0.1, 0.3))
106 eTinCone += visible_particle.pT();
107 }
108
109 // Apply reconstruction efficiency and simulate reco
110 int muon_id = 13;
111 if ( mu.hasAncestorWith(Cuts::pid == 15) || mu.hasAncestorWith(Cuts::pid == -15)) muon_id = 14;
112 const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(muon_id, mu) : 1.0;
113 const bool keep_muon = rand()/static_cast<double>(RAND_MAX) <= eff;
114
115 // Keep muon if pTCone30/pT < 0.15 and eTCone30/pT < 0.2 and reconstructed
116 if (keep_muon && pTinCone/mu.pT() <= 0.15 && eTinCone/mu.pT() < 0.2)
117 muon_candidates.push_back(mu);
118 }
119
120
121 // Electrons
122 Particles electron_candidates;
123 for (const Particle& e : apply<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
124 // Neglect electrons in crack regions
125 if (inRange(e.abseta(), 1.37, 1.52)) continue;
126
127 // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of electron itself)
128 double pTinCone = -e.pT();
129 for (const Particle& track : charged_tracks) {
130 if (deltaR(e.momentum(), track.momentum()) < 0.3) pTinCone += track.pT();
131 }
132
133 // Calculate eTCone30 variable (pT of all visible particles (except muons) within dR<0.3)
134 double eTinCone = 0.;
135 for (const Particle& visible_particle : visible_particles) {
136 if (visible_particle.abspid() != PID::MUON && inRange(deltaR(e.momentum(), visible_particle.momentum()), 0.1, 0.3))
137 eTinCone += visible_particle.pT();
138 }
139
140 // Apply reconstruction efficiency and simulate reco
141 int elec_id = 11;
142 if (e.hasAncestorWith(Cuts::pid == 15) || e.hasAncestorWith(Cuts::pid == -15)) elec_id = 12;
143 const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(elec_id, e) : 1.0;
144 const bool keep_elec = rand()/static_cast<double>(RAND_MAX) <= eff;
145
146 // Keep electron if pTCone30/pT < 0.13 and eTCone30/pT < 0.2 and reconstructed
147 if (keep_elec && pTinCone/e.pT() <= 0.13 && eTinCone/e.pT() < 0.2)
148 electron_candidates.push_back(e);
149 }
150
151
152 // Taus
153 /// @todo This could benefit from a tau finder projection
154 Particles tau_candidates;
155 for (const Particle& tau : apply<UnstableParticles>(event, "UFS").particlesByPt()) {
156 // Only pick taus out of all unstable particles
157 if (tau.abspid() != PID::TAU) continue;
158
159 // Check that tau has decayed into daughter particles
160 /// @todo Huh? Unstable taus with no decay vtx? Can use Particle.isStable()? But why in this situation?
161 if (tau.genParticle()->end_vertex() == 0) continue;
162
163 // Calculate visible tau pT from pT of tau neutrino in tau decay for pT and |eta| cuts
164 FourMomentum daughter_tau_neutrino_momentum = get_tau_neutrino_mom(tau);
165 Particle tau_vis = tau;
166 tau_vis.setMomentum(tau.momentum()-daughter_tau_neutrino_momentum);
167 // keep only taus in certain eta region and above 15 GeV of visible tau pT
168 if ( tau_vis.pT() <= 15.0*GeV || tau_vis.abseta() > 2.5) continue;
169
170 // Get prong number (number of tracks) in tau decay and check if tau decays leptonically
171 unsigned int nprong = 0;
172 bool lep_decaying_tau = false;
173 get_prong_number(tau.genParticle(), nprong, lep_decaying_tau);
174
175 // Apply reconstruction efficiency
176 int tau_id = 15;
177 if (nprong == 1) tau_id = 15;
178 else if (nprong == 3) tau_id = 16;
179
180 // Get fiducial lepton efficiency simulate reco efficiency
181 const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(tau_id, tau_vis) : 1.0;
182 const bool keep_tau = rand()/static_cast<double>(RAND_MAX) <= eff;
183
184 // Keep tau if nprong = 1, it decays hadronically, and it's reconstructed by the detector
185 if ( !lep_decaying_tau && nprong == 1 && keep_tau) tau_candidates.push_back(tau_vis);
186 }
187
188
189 // Jets (all anti-kt R=0.4 jets with pT > 25 GeV and eta < 4.9)
190 Jets jet_candidates = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 4.9);
191
192 // ETmiss
193 Particles vfs_particles = apply<VisibleFinalState>(event, "VFS").particles();
194 FourMomentum pTmiss;
195 for (const Particle& p : vfs_particles) pTmiss -= p.momentum();
196 double eTmiss = pTmiss.pT()/GeV;
197
198
199 //------------------
200 // Overlap removal
201
202 // electron - electron
203 Particles electron_candidates_2;
204 for (size_t ie = 0; ie < electron_candidates.size(); ++ie) {
205 const Particle & e = electron_candidates[ie];
206 bool away = true;
207 // If electron pair within dR < 0.1: remove electron with lower pT
208 for (size_t ie2=0; ie2 < electron_candidates_2.size(); ++ie2) {
209 if ( deltaR( e.momentum(), electron_candidates_2[ie2].momentum()) < 0.1 ) {
210 away = false;
211 break;
212 }
213 }
214 // If isolated keep it
215 if ( away )
216 electron_candidates_2.push_back( e );
217 }
218 // jet - electron
219 Jets recon_jets;
220 for (const Jet& jet : jet_candidates) {
221 bool away = true;
222 // if jet within dR < 0.2 of electron: remove jet
223 for (const Particle& e : electron_candidates_2) {
224 if (deltaR(e.momentum(), jet.momentum()) < 0.2) {
225 away = false;
226 break;
227 }
228 }
229 // jet - tau
230 if (away) {
231 // If jet within dR < 0.2 of tau: remove jet
232 for (const Particle& tau : tau_candidates) {
233 if (deltaR(tau.momentum(), jet.momentum()) < 0.2) {
234 away = false;
235 break;
236 }
237 }
238 }
239 // If isolated keep it
240 if ( away )
241 recon_jets.push_back( jet );
242 }
243
244
245 // electron - jet
246 Particles recon_leptons, recon_e;
247 for (size_t ie = 0; ie < electron_candidates_2.size(); ++ie) {
248 const Particle& e = electron_candidates_2[ie];
249 // If electron within 0.2 < dR < 0.4 from any jets: remove electron
250 bool away = true;
251 for (const Jet& jet : recon_jets) {
252 if (deltaR(e.momentum(), jet.momentum()) < 0.4) {
253 away = false;
254 break;
255 }
256 }
257 // electron - muon
258 // if electron within dR < 0.1 of a muon: remove electron
259 if (away) {
260 for (const Particle& mu : muon_candidates) {
261 if (deltaR(mu.momentum(), e.momentum()) < 0.1) {
262 away = false;
263 break;
264 }
265 }
266 }
267 // If isolated keep it
268 if (away) {
269 recon_e += e;
270 recon_leptons += e;
271 }
272 }
273
274
275 // tau - electron
276 Particles recon_tau;
277 for ( const Particle& tau : tau_candidates ) {
278 bool away = true;
279 // If tau within dR < 0.2 of an electron: remove tau
280 for ( const Particle& e : recon_e ) {
281 if (deltaR( tau.momentum(), e.momentum()) < 0.2) {
282 away = false;
283 break;
284 }
285 }
286 // tau - muon
287 // If tau within dR < 0.2 of a muon: remove tau
288 if (away) {
289 for (const Particle& mu : muon_candidates) {
290 if (deltaR(tau.momentum(), mu.momentum()) < 0.2) {
291 away = false;
292 break;
293 }
294 }
295 }
296 // If isolated keep it
297 if (away) recon_tau.push_back( tau );
298 }
299
300 // Muon - jet isolation
301 Particles recon_mu, trigger_mu;
302 // If muon within dR < 0.4 of a jet, remove muon
303 for (const Particle& mu : muon_candidates) {
304 bool away = true;
305 for (const Jet& jet : recon_jets) {
306 if ( deltaR( mu.momentum(), jet.momentum()) < 0.4 ) {
307 away = false;
308 break;
309 }
310 }
311 if (away) {
312 recon_mu.push_back( mu );
313 recon_leptons.push_back( mu );
314 if (mu.abseta() < 2.4) trigger_mu.push_back( mu );
315 }
316 }
317
318 // End overlap removal
319 //------------------
320
321
322 // Jet cleaning
323 if (rand()/static_cast<double>(RAND_MAX) <= 0.42) {
324 for (const Jet& jet : recon_jets) {
325 const double eta = jet.rapidity();
326 const double phi = jet.azimuthalAngle(MINUSPI_PLUSPI);
327 if (jet.pT() > 25*GeV && inRange(eta, -0.1, 1.5) && inRange(phi, -0.9, -0.5)) vetoEvent;
328 }
329 }
330
331
332 // Post-isolation event cuts
333 // Require at least 3 charged tracks in event
334 if (charged_tracks.size() < 3) vetoEvent;
335
336 // And at least one e/mu passing trigger
337 if (!( !recon_e .empty() && recon_e[0] .pT() > 25*GeV) &&
338 !( !trigger_mu.empty() && trigger_mu[0].pT() > 25*GeV) ) {
339 MSG_DEBUG("Hardest lepton fails trigger");
340 vetoEvent;
341 }
342
343 // And only accept events with at least 2 electrons and muons and at least 3 leptons in total
344 if (recon_mu.size() + recon_e.size() + recon_tau.size() < 3 || recon_leptons.size() < 2) vetoEvent;
345
346 // Sort leptons by decreasing pT
347 sortByPt(recon_leptons);
348 sortByPt(recon_tau);
349
350 // Calculate HTlep, fill lepton pT histograms & store chosen combination of 3 leptons
351 double HTlep = 0.;
352 Particles chosen_leptons;
353 if ( recon_leptons.size() > 2 ) {
354 _h_pt_1_3l->fill(recon_leptons[0].perp()/GeV);
355 _h_pt_2_3l->fill(recon_leptons[1].perp()/GeV);
356 _h_pt_3_3l->fill(recon_leptons[2].perp()/GeV);
357 HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + recon_leptons[2].pT())/GeV;
358 chosen_leptons.push_back( recon_leptons[0] );
359 chosen_leptons.push_back( recon_leptons[1] );
360 chosen_leptons.push_back( recon_leptons[2] );
361 }
362 else {
363 _h_pt_1_2ltau->fill(recon_leptons[0].perp()/GeV);
364 _h_pt_2_2ltau->fill(recon_leptons[1].perp()/GeV);
365 _h_pt_3_2ltau->fill(recon_tau[0].perp()/GeV);
366 HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + recon_tau[0].pT())/GeV ;
367 chosen_leptons.push_back( recon_leptons[0] );
368 chosen_leptons.push_back( recon_leptons[1] );
369 chosen_leptons.push_back( recon_tau[0] );
370 }
371
372 // Number of prompt e/mu and had taus
373 _h_e_n ->fill(recon_e.size());
374 _h_mu_n ->fill(recon_mu.size());
375 _h_tau_n->fill(recon_tau.size());
376
377 // Calculate HTjets
378 double HTjets = 0.;
379 for ( const Jet & jet : recon_jets )
380 HTjets += jet.perp()/GeV;
381
382 // Calculate meff
383 double meff = eTmiss + HTjets;
384 Particles all_leptons;
385 for ( const Particle & e : recon_e ) {
386 meff += e.perp()/GeV;
387 all_leptons.push_back( e );
388 }
389 for ( const Particle & mu : recon_mu ) {
390 meff += mu.perp()/GeV;
391 all_leptons.push_back( mu );
392 }
393 for ( const Particle & tau : recon_tau ) {
394 meff += tau.perp()/GeV;
395 all_leptons.push_back( tau );
396 }
397
398 // Fill histogram of kinematic variables
399 _h_HTlep_all ->fill(HTlep);
400 _h_HTjets_all->fill(HTjets);
401 _h_MET_all ->fill(eTmiss);
402 _h_Meff_all ->fill(meff);
403
404 // Determine signal region (3l/2ltau, onZ/offZ)
405 string basic_signal_region;
406 if ( recon_mu.size() + recon_e.size() > 2 )
407 basic_signal_region += "3l_";
408 else if ( (recon_mu.size() + recon_e.size() == 2) && (recon_tau.size() > 0))
409 basic_signal_region += "2ltau_";
410 // Is there an OSSF pair or a three lepton combination with an invariant mass close to the Z mass
411 int onZ = isonZ(chosen_leptons);
412 if (onZ == 1) basic_signal_region += "onZ";
413 else if (onZ == 0) basic_signal_region += "offZ";
414 // Check in which signal regions this event falls and adjust event counters
415 fillEventCountsPerSR(basic_signal_region, onZ, HTlep, eTmiss, HTjets, meff);
416 }
417
418
419 /// Normalise histograms etc., after the run
420 void finalize() {
421
422 // Normalize to an integrated luminosity of 1 fb-1
423 double norm = crossSection()/femtobarn/sumOfWeights();
424 string best_signal_region = "";
425 double ratio_best_SR = 0.;
426
427 // Loop over all signal regions and find signal region with best sensitivity (ratio signal events/visible cross-section)
428 for (size_t i = 0; i < _signal_regions.size(); i++) {
429 double signal_events = _eventCountsPerSR[_signal_regions[i]]->val() * norm;
430 // Use expected upper limits to find best signal region
431 double UL95 = getUpperLimit(_signal_regions[i], false);
432 double ratio = signal_events / UL95;
433 if (ratio > ratio_best_SR) {
434 best_signal_region = _signal_regions[i];
435 ratio_best_SR = ratio;
436 }
437 }
438
439 double signal_events_best_SR = _eventCountsPerSR[best_signal_region]->val() * norm;
440 double exp_UL_best_SR = getUpperLimit(best_signal_region, false);
441 double obs_UL_best_SR = getUpperLimit(best_signal_region, true);
442
443 // Print out result
444 cout << "----------------------------------------------------------------------------------------" << endl;
445 cout << "Best signal region: " << best_signal_region << endl;
446 cout << "Normalized number of signal events in this best signal region (per fb-1): " << signal_events_best_SR << endl;
447 cout << "Efficiency*Acceptance: " << _eventCountsPerSR[best_signal_region]->val()/sumOfWeights() << endl;
448 cout << "Cross-section [fb]: " << crossSection()/femtobarn << endl;
449 cout << "Expected visible cross-section (per fb-1): " << exp_UL_best_SR << endl;
450 cout << "Ratio (signal events / expected visible cross-section): " << ratio_best_SR << endl;
451 cout << "Observed visible cross-section (per fb-1): " << obs_UL_best_SR << endl;
452 cout << "Ratio (signal events / observed visible cross-section): " << signal_events_best_SR/obs_UL_best_SR << endl;
453 cout << "----------------------------------------------------------------------------------------" << endl;
454
455 cout << "Using the EXPECTED limits (visible cross-section) of the analysis: " << endl;
456 if (signal_events_best_SR > exp_UL_best_SR) {
457 cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl;
458 _h_excluded->fill(1);
459 }
460 else {
461 cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl;
462 _h_excluded->fill(0);
463 }
464 cout << "----------------------------------------------------------------------------------------" << endl;
465
466 cout << "Using the OBSERVED limits (visible cross-section) of the analysis: " << endl;
467 if (signal_events_best_SR > obs_UL_best_SR) {
468 cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl;
469 _h_excluded->fill(1);
470 }
471 else {
472 cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl;
473 _h_excluded->fill(0);
474 }
475 cout << "----------------------------------------------------------------------------------------" << endl;
476
477
478 // Normalize to cross section
479 if (norm != 0) {
480 scale(_h_HTlep_all, norm);
481 scale(_h_HTjets_all, norm);
482 scale(_h_MET_all, norm);
483 scale(_h_Meff_all, norm);
484
485 scale(_h_pt_1_3l, norm);
486 scale(_h_pt_2_3l, norm);
487 scale(_h_pt_3_3l, norm);
488 scale(_h_pt_1_2ltau, norm);
489 scale(_h_pt_2_2ltau, norm);
490 scale(_h_pt_3_2ltau, norm);
491
492 scale(_h_e_n, norm);
493 scale(_h_mu_n, norm);
494 scale(_h_tau_n, norm);
495
496 scale(_h_excluded, signal_events_best_SR);
497 }
498 }
499
500
501 /// Helper functions
502 /// @{
503
504 /// Function giving a list of all signal regions
505 vector<string> getSignalRegions() {
506
507 // List of basic signal regions
508 vector<string> basic_signal_regions;
509 basic_signal_regions.push_back("3l_offZ");
510 basic_signal_regions.push_back("3l_onZ");
511 basic_signal_regions.push_back("2ltau_offZ");
512 basic_signal_regions.push_back("2ltau_onZ");
513
514 // List of kinematic variables
515 vector<string> kinematic_variables;
516 kinematic_variables.push_back("HTlep");
517 kinematic_variables.push_back("METStrong");
518 kinematic_variables.push_back("METWeak");
519 kinematic_variables.push_back("Meff");
520 kinematic_variables.push_back("MeffStrong");
521
522 vector<string> signal_regions;
523 // Loop over all kinematic variables and basic signal regions
524 for (size_t i0 = 0; i0 < kinematic_variables.size(); i0++) {
525 for (size_t i1 = 0; i1 < basic_signal_regions.size(); i1++) {
526 // Is signal region onZ?
527 int onZ = (basic_signal_regions[i1].find("onZ") != string::npos) ? 1 : 0;
528 // Get cut values for this kinematic variable
529 vector<int> cut_values = getCutsPerSignalRegion(kinematic_variables[i0], onZ);
530 // Loop over all cut values
531 for (size_t i2 = 0; i2 < cut_values.size(); i2++) {
532 // push signal region into vector
533 signal_regions.push_back( (kinematic_variables[i0] + "_" + basic_signal_regions[i1] + "_cut_" + toString(i2)) );
534 }
535 }
536 }
537 return signal_regions;
538 }
539
540
541 /// Function giving all cut vales per kinematic variable (taking onZ for MET into account)
542 vector<int> getCutsPerSignalRegion(const string& signal_region, int onZ=0) {
543 vector<int> cutValues;
544
545 // Cut values for HTlep
546 if (signal_region.compare("HTlep") == 0) {
547 cutValues.push_back(0);
548 cutValues.push_back(100);
549 cutValues.push_back(150);
550 cutValues.push_back(200);
551 cutValues.push_back(300);
552 }
553 // Cut values for METStrong (HTjets > 100 GeV) and METWeak (HTjets < 100 GeV)
554 else if (signal_region.compare("METStrong") == 0 || signal_region.compare("METWeak") == 0) {
555 if (onZ == 0) cutValues.push_back(0);
556 else if (onZ == 1) cutValues.push_back(20);
557 cutValues.push_back(50);
558 cutValues.push_back(75);
559 }
560 // Cut values for Meff and MeffStrong (MET > 75 GeV)
561 if (signal_region.compare("Meff") == 0 || signal_region.compare("MeffStrong") == 0) {
562 cutValues.push_back(0);
563 cutValues.push_back(150);
564 cutValues.push_back(300);
565 cutValues.push_back(500);
566 }
567
568 return cutValues;
569 }
570
571
572 /// function fills map EventCountsPerSR by looping over all signal regions
573 /// and looking if the event falls into this signal region
574 void fillEventCountsPerSR(const string& basic_signal_region, int onZ,
575 double HTlep, double eTmiss,
576 double HTjets, double meff) {
577
578 // Get cut values for HTlep, loop over them and add event if cut is passed
579 vector<int> cut_values = getCutsPerSignalRegion("HTlep", onZ);
580 for (size_t i = 0; i < cut_values.size(); i++) {
581 if (HTlep > cut_values[i])
582 _eventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
583 }
584
585 // Get cut values for METStrong, loop over them and add event if cut is passed
586 cut_values = getCutsPerSignalRegion("METStrong", onZ);
587 for (size_t i = 0; i < cut_values.size(); i++) {
588 if (eTmiss > cut_values[i] && HTjets > 100.)
589 _eventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
590 }
591
592 // Get cut values for METWeak, loop over them and add event if cut is passed
593 cut_values = getCutsPerSignalRegion("METWeak", onZ);
594 for (size_t i = 0; i < cut_values.size(); i++) {
595 if (eTmiss > cut_values[i] && HTjets <= 100.)
596 _eventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
597 }
598
599 // Get cut values for Meff, loop over them and add event if cut is passed
600 cut_values = getCutsPerSignalRegion("Meff", onZ);
601 for (size_t i = 0; i < cut_values.size(); i++) {
602 if (meff > cut_values[i])
603 _eventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
604 }
605
606 // Get cut values for MeffStrong, loop over them and add event if cut is passed
607 cut_values = getCutsPerSignalRegion("MeffStrong", onZ);
608 for (size_t i = 0; i < cut_values.size(); i++) {
609 if (meff > cut_values[i] && eTmiss > 75.)
610 _eventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
611 }
612 }
613
614
615 /// Function returning 4-vector of daughter-particle if it is a tau neutrino
616 /// @todo Move to TauFinder and make less HepMC-ish
617 FourMomentum get_tau_neutrino_mom(const Particle& p) {
618 assert(p.abspid() == PID::TAU);
619 ConstGenVertexPtr dv = p.genParticle()->end_vertex();
620 assert(dv != nullptr);
621 for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
622 if (abs(pp->pdg_id()) == PID::NU_TAU) return FourMomentum(pp->momentum());
623 }
624 return FourMomentum();
625 }
626
627
628 /// Function calculating the prong number of taus
629 /// @todo Move to TauFinder and make less HepMC-ish
630 void get_prong_number(ConstGenParticlePtr p, unsigned int& nprong, bool& lep_decaying_tau) {
631 assert(p != nullptr);
632 //const int tau_barcode = p->barcode();
633 ConstGenVertexPtr dv = p->end_vertex();
634 assert(dv != nullptr);
635 for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
636 // If they have status 1 and are charged they will produce a track and the prong number is +1
637 if (pp->status() == 1 ) {
638 const int id = pp->pdg_id();
639 if (Rivet::PID::charge(id) != 0 ) ++nprong;
640 // Check if tau decays leptonically
641 // @todo Can a tau decay include a tau in its decay daughters?!
642 if ((abs(id) == PID::ELECTRON || abs(id) == PID::MUON || abs(id) == PID::TAU) && abs(p->pdg_id()) == PID::TAU) lep_decaying_tau = true;
643 }
644 // If the status of the daughter particle is 2 it is unstable and the further decays are checked
645 else if (pp->status() == 2 ) {
646 get_prong_number(pp, nprong, lep_decaying_tau);
647 }
648 }
649 }
650
651
652 /// Function giving fiducial lepton efficiency
653 double apply_reco_eff(int flavor, const Particle& p) {
654 float pt = p.pT()/GeV;
655 float eta = p.eta();
656
657 double eff = 0.;
658 //double err = 0.;
659
660 if (flavor == 11) { // weight prompt electron -- now including data/MC ID SF in eff.
661 //float rho = 0.820;
662 float p0 = 7.34; float p1 = 0.8977;
663 //float ep0= 0.5 ; float ep1= 0.0087;
664 eff = p1 - p0/pt;
665
666 //double err0 = ep0/pt; // d(eff)/dp0
667 //double err1 = ep1; // d(eff)/dp1
668 //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
669
670 double avgrate = 0.6867;
671 float wz_ele_eta[] = {0.588717,0.603674,0.666135,0.747493,0.762202,0.675051,0.751606,0.745569,0.665333,0.610432,0.592693,};
672 //float ewz_ele_eta[] ={0.00292902,0.002476,0.00241209,0.00182319,0.00194339,0.00299785,0.00197339,0.00182004,0.00241793,0.00245997,0.00290394,};
673 int ibin = 3;
674
675 if (eta >= -2.5 && eta < -2.0) ibin = 0;
676 if (eta >= -2.0 && eta < -1.5) ibin = 1;
677 if (eta >= -1.5 && eta < -1.0) ibin = 2;
678 if (eta >= -1.0 && eta < -0.5) ibin = 3;
679 if (eta >= -0.5 && eta < -0.1) ibin = 4;
680 if (eta >= -0.1 && eta < 0.1) ibin = 5;
681 if (eta >= 0.1 && eta < 0.5) ibin = 6;
682 if (eta >= 0.5 && eta < 1.0) ibin = 7;
683 if (eta >= 1.0 && eta < 1.5) ibin = 8;
684 if (eta >= 1.5 && eta < 2.0) ibin = 9;
685 if (eta >= 2.0 && eta < 2.5) ibin = 10;
686
687 double eff_eta = wz_ele_eta[ibin];
688 //double err_eta = ewz_ele_eta[ibin];
689
690 eff = (eff*eff_eta)/avgrate;
691 }
692
693 if (flavor == 12) { // weight electron from tau
694 //float rho = 0.884;
695 float p0 = 6.799; float p1 = 0.842;
696 //float ep0= 0.664; float ep1= 0.016;
697 eff = p1 - p0/pt;
698
699 //double err0 = ep0/pt; // d(eff)/dp0
700 //double err1 = ep1; // d(eff)/dp1
701 //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
702
703 double avgrate = 0.5319;
704 float wz_elet_eta[] = {0.468945,0.465953,0.489545,0.58709,0.59669,0.515829,0.59284,0.575828,0.498181,0.463536,0.481738,};
705 //float ewz_elet_eta[] ={0.00933795,0.00780868,0.00792679,0.00642083,0.00692652,0.0101568,0.00698452,0.00643524,0.0080002,0.00776238,0.0094699,};
706 int ibin = 3;
707
708 if (eta >= -2.5 && eta < -2.0) ibin = 0;
709 if (eta >= -2.0 && eta < -1.5) ibin = 1;
710 if (eta >= -1.5 && eta < -1.0) ibin = 2;
711 if (eta >= -1.0 && eta < -0.5) ibin = 3;
712 if (eta >= -0.5 && eta < -0.1) ibin = 4;
713 if (eta >= -0.1 && eta < 0.1) ibin = 5;
714 if (eta >= 0.1 && eta < 0.5) ibin = 6;
715 if (eta >= 0.5 && eta < 1.0) ibin = 7;
716 if (eta >= 1.0 && eta < 1.5) ibin = 8;
717 if (eta >= 1.5 && eta < 2.0) ibin = 9;
718 if (eta >= 2.0 && eta < 2.5) ibin = 10;
719
720 double eff_eta = wz_elet_eta[ibin];
721 //double err_eta = ewz_elet_eta[ibin];
722
723 eff = (eff*eff_eta)/avgrate;
724
725 }
726
727 if (flavor == 13) {// weight prompt muon
728
729 //if eta>0.1
730 float p0 = -18.21; float p1 = 14.83; float p2 = 0.9312;
731 //float ep0= 5.06; float ep1= 1.9; float ep2=0.00069;
732
733 if ( fabs(eta) < 0.1) {
734 p0 = 7.459; p1 = 2.615; p2 = 0.5138;
735 //ep0 = 10.4; ep1 = 4.934; ep2 = 0.0034;
736 }
737
738 double arg = ( pt-p0 )/( 2.*p1 ) ;
739 eff = 0.5 * p2 * (1.+erf(arg));
740 //err = 0.1*eff;
741 }
742
743 if (flavor == 14) {// weight muon from tau
744
745 if (fabs(eta) < 0.1) {
746 float p0 = -1.756; float p1 = 12.38; float p2 = 0.4441;
747 //float ep0= 10.39; float ep1= 7.9; float ep2=0.022;
748 double arg = ( pt-p0 )/( 2.*p1 ) ;
749 eff = 0.5 * p2 * (1.+erf(arg));
750 //err = 0.1*eff;
751 }
752 else {
753 float p0 = 2.102; float p1 = 0.8293;
754 //float ep0= 0.271; float ep1= 0.0083;
755 eff = p1 - p0/pt;
756 //double err0 = ep0/pt; // d(eff)/dp0
757 //double err1 = ep1; // d(eff)/dp1
758 //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
759 }
760 }
761
762 if (flavor == 15) {// weight hadronic tau 1p
763
764 float wz_tau1p[] = {0.0249278,0.146978,0.225049,0.229212,0.21519,0.206152,0.201559,0.197917,0.209249,0.228336,0.193548,};
765 //float ewz_tau1p[] ={0.00178577,0.00425252,0.00535052,0.00592126,0.00484684,0.00612941,0.00792099,0.0083006,0.0138307,0.015568,0.0501751,};
766 int ibin = 0;
767 if (pt > 15) ibin = 1;
768 if (pt > 20) ibin = 2;
769 if (pt > 25) ibin = 3;
770 if (pt > 30) ibin = 4;
771 if (pt > 40) ibin = 5;
772 if (pt > 50) ibin = 6;
773 if (pt > 60) ibin = 7;
774 if (pt > 80) ibin = 8;
775 if (pt > 100) ibin = 9;
776 if (pt > 200) ibin = 10;
777
778 eff = wz_tau1p[ibin];
779 //err = ewz_tau1p[ibin];
780
781
782 double avgrate = 0.1718;
783 float wz_tau1p_eta[] = {0.162132,0.176393,0.139619,0.178813,0.185144,0.210027,0.203937,0.178688,0.137034,0.164216,0.163713,};
784 //float ewz_tau1p_eta[] ={0.00706705,0.00617989,0.00506798,0.00525172,0.00581865,0.00865675,0.00599245,0.00529877,0.00506368,0.00617025,0.00726219,};
785
786 ibin = 3;
787 if (eta >= -2.5 && eta < -2.0) ibin = 0;
788 if (eta >= -2.0 && eta < -1.5) ibin = 1;
789 if (eta >= -1.5 && eta < -1.0) ibin = 2;
790 if (eta >= -1.0 && eta < -0.5) ibin = 3;
791 if (eta >= -0.5 && eta < -0.1) ibin = 4;
792 if (eta >= -0.1 && eta < 0.1) ibin = 5;
793 if (eta >= 0.1 && eta < 0.5) ibin = 6;
794 if (eta >= 0.5 && eta < 1.0) ibin = 7;
795 if (eta >= 1.0 && eta < 1.5) ibin = 8;
796 if (eta >= 1.5 && eta < 2.0) ibin = 9;
797 if (eta >= 2.0 && eta < 2.5) ibin = 10;
798
799 double eff_eta = wz_tau1p_eta[ibin];
800 //double err_eta = ewz_tau1p_eta[ibin];
801
802 eff = (eff*eff_eta)/avgrate;
803 }
804
805 if (flavor == 16) { //weight hadronic tau 3p
806
807 float wz_tau3p[] = {0.000587199,0.00247181,0.0013031,0.00280112,};
808 //float ewz_tau3p[] ={0.000415091,0.000617187,0.000582385,0.00197792,};
809
810 int ibin = 0;
811 if (pt > 15) ibin = 1;
812 if (pt > 20) ibin = 2;
813 if (pt > 40) ibin = 3;
814 if (pt > 80) ibin = 4;
815
816 eff = wz_tau3p[ibin];
817 //err = ewz_tau3p[ibin];
818 }
819
820 return eff;
821 }
822
823
824 /// Function giving observed upper limit (visible cross-section)
825 double getUpperLimit(const string& signal_region, bool observed) {
826 map<string,double> upperLimitsObserved;
827 upperLimitsObserved["HTlep_3l_offZ_cut_0"] = 11.;
828 upperLimitsObserved["HTlep_3l_offZ_cut_100"] = 8.7;
829 upperLimitsObserved["HTlep_3l_offZ_cut_150"] = 4.0;
830 upperLimitsObserved["HTlep_3l_offZ_cut_200"] = 4.4;
831 upperLimitsObserved["HTlep_3l_offZ_cut_300"] = 1.6;
832 upperLimitsObserved["HTlep_2ltau_offZ_cut_0"] = 25.;
833 upperLimitsObserved["HTlep_2ltau_offZ_cut_100"] = 14.;
834 upperLimitsObserved["HTlep_2ltau_offZ_cut_150"] = 6.1;
835 upperLimitsObserved["HTlep_2ltau_offZ_cut_200"] = 3.3;
836 upperLimitsObserved["HTlep_2ltau_offZ_cut_300"] = 1.2;
837 upperLimitsObserved["HTlep_3l_onZ_cut_0"] = 48.;
838 upperLimitsObserved["HTlep_3l_onZ_cut_100"] = 38.;
839 upperLimitsObserved["HTlep_3l_onZ_cut_150"] = 14.;
840 upperLimitsObserved["HTlep_3l_onZ_cut_200"] = 7.2;
841 upperLimitsObserved["HTlep_3l_onZ_cut_300"] = 4.5;
842 upperLimitsObserved["HTlep_2ltau_onZ_cut_0"] = 85.;
843 upperLimitsObserved["HTlep_2ltau_onZ_cut_100"] = 53.;
844 upperLimitsObserved["HTlep_2ltau_onZ_cut_150"] = 11.0;
845 upperLimitsObserved["HTlep_2ltau_onZ_cut_200"] = 5.2;
846 upperLimitsObserved["HTlep_2ltau_onZ_cut_300"] = 3.0;
847 upperLimitsObserved["METStrong_3l_offZ_cut_0"] = 2.6;
848 upperLimitsObserved["METStrong_3l_offZ_cut_50"] = 2.1;
849 upperLimitsObserved["METStrong_3l_offZ_cut_75"] = 2.1;
850 upperLimitsObserved["METStrong_2ltau_offZ_cut_0"] = 4.2;
851 upperLimitsObserved["METStrong_2ltau_offZ_cut_50"] = 3.1;
852 upperLimitsObserved["METStrong_2ltau_offZ_cut_75"] = 2.6;
853 upperLimitsObserved["METStrong_3l_onZ_cut_20"] = 11.0;
854 upperLimitsObserved["METStrong_3l_onZ_cut_50"] = 6.4;
855 upperLimitsObserved["METStrong_3l_onZ_cut_75"] = 5.1;
856 upperLimitsObserved["METStrong_2ltau_onZ_cut_20"] = 5.9;
857 upperLimitsObserved["METStrong_2ltau_onZ_cut_50"] = 3.4;
858 upperLimitsObserved["METStrong_2ltau_onZ_cut_75"] = 1.2;
859 upperLimitsObserved["METWeak_3l_offZ_cut_0"] = 11.;
860 upperLimitsObserved["METWeak_3l_offZ_cut_50"] = 5.3;
861 upperLimitsObserved["METWeak_3l_offZ_cut_75"] = 3.1;
862 upperLimitsObserved["METWeak_2ltau_offZ_cut_0"] = 23.;
863 upperLimitsObserved["METWeak_2ltau_offZ_cut_50"] = 4.3;
864 upperLimitsObserved["METWeak_2ltau_offZ_cut_75"] = 3.1;
865 upperLimitsObserved["METWeak_3l_onZ_cut_20"] = 41.;
866 upperLimitsObserved["METWeak_3l_onZ_cut_50"] = 16.;
867 upperLimitsObserved["METWeak_3l_onZ_cut_75"] = 8.0;
868 upperLimitsObserved["METWeak_2ltau_onZ_cut_20"] = 80.;
869 upperLimitsObserved["METWeak_2ltau_onZ_cut_50"] = 4.4;
870 upperLimitsObserved["METWeak_2ltau_onZ_cut_75"] = 1.8;
871 upperLimitsObserved["Meff_3l_offZ_cut_0"] = 11.;
872 upperLimitsObserved["Meff_3l_offZ_cut_150"] = 8.1;
873 upperLimitsObserved["Meff_3l_offZ_cut_300"] = 3.1;
874 upperLimitsObserved["Meff_3l_offZ_cut_500"] = 2.1;
875 upperLimitsObserved["Meff_2ltau_offZ_cut_0"] = 25.;
876 upperLimitsObserved["Meff_2ltau_offZ_cut_150"] = 12.;
877 upperLimitsObserved["Meff_2ltau_offZ_cut_300"] = 3.9;
878 upperLimitsObserved["Meff_2ltau_offZ_cut_500"] = 2.2;
879 upperLimitsObserved["Meff_3l_onZ_cut_0"] = 48.;
880 upperLimitsObserved["Meff_3l_onZ_cut_150"] = 37.;
881 upperLimitsObserved["Meff_3l_onZ_cut_300"] = 11.;
882 upperLimitsObserved["Meff_3l_onZ_cut_500"] = 4.8;
883 upperLimitsObserved["Meff_2ltau_onZ_cut_0"] = 85.;
884 upperLimitsObserved["Meff_2ltau_onZ_cut_150"] = 28.;
885 upperLimitsObserved["Meff_2ltau_onZ_cut_300"] = 5.9;
886 upperLimitsObserved["Meff_2ltau_onZ_cut_500"] = 1.9;
887 upperLimitsObserved["MeffStrong_3l_offZ_cut_0"] = 3.8;
888 upperLimitsObserved["MeffStrong_3l_offZ_cut_150"] = 3.8;
889 upperLimitsObserved["MeffStrong_3l_offZ_cut_300"] = 2.8;
890 upperLimitsObserved["MeffStrong_3l_offZ_cut_500"] = 2.1;
891 upperLimitsObserved["MeffStrong_2ltau_offZ_cut_0"] = 3.9;
892 upperLimitsObserved["MeffStrong_2ltau_offZ_cut_150"] = 4.0;
893 upperLimitsObserved["MeffStrong_2ltau_offZ_cut_300"] = 2.9;
894 upperLimitsObserved["MeffStrong_2ltau_offZ_cut_500"] = 1.5;
895 upperLimitsObserved["MeffStrong_3l_onZ_cut_0"] = 10.0;
896 upperLimitsObserved["MeffStrong_3l_onZ_cut_150"] = 10.0;
897 upperLimitsObserved["MeffStrong_3l_onZ_cut_300"] = 6.8;
898 upperLimitsObserved["MeffStrong_3l_onZ_cut_500"] = 3.9;
899 upperLimitsObserved["MeffStrong_2ltau_onZ_cut_0"] = 1.6;
900 upperLimitsObserved["MeffStrong_2ltau_onZ_cut_150"] = 1.4;
901 upperLimitsObserved["MeffStrong_2ltau_onZ_cut_300"] = 1.5;
902 upperLimitsObserved["MeffStrong_2ltau_onZ_cut_500"] = 0.9;
903
904 // Expected upper limits are also given but not used in this analysis
905 map<string,double> upperLimitsExpected;
906 upperLimitsExpected["HTlep_3l_offZ_cut_0"] = 11.;
907 upperLimitsExpected["HTlep_3l_offZ_cut_100"] = 8.5;
908 upperLimitsExpected["HTlep_3l_offZ_cut_150"] = 4.6;
909 upperLimitsExpected["HTlep_3l_offZ_cut_200"] = 3.6;
910 upperLimitsExpected["HTlep_3l_offZ_cut_300"] = 1.9;
911 upperLimitsExpected["HTlep_2ltau_offZ_cut_0"] = 23.;
912 upperLimitsExpected["HTlep_2ltau_offZ_cut_100"] = 14.;
913 upperLimitsExpected["HTlep_2ltau_offZ_cut_150"] = 6.4;
914 upperLimitsExpected["HTlep_2ltau_offZ_cut_200"] = 3.6;
915 upperLimitsExpected["HTlep_2ltau_offZ_cut_300"] = 1.5;
916 upperLimitsExpected["HTlep_3l_onZ_cut_0"] = 33.;
917 upperLimitsExpected["HTlep_3l_onZ_cut_100"] = 25.;
918 upperLimitsExpected["HTlep_3l_onZ_cut_150"] = 12.;
919 upperLimitsExpected["HTlep_3l_onZ_cut_200"] = 6.5;
920 upperLimitsExpected["HTlep_3l_onZ_cut_300"] = 3.1;
921 upperLimitsExpected["HTlep_2ltau_onZ_cut_0"] = 94.;
922 upperLimitsExpected["HTlep_2ltau_onZ_cut_100"] = 61.;
923 upperLimitsExpected["HTlep_2ltau_onZ_cut_150"] = 9.9;
924 upperLimitsExpected["HTlep_2ltau_onZ_cut_200"] = 4.5;
925 upperLimitsExpected["HTlep_2ltau_onZ_cut_300"] = 1.9;
926 upperLimitsExpected["METStrong_3l_offZ_cut_0"] = 3.1;
927 upperLimitsExpected["METStrong_3l_offZ_cut_50"] = 2.4;
928 upperLimitsExpected["METStrong_3l_offZ_cut_75"] = 2.3;
929 upperLimitsExpected["METStrong_2ltau_offZ_cut_0"] = 4.8;
930 upperLimitsExpected["METStrong_2ltau_offZ_cut_50"] = 3.3;
931 upperLimitsExpected["METStrong_2ltau_offZ_cut_75"] = 2.1;
932 upperLimitsExpected["METStrong_3l_onZ_cut_20"] = 8.7;
933 upperLimitsExpected["METStrong_3l_onZ_cut_50"] = 4.9;
934 upperLimitsExpected["METStrong_3l_onZ_cut_75"] = 3.8;
935 upperLimitsExpected["METStrong_2ltau_onZ_cut_20"] = 7.3;
936 upperLimitsExpected["METStrong_2ltau_onZ_cut_50"] = 2.8;
937 upperLimitsExpected["METStrong_2ltau_onZ_cut_75"] = 1.5;
938 upperLimitsExpected["METWeak_3l_offZ_cut_0"] = 10.;
939 upperLimitsExpected["METWeak_3l_offZ_cut_50"] = 4.7;
940 upperLimitsExpected["METWeak_3l_offZ_cut_75"] = 3.0;
941 upperLimitsExpected["METWeak_2ltau_offZ_cut_0"] = 21.;
942 upperLimitsExpected["METWeak_2ltau_offZ_cut_50"] = 4.0;
943 upperLimitsExpected["METWeak_2ltau_offZ_cut_75"] = 2.6;
944 upperLimitsExpected["METWeak_3l_onZ_cut_20"] = 30.;
945 upperLimitsExpected["METWeak_3l_onZ_cut_50"] = 10.;
946 upperLimitsExpected["METWeak_3l_onZ_cut_75"] = 5.4;
947 upperLimitsExpected["METWeak_2ltau_onZ_cut_20"] = 88.;
948 upperLimitsExpected["METWeak_2ltau_onZ_cut_50"] = 5.5;
949 upperLimitsExpected["METWeak_2ltau_onZ_cut_75"] = 2.2;
950 upperLimitsExpected["Meff_3l_offZ_cut_0"] = 11.;
951 upperLimitsExpected["Meff_3l_offZ_cut_150"] = 8.8;
952 upperLimitsExpected["Meff_3l_offZ_cut_300"] = 3.7;
953 upperLimitsExpected["Meff_3l_offZ_cut_500"] = 2.1;
954 upperLimitsExpected["Meff_2ltau_offZ_cut_0"] = 23.;
955 upperLimitsExpected["Meff_2ltau_offZ_cut_150"] = 13.;
956 upperLimitsExpected["Meff_2ltau_offZ_cut_300"] = 4.9;
957 upperLimitsExpected["Meff_2ltau_offZ_cut_500"] = 2.4;
958 upperLimitsExpected["Meff_3l_onZ_cut_0"] = 33.;
959 upperLimitsExpected["Meff_3l_onZ_cut_150"] = 25.;
960 upperLimitsExpected["Meff_3l_onZ_cut_300"] = 9.;
961 upperLimitsExpected["Meff_3l_onZ_cut_500"] = 3.9;
962 upperLimitsExpected["Meff_2ltau_onZ_cut_0"] = 94.;
963 upperLimitsExpected["Meff_2ltau_onZ_cut_150"] = 35.;
964 upperLimitsExpected["Meff_2ltau_onZ_cut_300"] = 6.8;
965 upperLimitsExpected["Meff_2ltau_onZ_cut_500"] = 2.5;
966 upperLimitsExpected["MeffStrong_3l_offZ_cut_0"] = 3.9;
967 upperLimitsExpected["MeffStrong_3l_offZ_cut_150"] = 3.9;
968 upperLimitsExpected["MeffStrong_3l_offZ_cut_300"] = 3.0;
969 upperLimitsExpected["MeffStrong_3l_offZ_cut_500"] = 2.0;
970 upperLimitsExpected["MeffStrong_2ltau_offZ_cut_0"] = 3.8;
971 upperLimitsExpected["MeffStrong_2ltau_offZ_cut_150"] = 3.9;
972 upperLimitsExpected["MeffStrong_2ltau_offZ_cut_300"] = 3.1;
973 upperLimitsExpected["MeffStrong_2ltau_offZ_cut_500"] = 1.6;
974 upperLimitsExpected["MeffStrong_3l_onZ_cut_0"] = 6.9;
975 upperLimitsExpected["MeffStrong_3l_onZ_cut_150"] = 7.1;
976 upperLimitsExpected["MeffStrong_3l_onZ_cut_300"] = 4.9;
977 upperLimitsExpected["MeffStrong_3l_onZ_cut_500"] = 3.0;
978 upperLimitsExpected["MeffStrong_2ltau_onZ_cut_0"] = 2.4;
979 upperLimitsExpected["MeffStrong_2ltau_onZ_cut_150"] = 2.5;
980 upperLimitsExpected["MeffStrong_2ltau_onZ_cut_300"] = 2.0;
981 upperLimitsExpected["MeffStrong_2ltau_onZ_cut_500"] = 1.1;
982
983 if (observed) return upperLimitsObserved[signal_region];
984 else return upperLimitsExpected[signal_region];
985 }
986
987
988 /// Function checking if there is an OSSF lepton pair or a combination of 3 leptons with an invariant mass close to the Z mass
989 /// @todo Should the reference Z mass be 91.2?
990 int isonZ (const Particles& particles) {
991 int onZ = 0;
992 double best_mass_2 = 999.;
993 double best_mass_3 = 999.;
994
995 // Loop over all 2 particle combinations to find invariant mass of OSSF pair closest to Z mass
996 for ( const Particle& p1 : particles ) {
997 for ( const Particle& p2 : particles ) {
998 double mass_difference_2_old = fabs(91.0 - best_mass_2);
999 double mass_difference_2_new = fabs(91.0 - (p1.momentum() + p2.momentum()).mass()/GeV);
1000
1001 // If particle combination is OSSF pair calculate mass difference to Z mass
1002 if ( (p1.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -169) ) {
1003
1004 // Get invariant mass closest to Z mass
1005 if (mass_difference_2_new < mass_difference_2_old)
1006 best_mass_2 = (p1.momentum() + p2.momentum()).mass()/GeV;
1007
1008 // In case there is an OSSF pair take also 3rd lepton into account (e.g. from FSR and photon to electron conversion)
1009 for ( const Particle & p3 : particles ) {
1010 double mass_difference_3_old = fabs(91.0 - best_mass_3);
1011 double mass_difference_3_new = fabs(91.0 - (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV);
1012 if (mass_difference_3_new < mass_difference_3_old)
1013 best_mass_3 = (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV;
1014 }
1015 }
1016 }
1017 }
1018
1019 // Pick the minimum invariant mass of the best OSSF pair combination and the best 3 lepton combination
1020 // If this mass is in a 20 GeV window around the Z mass, the event is classified as onZ
1021 double best_mass = min(best_mass_2, best_mass_3);
1022 if (fabs(91.0 - best_mass) < 20) onZ = 1;
1023 return onZ;
1024 }
1025
1026 /// @}
1027
1028
1029 private:
1030
1031 /// Histograms
1032 /// @{
1033 Histo1DPtr _h_HTlep_all, _h_HTjets_all, _h_MET_all, _h_Meff_all;
1034 Histo1DPtr _h_pt_1_3l, _h_pt_2_3l, _h_pt_3_3l, _h_pt_1_2ltau, _h_pt_2_2ltau, _h_pt_3_2ltau;
1035 Histo1DPtr _h_e_n, _h_mu_n, _h_tau_n;
1036 Histo1DPtr _h_excluded;
1037 /// @}
1038
1039 /// Fiducial efficiencies to model the effects of the ATLAS detector
1040 bool _use_fiducial_lepton_efficiency;
1041
1042 /// List of signal regions and event counts per signal region
1043 vector<string> _signal_regions;
1044
1045 map<string, CounterPtr> _eventCountsPerSR;
1046
1047 };
1048
1049
1050 RIVET_DECLARE_PLUGIN(ATLAS_2012_I1204447);
1051
1052}
|