1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/Sphericity.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6#include "Rivet/Projections/FinalState.hh"
  7#include "Rivet/Projections/Hemispheres.hh"
  8
  9namespace Rivet {
 10
 11  /// @brief Event shapes at 55.2
 12  class AMY_1990_I283337 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(AMY_1990_I283337);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      declare(Beam(), "Beams");
 25      const FinalState fs;
 26      declare(fs, "FS");
 27      const Sphericity sphere(fs);
 28      declare(sphere, "Sphericity");
 29      const Thrust thrust(fs);
 30      declare(thrust, "Thrust");
 31      declare(Hemispheres(sphere), "Hemispheres");
 32      // histograms
 33      book(_histRapidityT , 1, 1, 1);
 34      book(_histMajor     ,13, 1, 1);
 35
 36      book(_histScaledMom , 2, 1, 1);
 37      book(_histPl        , 3, 1, 1);
 38      book(_histPt        , 4, 1, 1);
 39      book(_histPt2       , 5, 1, 1);
 40      book(_histPtIn      , 6, 1, 1);
 41      book(_histPtOut     , 7, 1, 1);
 42      book(_histMeanPtIn2 , 8, 1, 1);
 43      book(_histMeanPtOut2, 9, 1, 1);
 44      book(_histNtheta    ,10, 1, 1);
 45      book(_histEtheta    ,11, 1, 1);
 46      book(_histThrust    ,12, 1, 1);
 47      book(_histMinor     ,14, 1, 1);
 48      book(_histOblateness,15, 1, 1);
 49      book(_histSphericity,16, 1, 1);
 50      book(_histAplanarity,17, 1, 1);
 51      book(_histQx        ,18, 1, 1);
 52      book(_histQ21       ,19, 1, 1);
 53      book(_histRhoLight  ,20, 1, 1);
 54      book(_histRhoHeavy  ,21, 1, 1);
 55      book(_histRhoDiff   ,22, 1, 1);
 56      book(_wSum,"TMP/wSum");
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      // First, veto on leptonic events by requiring at least 4 charged FS particles
 63      const FinalState& fs = apply<FinalState>(event, "FS");
 64      const size_t numParticles = fs.particles().size();
 65      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 66      if (numParticles < 2) {
 67        MSG_DEBUG("Failed leptonic event cut");
 68        vetoEvent;
 69      }
 70      MSG_DEBUG("Passed leptonic event cut");
 71      _wSum->fill();
 72
 73      // Get beams and average beam momentum
 74      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 75      const double meanBeamMom = ( beams.first.p3().mod() +
 76                                   beams.second.p3().mod() ) / 2.0;
 77      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 78
 79      // Thrusts
 80      MSG_DEBUG("Calculating thrust");
 81      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 82      _histThrust    ->fill(thrust.thrust()     );
 83      _histMajor     ->fill(thrust.thrustMajor());
 84      _histMinor     ->fill(thrust.thrustMinor());
 85      _histOblateness->fill(thrust.oblateness() );
 86      // Sphericities
 87      MSG_DEBUG("Calculating sphericity");
 88      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 89      _histSphericity->fill(sphericity.sphericity());
 90      _histAplanarity->fill(sphericity.aplanarity());
 91      _histQx        ->fill(sqrt(1./3.)*(sphericity.lambda1()-sphericity.lambda2()));
 92      _histQ21       ->fill(sphericity.lambda2()-sphericity.lambda3());
 93      // Hemispheres
 94      MSG_DEBUG("Calculating hemisphere variables");
 95      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 96      _histRhoHeavy->fill(hemi.scaledM2high());
 97      _histRhoLight->fill(hemi.scaledM2low() );
 98      _histRhoDiff ->fill(hemi.scaledM2diff());
 99      // single particle distributions
100      double pTIn2(0.),pTOut2(0.);
101      unsigned int nCharged(0);
102      for (const Particle& p : fs.particles()) {
103        // Get momentum and energy of each particle.
104        const Vector3 mom3 = p.p3();
105        const double energy = p.E();
106        const double mom = mom3.mod();
107        const double scaledMom = mom/meanBeamMom;
108        const double momT = dot(thrust.thrustAxis(), mom3);
109        const double momS = dot(sphericity.sphericityAxis(), mom3);
110        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
111        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
112        const double pT = sqrt(pow(pTinS, 2) + pow(pToutS, 2));
113
114        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
115        double angle = sphericity.sphericityAxis().angle(p.p3())/M_PI*180.;
116        if(angle>90.) angle=180.-angle;
117        if(PID::isCharged(p.pid())) {
118          _histScaledMom->fill(scaledMom);
119          _histRapidityT->fill(fabs(rapidityT));
120          _histPl       ->fill(fabs(momS)     );
121          _histPt       ->fill(pT             );
122          _histPt2      ->fill(sqr(pT)        );
123          _histPtIn     ->fill(fabs(pTinS)    );
124          _histPtOut    ->fill(fabs(pToutS)   );
125          pTIn2  += sqr(pTinS);
126          pTOut2 += sqr(pToutS);
127          _histNtheta->fill(angle);
128          ++nCharged;
129        }
130        _histEtheta->fill(angle,energy); 
131      }
132      if (nCharged) {
133        _histMeanPtIn2 ->fill( pTIn2/nCharged);
134        _histMeanPtOut2->fill(pTOut2/nCharged);
135      }
136    }
137
138
139    /// Normalise histograms etc., after the run
140    void finalize() {
141      // histograms
142      scale(_histRapidityT , 1./ *_wSum);
143      scale(_histScaledMom , 1./ *_wSum);
144      scale(_histPl        , 1./ *_wSum);
145      scale(_histPt        , 1./ *_wSum);
146      scale(_histPt2       , 1./ *_wSum);
147      scale(_histPtIn      , 1./ *_wSum);
148      scale(_histPtOut     , 1./ *_wSum);
149      scale(_histMeanPtIn2 , 1./ *_wSum);
150      scale(_histMeanPtOut2, 1./ *_wSum);
151      scale(_histNtheta    , 1./ *_wSum);
152      scale(_histEtheta    , 1./ *_wSum);
153      scale(_histThrust    , 1./ *_wSum);
154      scale(_histMajor     , 1./ *_wSum);
155      scale(_histMinor     , 1./ *_wSum);
156      scale(_histOblateness, 1./ *_wSum);
157      scale(_histSphericity, 1./ *_wSum);
158      scale(_histAplanarity, 1./ *_wSum);
159      scale(_histQx        , 1./ *_wSum);
160      scale(_histQ21       , 1./ *_wSum);
161      scale(_histRhoLight  , 1./ *_wSum);
162      scale(_histRhoHeavy  , 1./ *_wSum);
163      scale(_histRhoDiff   , 1./ *_wSum);
164    }
165
166    /// @}
167
168
169    /// @name Histograms
170    /// @{
171    Histo1DPtr _histRapidityT, _histScaledMom, _histPl, _histPt, _histPt2, _histPtIn, _histPtOut,
172      _histMeanPtIn2, _histMeanPtOut2, _histNtheta, _histEtheta, _histThrust, _histMajor, _histMinor,
173      _histOblateness, _histSphericity, _histAplanarity, _histQx, _histQ21, _histRhoLight,
174      _histRhoHeavy, _histRhoDiff;
175    CounterPtr _wSum;
176    /// @}
177
178
179  };
180
181
182  RIVET_DECLARE_PLUGIN(AMY_1990_I283337);
183
184
185}