Rivet analyses referenceATLAS_2019_I1725190Dilepton mass spectrum in 13 TeV pp collisions with 139/fb Run 2 datasetExperiment: ATLAS (LHC) Inspire ID: 1725190 Status: VALIDATED SINGLEWEIGHT NOREENTRY Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
A search for high-mass dielectron and dimuon resonances in the mass range of 250 GeV to 6 TeV. Functional forms have been fitted to the components of Crystal Ball + Gaussian dilepton invariant-mass distribution smearing functions for electrons and muons separately, which are encoded in this analysis to allow particle-level MC to be compared to the reco-level experimental data. Source code: ATLAS_2019_I1725190.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/PromptFinalState.hh"
6#include "Rivet/Projections/FastJets.hh"
7#include "Rivet/Projections/Smearing.hh"
8#include "Rivet/Tools/Random.hh"
9
10namespace Rivet {
11
12
13 /// @brief Dilepton high-mass resonance search
14 ///
15 /// @todo Use the proper smearing system rather than hand-rolled sampling
16 class ATLAS_2019_I1725190 : public Analysis {
17 public:
18
19 /// Constructor
20 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1725190);
21
22
23 /// @name Analysis methods
24 /// @{
25
26 /// Book histograms and initialise projections before the run
27 void init() {
28
29 PromptFinalState electrons(Cuts::abspid == PID::ELECTRON && Cuts::Et > 30*GeV &&
30 (Cuts::abseta < 2.47 && !Cuts::absetaIn(1.37, 1.52)));
31 SmearedParticles recoelectrons(electrons, PARTICLE_EFF_ONE); //ELECTRON_EFF_ATLAS_RUN2_MEDIUM);
32 declare(recoelectrons, "Elecs");
33
34 PromptFinalState muons(Cuts::abspid == PID::MUON && Cuts::pT > 30*GeV && Cuts::abseta < 2.5);
35 SmearedParticles recomuons(muons, PARTICLE_EFF_ONE);
36 // [](const Particle& m) -> double {
37 // if (m.pT() < 1*TeV) return 0.69;
38 // if (m.pT() > 2.5*TeV) return 0.57;
39 // return 0.69 - 0.12*(m.pT() - 1*TeV)/(2.5*TeV - 1*TeV);
40 // });
41 declare(recomuons, "Muons");
42
43 // Book histograms
44 book(_h_mee, 1, 1, 1);
45 book(_h_mmm, 2, 1, 1);
46 }
47
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51
52 // Get leptons
53 const Particles elecs = apply<ParticleFinder>(event, "Elecs").particlesByPt();
54 const Particles muons = apply<ParticleFinder>(event, "Muons").particlesByPt();
55
56 // MSG_INFO("Num e, mu = " << elecs.size() << ", " << muons.size());
57 // for (const Particle& e : elecs) MSG_INFO(e.mom());
58 // for (const Particle& m : muons) MSG_INFO(m.mom());
59
60 // Isolation
61 /// @todo Can't be done from provided information?
62 // Particles isoelecs, isomuons;
63 // for (const Particle& e : elecs) {
64 // }
65
66 // Choose the highest-pT lepton pair, preferring electrons
67 if (elecs.size() < 2 && muons.size() < 2) vetoEvent;
68 const Particle l1 = (elecs.size() >= 2) ? elecs[0] : muons[0];
69 const Particle l2 = (elecs.size() >= 2) ? elecs[1] : muons[1];
70
71 // Require opposite sign for muons only
72 const bool mumu = (l1.abspid() == PID::MUON);
73 if (mumu && l1.pid()*l2.pid() > 0) vetoEvent;
74
75 // Get the true dilepton pair
76 const FourMomentum pll = l1.mom() + l2.mom();
77 const double mll = pll.mass();
78
79 // Make sure we're in a region where the smearing and efficiencies are well-behaved
80 if (mll < 200*GeV) vetoEvent;
81
82 // Apply dilepton efficiency curves (as function of mX ~ mll)
83 const double eff = mumu ?
84 (0.54 - (mll - 200*GeV)/(6000*GeV - 200*GeV) * (0.54 - 0.38)) :
85 (0.74 - 0.04*exp(-(mll-200*GeV)/100*GeV) - 0.08*exp(-(mll-200*GeV)/1000*GeV));
86 if (rand01() > eff) vetoEvent;
87
88 // Smear the dilepton mass with a CB + Gauss function
89 double muCB, sigCB, alpCB, nCB, muG, sigG, kappa;
90 if (!mumu) {
91 const double lnm = log(mll);
92 static const vector<double> pmuCB = {0.13287, -0.410663, -0.0126743, 2.9547e-6};
93 muCB = pmuCB[0] + pmuCB[1]/lnm + pmuCB[2]*lnm + pmuCB[3]*pow(lnm, 4);
94 static const vector<double> psigCB = {0.0136624, 0.230678, 1.73254};
95 sigCB = sqrt(pow(psigCB[0],2) + pow(psigCB[1],2)/mll + pow(psigCB[2]/mll, 2));
96 alpCB = 1.59112;
97 static const vector<double> pnCB = {1.13055, 0.76705, 0.00298312};
98 nCB = pnCB[0] + pnCB[1]*exp(-pnCB[2]*mll);
99 static const vector<double> pmuG = {-0.00402708, 0.814172, -3.94281e-7, 7.97076e-6, -87.6397, -1.64806e-11};
100 muG = pmuG[0] + pmuG[1]/mll + pmuG[2]*mll + pmuG[3]*pow(lnm,3) + pmuG[4]/sqr(mll) + pmuG[5]*sqr(mll);
101 static const vector<double> psigG = {0.00553858, 0.140909, 0.644418};
102 sigG = sqrt(sqr(psigG[0]) + sqr(psigG[1])/mll + sqr(psigG[2]/mll));
103 static const vector<double> pkappa = {0.347003, 0.135768, 0.00372497, -2.2822e-5, 5.06351e-13};
104 kappa = pkappa[0] + pkappa[1]*exp(-pkappa[2]*mll) + pkappa[3]*mll + pkappa[4]*pow(mll,3);
105 } else {
106 static const vector<double> pmuCB = {-0.0891397, 10.6169, -951.712, 74775.3, 5.60192e-5, -1.58827e-9, -3.81706e-13};
107 muCB = pmuCB[0] + pmuCB[1]/mll + pmuCB[2]/sqr(mll) + pmuCB[3]/pow(mll,3) + pmuCB[4]*mll + pmuCB[5]*sqr(mll) + pmuCB[6]*pow(mll,3);
108 static const vector<double> psigCB = {0.0836349, -8.98476, 491.19, 5.18068e-5, -3.45042e-10};
109 sigCB = psigCB[0] + psigCB[1]/mll + psigCB[2]/sqr(mll) + psigCB[3]*mll + psigCB[4]*sqr(mll);
110 static const vector<double> palpCB = {0.512577, 252.922, -79337.4, 7.31863e6, 0.000237883};
111 alpCB = palpCB[0] + palpCB[1]/mll + palpCB[2]/sqr(mll) + palpCB[3]/pow(mll,3) + palpCB[4]*mll;
112 static const vector<double> pnCB = {1.13055, 0.76705, 0.00298312};
113 nCB = 6.08818;
114 static const vector<double> pmuG = {-0.00410659, 2.82352e-6, -6.264e-9, 1.25547e-12, -9.94431e-17};
115 muG = pmuG[0] + pmuG[1]*mll + pmuG[2]*sqr(mll) + pmuG[3]*pow(mll,3) + pmuG[4]*pow(mll,4);
116 static const vector<double> psigG = {0.0214264, -0.795058, 15.4726, 3.38205e-5, -1.64984e-9};
117 sigG = psigG[0] + psigG[1]/mll + psigG[2]/sqr(mll) + psigG[3]*mll + psigG[4]*sqr(mll);
118 static const vector<double> pkappa = {0.235477, -31.2227, 3447.34, 4.54408e-5, -3.25374e-9};
119 kappa = pkappa[0] + pkappa[1]/mll + pkappa[2]/sqr(mll) + pkappa[3]*mll + pkappa[4]*sqr(mll);
120 }
121
122 // Do the smearing
123 double mll_scale = -1;
124 do {
125 mll_scale = (rand01() > kappa) ? randnorm(muG, sigG) : randcrystalball(alpCB, nCB, muCB, sigCB);
126 } while (fabs(mll_scale) > 0.5);
127 const double mll_reco = (1 + mll_scale) * mll;
128
129 // Require high-Mll to avoid the Z peak
130 if (mll_reco < 225*GeV) vetoEvent;
131
132 // Fill Mll histograms
133 (mumu ? _h_mmm : _h_mee)->fill(mll_reco/GeV);
134 }
135
136
137 /// Normalise histograms etc., after the run.
138 // The plot is expressed as number of events in a 10 GeV bin.
139 // Multiply by luminosity*cross section (in same units!) to get number
140 // of events a 1 GeV bin, then by 10 to get number of events in a 10 GeV.
141 void finalize() {
142 scale(_h_mee, 10.*crossSection()/picobarn*luminosity()/sumOfWeights());
143 scale(_h_mmm, 10.*crossSection()/picobarn*luminosity()/sumOfWeights());
144 }
145
146 /// @}
147
148
149 /// @name Histograms
150 /// @{
151 Histo1DPtr _h_mee, _h_mmm;
152 /// @}
153
154
155 };
156
157
158 RIVET_DECLARE_PLUGIN(ATLAS_2019_I1725190);
159
160}
|