Rivet analyses referenceH1_2005_I676166Measurement of beauty production at HERA using events with muons and jets (H1)Experiment: H1 (HERA) Inspire ID: 676166 Status: VALIDATED Authors:
Beam energies: (27.6, 920.0); (920.0, 27.6) GeV Run details:
A measurement of the beauty production cross section in $ep$ collisions at a centre-of-mass energy of 319 GeV is presented. The data were collected with the H1 detector at the HERA collider in the years 1999-2000. Events are selected by requiring the presence of jets and muons in the final state. Both the long lifetime and the large mass of b-flavoured hadrons are exploited to identify events containing beauty quarks. Differential cross sections are measured in photoproduction, with photon virtualities $Q^2 < 1$ GeV$^2$, and in deep inelastic scattering, where $2 < Q^2 < 100$ GeV$^2$. Source code: H1_2005_I676166.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/DISKinematics.hh"
6#include "Rivet/Projections/DISLepton.hh"
7#include "Rivet/Projections/DISFinalState.hh"
8
9
10namespace Rivet {
11
12
13 /// @brief Measurement of beauty production at HERA using events with muons and jets
14 class H1_2005_I676166 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2005_I676166);
19
20
21 /// @name Analysis methods
22 /// @{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 // Initialise and register projections
28 declare(FinalState(Cuts::abspid == PID::MUON), "muons");
29 declare(DISKinematics(), "Kinematics");
30 const DISLepton dl;
31 //declare(dl, "Lepton");
32 declare(dl.remainingFinalState(), "FS");
33
34 // DIS events: final state particles boosted to Breit frame then clustered
35 // using FastJet KT algorithm with jet radius parameter 1 , pT recombination scheme
36 const DISFinalState DISfs(DISFrame::BREIT);
37 declare(FastJets(DISfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::pt_scheme, 1.0,
38 JetMuons::ALL, JetInvisibles::NONE, nullptr), "DISjets");
39
40 // Photoproduction events: final state particles in lab frame then clustered
41 // using FastJet KT algorithm with jet radius parameter 1 , pT recombination scheme
42 const DISFinalState PHOfs(DISFrame::LAB);
43 declare(FastJets(PHOfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::pt_scheme, 1.0,
44 JetMuons::ALL, JetInvisibles::NONE, nullptr), "PHOjets");
45
46 // Photoproduction (Table 4)
47 book(_h["PHO_eta_mu"], 1, 1, 1);
48 book(_h["PHO_pT_mu"], 2, 1, 1);
49 book(_h["PHO_pT_jet"], 3, 1, 1);
50 book(_h["PHO_x_obs_gamma"], 4, 1, 1);
51
52 // Electroproduction (Table 5)
53 book(_h["DIS_Q^2"], 5, 1, 1);
54 book(_h["DIS_x"], 6, 1, 1);
55 book(_h["DIS_eta_mu"], 7, 1, 1);
56 book(_h["DIS_pT_mu"], 8, 1, 1);
57 book(_h["DIS_pT_jet_Breit"], 9, 1, 1);
58
59 isDIS = false;
60 nPHO = 0;
61 nDIS = 0;
62 nVeto0 = 0;
63 nVeto1 = 0;
64 nVeto2 = 0;
65 nVeto3 = 0;
66 nVeto4 = 0;
67 nVeto5 = 0;
68 nVeto6 = 0;
69 nVeto7 = 0;
70
71 }
72
73
74 /// Perform the per-event analysis
75 void analyze(const Event& event) {
76
77 isDIS = false;
78
79 //First Lorentz invariant quantities in Lab frame
80 DISKinematics dis = apply<DISKinematics>(event, "Kinematics");
81 //const DISLepton& dl = apply<DISLepton>(event,"Lepton");
82 double Q2 = dis.Q2()/GeV2;
83 double x = dis.x();
84 double y = dis.y();
85
86 //after x_g
87 //calculated from mass of
88
89 // Separate into DIS and PHO regimes else veto
90 if (Q2 < 1 && inRange(y, 0.2, 0.8)) {
91 isDIS = false;
92 ++nPHO;
93 } else if (inRange(Q2, 2.0, 100) && inRange(y, 0.1, 0.7)) {
94 isDIS = true;
95 ++nDIS;
96 } else {
97 vetoEvent;
98 }
99 ++nVeto0;
100
101 // Extract the particles other than the lepton
102 Particles particles = apply<FinalState>(event, "FS").particles();
103
104 // Get Lorentz transforms for Breit Boost and Lab Boost
105 const LorentzTransform breitboost = dis.boostBreit();
106 const LorentzTransform hcmboost = dis.boostHCM(); // Hadron cm system
107 const LorentzTransform labboost = breitboost.inverse();
108
109
110 // If DIS regime, cluster in Breit frame. If PHO regime, cluster in lab frame
111 // Apply jet cuts for relevant regime
112 Jets DISjets, PHOjets;
113 Jets cutJetsDIS, cutJetsPHO;
114 Jet pjet1, pjet2;
115
116 if (isDIS) {
117
118 DISjets = apply<FastJets>(event, "DISjets").jetsByPt(Cuts::pT > 6*GeV);
119 // Cut on Pseudorapidity in lab frame
120 for (vector<int>::size_type i=0; i<DISjets.size(); i++){
121 const FourMomentum LabMom = labboost.transform(DISjets[i]);
122 // eta cut is in lab frame
123 double etaDISJet =abs(LabMom.eta());
124 if(etaDISJet < 2.5 ){
125 cutJetsDIS.push_back(DISjets[i]);
126 }
127 }
128
129 // Apply number of jet cuts: in DIS at least one jet is needed
130 if (cutJetsDIS.size()<1) vetoEvent;
131 ++nVeto1;
132
133 } else {
134
135 PHOjets = apply<FastJets>(event, "PHOjets").jetsByPt(Cuts::pT > 6*GeV);
136
137 // Cut on Pseudorapidity
138 for(vector<int>::size_type i=0; i<PHOjets.size(); i++){
139 double etaPHOJet = PHOjets[i].eta();
140 if(abs(etaPHOJet) < 2.5 ){
141 cutJetsPHO.push_back(PHOjets[i]);
142 }
143 }
144
145 // Apply number of jet cuts: in phtotproduction at least 2 jets are needed
146
147 if(cutJetsPHO.size()<2) vetoEvent;
148 ++nVeto2;
149
150 // Ensure 1st (2nd) hardest jets have pT > 7 (6) GeV
151 pjet1 = cutJetsPHO[0];
152 pjet2 = cutJetsPHO[1];
153
154 if(!(pjet1.pT()>7 || pjet2.pT()>7) ){
155 vetoEvent;
156 }
157
158 ++nVeto3;
159
160 }
161
162 // Apply muon cuts in relevant regime
163 const Particles& all_muons = apply<FinalState>(event, "muons").particlesByPt();
164 Particles DIS_muons_cut;
165 Particles PHO_muons_cut;
166
167 // Apply muon eta / pT cuts
168 const Particles DIS_muons = select(all_muons, [](const Particle& m) {
169 return m.eta() > -0.75 && m.eta() < 1.15 && m.pT() > 2.5*GeV; });
170 const Particles PHO_muons = select(all_muons, [](const Particle& m) {
171 return m.eta() > -0.55 && m.eta() < 1.10 && m.pT() > 2.5*GeV; });
172
173 if (isDIS) {
174 // Veto event events with no muons
175 if (DIS_muons.size() == 0) vetoEvent;
176 ++nVeto4;
177 }
178 else {
179 // Veto event events with no muons
180 if (PHO_muons.size() == 0) vetoEvent;
181 ++nVeto5;
182 }
183
184 // Calculate fraction of photon energy entering the hard interaction observable
185 double sumJet1 = 0;
186 double sumJet2 = 0;
187 double sumAllHad = 0;
188 double x_gamma;
189
190 if (!isDIS) {
191 //const Particles& hadfs = apply<HadronicFinalState>(event, "HFS").particles();
192 //for (const Particle& h : hadfs) {
193 for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
194 Particle& h = particles[ip1];
195 FourMomentum hcmMom = hcmboost.transform(h.momentum());
196 // Need to change sign: by default hcmboost has gamma* in +z dir,
197 // and p in -z dir, but here we need: gamma* -in -z and proton in +z.
198 sumAllHad += (hcmMom.E() + hcmMom.pz());
199 if ( deltaR(h, pjet1) <= 1.0 ) {
200 sumJet1 += (hcmMom.E() + hcmMom.pz());
201 }
202 else if ( deltaR(h, pjet2) <= 1.0 ) {
203 sumJet2 += (hcmMom.E() + hcmMom.pz());
204 }
205 }
206
207 //sumJet1 = hcmboost.transform(pjet1).E() + hcmboost.transform(pjet1).pz();
208 //sumJet2 = hcmboost.transform(pjet2).E() + hcmboost.transform(pjet2).pz();
209
210 x_gamma = (sumJet1 + sumJet2)/sumAllHad;
211
212 }
213
214
215 // Fill histos
216 if (isDIS ) {
217
218 for (const Particle& m : DIS_muons) {
219 _h["DIS_pT_mu"] -> fill(m.pT());
220 _h["DIS_eta_mu"] -> fill(m.eta());
221 }
222
223 _h["DIS_pT_jet_Breit"] -> fill(cutJetsDIS[0].pT()/GeV);
224 _h["DIS_Q^2"] -> fill(Q2);
225 _h["DIS_x"] -> fill(log10(x));
226
227 } else {
228
229 for (const Particle& m : PHO_muons) {
230 _h["PHO_pT_mu"] -> fill(m.pT());
231 _h["PHO_eta_mu"] -> fill(m.eta());
232 }
233
234 _h["PHO_pT_jet"] -> fill(cutJetsPHO[0].pT());
235 _h["PHO_x_obs_gamma"] -> fill(x_gamma);
236
237 }
238
239 }
240
241
242 /// Normalise histograms etc., after the run
243 void finalize() {
244
245 double normpb = crossSection()/picobarn/sumW();
246
247 scale(_h["PHO_eta_mu"], normpb);
248 scale(_h["PHO_pT_mu"], normpb);
249 scale(_h["PHO_pT_jet"], normpb);
250 scale(_h["PHO_x_obs_gamma"], normpb);
251
252 scale(_h["DIS_Q^2"], normpb);
253 scale(_h["DIS_x"], normpb);
254 scale(_h["DIS_eta_mu"], normpb);
255 scale(_h["DIS_pT_mu"], normpb);
256 scale(_h["DIS_pT_jet_Breit"], normpb);
257
258 MSG_DEBUG("Events passing Q2/y cuts = " << nVeto0 );
259 MSG_DEBUG("PHO events = " << nPHO );
260 MSG_DEBUG("DIS events = " << nDIS );
261 MSG_DEBUG("DIS Events passing number of jets cuts= " << nVeto1 );
262 MSG_DEBUG("PHO Events passing number of jets cuts= " << nVeto2 );
263 MSG_DEBUG("PHO Events passing pT jet cuts= " << nVeto3 );
264 MSG_DEBUG("DIS Events passing one muon cut= " << nVeto4 );
265 MSG_DEBUG("PHO Events passing one muon cut= " << nVeto5 );
266 MSG_DEBUG("DIS Events passing muon in jet cut = " << nVeto6 );
267 MSG_DEBUG("PHO Events passing muon in jet cut = " << nVeto7 );
268
269 }
270
271 /// @}
272
273
274 private:
275
276 /// @name Histograms
277 /// @{
278 map<string, Histo1DPtr> _h;
279 map<string, Profile1DPtr> _p;
280 map<string, CounterPtr> _c;
281 /// @}
282
283 bool isDIS;
284 int nPHO, nDIS;
285 int nVeto0, nVeto1, nVeto2, nVeto3, nVeto4, nVeto5, nVeto6, nVeto7;
286
287 };
288
289
290 RIVET_DECLARE_PLUGIN(H1_2005_I676166);
291
292}
|