1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Charged particle spectra for a range of energies
 10  class TASSO_1982_I177174 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1982_I177174);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(Beam(), "Beams");
 24      declare(ChargedFinalState(), "FS");
 25      if(isCompatibleWithSqrtS(12*GeV)) {
 26	book(_h_x2,2,1,1);
 27	book(_h_x3,3,1,1);
 28      }
 29      else if (isCompatibleWithSqrtS(14*GeV)) {
 30	book(_h_x1,1,1,1);
 31	book(_h_x2,2,1,2);
 32	book(_h_x3,3,1,2);
 33      }
 34      else if (isCompatibleWithSqrtS(22*GeV)) {
 35	book(_h_x1,1,1,2);
 36	book(_h_x2,2,1,3);
 37	book(_h_x3,3,1,3);
 38      }
 39      else if (isCompatibleWithSqrtS(25*GeV)) {
 40	book(_h_x2,2,1,4);
 41	book(_h_x3,3,1,4);
 42      }
 43      else if (isCompatibleWithSqrtS(30*GeV)) {
 44	book(_h_x2,2,1,5);
 45	book(_h_x3,3,1,5);
 46      }
 47      else if (isCompatibleWithSqrtS(34*GeV)) {
 48	book(_h_x2,2,1,6);
 49	book(_h_x3,3,1,6);
 50      }
 51      else if (isCompatibleWithSqrtS(35*GeV)) {
 52        book(_h_x2, 2,1,7);
 53        book(_h_x3, 3,1,7);
 54      }
 55
 56      if(inRange(sqrtS()/GeV, 29.9,36.7))
 57        book(_h_x1, 1,1,3);
 58
 59      if(_h_x1==Histo1DPtr() && _h_x2==Histo1DPtr() && _h_x3==Histo1DPtr())
 60      	MSG_ERROR("Beam energy not supported!");
 61
 62    }
 63
 64
 65    /// Perform the per-event analysis
 66    void analyze(const Event& event) {
 67      // First, veto on leptonic events by requiring at least 4 charged FS particles
 68      const ChargedFinalState& fs = apply<ChargedFinalState>(event, "FS");
 69      const size_t numParticles = fs.particles().size();
 70
 71      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 72      if (numParticles < 2) {
 73        MSG_DEBUG("Failed leptonic event cut");
 74        vetoEvent;
 75      }
 76      MSG_DEBUG("Passed leptonic event cut");
 77
 78      // Get beams and average beam momentum
 79      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 80      const double meanBeamMom = ( beams.first.p3().mod() +
 81                                   beams.second.p3().mod() ) / 2.0;
 82      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 83
 84      for (const Particle& p : fs.particles()) {
 85	double xp = p.p3().mod()/meanBeamMom;
 86	if(_h_x1!=Histo1DPtr()) _h_x1->fill(xp);
 87	if(_h_x2!=Histo1DPtr()) _h_x2->fill(xp);
 88	if(_h_x3!=Histo1DPtr()) _h_x3->fill(xp);
 89      }
 90    }
 91
 92
 93    /// Normalise histograms etc., after the run
 94    void finalize() {
 95
 96      scale(_h_x1, crossSection()/microbarn/sumOfWeights()*sqr(sqrtS()));
 97      scale(_h_x2, crossSection()/microbarn/sumOfWeights()*sqr(sqrtS()));
 98      scale(_h_x3, 1./sumOfWeights());
 99
100    }
101
102    /// @}
103
104
105    /// @name Histograms
106    /// @{
107    Histo1DPtr _h_x1,_h_x2,_h_x3;
108    /// @}
109
110
111  };
112
113
114  RIVET_DECLARE_PLUGIN(TASSO_1982_I177174);
115
116
117}