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