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)
 75	    ++n_pos_05;
 76	  else
 77	    ++n_neg_05;
 78	}
 79	if(abs(rapidityT)<1.0) {
 80	  ++n_all_10;
 81	  if(rapidityT>0)
 82	    ++n_pos_10;
 83	  else
 84	    ++n_neg_10;
 85	}
 86	if(abs(rapidityT)<1.5) {
 87	  ++n_all_15;
 88	  if(rapidityT>0)
 89	    ++n_pos_15;
 90	  else
 91	    ++n_neg_15;
 92	}
 93	if(abs(rapidityT)<2.0) {
 94	  ++n_all_20;
 95	  if(rapidityT>0)
 96	    ++n_pos_20;
 97	  else
 98	    ++n_neg_20;
 99	}
100	if(abs(rapidityT)<3.0) {
101	  if(rapidityT>0)
102	    ++n_pos_30;
103	  else
104	    ++n_neg_30;
105	}
106	if(abs(rapidityT)<4.0) {
107	  if(rapidityT>0)
108	    ++n_pos_40;
109	  else
110	    ++n_neg_40;
111	}
112	if(abs(rapidityT)<5.0) {
113	  if(rapidityT>0)
114	    ++n_pos_50;
115	  else
116	    ++n_neg_50;
117	}
118	if(rapidityT>0)
119	  ++n_pos_all;
120	else
121	  ++n_neg_all;
122      }
123      _h_all_05 ->fill(n_all_05 );
124      _h_all_10 ->fill(n_all_10 );
125      _h_all_15 ->fill(n_all_15 );
126      _h_all_20 ->fill(n_all_20 );
127      _h_all_all->fill(n_all_all);
128      _h_hemi_05 ->fill(n_pos_05 );
129      _h_hemi_10 ->fill(n_pos_10 );
130      _h_hemi_15 ->fill(n_pos_15 );
131      _h_hemi_20 ->fill(n_pos_20 );
132      _h_hemi_30 ->fill(n_pos_30 );
133      _h_hemi_40 ->fill(n_pos_40 );
134      _h_hemi_50 ->fill(n_pos_50 );
135      _h_hemi_all->fill(n_pos_all);
136      _h_hemi_05 ->fill(n_neg_05 );
137      _h_hemi_10 ->fill(n_neg_10 );
138      _h_hemi_15 ->fill(n_neg_15 );
139      _h_hemi_20 ->fill(n_neg_20 );
140      _h_hemi_30 ->fill(n_neg_30 );
141      _h_hemi_40 ->fill(n_neg_40 );
142      _h_hemi_50 ->fill(n_neg_50 );
143      _h_hemi_all->fill(n_neg_all);
144    }
145
146
147    /// Normalise histograms etc., after the run
148    void finalize() {
149      normalize( _h_all_05  , 1000.);
150      normalize( _h_all_10  , 1000.);
151      normalize( _h_all_15  , 1000.);
152      normalize( _h_all_20  , 1000.);
153      normalize( _h_all_all , 2000.);
154      normalize( _h_hemi_05 , 1000.);
155      normalize( _h_hemi_10 , 1000.);
156      normalize( _h_hemi_15 , 1000.);
157      normalize( _h_hemi_20 , 1000.);
158      normalize( _h_hemi_30 , 1000.);
159      normalize( _h_hemi_40 , 1000.);
160      normalize( _h_hemi_50 , 1000.);
161      normalize( _h_hemi_all, 1000.);
162    }
163
164    //@}
165
166
167    /// @name Histograms
168    //@{
169    Histo1DPtr _h_all_05, _h_all_10, _h_all_15, _h_all_20, _h_all_all;
170    Histo1DPtr _h_hemi_05, _h_hemi_10, _h_hemi_15, _h_hemi_20,
171      _h_hemi_30, _h_hemi_40, _h_hemi_50, _h_hemi_all;
172    //@}
173
174
175  };
176
177
178  // The hook for the plugin system
179  RIVET_DECLARE_PLUGIN(DELPHI_1991_I324035);
180
181
182}