rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_2005_I676166

Measurement of beauty production at HERA using events with muons and jets (H1)
Experiment: H1 (HERA)
Inspire ID: 676166
Status: VALIDATED
Authors:
  • Danielle Wislon
  • Hannes Jung <
References:
  • Eur.Phys.J.C 41 (2005) 453
  • DOI:10.1140/epjc/s2005-02267-0
  • arXiv: hep-ex/0502010
  • DESY-05-004
Beams: e+ p+, p+ e+
Beam energies: (27.6, 920.0); (920.0, 27.6) GeV
Run details:
  • validation with 3000001 events generated RAPGAP33, 26.7x820 GeV, Q2>2, 0.01 < y <0.95, IPRO=12, pt2cut=9, mixing with O(alphas) process, Nf_QCDC=4, NFLA=5

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}