Rivet analyses

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

namespace Rivet {


  /// @brief eta' -> eta pi0pi0
  class BESIII_2022_I2105430 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2105430);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // Initialise and register projections
      UnstableParticles ufs = UnstableParticles(Cuts::pid == PID::ETAPRIME);
      declare(ufs, "UFS");
      DecayedParticles ETA(ufs);
      ETA.addStable(PID::PI0);
      ETA.addStable(PID::K0S);
      ETA.addStable(PID::ETA);
      declare(ETA, "ETA");
      // histograms
      for (unsigned int ix = 0; ix < 4; ++ix) {
        if (ix < 2)
          book(_h[ix], 1 + ix, 1, 1);
        else
          book(_h[ix], "TMP/h_" + toString(ix + 1), refData(1 + ix, 1, 1));
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      static const map<PdgId, unsigned int>& mode = {{221, 1}, {111, 2}};
      DecayedParticles ETA = apply<DecayedParticles>(event, "ETA");
      // loop over particles
      for (unsigned int ix = 0; ix < ETA.decaying().size(); ++ix) {
        // select right decay mode
        if (!ETA.modeMatches(ix, 3, mode)) continue;
        const Particle& eta = ETA.decayProducts()[ix].at(221)[0];
        const Particles& pi0 = ETA.decayProducts()[ix].at(111);
        double s1 = (eta.momentum() + pi0[0].momentum()).mass2();
        double s2 = (eta.momentum() + pi0[1].momentum()).mass2();
        double s3 = (pi0[0].momentum() + pi0[1].momentum()).mass2();
        double mOut = eta.mass() + pi0[0].mass() + pi0[1].mass();
        double Q = ETA.decaying()[ix].mass() - mOut;
        double X = sqrt(3.) / 2. / ETA.decaying()[ix].mass() / Q * fabs(s1 - s2);
        double Y = mOut / 2. / Q / pi0[0].mass() / ETA.decaying()[ix].mass()
                * (sqr(ETA.decaying()[ix].mass() - eta.mass()) - s3)
            - 1.;
        _h[0]->fill(X);
        _h[1]->fill(Y);
        _h[2]->fill(s1);
        _h[2]->fill(s2);
        _h[3]->fill(s3);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // phase space volumes in the different bins (PR has python program which computed these)
      vector<double> phsp[2] = {
          {1.29921180560319,   1.9664163781008512, 2.437411090783833,  2.8167198110075895, 3.1384424839069816,
           3.4191598743922667, 3.6684846236204653, 3.892603483754069,  4.0957879360431315, 4.281140528359467,
           4.451003067955206,  4.607197285803388,  4.751175204530485,  4.88411754014731,   5.007000581045515,
           5.120643099873668,  5.225740150956019,  5.322887984241281,  5.412602779591596,  5.4953349814658825,
           5.571480436748339,  5.641389167221586,  5.705372363331301,  5.763708020726753,  5.8166455273826925,
           5.864409429427267,  5.907202547019624,  5.945208570543684,  5.978594237253912,  6.007511166123523,
           6.032097411815989,  6.0524787859064215, 6.068769983650643,  6.081075546972378,  6.08949068836251,
           6.094101995655494,  6.094988033869798,  6.092219857249216,  6.085861442152403,  6.075970049380721,
           6.062596522809926,  6.045785529721536,  6.025575746951852,  6.00199999583944,   5.975085327911516,
           5.944853062268092,  5.911318774664203,  5.874492237320916,  5.834377307479565,  5.790971761612199,
           5.744267070970989,  5.6942481127479825, 5.640892809460393,  5.584171687195311,  5.524047340938203,
           5.460473792237251,  5.39339572074258,   5.32274754646675,   5.248452333616691,  5.1704204791023605,
           5.088548138714816,  5.002715330614686,  4.91278363793704,   4.818593408199339,  4.719960314146789,
           4.616671094724371,  4.50847822999662,   4.395093210678961,  4.276177926624657,  4.151333494932917,
           4.02008553685901,   3.881864423721346,  3.735978221227386,  3.5815747385987744, 3.4175867864493514,
           3.2426505513681274, 3.0549789164129533, 2.8521549298248745, 2.6307734449059827, 2.3857664628727737,
           2.10898208394632,   1.7856475571555508, 1.382647074122462,  0.7594829443647874},
          {2.3997722961635115, 4.231978624872276,  5.395770701813711,  6.303248399189668,  7.05746502255801,
           7.705360260187261,  8.273295064456326,  8.777930815906547,  9.230644528844046,  9.63964259335122,
           10.011092003070875, 10.349776928006058, 10.659504191013701, 10.943366292678283, 11.203918989983865,
           11.443305233235264, 11.66334412312612,  11.865596310404976, 12.051413085183789, 12.221973898377072,
           12.378315503623586, 12.521354914670614, 12.65190772145681,  12.770702870483907, 12.878394715019976,
           12.975572931081437, 13.062770746227885, 13.140471820760519, 13.209116042291283, 13.269104436365959,
           13.32080335209865,  13.364548048606796, 13.400645782615129, 13.429378477926482, 13.45100504209833,
           13.465763383568648, 13.473872172872744, 13.475532383910185, 13.470928645028556, 13.46023042466219,
           13.443593072150625, 13.421158730971719, 13.393057138809006, 13.35940632651546,  13.320313226043833,
           13.275874195713712, 13.22617546971762,  13.171293537484289, 13.11129545737622,  13.046239108166182,
           12.9761733807838,   12.901138311921745, 12.821165160216768, 12.736276424849162, 12.646485805511595,
           12.551798101758047, 12.45220904872841,  12.347705085120998, 12.238263048016853, 12.123849787699887,
           12.004421693910436, 11.879924122946402, 11.750290712597788, 11.61544256895382,  11.475287305511241,
           11.329717910545805, 11.178611413135048, 11.021827311202296, 10.859205716037735, 10.690565156327159,
           10.515699969937556, 10.334377192397056, 10.146332825512,    9.951267335531437,  9.748840184294515,
           9.538663133883416,  9.320291977999595,  9.093216230262732,  8.856846123359322,  8.610496015674501,
           8.3533629186997,    8.084498274135298,  7.802770195530031,  7.506811917146632,  7.194949741305242,
           6.865099532725173,  6.5146131262832965, 6.140041338435136,  5.736750353635648,  5.298262152423473,
           4.81502775442034,   4.271887046472109,  3.6419114500505065, 2.8669813619285884}};
      // normalize the histograms
      for (unsigned int ix = 0; ix < 4; ++ix) normalize(_h[ix]);
      // last two plots  convert to scatter and normalize to phase space volume in bin
      for (unsigned int ix = 0; ix < 2; ++ix) {
        // convert to scatter
        Estimate1DPtr tmp;
        book(tmp, 3 + ix, 1, 1);
        barchart(_h[ix + 2], tmp);
        double step[2] = {0.0025, 0.001};
        // divide by phase space volume
        for (unsigned int ip = 0; ip < tmp->numBins(); ++ip) {
          tmp->bin(ip + 1).scale(1. / phsp[ix][ip] / step[ix]);
        }
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h[4];
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BESIII_2022_I2105430);

}