1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/Thrust.hh"
  7#include "Rivet/Projections/Sphericity.hh"
  8#include "Rivet/Projections/Hemispheres.hh"
  9#include "Rivet/Projections/ParisiTensor.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief event shapes at 133, 161 172 and 183 GeV
 15  class DELPHI_1999_I499183 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1999_I499183);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      // Initialise and register projections
 29      declare(Beam(), "Beams");
 30      const FinalState fs;
 31      declare(fs, "FS");
 32      const Thrust thrust(fs);
 33      declare(thrust, "Thrust");
 34      declare(Sphericity(fs), "Sphericity");
 35      declare(ParisiTensor(fs), "Parisi");
 36      declare(Hemispheres(thrust), "Hemispheres");
 37
 38      // Book histograms
 39      unsigned int offset = 0;
 40      int offset2 = 0;
 41
 42      if (isCompatibleWithSqrtS(133*GeV)) {
 43	offset  = 0;
 44	offset2 = 1;
 45      }
 46      else if (isCompatibleWithSqrtS(161*GeV)) {
 47	offset  = 0;
 48	offset2 = 2;
 49      }
 50      else if (isCompatibleWithSqrtS(172*GeV)) {
 51	offset  = 0;
 52	offset2 = 3;
 53      }
 54      else if (isCompatibleWithSqrtS(183*GeV)) {
 55	offset  = 1;
 56	offset2 = 1;
 57      }
 58
 59      book(_h_thrust          , 13+offset, 1, offset2);
 60      book(_h_major           , 15+offset, 1, offset2);
 61      book(_h_minor           , 17+offset, 1, offset2);
 62      book(_h_oblateness      , 19+offset, 1, offset2);
 63      book(_h_sphericity      , 21+offset, 1, offset2);
 64      book(_h_planarity       , 23+offset, 1, offset2);
 65      book(_h_aplanarity      , 25+offset, 1, offset2);
 66      book(_h_heavy_jet_mass  , 27+offset, 1, offset2);
 67      book(_h_light_jet_mass  , 29+offset, 1, offset2);
 68      book(_h_diff_jet_mass   , 31+offset, 1, offset2);
 69      book(_h_wide_broading   , 33+offset, 1, offset2);
 70      book(_h_narrow_broading , 35+offset, 1, offset2);
 71      book(_h_total_broading  , 37+offset, 1, offset2);
 72      book(_h_diff_broading   , 39+offset, 1, offset2);
 73      book(_h_CParam          , 41+offset, 1, offset2);
 74      book(_h_DParam          , 43+offset, 1, offset2);
 75    }
 76
 77
 78    /// Perform the per-event analysis
 79    void analyze(const Event& event) {
 80
 81      // Get beams and average beam momentum
 82      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 83      const double meanBeamMom = ( beams.first.p3().mod() +
 84                                   beams.second.p3().mod() ) / 2.0;
 85      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 86
 87      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 88      // thrust related observables
 89      _h_thrust    ->fill(1.-thrust.thrust()  );
 90      _h_major     ->fill(thrust.thrustMajor());
 91      _h_minor     ->fill(thrust.thrustMinor());
 92      _h_oblateness->fill(thrust.oblateness() );
 93
 94      // sphericity related
 95      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 96      _h_sphericity->fill(sphericity.sphericity());
 97      _h_planarity ->fill(sphericity.planarity() );
 98      _h_aplanarity->fill(sphericity.aplanarity());
 99      // hemisphere related
100      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
101      // standard jet masses
102      _h_heavy_jet_mass->fill(hemi.scaledM2high());
103      _h_light_jet_mass->fill(hemi.scaledM2low() );
104      _h_diff_jet_mass ->fill(hemi.scaledM2diff());
105      // jet broadening
106      _h_wide_broading  ->fill(hemi.Bmax() );
107      _h_narrow_broading->fill(hemi.Bmin() );
108      _h_total_broading ->fill(hemi.Bsum() );
109      _h_diff_broading  ->fill(hemi.Bdiff());
110      MSG_DEBUG("Calculating Parisi params");
111      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
112      _h_CParam->fill(parisi.C());
113      _h_DParam->fill(parisi.D());
114    }
115
116
117    /// Normalise histograms etc., after the run
118    void finalize() {
119
120      normalize(_h_thrust          );
121      normalize(_h_major           );
122      normalize(_h_minor           );
123      normalize(_h_sphericity      );
124      normalize(_h_planarity       );
125      normalize(_h_aplanarity       );
126      normalize(_h_oblateness      );
127      normalize(_h_heavy_jet_mass  );
128      normalize(_h_light_jet_mass  );
129      normalize(_h_diff_jet_mass   );
130      normalize(_h_wide_broading   );
131      normalize(_h_narrow_broading );
132      normalize(_h_total_broading  );
133      normalize(_h_diff_broading   );
134      normalize(_h_CParam   );
135      normalize(_h_DParam   );
136
137    }
138
139    /// @}
140
141
142    /// @name Histograms
143    /// @{
144    Histo1DPtr _h_thrust,_h_major,_h_minor;
145    Histo1DPtr _h_sphericity,_h_planarity,_h_aplanarity,_h_oblateness;
146    Histo1DPtr _h_heavy_jet_mass,_h_light_jet_mass,_h_diff_jet_mass;
147    Histo1DPtr _h_wide_broading,_h_narrow_broading,_h_total_broading,_h_diff_broading;
148    Histo1DPtr _h_CParam,_h_DParam;
149    /// @}
150
151
152  };
153
154
155  RIVET_DECLARE_PLUGIN(DELPHI_1999_I499183);
156
157
158}