1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief pi, K and p spectra at 34 and 44 GeV
 11  class TASSO_1989_I267755 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1989_I267755);
 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      declare(Beam(), "Beams");
 26      declare(ChargedFinalState(), "FS");
 27      declare(UnstableParticles(), "UFS");
 28
 29      // Book histograms
 30      // Book histograms
 31      _iHist=-1;
 32      if (isCompatibleWithSqrtS(34*GeV)) {
 33        _iHist = 0;
 34      }
 35      else if (isCompatibleWithSqrtS(44*GeV)) {
 36        _iHist = 1;
 37      }
 38      else
 39        MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
 40          << " doesn't match any available analysis energy .");
 41
 42      book(_h_x_pi, 3*_iHist+7,1,1);
 43      book(_h_x_K , 3*_iHist+8,1,1);
 44      book(_h_x_p , 3*_iHist+9,1,1);
 45      if(_iHist==1) book(_h_x_pi0,13,1,1);
 46      book(_n_pi,"TMP/n_pi",refData(3*_iHist+1,1,1));
 47      book(_d_pi,"TMP/d_pi",refData(3*_iHist+1,1,1));
 48      book(_n_K ,"TMP/n_K" ,refData(3*_iHist+2,1,1));
 49      book(_d_K ,"TMP/d_K" ,refData(3*_iHist+2,1,1));
 50      book(_n_p ,"TMP/n_p" ,refData(3*_iHist+3,1,1));
 51      book(_d_p ,"TMP/d_p" ,refData(3*_iHist+3,1,1));
 52    }
 53
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      // First, veto on leptonic events by requiring at least 4 charged FS particles
 58      const ChargedFinalState& fs = apply<ChargedFinalState>(event, "FS");
 59      const size_t numParticles = fs.particles().size();
 60
 61      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 62      if (numParticles < 2) {
 63        MSG_DEBUG("Failed leptonic event cut");
 64        vetoEvent;
 65      }
 66      MSG_DEBUG("Passed leptonic event cut");
 67
 68      // Get beams and average beam momentum
 69      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 70      const double meanBeamMom = ( beams.first.p3().mod() +
 71                                   beams.second.p3().mod() ) / 2.0;
 72      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 73
 74      for (const Particle& p : fs.particles()) {
 75      	double xP = p.p3().mod()/meanBeamMom;
 76      	_d_pi->fill(xP);
 77      	_d_K ->fill(xP);
 78      	_d_p ->fill(xP);
 79      	if(abs(p.pid())==211) {
 80      	  _h_x_pi->fill(xP);
 81      	  _n_pi  ->fill(xP);
 82      	}
 83      	else if(abs(p.pid())==321) {
 84      	  _h_x_K->fill(xP);
 85      	  _n_K  ->fill(xP);
 86      	}
 87      	else if(abs(p.pid())==2212) {
 88      	  _h_x_p->fill(xP);
 89      	  _n_p  ->fill(xP);
 90      	}
 91      }
 92      if(_h_x_pi0) {
 93      	for(const Particle & p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==111)) {
 94      	  double xP = p.p3().mod()/meanBeamMom;
 95      	  _h_x_pi0->fill(xP);
 96      	}
 97      }
 98    }
 99
100
101    /// Normalise histograms etc., after the run
102    void finalize() {
103      scale(_h_x_pi , 1./sumOfWeights());
104      scale(_h_x_K  , 1./sumOfWeights());
105      scale(_h_x_p  , 1./sumOfWeights());
106      if(_h_x_pi0) scale(_h_x_pi0, 1./sumOfWeights());
107      Estimate1DPtr temp1,temp2,temp3;
108      book(temp1,3*_iHist+1,1,1);
109      book(temp2,3*_iHist+2,1,1);
110      book(temp3,3*_iHist+3,1,1);
111
112      divide(_n_pi,_d_pi, temp1);
113      divide(_n_K ,_d_K , temp2);
114      divide(_n_p ,_d_p , temp3);
115    }
116
117    /// @}
118
119
120    /// @name Histograms
121    /// @{
122    Histo1DPtr _h_p_pi, _h_x_pi, _h_p_K, _h_x_K, _h_p_p, _h_x_p, _h_x_pi0;
123    Histo1DPtr _n_pi,_d_pi,_n_K,_d_K,_n_p,_d_p;
124    int _iHist;
125    /// @}
126
127
128  };
129
130
131  RIVET_DECLARE_PLUGIN(TASSO_1989_I267755);
132
133
134}