1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/FastJets.hh"
 5#include "Rivet/Projections/IdentifiedFinalState.hh"
 6#include "Rivet/Projections/MissingMomentum.hh"
 7#include "Rivet/Projections/DileptonFinder.hh"
 8
 9namespace Rivet {
10
11
12  /// Drell Yan production at low masses at $\sqrt{s} = 200$ GeV
13  class PHENIX_2019_I1672015 : public Analysis {
14  public:
15
16    /// Constructor
17    RIVET_DEFAULT_ANALYSIS_CTOR(PHENIX_2019_I1672015);
18
19    /// @name Analysis methods
20    /// @{
21
22    /// Book histograms and initialise projections before the run
23    void init() {
24
25      // Initialise and register projections
26      Cut cut = Cuts::etaIn(-10.,10.);
27      DileptonFinder zfinder(91.2*GeV, 0.0, cut && Cuts::abspid == PID::MUON,
28                             Cuts::massIn(4.0*GeV, 100.0*GeV), LeptonOrigin::PROMPT);
29      declare(zfinder, "DileptonFinder");
30
31      // Book histograms
32      book(_h_pT ,1, 1, 1);
33      book(_h_mass,2, 1, 1);
34      int Nbin = 50;
35      book(_h_m_DiMuon,"DiMuon_mass",Nbin,0.0,30.0);
36      book(_h_pT_DiMuon,"DiMuon_pT",Nbin,0.0,20.0);
37      book(_h_y_DiMuon,"DiMuon_y",Nbin,-8.0, 8.0);
38
39    }
40
41
42    /// Perform the per-event analysis
43    void analyze(const Event& event) {
44
45      double sqrts_tol = 10. ;
46      if (!isCompatibleWithSqrtS(200*GeV, sqrts_tol)) {
47        MSG_ERROR("Incorrect beam energy used: " << sqrtS()/GeV);
48        throw Error("Unexpected sqrtS ! Only 200 GeV is supported");
49      }
50      const DileptonFinder& zfinder = apply<DileptonFinder>(event, "DileptonFinder");
51      if (zfinder.particles().size() < 1) vetoEvent;
52      double mass = zfinder.bosons()[0].momentum().mass()/GeV;
53      double pt_DY   = zfinder.bosons()[0].momentum().pT()/GeV;
54      double y_DY    = abs(zfinder.bosons()[0].momentum().rapidity());
55
56      _h_m_DiMuon->fill(mass/GeV);
57      _h_pT_DiMuon->fill(pt_DY);
58      _h_y_DiMuon ->fill(y_DY);
59      if (inRange( y_DY, 1.2, 2.2) && inRange(mass, 4.8, 8.2) && pt_DY > 0)
60        _h_pT->fill(pt_DY, 1./2./pt_DY/M_PI);
61      if ( y_DY >= 1.2 && y_DY <= 2.2)
62        _h_mass->fill(mass);
63
64    }
65
66
67    /// Normalise histograms etc., after the run
68    void finalize() {
69
70      // Divide by 2 because of rapidity range
71      normalize(_h_m_DiMuon);
72      normalize(_h_pT_DiMuon);
73      normalize(_h_y_DiMuon);
74      scale(_h_pT, crossSection()/picobarn/(sumOfWeights())/2.);
75      scale(_h_mass, crossSection()/picobarn/(sumOfWeights())/2.);
76
77    }
78
79    /// @}
80
81
82    /// @name Histograms
83    /// @{
84    Histo1DPtr _h_pT, _h_y, _h_mass, _h_ZZZZ;
85    Histo1DPtr _h_m_DiMuon;
86    Histo1DPtr _h_pT_DiMuon;
87    Histo1DPtr _h_y_DiMuon;
88    Histo1DPtr _h_xF_DiMuon;
89    /// @}
90
91  };
92
93
94  RIVET_DECLARE_PLUGIN(PHENIX_2019_I1672015);
95
96}