1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// STAR identified hadron spectra in pp at 200 GeV
 10  class STAR_2006_I709170 : public Analysis {
 11  public:
 12
 13    RIVET_DEFAULT_ANALYSIS_CTOR(STAR_2006_I709170);
 14
 15
 16    /// Book projections and histograms
 17    void init() {
 18      ChargedFinalState bbc1(Cuts::etaIn(-5.0,-3.3)); // beam-beam-counter trigger
 19      ChargedFinalState bbc2(Cuts::etaIn( 3.3, 5.0)); // beam-beam-counter trigger
 20      declare(bbc1, "BBC1");
 21      declare(bbc2, "BBC2");
 22
 23      IdentifiedFinalState pionfs(Cuts::abseta < 2.5 && Cuts::pT > 0.3*GeV);
 24      IdentifiedFinalState protonfs(Cuts::abseta < 2.5 && Cuts::pT > 0.4*GeV);
 25      pionfs.acceptIdPair(PID::PIPLUS);
 26      protonfs.acceptIdPair(PID::PROTON);
 27      declare(pionfs, "PionFS");
 28      declare(protonfs, "ProtonFS");
 29
 30      book(_h_pT_piplus     ,2, 1, 1); // full range pion binning
 31      book(_h_pT_piminus    ,7, 1, 1); // full range pion binning
 32      book(_tmp_pT_piplus   ,"TMP/pT_piplus" , refData(12, 1, 1)); // pi histo compatible with more restricted proton binning
 33      book(_tmp_pT_piminus  ,"TMP/pT_piminus", refData(12, 1, 1)); // pi histo compatible with more restricted proton binning
 34      book(_h_pT_proton     ,12, 1, 1);
 35      book(_h_pT_antiproton ,17, 1, 1);
 36
 37      book(_sumWeightSelected, "_sumWeightSelected");
 38    }
 39
 40
 41    /// Do the analysis
 42    void analyze(const Event& event) {
 43      const ChargedFinalState& bbc1 = apply<ChargedFinalState>(event, "BBC1");
 44      const ChargedFinalState& bbc2 = apply<ChargedFinalState>(event, "BBC2");
 45      if (bbc1.size() < 1 || bbc2.size() < 1) {
 46        MSG_DEBUG("Failed beam-beam-counter trigger");
 47        vetoEvent;
 48      }
 49
 50      const IdentifiedFinalState& pionfs = apply<IdentifiedFinalState>(event, "PionFS");
 51      for (const Particle& p : pionfs.particles()) {
 52        if (p.absrap() < 0.5) {
 53          /// @todo Use a binned counter to avoid this bin width cancellation hack
 54          const double pT = p.pT() / GeV;
 55          ((p.pid() > 0) ? _h_pT_piplus : _h_pT_piminus)->fill(pT, 1.0/pT);
 56          ((p.pid() > 0) ? _tmp_pT_piplus : _tmp_pT_piminus)->fill(pT, 1.0/pT);
 57        }
 58      }
 59
 60      const IdentifiedFinalState& protonfs = apply<IdentifiedFinalState>(event, "ProtonFS");
 61      for (const Particle& p : protonfs.particles()) {
 62        if (p.absrap() < 0.5) {
 63          /// @todo Use a binned counter to avoid this bin width cancellation hack
 64          const double pT = p.pT() / GeV;
 65          ((p.pid() > 0) ? _h_pT_proton : _h_pT_antiproton)->fill(pT, 1.0/pT);
 66        }
 67      }
 68      _sumWeightSelected->fill();
 69    }
 70
 71
 72    /// Finalize
 73    void finalize() {
 74      
 75      Estimate1DPtr _s_piminus_piplus;
 76      book(_s_piminus_piplus,"/TMP/s_piminus_piplus",refData(2,1,1).xEdges());
 77      divide(_h_pT_piminus, _h_pT_piplus, _s_piminus_piplus);
 78      BinnedEstimatePtr<string> _r_piminus_piplus;
 79      book(_r_piminus_piplus,23,1,2);
 80      for(const auto & b : _s_piminus_piplus->bins()) {
 81        const size_t idx = b.index();
 82        _r_piminus_piplus->bin(idx).setVal(b.val());
 83        for(auto & item : b.sources())
 84          _r_piminus_piplus->bin(idx).setErr(b.err(item),item);
 85      }
 86      Estimate1DPtr _s_antipr_pr;
 87      book(_s_antipr_pr,"/TMP/s_antipr_pr",refData(12,1,1).xEdges());
 88      divide(_h_pT_antiproton, _h_pT_proton, _s_antipr_pr);
 89      BinnedEstimatePtr<string> _r_antipr_pr;
 90      book(_r_antipr_pr,24,1,2);
 91      for(const auto & b : _s_antipr_pr->bins()) {
 92        const size_t idx = b.index();
 93        _r_antipr_pr->bin(idx).setVal(b.val());
 94        for(auto & item : b.sources())
 95          _r_antipr_pr->bin(idx).setErr(b.err(item),item);
 96      }
 97      
 98      Estimate1DPtr _s_pr_piplus;
 99      book(_s_pr_piplus,"/TMP/s_pr_piplus",refData(12,1,1).xEdges());
100      divide(_h_pT_proton, _tmp_pT_piplus, _s_pr_piplus);
101      BinnedEstimatePtr<string> _r_pr_piplus;
102      book(_r_pr_piplus,25,1,2);
103      for(const auto & b : _s_pr_piplus->bins()) {
104        const size_t idx = b.index();
105        _r_pr_piplus->bin(idx).setVal(b.val());
106        for(auto & item : b.sources())
107          _r_pr_piplus->bin(idx).setErr(b.err(item),item);
108      }
109      
110      Estimate1DPtr _s_antipr_piminus;
111      book(_s_antipr_piminus,"/TMP/s_antipr_piminus",refData(12,1,1).xEdges());
112      divide(_h_pT_antiproton, _tmp_pT_piminus, _s_antipr_piminus);
113      BinnedEstimatePtr<string> _r_antipr_piminus;
114      book(_r_antipr_piminus,26,1,2);
115      for(const auto & b : _s_antipr_piminus->bins()) {
116        const size_t idx = b.index();
117        _r_antipr_piminus->bin(idx).setVal(b.val());
118        for(auto & item : b.sources())
119          _r_antipr_piminus->bin(idx).setErr(b.err(item),item);
120      }
121      
122      const double factor = ((1/(2*M_PI)) / _sumWeightSelected->val());
123      scale(_h_pT_piplus,     factor);
124      scale(_h_pT_piminus,    factor);
125      scale(_h_pT_proton,     factor);
126      scale(_h_pT_antiproton, factor);
127    }
128
129
130  private:
131
132    /// @{
133    CounterPtr _sumWeightSelected;
134    Histo1DPtr _h_pT_piplus, _h_pT_piminus, _h_pT_proton, _h_pT_antiproton;
135    Histo1DPtr _tmp_pT_piplus, _tmp_pT_piminus;
136    /// @}
137
138  };
139
140
141
142  RIVET_DECLARE_ALIASED_PLUGIN(STAR_2006_I709170, STAR_2006_S6500200);
143
144}