1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/Beam.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief K0 and Lambda spectra at 14, 22 and 34 GeV.
 10  class TASSO_1985_I205119 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1985_I205119);
 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      declare(Beam(), "Beams");
 25      declare(UnstableParticles(), "UFS");
 26
 27      // Book histograms
 28      sqs = 1.;
 29      if(isCompatibleWithSqrtS(14*GeV)) {
 30	book(_h_kaon_x  ,  1,1,1);
 31	book(_h_lambda_x,  4,1,1);
 32	book(_h_kaon_p  ,  7,1,1);
 33	book(_h_lambda_p, 10,1,1);
 34	sqs = 14.;
 35      }
 36      else if (isCompatibleWithSqrtS(22*GeV)) {
 37	book(_h_kaon_x  ,  2,1,1);
 38	book(_h_lambda_x,  5,1,1);
 39	book(_h_kaon_p  ,  8,1,1);
 40	book(_h_lambda_p, 11,1,1);
 41	sqs = 22.;
 42      }
 43      else if (isCompatibleWithSqrtS(34*GeV)) {
 44	book(_h_kaon_x  , 3,1,1);
 45	book(_h_lambda_x, 6,1,1);
 46	book(_h_kaon_p  , 9,1,1);
 47	book(_h_lambda_p,12,1,1);
 48	sqs = 34.;
 49      }
 50      else
 51	MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
 52                    << " doesn't match any available analysis energy .");
 53    }
 54
 55
 56    /// Perform the per-event analysis
 57    void analyze(const Event& event) {
 58      // Get beams and average beam momentum
 59      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 60      const double meanBeamMom = ( beams.first.p3().mod() +
 61                                   beams.second.p3().mod() ) / 2.0;
 62      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 63
 64      for (const Particle& p : apply<UnstableParticles>(event, "UFS").
 65	       particles(Cuts::abspid==PID::LAMBDA or Cuts::pid==130 or Cuts::pid==310)) {
 66	double xE = p.E()/meanBeamMom;
 67	double modp = p.p3().mod();
 68	double beta = modp/p.E();
 69	if(abs(p.pid())==PID::LAMBDA) {
 70	  _h_lambda_x->fill(xE,1./beta);
 71	  _h_lambda_p->fill(modp,1.);
 72	}
 73	else {
 74	  _h_kaon_x->fill(xE,1./beta);
 75	  _h_kaon_p->fill(modp,1.);
 76	}
 77      }
 78    }
 79
 80
 81    /// Normalise histograms etc., after the run
 82    void finalize() {
 83
 84      scale(_h_kaon_x  , sqr(sqs)*crossSection()/microbarn/sumOfWeights());
 85      scale(_h_lambda_x, sqr(sqs)*crossSection()/microbarn/sumOfWeights());
 86      scale(_h_kaon_p  , crossSection()/nanobarn/sumOfWeights());
 87      scale(_h_lambda_p, crossSection()/nanobarn/sumOfWeights());
 88    }
 89
 90    /// @}
 91
 92
 93    /// @name Histograms
 94    /// @{
 95    Histo1DPtr _h_kaon_x, _h_lambda_x, _h_kaon_p, _h_lambda_p;
 96    double sqs;
 97    /// @}
 98
 99
100  };
101
102
103  RIVET_DECLARE_PLUGIN(TASSO_1985_I205119);
104
105
106}