1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/DileptonFinder.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// Drell-Yan dimuon absolute cross-sections in 800 GeV pp and pd collisions
 12  class NUSEA_2003_I613362 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(NUSEA_2003_I613362);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Projections
 26      DileptonFinder zfinder(91.2*GeV, 0.0, Cuts::abseta < 10. && Cuts::abspid == PID::MUON,
 27                                       Cuts::massIn(4.0*GeV, 100.0*GeV));
 28      declare(zfinder, "DileptonFinder");
 29
 30      // Booking histograms
 31      // hydrogen d01-d16
 32      book(_hist_M_xF, {-0.05, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8});
 33      for (auto& b : _hist_M_xF->bins()) book(b, b.index(), 1, 1);
 34
 35      // deuterium d17-d32
 36
 37      // hydrogen d40
 38      book(_hist_pT_M, {4.2, 5.2, 6.2, 7.2, 8.7, 10.85, 12.85});
 39      _hist_pT_M->maskBin(5); int idx = 0;
 40      for (auto& b : _hist_pT_M->bins()) book(b, 40, 1, ++idx);
 41
 42    }
 43
 44
 45    /// Perform the per-event analysis
 46    void analyze(const Event& event) {
 47
 48      const double sqrts_tol = 10.;
 49      if (!isCompatibleWithSqrtS(38.8*GeV, sqrts_tol)) {
 50        MSG_ERROR("Incorrect beam energy used: " << sqrtS()/GeV);
 51        throw Error("Unexpected sqrtS ! Only 38.8 GeV is supported");
 52      }
 53
 54      // Muons
 55      const DileptonFinder& zfinder = apply<DileptonFinder>(event, "DileptonFinder");
 56      if (zfinder.particles().size() <= 0) vetoEvent;
 57
 58      double Zmass = zfinder.bosons()[0].momentum().mass()/GeV;
 59      double Zpt   = zfinder.bosons()[0].momentum().pT()/GeV;
 60      double Zpl   = zfinder.bosons()[0].momentum().pz()/GeV;
 61      double ZE    = zfinder.bosons()[0].momentum().E();
 62      double xF = 2.*Zpl/sqrtS();
 63
 64      // Filling dimuon mass in bins of xF
 65      _hist_M_xF->fill(xF, Zmass/GeV, Zmass*sqr(Zmass));
 66
 67      // Filling pT in bins of Zmass
 68      if ( xF > -0.05 && xF <= 0.15 ) {
 69        // Include here all factors which are run-dependent for later scaling
 70        if (Zpt > 0) _hist_pT_M->fill(Zmass,Zpt, 1./2./Zpt*2.*ZE/sqrtS());
 71      }
 72
 73      MSG_DEBUG("Dimuon pT = "<< Zpt<<"   Dimuon E = ");
 74      MSG_DEBUG("DiMuon mass " << Zmass/GeV);
 75      MSG_DEBUG("DiMuon pT "<<Zpt);
 76    }
 77
 78
 79    /// Normalise histograms etc., after the run
 80    void finalize() {
 81      // xf bin width = 0.2, x-section in picobarn
 82      const double scalefactor=crossSection()/picobarn/(sumOfWeights() * M_PI *0.2 );
 83      scale(_hist_pT_M, scalefactor);
 84
 85      // x-section is quoted in nanobarn
 86      scale(_hist_M_xF, crossSection()/nanobarn/sumOfWeights());
 87
 88      divByGroupWidth({_hist_pT_M, _hist_M_xF});
 89    }
 90
 91    /// @}
 92
 93
 94  private:
 95
 96    /// @name Histograms
 97    /// @{
 98    Histo1DGroupPtr _hist_pT_M, _hist_M_xF;
 99    /// @}
100
101  };
102
103
104  RIVET_DECLARE_PLUGIN(NUSEA_2003_I613362);
105
106}