Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {


  /// @brief charm spectra at 10.55
  class CLEO_1988_I22954 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1988_I22954);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare(UnstableParticles(Cuts::abspid == 411 || Cuts::abspid == 421 || Cuts::abspid == 413
                                || Cuts::abspid == 423 || Cuts::abspid == 431 || Cuts::abspid == 4122),
              "UFS");
      // histos
      for (unsigned int ix = 0; ix < 4; ++ix) {
        book(_h_spect[ix], 1 + 2 * ix, 1, 1);
        book(_h_scale[ix], 2 + 2 * ix, 1, 1);
        book(_h_sigma[ix], 14, 1, 1 + ix);
        if (ix == 3) continue;
        book(_h_spect[4 + ix], 9 + ix, 1, 1);
      }
      book(_h_scale[4], 12, 1, 1);
      for (unsigned int ix = 0; ix < 2; ++ix) {
        for (unsigned int iy = 0; iy < 7; ++iy) {
          book(_h_sigma_br[ix][iy], 13, 1 + ix, 1 + iy);
        }
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      for (const Particle& p : ufs.particles()) {
        double Emax = 0.5 * sqrtS();
        double pmax = sqrt(sqr(Emax) - p.mass2());
        double xp = (p.E() + p.p3().mod()) / (Emax + pmax);
        if (p.abspid() == 413) {
          _h_spect[0]->fill(xp);
          _h_scale[0]->fill(xp);
          _h_spect[1]->fill(xp);
          _h_scale[1]->fill(xp);
          _h_sigma[2]->fill("10.55"s);
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][2]->fill(xp);
          }
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][3]->fill(xp);
          }
        }
        else if (p.abspid() == 421) {
          _h_spect[2]->fill(xp);
          _h_scale[2]->fill(xp);
          _h_sigma[0]->fill("10.55"s);
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][0]->fill(xp);
          }
        }
        else if (p.abspid() == 411) {
          _h_spect[3]->fill(xp);
          _h_scale[3]->fill(xp);
          _h_sigma[1]->fill("10.55"s);
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][1]->fill(xp);
          }
        }
        else if (p.abspid() == 431) {
          _h_spect[4]->fill(xp);
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][5]->fill(xp);
          }
        }
        else if (p.abspid() == 4122) {
          _h_spect[5]->fill(xp);
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][6]->fill(xp);
          }
        }
        else if (p.abspid() == 423) {
          _h_spect[6]->fill(xp);
          _h_scale[4]->fill(xp);
          _h_sigma[3]->fill("10.55"s);
          for (unsigned int ix = 0; ix < 2; ++ix) {
            _h_sigma_br[ix][4]->fill(xp);
          }
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      const vector<double> br = {0.677 * 0.03947, 0.677 * 0.0822, 0.03947, 0.0938, 0.045, 0.0628, 0.03947};
      const double fact1 = crossSection() / picobarn / sumOfWeights();
      const double fact2 = crossSection() / nanobarn / sumOfWeights() * sqr(sqrtS());
      for (unsigned int ix = 0; ix < 7; ++ix) {
        scale(_h_spect[ix], fact1 * br[ix]);
        if (ix > 4) continue;
        scale(_h_scale[ix], fact2);
        if (ix == 4) continue;
        scale(_h_sigma[ix], fact1 * 1e-3);
      }
      const double width[2] = {1., 0.5};
      for (unsigned int ix = 0; ix < 2; ++ix) {
        scale(_h_sigma_br[ix][0], width[ix] * fact1 * br[2]);
        scale(_h_sigma_br[ix][1], width[ix] * fact1 * br[3]);
        scale(_h_sigma_br[ix][2], width[ix] * fact1 * br[0]);
        scale(_h_sigma_br[ix][3], width[ix] * fact1 * br[1]);
        scale(_h_sigma_br[ix][4], width[ix] * fact1 * br[6]);
        scale(_h_sigma_br[ix][5], width[ix] * fact1 * br[4]);
        scale(_h_sigma_br[ix][6], width[ix] * fact1 * br[5]);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_spect[7], _h_scale[5], _h_sigma_br[2][7];
    BinnedHistoPtr<string> _h_sigma[4];
    /// @}
  };


  RIVET_DECLARE_PLUGIN(CLEO_1988_I22954);

}