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/UnstableParticles.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Measurement of D*+- meson production in deep inelastic scattering at HERA (H1)
 12  class H1_2002_I561885 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2002_I561885);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25
 26      declare(DISKinematics(), "Kinematics");
 27      declare(UnstableParticles(), "Dstars");
 28      //Cuts::abspid == PID::DSTARPLUS
 29
 30      // Initialise and register projections
 31
 32
 33      Histo1DPtr dummy; //Introducing
 34
 35
 36      // Book histograms
 37      book(_h["W_GeV"], 2, 1, 1);
 38      book(_h["p_tD*"], 3, 1, 1);
 39      book(_h["logx"], 4, 1, 1);
 40      book(_h["etaD*"], 5, 1, 1);
 41      book(_h["Q2"], 6, 1, 1);
 42      book(_h["Z_D*"], 7, 1, 1);
 43
 44      book(_h_Q2eta, {-1.5, -0.5, 0.5, 1.5}, {"d08-x01-y01", "d08-x01-y02", "d08-x01-y03"});
 45      book(_h_Q2pt, {1.5, 4., 10.}, {"d09-x01-y01", "d09-x01-y02"});
 46      book(_h_eta1, {0., 0.25, 0.5, 1.}, {"d10-x01-y01", "d10-x01-y02", "d10-x01-y03"});
 47      book(_h_eta2, {1.5, 2.5, 4., 10.}, {"d11-x01-y01", "d11-x01-y02", "d11-x01-y03"});
 48      book(_h_zD1, {1.5, 2.5, 4., 10.}, {"d12-x01-y01", "d12-x01-y02", "d12-x01-y03"});
 49
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
 57
 58       // Q2 and inelasticity cuts
 59      if (!inRange(kin.Q2(), 1.0*GeV2, 100*GeV2)) vetoEvent;
 60      if (!inRange(kin.y(), 0.05, 0.7)) vetoEvent;
 61
 62        // D* reconstruction
 63      // const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV && Cuts::abseta < 1.5);
 64      const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV );
 65      const Particles dstars = select(unstables, [](const Particle& p){ return p.abspid() == PID::DSTARPLUS; });
 66      if (dstars.empty()) vetoEvent;
 67      // MSG_DEBUG("#D* = " << dstars.size());
 68      //const Particle& dstar = dstars.front();
 69      for (const Particle& dstar : dstars){
 70        const double zD = (dstar.E() - dstar.pz()) / (2*kin.beamLepton().E()*kin.y());
 71
 72
 73        // Single-differential histograms
 74        _h["p_tD*"]->fill(dstar.pT()/GeV);
 75        _h["etaD*"]->fill(dstar.eta());
 76        _h["Z_D*"]->fill(zD/GeV);
 77        _h["Q2"]->fill(kin.Q2()/GeV2);
 78        _h["W_GeV"]->fill(sqrt(kin.W2())/GeV);
 79        _h["logx"]->fill(log10(kin.x()));
 80
 81        // Double-differential (y,Q2) histograms
 82        _h_Q2eta->fill(dstar.eta(),kin.Q2());
 83        _h_Q2pt->fill(dstar.pT(),kin.Q2());
 84        _h_eta1->fill(zD, dstar.eta());
 85        _h_eta2->fill(dstar.pT(), dstar.eta());
 86        _h_zD1->fill(dstar.pT(), zD/GeV);
 87
 88      }
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94
 95
 96     const double sf = crossSection()/nanobarn/sumOfWeights();
 97      scale(_h["p_tD*"], sf);
 98      scale(_h["etaD*"], sf);
 99      scale(_h["Z_D*"], sf);
100      scale(_h["Q2"], sf);
101      scale(_h["W_GeV"], sf);
102      scale(_h["logx"], sf);
103
104      scale(_h_Q2eta, sf);
105      scale(_h_Q2pt, sf);
106      scale(_h_eta1, sf);
107      scale(_h_eta2, sf);
108      scale(_h_zD1, sf);
109
110
111    }
112
113    ///@}
114
115
116    /// @name Histograms
117    ///@{
118    map<string, Histo1DPtr> _h;
119    map<string, Profile1DPtr> _p;
120    map<string, CounterPtr> _c;
121    Histo1DGroupPtr _h_eta1,_h_eta2,_h_zD1, _h_Q2eta, _h_Q2pt;
122    ///@}
123
124
125  };
126
127
128  RIVET_DECLARE_PLUGIN(H1_2002_I561885);
129
130}