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      
 47      book(_h_eta1, {0., 0.25, 0.5, 1.}, {"d10-x01-y01", "d10-x01-y02", "d10-x01-y03"});
 48      book(_h_eta2, {1.5, 2.5, 4., 10.}, {"d11-x01-y01", "d11-x01-y02", "d11-x01-y03"});
 49      book(_h_zD1, {1.5, 2.5, 4., 10.}, {"d12-x01-y01", "d12-x01-y02", "d12-x01-y03"});
 50      _axes["Q2"]    = YODA::Axis<double>{1., 2.371, 4.218, 10.0, 17.78, 31.62,100.0};
 51      _axes["Q2eta"] = YODA::Axis<double>{1.0, 9.0, 20.0,100.0};
 52      _axes["Q2pt"]  = YODA::Axis<double>{1.0, 3.0, 8.0, 20.0, 100.0};
 53    }
 54
 55
 56    /// Perform the per-event analysis
 57    void analyze(const Event& event) {
 58      if (_edges.empty()) {
 59        _edges["Q2"   ] = _h_Q2->xEdges();
 60        _edges["Q2eta"] = _h_Q2eta->bins()[1]->xEdges();
 61        _edges["Q2pt" ] = _h_Q2pt ->bins()[1]->xEdges();
 62      }
 63      const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
 64
 65       // Q2 and inelasticity cuts
 66      if (!inRange(kin.Q2(), 1.0*GeV2, 100*GeV2)) vetoEvent;
 67      if (!inRange(kin.y(), 0.05, 0.7)) vetoEvent;
 68
 69        // D* reconstruction
 70      // const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV && Cuts::abseta < 1.5);
 71      const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV );
 72      const Particles dstars = select(unstables, [](const Particle& p){ return p.abspid() == PID::DSTARPLUS; });
 73      if (dstars.empty()) vetoEvent;
 74      // MSG_DEBUG("#D* = " << dstars.size());
 75      //const Particle& dstar = dstars.front();
 76      for (const Particle& dstar : dstars){
 77        const double zD = (dstar.E() - dstar.pz()) / (2*kin.beamLepton().E()*kin.y());
 78
 79
 80        // Single-differential histograms
 81        _h["p_tD*"]->fill(dstar.pT()/GeV);
 82        _h["etaD*"]->fill(dstar.eta());
 83        _h["Z_D*"]->fill(zD/GeV);
 84        
 85        size_t idx = _axes["Q2"].index(kin.Q2()/GeV2);
 86        string edge = "OTHER";
 87        if(idx && idx <= _edges["Q2"].size()) edge=_edges["Q2"][idx-1];
 88        _h_Q2->fill(edge);
 89        
 90        _h["W_GeV"]->fill(sqrt(kin.W2())/GeV);
 91        _h["logx"]->fill(log10(kin.x()));
 92
 93        // Double-differential (y,Q2) histograms
 94
 95        idx = _axes["Q2eta"].index(kin.Q2()/GeV2);
 96        edge = "OTHER";
 97        if(idx && idx <= _edges["Q2eta"].size()) edge=_edges["Q2eta"][idx-1];
 98        _h_Q2eta->fill(dstar.eta(),edge);
 99        idx = _axes["Q2pt"].index(kin.Q2()/GeV2);
100        edge = "OTHER";
101        if(idx && idx <= _edges["Q2pt"].size()) edge=_edges["Q2pt"][idx-1];
102        _h_Q2pt->fill(dstar.pT(),edge);
103        _h_eta1->fill(zD, dstar.eta());
104        _h_eta2->fill(dstar.pT(), dstar.eta());
105        _h_zD1->fill(dstar.pT(), zD/GeV);
106
107      }
108    }
109
110
111    /// Normalise histograms etc., after the run
112    void finalize() {
113
114
115     const double sf = crossSection()/nanobarn/sumOfWeights();
116      scale(_h["p_tD*"], sf);
117      scale(_h["etaD*"], sf);
118      scale(_h["Z_D*"], sf);
119      scale(_h_Q2     , sf);
120      for(auto & b : _h_Q2->bins()) {
121        const size_t idx=b.index();
122        b.scaleW(1./_axes["Q2"].width(idx));
123      }
124      scale(_h["W_GeV"], sf);
125      scale(_h["logx"], sf);
126
127      scale(_h_Q2eta, sf);
128      for(auto & histo : _h_Q2eta->bins()) {
129        for(auto & b : histo->bins()) {
130          const size_t idx=b.index();
131          b.scaleW(1./_axes["Q2eta"].width(idx));
132        }
133      }
134      scale(_h_Q2pt, sf);
135      for(auto & histo : _h_Q2pt->bins()) {
136        for(auto & b : histo->bins()) {
137          const size_t idx=b.index();
138          b.scaleW(1./_axes["Q2pt"].width(idx));
139        }
140      }
141      scale(_h_eta1, sf);
142      scale(_h_eta2, sf);
143      scale(_h_zD1, sf);
144      // distributionsa are really double differential divide by group width
145      _h_Q2eta->divByGroupWidth();
146      _h_Q2pt ->divByGroupWidth();
147      _h_eta1 ->divByGroupWidth();
148      _h_eta2 ->divByGroupWidth();
149      _h_zD1  ->divByGroupWidth();
150    }
151
152    ///@}
153
154
155    /// @name Histograms
156    ///@{
157    map<string, Histo1DPtr> _h;
158    map<string, Profile1DPtr> _p;
159    map<string, CounterPtr> _c;
160    HistoGroupPtr<double,string> _h_Q2eta, _h_Q2pt;
161    Histo1DGroupPtr _h_eta1,_h_eta2,_h_zD1;
162    BinnedHistoPtr<string> _h_Q2;
163    map<string,YODA::Axis<double> > _axes;
164    map<string,vector<string> > _edges;
165    ///@}
166
167
168  };
169
170
171  RIVET_DECLARE_PLUGIN(H1_2002_I561885);
172
173}