Rivet analyses referenceATLAS_2017_I1598613BB to Jpsi plus mu at 8 TeVExperiment: ATLAS (LHC) Inspire ID: 1598613 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
A measurement of $b$-hadron pair production is presented, based on a dataset corresponding to an integrated luminosity of 11.4 fb${}^{-1}$ of proton--proton collisions recorded at $\sqrt{s}=8$ TeV with the ATLAS detector at the LHC. Events are selected in which both a $b$-hadron $\to J/\psi(\to\mu\mu)+X$ and $b$-hadron $\to \mu + X$ were identified, and results are presented in a fiducial volume defined by kinematic requirements on three muons based on those used in the analysis. The fiducial cross section is measured to be $17.7\pm 0.1$(stat.)$\pm2.0$(syst.) nb. A number of normalised differential cross sections are i also measured, and compared to predictions from the Pythia8, Herwig++, MadGraph5\_aMC@NLO+Pythia8 and Sherpa event generators, providing constraints on heavy flavour production. Source code: ATLAS_2017_I1598613.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/HeavyHadrons.hh"
4#include "Rivet/Projections/DressedLeptons.hh"
5#include "Rivet/Projections/FinalState.hh"
6#include "Rivet/Projections/PromptFinalState.hh"
7#include "Rivet/Projections/IdentifiedFinalState.hh"
8
9namespace Rivet {
10
11
12 class ATLAS_2017_I1598613 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1598613);
17
18 struct HistoHandler {
19 Histo1DPtr histo;
20 Scatter2DPtr scatter;
21 unsigned int d, x, y;
22
23 HistoHandler() {}
24
25 void fill(double value) {
26 histo->fill(value);
27 }
28 };
29
30
31 /// Book histograms and initialise projections before the run
32 void init() {
33
34 // default to widest cut, electrons and muons.
35 _mode = 0;
36 if ( getOption("BMODE") == "BB" ) _mode = 1;
37
38 // Get the particles needed for each running mode:
39 if (_mode == 0) {
40 // Get photons to dress leptons
41 FinalState photons(Cuts::abspid == PID::PHOTON);
42 FinalState muons(Cuts::abspid == PID::MUON);
43 Cut eta_lep = Cuts::abseta < 2.5;
44 DressedLeptons dressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT >= 6*GeV, true);
45 declare(dressedmuons, "dressedmuons");
46 } else {
47 declare(HeavyHadrons(Cuts::absrap < 2.4 && Cuts::pT > 15.5*GeV), "BHadrons");
48 }
49
50 //Book the histograms:
51 bookHandler(_h["dR"], 1);
52 bookHandler(_h["highpT_dR"], 4);
53 bookHandler(_h["lowpT_dR"], 7);
54 bookHandler(_h["dPhi"], 10);
55 bookHandler(_h["dy"], 13);
56 bookHandler(_h["MopT"], 16);
57 bookHandler(_h["pToM"], 19);
58 bookHandler(_h["pT"], 22);
59 bookHandler(_h["M"], 25);
60 bookHandler(_h["yboost"], 29);
61 }
62
63
64 void bookHandler(HistoHandler& handler, unsigned int id_xsec) {
65 if (_mode) {
66 book(handler.histo, "_aux_hist" + toString(id_xsec), refData(id_xsec, 1, 1));
67 book(handler.scatter, id_xsec, 1, 1, true);
68 handler.d = id_xsec + 1; // transfer function
69 handler.x = 1; handler.y = 1;
70 }
71 else book(handler.histo, id_xsec, 1, 1);
72 }
73
74
75 /// Perform the per-event analysis
76 void analyze(const Event& event) {
77
78 if (_mode == 1) { // make the 2-b-hadron-level plots
79 const Particles& bHadrons = apply<HeavyHadrons>(event, "BHadrons").bHadrons();
80 if (bHadrons.size() > 1) {
81 sortBy(bHadrons, cmpMomByPt);
82
83 float dphiBB = deltaPhi(bHadrons[0], bHadrons[1]);
84 float dRBB = deltaR(bHadrons[0], bHadrons[1], RAPIDITY);
85 float dyBB = fabs(bHadrons[0].rapidity() - bHadrons[1].rapidity());
86 float yboostBB = 0.5*fabs(bHadrons[0].rapidity() + bHadrons[1].rapidity());
87 FourMomentum systemBB = bHadrons[0].momentum() + bHadrons[1].momentum();
88 // Due to the additional particles produced in the decay,
89 // the 3 muons carry only a fraction of the momentum of the b-hadrons,
90 // scale down mass and pT to match 3-muon-level more closely
91 float MBB = systemBB.mass()/1.75;
92 float pTBB = systemBB.pT()/1.75;
93
94 _h["dPhi"].fill(dphiBB);
95 _h["dy"].fill(dyBB);
96 _h["yboost"].fill(yboostBB);
97 _h["dR"].fill(dRBB);
98 _h["M"].fill(MBB/GeV);
99 _h["pT"].fill(pTBB/GeV);
100 _h["MopT"].fill(MBB/pTBB);
101 _h["pToM"].fill(pTBB/MBB);
102
103 if (pTBB >= 20*GeV) _h["highpT_dR"].fill(dRBB);
104 else _h["lowpT_dR"].fill(dRBB);
105 }
106 }
107
108
109 if (_mode == 0) { // the 3-muon-level analysis
110
111 // First, simply check that we have enough muons
112 const vector<DressedLepton> muons = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
113 if (muons.size() < 3) vetoEvent;
114
115 // Not sure if this is going to work, but ..
116 vector<DressedLepton> Jpsi_muons, third_muons;
117 for (DressedLepton mu : muons) {
118 if (mu.constituentLepton().fromBottom() && mu.constituentLepton().hasAncestor(PID::JPSI)) {
119 Jpsi_muons.push_back(mu);
120 }
121 else if (mu.constituentLepton().fromBottom()) {
122 third_muons.push_back(mu);
123 }
124 }
125
126 // Veto events without enough muons:
127 if (Jpsi_muons.size() < 2) vetoEvent;
128
129 // At this point, we must always have a Jpsi. So get its 4-vector:
130 FourMomentum Jpsi = Jpsi_muons[0].momentum() + Jpsi_muons[1].momentum();
131
132 // If there is more than one J/psi, take the one closest to PDG mass,
133 // and push all the other muons back to the 3rd muon list
134 size_t mu1 = 0, mu2 = 1;
135 if (Jpsi_muons.size() > 2) {
136 for (size_t i = 0; i < Jpsi_muons.size(); ++i) {
137 for (size_t j = i; j < Jpsi_muons.size(); ++j) {
138 FourMomentum temp = Jpsi_muons[i].momentum() + Jpsi_muons[j].momentum();
139 if (fabs(temp.mass() - 3.096) < fabs(Jpsi.mass() - 3.096)) {
140 Jpsi = temp;
141 mu1 = i;
142 mu2 = j;
143 }
144 }
145 }
146
147 for (size_t i = 0; i < Jpsi_muons.size(); ++i) {
148 if (i == mu1 || i == mu2) continue;
149 third_muons.push_back(Jpsi_muons[i]);
150 }
151 }
152
153 // There *is* more than one Jpsi:
154 if (Jpsi_muons[mu1].abseta() >= 2.3) vetoEvent;
155 if (Jpsi_muons[mu2].abseta() >= 2.3) vetoEvent;
156
157 // We should now have the full list of 3rd muons to consider. Make sure we have one:
158 if (third_muons.empty()) vetoEvent;
159
160 // Sort the third muons by pT and pick highest one
161 std::sort(third_muons.begin(), third_muons.end(), [](const DressedLepton &l1, const DressedLepton &l2) {
162 return (l1.pT() > l2.pT());
163 });
164 FourMomentum third_mu = third_muons[0].momentum();
165
166 // Finally, make the plots!
167 float dphi = deltaPhi(Jpsi, third_mu);
168 float dR = deltaR(Jpsi, third_mu, RAPIDITY);
169 float dy = fabs(Jpsi.rapidity() - third_mu.rapidity());
170 float yboost = 0.5*fabs(Jpsi.rapidity() + third_mu.rapidity());
171 FourMomentum system = Jpsi + third_mu;
172 float M = system.mass();
173 float pT = system.pT();
174
175 _h["dPhi"].fill(dphi);
176 _h["dy"].fill(dy);
177 _h["yboost"].fill(yboost);
178 _h["dR"].fill(dR);
179 if (pT >= 20*GeV) _h["highpT_dR"].fill(dR);
180 else _h["lowpT_dR"].fill(dR);
181
182 _h["M"].fill(M);
183 _h["pT"].fill(pT);
184 _h["MopT"].fill(M/pT);
185 _h["pToM"].fill(pT/M);
186
187 } //< end of muon analysis.
188 }
189
190
191 /// Normalise histograms etc., after the run
192 void finalize() {
193 for (map<string, HistoHandler>::iterator hit = _h.begin(); hit != _h.end(); ++hit) {
194 normalize(hit->second.histo);
195 if (_mode == 1) applyTransferFnAndNorm(hit->second);
196 }
197 }
198
199
200 void applyTransferFnAndNorm(HistoHandler &handler) { ///< @todo Pass as const reference?
201 // Load transfer function from reference data file
202 const YODA::Scatter2D& myTransferFn = refData(handler.d, handler.x, handler.y);
203 double area = 0.0;
204 for (size_t i = 0; i < handler.scatter->numPoints(); ++i) {
205 const Point2D& f = myTransferFn.point(i);
206 Point2D& p = handler.scatter->point(i);
207 const HistoBin1D& b = handler.histo->bin(i);
208 double newy;
209 try {
210 newy = b.height();
211 } catch (const Exception&) { // LowStatsError or WeightError
212 newy = 0;
213 }
214 double newey;
215 try {
216 newey = b.heightErr();
217 } catch (const Exception&) { // LowStatsError or WeightError
218 newey = 0;
219 }
220 // apply transfer function here
221 newy *= f.y(); newey *= f.y();
222 double rp = safediv(newey, newy);
223 double rf = safediv(f.yErrAvg(), f.y());
224 newey = newy * sqrt(rp*rp + rf*rf);
225 // set new values
226 p.setY(newy);
227 p.setYErrMinus(newey);
228 p.setYErrPlus(newey);
229 area += newy * (p.xMax() - p.xMin());
230 }
231 if (area > 0.) { // normalise to unity
232 for (size_t i = 0; i < handler.scatter->numPoints(); ++i)
233 handler.scatter->point(i).scaleY(1.0 / area);
234 }
235 }
236
237
238 protected:
239
240 /// Analysis-mode switch
241 size_t _mode;
242
243 /// Histograms
244 map<string, HistoHandler> _h;
245
246 };
247
248
249 // Hooks for the plugin system
250 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1598613);
251
252}
|