1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/Thrust.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/UnstableParticles.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Add a short analysis description here
 12  class HRS_1990_I280958 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1990_I280958);
 17
 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      declare(Beam(), "Beams");
 27      declare(UnstableParticles(), "UFS");
 28      const ChargedFinalState cfs;
 29      declare(cfs, "CFS");
 30      declare(Thrust(cfs), "Thrust");
 31
 32      // Book histograms
 33      book(_h["X"], 3, 1, 1);
 34      book(_h["rap_all"], 4, 1, 1);
 35      book(_h["rap_charm"], 5, 1, 1);
 36      book(_h["rap_light"], 6, 1, 1);
 37      book(_wLight,"TMP/wLight");
 38      book(_wCharm,"TMP/wCharm");
 39
 40      _axes["X"] = YODA::Axis<double>({0.0, 0.1, 0.14, 0.18, 0.22, 0.26, 0.30,
 41                                       0.34, 0.38, 0.42, 0.46, 0.5, 0.6, 0.7});
 42      _axes["rap_all"] = YODA::Axis<double>({0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75,
 43                                             2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75});
 44      _axes["rap_charm"] = YODA::Axis<double>({-3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5,
 45                                               0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5});
 46      _axes["rap_light"] = _axes["rap_charm"];
 47
 48    }
 49
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53      if (_edges.empty()) {
 54        for (const auto& item : _h) {
 55          _edges[item.first] = item.second->xEdges();
 56        }
 57      }
 58      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 59      int nch = cfs.particles().size();
 60      if (nch<5) vetoEvent;
 61      // Get beams and average beam momentum
 62      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 63      const double meanBeamMom = ( beams.first.p3().mod() +
 64                                   beams.second.p3().mod() ) / 2.0;
 65      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 66      // get the thrust axes
 67      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 68      const Vector3 & axis = thrust.thrustAxis();
 69      // unstable particles
 70      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 71      Particle pTag;
 72      // get the tags
 73      Particles Dstar = ufs.particles(Cuts::abspid==413);
 74      bool charmTagged = !Dstar.empty();
 75      if (charmTagged) {
 76        pTag = Dstar[0];
 77        for (const Particle & p : Dstar) {
 78          if (p.E()>pTag.E()) pTag=p;
 79        }
 80        _wCharm->fill();
 81      }
 82      bool lightTagged = false;
 83      if (!charmTagged) {
 84        for (const Particle& p : cfs.particles()) {
 85          if (p.p3().mod()>9.43*GeV) {
 86            pTag=p;
 87            lightTagged=true;
 88            _wLight->fill();
 89            break;
 90          }
 91        }
 92      }
 93      // sign of hemispheres if tagged
 94      double sign=1.;
 95      if(charmTagged || lightTagged) {
 96        if(dot(axis,pTag.p3())<0.) sign=-1.;
 97      }
 98      // now loop over the kaons
 99      for (const Particle & p : ufs.particles(Cuts::pid==130 || Cuts::pid==310)) {
100         double xE = p.E()/meanBeamMom;
101         const double energy = p.E();
102         const double momT = dot(axis, p.p3());
103         discfill("X", xE);
104         double rap = 0.5 * std::log((energy + momT) / (energy - momT));
105         discfill("rap_all", fabs(rap));
106         rap *= sign;
107         if (charmTagged) {
108           discfill("rap_charm", rap);
109         }
110         else if (lightTagged) {
111           discfill("rap_light", rap);
112         }
113      }
114    }
115
116    void discfill(const string& name, const double value) {
117      string edge = "OTHER";
118      const size_t idx = _axes[name].index(value);
119      if (idx && idx <= _edges[name].size())  edge = _edges[name][idx-1];
120      _h[name]->fill(edge);
121    }
122
123
124    /// Normalise histograms etc., after the run
125    void finalize() {
126      scale(_h["X"], crossSection()*sqr(sqrtS())/nanobarn/sumOfWeights());
127      scale(_h["rap_all"], 1./sumOfWeights());
128      scale(_h["rap_light"], 1./ *_wLight);
129      scale(_h["rap_charm"], 1./ *_wCharm);
130      for( auto & hist : _h) {
131        for(auto & b: hist.second->bins()) {
132          const size_t idx = b.index();
133          b.scaleW(1./_axes[hist.first].width(idx));
134        }
135      }
136    }
137
138    /// @}
139
140
141    /// @name Histograms
142    /// @{
143    CounterPtr _wLight, _wCharm;
144    map<string, BinnedHistoPtr<string>> _h;
145    map<string, YODA::Axis<double>> _axes;
146    map<string, vector<string>> _edges;
147    /// @}
148
149
150  };
151
152
153  RIVET_DECLARE_PLUGIN(HRS_1990_I280958);
154
155
156}