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        _h_Q2eta.add(-1.5,-0.5, book(dummy,8,1,1)); 
 45        _h_Q2eta.add(-0.5,0.5, book(dummy,8,1,2)); 
 46        _h_Q2eta.add(0.5,1.5, book(dummy,8,1,3)); 
 47
 48
 49        _h_Q2pt.add(1.5,4., book(dummy,9,1,1)); 
 50        _h_Q2pt.add(4.,10., book(dummy,9,1,2)); 
 51
 52
 53        _h_eta1.add(0,0.25, book(dummy,10,1,1));
 54        _h_eta1.add(0.25,0.50, book(dummy,10,1,2));
 55        _h_eta1.add(0.5,1, book(dummy,10,1,3));
 56         
 57         
 58        _h_eta2.add(1.5,2.5, book(dummy,11,1,1));
 59        _h_eta2.add(2.5,4.0, book(dummy,11,1,2));
 60        _h_eta2.add(4.0,10.0,book(dummy,11,1,3));
 61         
 62         
 63        _h_zD1.add(1.5,2.5, book(dummy,12,1,1));
 64        _h_zD1.add(2.5,4.0, book(dummy,12,1,2));
 65        _h_zD1.add(4.0,10.0, book(dummy,12,1,3));
 66                 
 67    }
 68
 69
 70    /// Perform the per-event analysis
 71    void analyze(const Event& event) {
 72          
 73      const DISKinematics& kin = applyProjection<DISKinematics>(event, "Kinematics");
 74      
 75       // Q2 and inelasticity cuts
 76      if (!inRange(kin.Q2(), 1.0*GeV2, 100*GeV2)) vetoEvent;
 77      if (!inRange(kin.y(), 0.05, 0.7)) vetoEvent;
 78      
 79        // D* reconstruction
 80      // const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV && Cuts::abseta < 1.5);
 81      const Particles unstables = apply<ParticleFinder>(event, "Dstars").particles(Cuts::pT > 1.5*GeV );
 82      const Particles dstars = filter_select(unstables, [](const Particle& p){ return p.abspid() == PID::DSTARPLUS; });
 83      if (dstars.empty()) vetoEvent;
 84      // MSG_DEBUG("#D* = " << dstars.size());
 85      //const Particle& dstar = dstars.front();
 86      for(const Particle& dstar : dstars){
 87      const double zD = (dstar.E() - dstar.pz()) / (2*kin.beamLepton().E()*kin.y());
 88  
 89  
 90  // Single-differential histograms
 91    //     cout << " ' D* found "<< dstar.pid() << " " <<  kin.Q2() << " " << kin.x() << " " << dstar.pT() << " " << dstar.eta() << " " << dstar.rapidity() << endl;
 92        _h["p_tD*"]->fill(dstar.pT()/GeV);
 93        _h["etaD*"]->fill(dstar.eta());
 94        _h["Z_D*"]->fill(zD/GeV);
 95        _h["Q2"]->fill(kin.Q2()/GeV2);
 96        _h["W_GeV"]->fill(sqrt(kin.W2())/GeV);
 97        _h["logx"]->fill(log10(kin.x()));
 98  // Double-differential (y,Q2) histograms
 99
100        _h_Q2eta.fill(dstar.eta(),kin.Q2());
101        _h_Q2pt.fill(dstar.pT(),kin.Q2());
102
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      scale(_h["W_GeV"], sf);
121      scale(_h["logx"], sf);
122      
123      _h_Q2eta.scale(sf, this );
124      _h_Q2pt.scale(sf, this );
125      _h_eta1.scale(sf, this );
126      _h_eta2.scale(sf, this );
127      _h_zD1.scale(sf, this );
128        
129
130    }
131
132    ///@}
133
134
135    /// @name Histograms
136    ///@{
137    map<string, Histo1DPtr> _h;
138    map<string, Profile1DPtr> _p;
139    map<string, CounterPtr> _c;
140    BinnedHistogram _h_eta1,_h_eta2,_h_zD1, _h_Q2eta, _h_Q2pt;
141    ///@}
142
143
144  };
145
146
147  RIVET_DECLARE_PLUGIN(H1_2002_I561885);
148
149}