1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/DileptonFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Dimuon production in proton-copper collisions at 38.8 GeV
 11  class E605_1991_I302822 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(E605_1991_I302822);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      Cut cut = Cuts::abseta < 10.;
 26      DileptonFinder zfinder(91.2*GeV, 0.0, cut && Cuts::abspid == PID::MUON,
 27                             Cuts::massIn(4.0*GeV, 100.0*GeV), LeptonOrigin::PROMPT);
 28      declare(zfinder, "DileptonFinder");
 29
 30      // Book histograms in mass ranges (measurement is not normalised to mass range)
 31      book(_h["pT_M_78"],   17, 1, 1);
 32      book(_h["pT_M_89"],   18, 1, 1);
 33      book(_h["pT_M_1011"], 19, 1, 1);
 34      book(_h["pT_M_1113"], 20, 1, 1);
 35      book(_h["pT_M_1318"], 21, 1, 1);
 36
 37      int Nbin = 50;
 38      book(_h_m_DiMuon, "DiMuon_mass",Nbin, 0.0,30.0);
 39      book(_h_pT_DiMuon,"DiMuon_pT",  Nbin, 0.0,20.0);
 40      book(_h_y_DiMuon, "DiMuon_y",   Nbin,-8.0, 8.0);
 41      book(_h_xF_DiMuon,"DiMuon_xF",  Nbin,-1.5, 1.5);
 42
 43      _axis = YODA::Axis<double>(18, 0.0, 3.6);
 44    }
 45
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49
 50      if (_edges.empty()) {
 51        for (string tag : vector<string>{"pT_M_78", "pT_M_89", "pT_M_1011", "pT_M_1113", "pT_M_1318"}) {
 52          _edges.insert({tag, _h[tag]->xEdges()});
 53        }
 54      }
 55
 56      const double sqrts_tol = 10. ;
 57      if (!isCompatibleWithSqrtS(38.8*GeV, sqrts_tol)) {
 58        MSG_ERROR("Incorrect beam energy used: " << sqrtS()/GeV);
 59        throw Error("Unexpected sqrtS ! Only 38.8 GeV is supported");
 60      }
 61
 62      const DileptonFinder& zfinder = apply<DileptonFinder>(event, "DileptonFinder");
 63      if (zfinder.particles().size() >= 1) {
 64
 65        double Zmass = zfinder.bosons()[0].momentum().mass()/GeV;
 66        double Zpt   = zfinder.bosons()[0].momentum().pT()/GeV;
 67        double Zpl   = zfinder.bosons()[0].momentum().pz()/GeV;
 68        double Zy    = zfinder.bosons()[0].momentum().rapidity();
 69        double ZE    = zfinder.bosons()[0].momentum().E();
 70
 71        double xf = 2.*Zpl/sqrtS();
 72        _h_xF_DiMuon->fill(xf);
 73        _h_m_DiMuon->fill(Zmass/GeV);
 74        _h_pT_DiMuon->fill(Zpt);
 75        _h_y_DiMuon ->fill(Zy);
 76
 77        if (xf > -0.1 && xf < 0.2 && Zpt > 0) {
 78          if (Zmass > 7. && Zmass < 8.)      discfill("pT_M_78",   Zpt, 1./2./Zpt *2.*ZE/sqrtS());
 79          if (Zmass > 8. && Zmass < 9.)      discfill("pT_M_89",   Zpt, 1./2./Zpt *2.*ZE/sqrtS());
 80          if (Zmass > 10.5 && Zmass < 11.5)  discfill("pT_M_1011", Zpt, 1./2./Zpt *2.*ZE/sqrtS());
 81          if (Zmass > 11.5 && Zmass < 13.5)  discfill("pT_M_1113", Zpt, 1./2./Zpt *2.*ZE/sqrtS());
 82          if (Zmass > 13.5 && Zmass < 18.)   discfill("pT_M_1318", Zpt, 1./2./Zpt *2.*ZE/sqrtS());
 83        }
 84      }
 85
 86    }
 87
 88
 89    /// Normalise histograms etc., after the run
 90    void finalize() {
 91      MSG_DEBUG(" Generator cross section [pb] " << crossSection()/picobarn);
 92
 93      normalize(_h_m_DiMuon);
 94      normalize(_h_pT_DiMuon);
 95      normalize(_h_xF_DiMuon);
 96      normalize(_h_y_DiMuon);
 97
 98      // normalisation of xF bin width = 0.3
 99      const double scalefactor=crossSection()/picobarn/(sumOfWeights() * M_PI *0.3 );
100      scale(_h, scalefactor/0.2); // divide by bin width
101    }
102
103    /// @}
104
105    void discfill(const string& tag, const double value, const double weight) {
106      string edge = "OTHER";
107      const size_t idx = _axis.index(value);
108      if (idx && idx <= _edges.at(tag).size()) {
109        edge = _edges.at(tag)[idx-1];
110      }
111      _h[tag]->fill(edge, weight);
112    }
113
114
115  private:
116
117    /// @name Histograms
118    /// @{
119    Histo1DPtr _h_m_DiMuon;
120    Histo1DPtr _h_pT_DiMuon;
121    Histo1DPtr _h_y_DiMuon;
122    Histo1DPtr _h_xF_DiMuon;
123    map<string, BinnedHistoPtr<string>> _h;
124    map<string, vector<string>> _edges;
125    YODA::Axis<double> _axis;
126    /// @}
127
128  };
129
130
131  RIVET_DECLARE_PLUGIN(E605_1991_I302822);
132
133}