1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Charged particle multiplicities and distributions
 11  class TASSO_1980_I143691 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1980_I143691);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      ChargedFinalState cfs;
 26      declare(cfs, "CFS");
 27      declare(Thrust(cfs), "Thrust");
 28
 29
 30      // Book histograms
 31      book(_mult, 1,1,1);
 32      unsigned int iloc(0);
 33      sqs = 1.0;
 34      if(isCompatibleWithSqrtS(13*GeV)) {
 35	iloc = 1;
 36	sqs = 13.0;
 37        sloc="13.0";
 38      }
 39      else if(isCompatibleWithSqrtS(17*GeV)) {
 40	iloc = 2;
 41        sqs = 17.0;
 42        sloc=  "17.0";
 43      }
 44      else if (isCompatibleWithSqrtS(22*GeV)) {
 45	iloc = 2;
 46        sloc= "22.0";
 47	sqs = 22.;
 48      }
 49      else if(isCompatibleWithSqrtS(27.6*GeV)) {
 50	iloc = 3;
 51	sqs = 27.6;
 52        sloc="27.6";
 53      }
 54      else if (isCompatibleWithSqrtS(30.3*GeV)) {
 55	iloc = 3;
 56	sqs = 30.3;
 57        sloc= "30.3";
 58      }
 59      else if (isCompatibleWithSqrtS(31.2*GeV)) {
 60	iloc = 3;
 61	sqs = 31.2;
 62        sloc= "31.2";
 63      }
 64      else
 65	MSG_ERROR("Beam energy not supported!");
 66
 67      book(_h_rap,iloc+1,1,1);
 68      book(_h_x  ,iloc+4,1,1);
 69    }
 70
 71
 72    /// Perform the per-event analysis
 73    void analyze(const Event& event) {
 74      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 75      // thrust
 76      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 77      Vector3 axis=thrust.thrustAxis();
 78      _mult->fill(sloc,cfs.particles().size());
 79      for (const Particle& p : cfs.particles()) {
 80	const Vector3 mom3 = p.p3();
 81	double pp = mom3.mod();
 82	double xP = 2.*pp/sqs;
 83	_h_x->fill(xP);
 84        const double mom = dot(axis, mom3);
 85	const double rap = 0.5 * log((p.E() + mom) /
 86				     (p.E() - mom));
 87	_h_rap->fill(fabs(rap));
 88      }
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94
 95      scale(_h_rap, 1./sumOfWeights());
 96      scale(_h_x  , crossSection()*sqr(sqs)/sumOfWeights()/microbarn);
 97    }
 98
 99    /// @}
100
101
102    /// @name Histograms
103    /// @{
104    Histo1DPtr _h_rap, _h_x;
105    BinnedProfilePtr<string> _mult;
106    double sqs;
107    string sloc;
108    /// @}
109
110
111  };
112
113
114  RIVET_DECLARE_PLUGIN(TASSO_1980_I143691);
115
116
117}