1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/ZFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Measurement of dimuon continuum in proton-nucleus collisions
 11  class E288_1981_I153009 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    E288_1981_I153009()
 16      : Analysis("E288_1981_I153009")
 17    {   }
 18
 19
 20    /// @name Analysis methods
 21    //@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27      const FinalState fs;
 28      declare(fs, "FS");
 29      Cut cut = Cuts::etaIn(-15.,15.);
 30      ZFinder zfinder(fs, cut, PID::MUON, 3.5*GeV, 30.0*GeV, 0.1, ZFinder::ClusterPhotons::NONE );
 31      declare(zfinder, "ZFinder");
 32
 33      // Book histograms
 34      // 400 GeV and y = 0.03
 35      Histo1DPtr dummy;
 36      _hist_pT_M_400.add(5., 6., book(dummy,9, 1, 1));
 37      _hist_pT_M_400.add(6., 7., book(dummy,9, 1, 2));
 38      _hist_pT_M_400.add(7., 8., book(dummy,9, 1, 3));
 39      _hist_pT_M_400.add(8., 9., book(dummy,9, 1, 4));
 40      _hist_pT_M_400.add(9., 10., book(dummy,9, 1, 5));
 41      _hist_pT_M_400.add(10.,11., book(dummy,9, 1, 6));
 42      _hist_pT_M_400.add(11.,12., book(dummy,9, 1, 7));
 43      _hist_pT_M_400.add(12.,13., book(dummy,9, 1, 8));
 44      _hist_pT_M_400.add(13.,14., book(dummy,9, 1, 9));
 45
 46      int Nbin = 50;
 47      book(_h_m_DiMuon ,"DiMuon_mass",Nbin,0.0,30.0);
 48      book(_h_pT_DiMuon,"DiMuon_pT",Nbin,0.0,20.0);
 49      book(_h_y_DiMuon,"DiMuon_y",Nbin,-8.0, 8.0);
 50      book(_h_xF_DiMuon,"DiMuon_xF",Nbin, -1.5,  1.5);
 51    }
 52
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56
 57      const double sqrts_tol = 10. ;
 58      if (!isCompatibleWithSqrtS(27.4, sqrts_tol)) {
 59        MSG_ERROR("Incorrect beam energy used: " << sqrtS()/GeV);
 60        throw Error("Unexpected sqrtS ! Only 27.4 GeV is supported");
 61      }
 62
 63      const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder");
 64      if (zfinder.particles().size() >= 1) {
 65
 66        double Zmass = zfinder.bosons()[0].momentum().mass()/GeV;
 67        double Zpt   = zfinder.bosons()[0].momentum().pT()/GeV;
 68        double Zpl   = zfinder.bosons()[0].momentum().pz()/GeV;
 69        double Zy    = zfinder.bosons()[0].momentum().rapidity();
 70        //double ZE    = zfinder.bosons()[0].momentum().E();
 71
 72        double xf = 2.*Zpl/sqrtS() ;
 73        _h_xF_DiMuon->fill(xf);
 74        _h_m_DiMuon->fill(Zmass/GeV);
 75        _h_pT_DiMuon->fill(Zpt);
 76        _h_y_DiMuon ->fill(Zy);
 77        double Zymin = -1.0;
 78        double Zymax = 1.03;
 79        double Z_y_width = Zymax - Zymin ;
 80        if ( Zy > Zymin && Zy < Zymax ) {
 81          // Edsigma^3/dp^3 = 2E/(pi*sqrts)dsigma/dx_F/dq_T^2 = 1/pi dsigma/dy/dq_T^2
 82          // normalisation of Zy bin width = Zwidth
 83          if (Zpt > 0) _hist_pT_M_400.fill(Zmass,Zpt, 1./2./Zpt/Z_y_width);
 84        }
 85      }
 86
 87    }
 88
 89
 90    /// Normalise histograms etc., after the run
 91    void finalize() {
 92      MSG_DEBUG("Generator cross section [pb] = " << crossSection()/picobarn);
 93      _hist_pT_M_400.scale(crossSection()/femtobarn/(sumOfWeights() * M_PI), this);
 94    }
 95
 96    //@}
 97
 98
 99    /// @name Histograms
100    //@{
101    BinnedHistogram _hist_pT_M_400,_hist_pT_M_300,_hist_pT_M_200;
102    Histo1DPtr _h_m_DiMuon ;
103    Histo1DPtr _h_pT_DiMuon;
104    Histo1DPtr _h_y_DiMuon;
105    Histo1DPtr _h_xF_DiMuon;
106    Histo1DPtr _h_XXXX, _h_YYYY, _h_ZZZZ;
107    Profile1DPtr _p_AAAA;
108    CounterPtr _c_BBBB;
109    //@}
110
111  };
112
113
114  RIVET_DECLARE_PLUGIN(E288_1981_I153009);
115
116}