rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_2007_I736052

Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)
Experiment: H1 (HERA)
Inspire ID: 736052
Status: VALIDATED
Authors:
  • Maksim Davydov
  • Hannes Jung
References:
  • Eur.Phys.J.C 51 (2007) 271
  • DOI: 10.1140/epjc/s10052-007-0296-5
  • arXiv: hep-ex/0701023
  • DESY-06-240
Beams: e+ p+, p+ e+
Beam energies: (27.5, 920.0); (920.0, 27.5) GeV
    No run details listed

Inclusive $D^*$ production is measured in deep-inelastic $ep$ scattering at HERA with the H1 detector. In addition, the production of dijets in events with a $D^*$ meson is investigated. The analysis covers values of photon virtuality $2< Q^2 \leq 100$ GeV$^2$ and of inelasticity $0.05 \leq y \leq 0.7$. Differential cross sections are measured as a function of $Q^2$ and $x$ and of various $D^*$ meson and jet observables.

Source code: H1_2007_I736052.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/DISFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/DISKinematics.hh"
  7#include "Rivet/Projections/DISLepton.hh"
  8#include "Rivet/Projections/UnstableParticles.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)
 14  class H1_2007_I736052 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2007_I736052);
 19
 20
 21    /// @name Analysis methods
 22    ///@{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27        declare(DISLepton(), "Lepton");
 28        declare(DISKinematics(), "Kinematics");
 29
 30        const FinalState fs;
 31        declare(fs, "FS");
 32        const UnstableParticles ufs;
 33        declare(ufs, "UFS");
 34
 35        double jet_radius = 1.0;
 36        const DISFinalState DISfs(DISFrame::BREIT);
 37        declare(FastJets(DISfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme, jet_radius), "jets_fs");
 38
 39        Histo1DPtr tmp;
 40        book(_h["411"], 4, 1, 1);
 41        book(_h["511"], 5, 1, 1);
 42        book(_h["611"], 6, 1, 1);
 43        book(_h["711"], 7, 1, 1);
 44        book(_h["811"], 8, 1, 1);
 45        book(_h["911"], 9, 1, 1);
 46        book(_h_binned["Q2xbj"], {2., 4.22, 10., 17.8, 31.6, 100.},
 47                                 {"d10-x01-y01", "d11-x01-y01", "d12-x01-y01",
 48                                  "d13-x01-y01", "d14-x01-y01"});
 49        book(_h["1511"], 15, 1, 1);
 50        book(_h["1611"], 16, 1, 1);
 51        book(_h["1711"], 17, 1, 1);
 52        book(_h["1811"], 18, 1, 1);
 53        book(_h["1911"], 19, 1, 1);
 54        book(_h["2011"], 20, 1, 1);
 55        book(_h["2111"], 21, 1, 1);
 56        book(_h["2211"], 22, 1, 1);
 57
 58        book(_h_binned["Q2phi"], {2., 10., 100.}, {"d23-x01-y01", "d24-x01-y01"});
 59        book(_h["2511"], 25, 1, 1);
 60        book(_h["2611"], 26, 1, 1);
 61        book(_h["2711"], 27, 1, 1);
 62        book(_h["2811"], 28, 1, 1);
 63        book(_h_binned["Q2xgam"], {2., 5., 10., 100.}, {"d29-x01-y01", "d29-x01-y02", "d29-x01-y03"});
 64        book(_h["3011"], 30, 1, 1);
 65        book(_h_binned["Q2xglue"], {2., 5., 10., 100.}, {"d31-x01-y01", "d31-x01-y02", "d31-x01-y03"});
 66    }
 67
 68
 69    /// Perform the per-event analysis
 70    void analyze(const Event& event) {
 71
 72        const FinalState& fs = apply<FinalState>(event, "FS");
 73        const size_t numParticles = fs.particles().size();
 74
 75        Jets jets_fs = apply<JetFinder>(event, "jets_fs").jetsByPt(); // Jets with cut on eta
 76        double jet_radius = 1.0;
 77
 78        const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
 79
 80        if (numParticles < 2){
 81            MSG_DEBUG("Failed leptonic event cut");
 82            vetoEvent;
 83        }
 84
 85        Particles Dstar;
 86        for(const Particle& p : select(ufs.particles(), Cuts::pT > 1.5*GeV and Cuts::pT < 15*GeV and Cuts::abseta < 1.5 and Cuts::abspid==413)) {
 87            Dstar.push_back(p);
 88        }
 89        if(Dstar.size() == 0){ // Cut on Dstar
 90            MSG_DEBUG("Failed Dstar cut");
 91            vetoEvent;
 92        }
 93
 94        const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 95        const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 96
 97
 98        double Q2 = dk.Q2();
 99        double y  = dk.y();
100
101        if(y < 0.05 or y > 0.7 or Q2 < 2 or Q2 > 100){ // Cut on event kinematics
102            MSG_DEBUG("Failed kinematics cut");
103            vetoEvent;
104        }
105
106        // Extract the particles other than the lepton
107        Particles particles;
108        particles.reserve(fs.particles().size());
109        ConstGenParticlePtr dislepGP = dl.out().genParticle();
110        for (const Particle& p : fs.particles()) {
111           ConstGenParticlePtr loopGP = p.genParticle();
112           if (loopGP == dislepGP) continue;
113           particles.push_back(p);
114        }
115
116
117
118        const LorentzTransform hcmboost = dk.boostHCM(); // Hadron cm system
119        const LorentzTransform breitboost = dk.boostBreit(); //Breit system
120        const LorentzTransform labboost = breitboost.inverse(); //Labsystem from Breit
121
122        double xbj  = dk.x();
123        double W2 = dk.W2();
124        double pT_cm(0);
125
126        _h["411"] -> fill(Q2);
127        _h["511"] -> fill(xbj);
128        _h["611"] -> fill(std::sqrt(W2));
129
130        _h_binned["Q2xbj"]->fill(Q2, xbj);
131
132        for(const Particle& p : Dstar){
133            _h["711"] -> fill(p.pT());
134            _h["811"] -> fill(p.eta());
135            _h["911"] -> fill((p.E() - p.pz())/(2*y*dk.beamLepton().E()));
136
137            const FourMomentum hcmMom = hcmboost.transform(p.momentum());
138            if(pT_cm < hcmMom.pT()) pT_cm = hcmMom.pT();
139
140            if(hcmMom.pT() > 2){              //Cut for 1711 and 1811 histograms
141                _h["1711"] -> fill(p.pT());
142                _h["1811"] -> fill(p.eta());
143            }
144        }
145
146        if(pT_cm > 2){
147            _h["1511"] -> fill(Q2);
148            _h["1611"] -> fill(xbj);
149        }
150
151        bool two_jets_with_cut(false);
152
153        Jets jet_cut;
154
155        for(const Jet& j : jets_fs){
156            double etalab = labboost.transform(j.momentum()).eta() ;
157            if(j.momentum().Et() > 3 and etalab > -1 and etalab < 2.5){  // Cut on jet energy in Breit sys.
158                jet_cut.push_back(j);                            }
159        }
160
161        if(jet_cut.size() >= 2 && jet_cut[0].momentum().Et() > 4 ) two_jets_with_cut = true;
162
163        double delta_phi;
164        FourMomentum jet1;
165        FourMomentum jet2;
166
167        Jets jet_DJ, jet_OJ;
168        bool found_DJ(false), found_OJ(false);
169        //cout << " check D and other jets " << endl;
170        if (two_jets_with_cut){
171            jet1 = jet_cut[0].momentum(); // momentum of jet #1 in Breit sys.
172            jet2 = jet_cut[1].momentum(); // momentum of jet #2 in Breit sys.
173
174            delta_phi = deltaPhi(jet1,jet2)/degree ;
175            _h["1911"] -> fill(Q2);
176            _h["2011"] -> fill(xbj);
177            _h["2111"] -> fill(jet_cut[0].momentum().Et()/GeV);
178            _h["2211"] -> fill(FourMomentum(jet1+jet2).mass()/GeV);  // Jets invariant mass in Breit sys.
179            //cout << " dphi " << delta_phi <<" " << deltaPhi(jet_cut[0].momentum(),jet_cut[1].momentum())/degree <<  endl;
180            _h_binned["Q2phi"]->fill(Q2, delta_phi);
181            for (const Jet& jet : jet_cut) {
182                for(const Particle & p : Dstar) {
183                    if(deltaR(breitboost.transform(p.momentum()), jet.momentum()) < jet_radius  ) {
184                        jet_DJ.push_back(jet);
185                        found_DJ = true;
186                    }
187                    else {
188                    jet_OJ.push_back(jet);
189                    found_OJ = true;
190                    }
191                }
192
193               if( found_DJ &&  found_OJ ) break ;
194            }
195            if(jet_DJ.size()>0 && jet_OJ.size()>0 ) {
196                double eta_DJ_breit = jet_DJ[0].momentum().eta();
197                double eta_OJ_breit = jet_OJ[0].momentum().eta();
198
199                _h["2511"] -> fill(eta_DJ_breit);
200                _h["2611"] -> fill(eta_OJ_breit);
201                _h["2711"] -> fill(abs(eta_DJ_breit - eta_OJ_breit));
202
203                double x_gamma, x_gluon;
204                double E_star_p_z_star(0); // E() - pz() sum for for x_gamma in \gamma p cm
205                double E_star_p_z_had(0);
206
207                // for jets in had cms: boost first back to lab and then to hcm
208               Jet jet_DJ_hcm, jet_OJ_hcm;
209               jet_DJ_hcm = hcmboost.transform(labboost(jet_DJ[0].momentum()));
210               jet_OJ_hcm = hcmboost.transform(labboost(jet_OJ[0].momentum()));
211
212               // Need to change sign: by default hcmboost has gamma* in +z dir,
213               // and p in -z dir, but here we need: gamma* -in -z and proton in +z.
214               // so - -> +: watch out for E-pz -> E+pz  and exp(eta_OJ_hcm) -> exp(-eta_OJ_hcm)
215               // this was already noted in  H1_2007_I746380
216
217               double Et_DJ_hcm = jet_DJ_hcm.Et();
218               double eta_DJ_hcm = jet_DJ_hcm.eta();
219               double Et_OJ_hcm = jet_OJ_hcm.Et();
220               double eta_OJ_hcm = jet_OJ_hcm.eta();
221
222                //observed fraction of the photon momentum carried by the parton involved in the hard subprocess
223
224                E_star_p_z_star = jet_DJ_hcm.momentum().E() + jet_DJ_hcm.momentum().pz() + jet_OJ_hcm.momentum().E() + jet_OJ_hcm.momentum().pz();
225
226                for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
227                  const Particle& p = particles[ip1];
228                  E_star_p_z_had += (hcmboost.transform(p.momentum())).E() + (hcmboost.transform(p.momentum())).pz();
229                }
230                x_gamma = E_star_p_z_star/E_star_p_z_had;
231                x_gluon = (Et_OJ_hcm*exp(-eta_OJ_hcm) + Et_DJ_hcm*exp(-eta_DJ_hcm))/(2.*hcmboost.transform(dk.beamHadron().momentum()).E());
232                //observed fraction of the proton momentum carried by the gluon
233
234                _h["2811"]->fill(x_gamma);
235                _h_binned["Q2xgam"]->fill(Q2, x_gamma);
236
237                _h["3011"]->fill(log10(x_gluon));
238                _h_binned["Q2xglue"]->fill(Q2, log10(x_gluon));
239            }
240        }
241    }
242
243
244    /// Normalise histograms etc., after the run
245    void finalize() {
246
247        const double norm = crossSection()/nanobarn/sumW();
248        scale(_h, norm);
249        scale(_h_binned["Q2xbj"], norm);
250        scale(_h_binned["Q2phi"], norm*180./M_PI);
251        scale(_h_binned["Q2xgam"], norm);
252        scale(_h_binned["Q2xglue"], norm);
253
254        // new scaling needed, since x bins are in log10(x)
255        for (auto& b : _h["3011"]->bins()) {
256           const double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
257           // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
258           const double factor = pow(10,b.xMid())/2.3026;
259           b.scaleW(scale_new/factor) ;
260        }
261
262        // new scaling needed, since x bins are in log10(x)
263        for (auto& histo : _h_binned["Q2xglue"]->bins()) {
264           for (auto& b : histo->bins()) {
265              const double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
266              // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
267              const double factor = pow(10,b.xMid())/2.3026;
268              b.scaleW(scale_new/factor) ;
269           }
270        }
271
272    }
273
274    ///@}
275
276
277    /// @name Histograms
278    ///@{
279    map<string, Histo1DPtr> _h;
280    map<string, Histo1DGroupPtr> _h_binned;
281    ///@}
282
283
284  };
285
286
287  RIVET_DECLARE_PLUGIN(H1_2007_I736052);
288
289}