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