1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Thrust.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Charged particle multiplicities in different regions
 10  class DELPHI_1991_I324035 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1991_I324035);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      const ChargedFinalState cfs;
 25      declare(cfs, "FS");
 26      const Thrust thrust(cfs);
 27      declare(thrust, "Thrust");
 28
 29      // Book histograms
 30      book(_h_all_05  ,  1, 1, 1);
 31      book(_h_all_10  ,  2, 1, 1);
 32      book(_h_all_15  ,  3, 1, 1);
 33      book(_h_all_20  ,  4, 1, 1);
 34      book(_h_all_all ,  5, 1, 1);
 35      book(_h_hemi_05 ,  6, 1, 1);
 36      book(_h_hemi_10 ,  7, 1, 1);
 37      book(_h_hemi_15 ,  8, 1, 1);
 38      book(_h_hemi_20 ,  9, 1, 1);
 39      book(_h_hemi_30 , 10, 1, 1);
 40      book(_h_hemi_40 , 11, 1, 1);
 41      book(_h_hemi_50 , 12, 1, 1);
 42      book(_h_hemi_all, 13, 1, 1);
 43    }
 44
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      // First, veto on leptonic events by requiring at least 4 charged FS particles
 49      const FinalState& fs = apply<FinalState>(event, "FS");
 50      const size_t numParticles = fs.particles().size();
 51      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 52      if (numParticles < 2) {
 53        MSG_DEBUG("Failed leptonic event cut");
 54        vetoEvent;
 55      }
 56      MSG_DEBUG("Passed leptonic event cut");
 57
 58      // Thrusts
 59      MSG_DEBUG("Calculating thrust");
 60      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 61      Vector3 axis=thrust.thrustAxis();
 62
 63      unsigned int n_all_05(0),n_all_10(0),n_all_15(0),n_all_20(0),n_all_all(0);
 64      unsigned int n_pos_05(0),n_pos_10(0),n_pos_15(0),n_pos_20(0),n_pos_30(0),n_pos_40(0),n_pos_50(0),n_pos_all(0);
 65      unsigned int n_neg_05(0),n_neg_10(0),n_neg_15(0),n_neg_20(0),n_neg_30(0),n_neg_40(0),n_neg_50(0),n_neg_all(0);
 66      for (const Particle &p : fs.particles()) {
 67        const Vector3 mom3 = p.p3();
 68        const double energy = p.E();
 69        const double momT = dot(axis, mom3);
 70        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
 71        ++n_all_all;
 72        if (abs(rapidityT)<0.5) {
 73          ++n_all_05;
 74          if(rapidityT>0)  ++n_pos_05;
 75          else             ++n_neg_05;
 76        }
 77        if (abs(rapidityT)<1.0) {
 78          ++n_all_10;
 79          if(rapidityT>0)  ++n_pos_10;
 80          else             ++n_neg_10;
 81        }
 82        if (abs(rapidityT)<1.5) {
 83          ++n_all_15;
 84          if(rapidityT>0)  ++n_pos_15;
 85          else             ++n_neg_15;
 86        }
 87        if (abs(rapidityT)<2.0) {
 88          ++n_all_20;
 89          if(rapidityT>0)  ++n_pos_20;
 90          else             ++n_neg_20;
 91        }
 92        if (abs(rapidityT)<3.0) {
 93          if(rapidityT>0)  ++n_pos_30;
 94          else             ++n_neg_30;
 95        }
 96        if (abs(rapidityT)<4.0) {
 97          if(rapidityT>0)  ++n_pos_40;
 98          else             ++n_neg_40;
 99        }
100        if(abs(rapidityT)<5.0) {
101          if(rapidityT>0)  ++n_pos_50;
102          else             ++n_neg_50;
103        }
104        if(rapidityT>0)  ++n_pos_all;
105        else             ++n_neg_all;
106      }
107      _h_all_05 ->fill(n_all_05 );
108      _h_all_10 ->fill(n_all_10 );
109      _h_all_15 ->fill(n_all_15 );
110      _h_all_20 ->fill(n_all_20 );
111      _h_all_all->fill(n_all_all);
112      _h_hemi_05 ->fill(n_pos_05 );
113      _h_hemi_10 ->fill(n_pos_10 );
114      _h_hemi_15 ->fill(n_pos_15 );
115      _h_hemi_20 ->fill(n_pos_20 );
116      _h_hemi_30 ->fill(n_pos_30 );
117      _h_hemi_40 ->fill(n_pos_40 );
118      _h_hemi_50 ->fill(n_pos_50 );
119      _h_hemi_all->fill(n_pos_all);
120      _h_hemi_05 ->fill(n_neg_05 );
121      _h_hemi_10 ->fill(n_neg_10 );
122      _h_hemi_15 ->fill(n_neg_15 );
123      _h_hemi_20 ->fill(n_neg_20 );
124      _h_hemi_30 ->fill(n_neg_30 );
125      _h_hemi_40 ->fill(n_neg_40 );
126      _h_hemi_50 ->fill(n_neg_50 );
127      _h_hemi_all->fill(n_neg_all);
128    }
129
130
131    /// Normalise histograms etc., after the run
132    void finalize() {
133      normalize( _h_all_05  , 1000.);
134      normalize( _h_all_10  , 1000.);
135      normalize( _h_all_15  , 1000.);
136      normalize( _h_all_20  , 1000.);
137      normalize( _h_all_all , 1000.);
138      normalize( _h_hemi_05 , 1000.);
139      normalize( _h_hemi_10 , 1000.);
140      normalize( _h_hemi_15 , 1000.);
141      normalize( _h_hemi_20 , 1000.);
142      normalize( _h_hemi_30 , 1000.);
143      normalize( _h_hemi_40 , 1000.);
144      normalize( _h_hemi_50 , 1000.);
145      normalize( _h_hemi_all, 1000.);
146    }
147
148    /// @}
149
150
151    /// @name Histograms
152    /// @{
153    BinnedHistoPtr<int> _h_all_05, _h_all_10, _h_all_15, _h_all_20, _h_all_all;
154    BinnedHistoPtr<int> _h_hemi_05, _h_hemi_10, _h_hemi_15, _h_hemi_20;
155    BinnedHistoPtr<int> _h_hemi_30, _h_hemi_40, _h_hemi_50, _h_hemi_all;
156    /// @}
157
158
159  };
160
161
162  RIVET_DECLARE_PLUGIN(DELPHI_1991_I324035);
163
164
165}