1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/Logging.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Math/Units.hh"
  6#include <vector>
  7
  8using namespace std;
  9
 10namespace Rivet {
 11
 12
 13  class  LHCB_2015_I1333223 : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    /// @{
 18
 19    /// Constructor
 20     LHCB_2015_I1333223()
 21      : Analysis("LHCB_2015_I1333223")
 22    {    }
 23
 24    /// @}
 25
 26
 27  public:
 28
 29    /// @name Analysis methods
 30    /// @{
 31
 32    /// Book histograms and initialise projections before the run
 33    void init() {
 34      // Charged particles
 35      declare(ChargedFinalState(Cuts::eta> 2.0 && Cuts::eta <4.5 && Cuts::pT >0.2*GeV), "CFS");
 36      // Reproducing only measurement for prompt charged particles
 37      book(_hInelasticXs ,1, 1, 1);
 38    }
 39
 40
 41    /// Perform the per-event analysis
 42    void analyze(const Event& event) {
 43      const ChargedFinalState   &cfs    = apply<ChargedFinalState> (event, "CFS");
 44
 45      // eliminate non-inelastic events and empty events in LHCb
 46      if (cfs.particles().size() == 0) vetoEvent;
 47
 48      // See if this event has at least one prompt particle
 49      for (const Particle &myp : cfs.particles()){
 50          double dPV = getPVDCA(myp);
 51          // if IP > 200 microns the particle is not considered prompt
 52          if ((dPV < 0.) || (dPV > 0.2 * millimeter)) {
 53            MSG_DEBUG(" Vetoing " << myp.pid() << " at " << dPV);
 54            continue;
 55          }
 56          // histo gets filled only for inelastic events (at least one prompt charged particle)
 57          _hInelasticXs->fill(sqrtS());
 58          break;
 59      } //end loop on particles
 60
 61    }
 62
 63
 64    /// Normalise histograms etc., after the run
 65    void finalize() {
 66      scale(_hInelasticXs, crossSection()/sumOfWeights()/millibarn);
 67    }
 68
 69    /// @}
 70
 71
 72  private:
 73
 74    /// Compute distance of closest approach in z range for one particle.
 75    /// Returns -1 if unable to compute the DCA to PV.
 76    double getPVDCA(const Particle& p) {
 77      ConstGenVertexPtr vtx = p.genParticle()->production_vertex();
 78      if ( 0 == vtx ) return -1.;
 79
 80      // Unit vector of particle's MOMENTUM three vector
 81      const Vector3 u = p.momentum().p3().unit();
 82
 83      // The interaction point is always at (0, 0,0,0) hence the
 84      // vector pointing from the PV to the particle production vertex is:
 85      Vector3 d(vtx->position().x(), vtx->position().y(), vtx->position().z());
 86
 87      // Subtract projection of d onto u from d
 88      double proj = d.dot(u);
 89      d -= (u * proj);
 90
 91      // d should be orthogonal to u and it's length give the distance of
 92      // closest approach
 93      return d.mod();
 94    }
 95
 96
 97    /// @name Histograms
 98    /// @{
 99    Histo1DPtr _hInelasticXs;
100    /// @}
101    //
102  };
103
104
105
106  RIVET_DECLARE_PLUGIN(LHCB_2015_I1333223);
107
108}