1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/Thrust.hh"
  5#include "Rivet/Projections/Sphericity.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief event shapes at 35 GeV
 12  class TASSO_1988_I263859 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1988_I263859);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      const ChargedFinalState cfs;
 25      declare(cfs, "CFS");
 26
 27      // Thrust and sphericity
 28      declare(Beam(), "Beams");
 29      declare(Thrust(cfs), "Thrust");
 30      declare(Sphericity(cfs), "Sphericity");
 31
 32      // Book histograms
 33      book(_h_sphericity,  1, 1, 1);
 34      book(_h_aplanarity,  2, 1, 1);
 35      book(_h_thrust    ,  3, 1, 1);
 36      book(_h_pTin2     ,  4, 1, 1);
 37      book(_h_pTout2    ,  5, 1, 1);
 38      book(_h_ncharged  ,  6, 1, 1);
 39      book(_h_pTin      ,  7, 1, 1);
 40      book(_h_pTout     ,  8, 1, 1);
 41      book(_h_pT        ,  9, 1, 1);
 42      book(_h_x         , 10, 1, 1);
 43      book(_h_rap       , 11, 1, 1);
 44
 45    }
 46
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 51
 52      // Get beams and average beam momentum
 53      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 54      const double meanBeamMom = ( beams.first.p3().mod() +
 55				   beams.second.p3().mod() ) / 2.0;
 56      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 57      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 58
 59      _h_sphericity->fill(sphericity.sphericity());
 60      _h_aplanarity->fill(sphericity.aplanarity());
 61      _h_thrust    ->fill(thrust.thrust());
 62
 63      double pTin2sum(0.), pTout2sum(0.);
 64      for (const Particle& p : cfs.particles()) {
 65        // Get momentum and energy of each particle.
 66        const Vector3 mom3 = p.p3();
 67        const double energy = p.E();
 68        // Scaled momenta.
 69        const double mom = mom3.mod();
 70        const double scaledMom = mom/meanBeamMom;
 71        _h_x->fill(scaledMom);
 72
 73        const double momS = dot(sphericity.sphericityAxis(), mom3);
 74        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
 75        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
 76        const double pT = sqrt(pow(pTinS, 2) + pow(pToutS, 2));
 77        const double rapidityS = 0.5 * std::log((energy + momS) / (energy - momS));
 78
 79        pTin2sum  += sqr(pTinS);
 80        pTout2sum += sqr(pToutS);
 81
 82        _h_pTin ->fill(abs(pTinS)/GeV );
 83        _h_pTout->fill(abs(pToutS)/GeV);
 84        _h_pT   ->fill(pT/GeV         );
 85        _h_rap  ->fill(abs(rapidityS) );
 86
 87      }
 88      unsigned int nCharged = cfs.particles().size();
 89      _h_ncharged->fill(nCharged);
 90      if (nCharged) {
 91        _h_pTin2 ->fill(pTin2sum /nCharged);
 92        _h_pTout2->fill(pTout2sum/nCharged);
 93      }
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      scale(_h_sphericity , 1./sumOfWeights());
100      scale(_h_aplanarity , 1./sumOfWeights());
101      scale(_h_thrust     , 1./sumOfWeights());
102      scale(_h_pTin2      , 1./sumOfWeights());
103      scale(_h_pTout2     , 1./sumOfWeights());
104      scale(_h_ncharged   , 2000./sumOfWeights());
105      scale(_h_pTin       , 1./sumOfWeights());
106      scale(_h_pTout      , 1./sumOfWeights());
107      scale(_h_pT         , 1./sumOfWeights());
108      scale(_h_x          , 1./sumOfWeights());
109      scale(_h_rap        , 1./sumOfWeights());
110    }
111
112    /// @}
113
114
115    /// @name Histograms
116    /// @{
117    Histo1DPtr _h_sphericity, _h_aplanarity, _h_thrust,
118      _h_pTin2, _h_pTout2, _h_ncharged,
119      _h_pTin, _h_pTout, _h_pT, _h_x, _h_rap;
120    /// @}
121
122  };
123
124
125  RIVET_DECLARE_PLUGIN(TASSO_1988_I263859);
126
127}