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/FinalState.hh"
  7#include "Rivet/Projections/ChargedFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Average pT w.r.t. thrust and sphericity axes
 13  class PLUTO_1983_I191161 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1983_I191161);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27      const FinalState fs;
 28      declare(fs, "FS");
 29      declare(ChargedFinalState(), "CFS");
 30      // Thrust
 31      const Thrust thrust(fs);
 32      declare(thrust, "Thrust");
 33      const Sphericity sphericity(fs);
 34      declare(sphericity, "Sphericity");
 35
 36      sqs = 1.0;
 37      if (isCompatibleWithSqrtS(7.7))       sqs = "7.7";
 38      else if (isCompatibleWithSqrtS(9.4))  sqs = "9.4";
 39      else if (isCompatibleWithSqrtS(12.))  sqs = "12.0";
 40      else if (isCompatibleWithSqrtS(13.))  sqs = "13.0";
 41      else if (isCompatibleWithSqrtS(17.))  sqs = "17.0";
 42      else if (isCompatibleWithSqrtS(22.))  sqs = "22.0";
 43      else if (isCompatibleWithSqrtS(27.6)) sqs = "27.6";
 44      else if (isCompatibleWithSqrtS(30.8)) sqs = "30.0 - 31.6";
 45      else MSG_ERROR("Beam energy " << sqrtS() << " not supported!");
 46      book(_p_thrust_pt      , 1,1,1);
 47      book(_p_thrust_pt2     , 1,1,2);
 48      book(_p_thrust_sum_pt  , 1,1,3);
 49      book( _p_thrust_sum_pt2, 1,1,4);
 50      book(_p_sphere_pt      , 2,1,1);
 51      book( _p_sphere_pt2    , 2,1,2);
 52      book( _p_sphere_sum_pt , 2,1,3);
 53      book( _p_sphere_sum_pt2, 2,1,4);
 54    }
 55
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59      // at least 4 charged particles
 60      if (apply<ChargedFinalState>(event, "CFS").particles().size()<4) vetoEvent;
 61      // Sphericities
 62      MSG_DEBUG("Calculating sphericity");
 63      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 64      MSG_DEBUG("Calculating thrust");
 65      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 66      const FinalState & fs = apply<FinalState>(event, "FS");
 67      // remove pi0->gammagamma decay products and replace with the pi0s
 68      // needed to get the average pTs right
 69      Particles fsParticles;
 70      set<ConstGenParticlePtr> pi0;
 71      for(const Particle & p : fs.particles()) {
 72        if ((p.pid()!=PID::PHOTON && p.abspid()!=PID::ELECTRON)|| p.parents().empty() || p.parents()[0].pid()!=PID::PI0)
 73          fsParticles.push_back(p);
 74        else {
 75          if (pi0.find(p.parents()[0].genParticle())==pi0.end())
 76            fsParticles.push_back(p.parents()[0]);
 77        }
 78      }
 79      double pT_T_sum(0.),pT2_T_sum(0.);
 80      double pT_S_sum(0.),pT2_S_sum(0.);
 81      for(const Particle & p : fsParticles) {
 82        const Vector3 mom3 = p.p3();
 83        const double pTinT = dot(mom3, thrust.thrustMajorAxis());
 84        const double pToutT = dot(mom3, thrust.thrustMinorAxis());
 85        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
 86        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
 87        const double pT2_T = sqr(pTinT) + sqr(pToutT);
 88        const double pT2_S = sqr(pTinS) + sqr(pToutS);
 89        const double pT_T  = sqrt(pT2_T);
 90        const double pT_S  = sqrt(pT2_S);
 91        pT_T_sum  += sqrt(pT2_T);
 92        pT2_T_sum +=      pT2_T ;
 93        pT_S_sum  += sqrt(pT2_S);
 94        pT2_S_sum +=      pT2_S ;
 95        _p_thrust_pt ->fill(sqs,pT_T /MeV         );
 96        _p_thrust_pt2->fill(sqs,pT2_T/1e3/sqr(MeV));
 97        _p_sphere_pt ->fill(sqs,pT_S /MeV         );
 98        _p_sphere_pt2->fill(sqs,pT2_S/1e3/sqr(MeV));
 99      }
100      _p_thrust_sum_pt ->fill(sqs,pT_T_sum /GeV);
101      _p_thrust_sum_pt2->fill(sqs,pT2_T_sum/GeV);
102      _p_sphere_sum_pt ->fill(sqs,pT_S_sum /GeV);
103      _p_sphere_sum_pt2->fill(sqs,pT2_S_sum/GeV);
104    }
105
106
107    /// Normalise histograms etc., after the run
108    void finalize() {
109    }
110
111    /// @}
112
113
114    /// @name Histograms
115    /// @{
116    BinnedProfilePtr<string> _p_thrust_pt, _p_thrust_pt2, _p_thrust_sum_pt, _p_thrust_sum_pt2;
117    BinnedProfilePtr<string> _p_sphere_pt, _p_sphere_pt2, _p_sphere_sum_pt, _p_sphere_sum_pt2;
118    string sqs;
119    /// @}
120
121  };
122
123
124  RIVET_DECLARE_PLUGIN(PLUTO_1983_I191161);
125
126}