Rivet analyses

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

namespace Rivet {


  /// @brief e+ e- -> double charmonium
  class BELLE_2004_I653673 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2004_I653673);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare("FS", FinalState());
      declare("UFS",
              UnstableParticles(Cuts::pid == 443 || Cuts::pid == 100443 || Cuts::pid == 441
                                || Cuts::pid == 10441 || Cuts::pid == 100441));
      // histograms
      for (unsigned int ix = 0; ix < 3; ++ix) {
        for (unsigned int iy = 0; iy < 3; ++iy) {
          book(_p[ix][iy], "TMP/p_" + toString(ix + 1) + "_" + toString(iy + 1),
               refData<YODA::BinnedEstimate<string>>(3, 1 + ix, 1 + iy));
          if (ix == 2) continue;
          book(_h_sigma[ix][iy], 1 + ix, 1, 1 + 2 * iy);
          book(_h_angle[ix][iy], 4, 1 + ix, 1 + iy);
        }
      }
    }

    void findChildren(const Particle& p, map<long, int>& nRes, int& ncount, unsigned int& nCharged) {
      for (const Particle& child : p.children()) {
        if (child.children().empty()) {
          --nRes[child.pid()];
          --ncount;
          if (PID::isCharged(p.pid())) ++nCharged;
        }
        else {
          findChildren(child, nRes, ncount, nCharged);
        }
      }
    }

    double helicityAngle(const Particle& p) const {
      if (p.children().size() != 2) return 10.;
      if (p.children()[0].abspid() != PID::MUON || p.children()[0].pid() != -p.children()[1].pid()) {
        return 10.;
      }
      Vector3 axis = p.p3().unit();
      const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
      for (const Particle& child : p.children()) {
        if (child.pid() != PID::MUON) continue;
        return axis.dot(boost.transform(child.momentum()).p3().unit());
      }
      return 10;
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // final state particles
      const FinalState& fs = apply<FinalState>(event, "FS");
      map<long, int> nCount;
      int ntotal(0);
      for (const Particle& p : fs.particles()) {
        nCount[p.pid()] += 1;
        ++ntotal;
      }
      // loop over J/psi and psi(2S)
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      bool matched = false;
      for (const Particle& p : ufs.particles(Cuts::pid == 443 || Cuts::pid == 100443)) {
        if (p.children().empty()) continue;
        map<long, int> nRes = nCount;
        int ncount = ntotal;
        unsigned int nCharged = 0;
        findChildren(p, nRes, ncount, nCharged);
        // eta_c, chi_c0, eta_c(2S)
        for (const Particle& p2 :
             ufs.particles(Cuts::pid == 441 || Cuts::pid == 10441 || Cuts::pid == 100441)) {
          map<long, int> nResB = nRes;
          int ncountB = ncount;
          unsigned int nChargedB = 0;
          findChildren(p2, nResB, ncountB, nChargedB);
          if (ncountB != 0) continue;
          matched = true;
          for (const auto& val : nResB) {
            if (val.second != 0) {
              matched = false;
              break;
            }
          }
          if (matched) {
            unsigned int ipsi = p.pid() / 100000;
            unsigned int ieta = p2.pid() / 10000;
            if (ieta > 1) ieta = 2;
            // fill the cross sections
            if ((ipsi == 0 && nChargedB > 2) || (ipsi == 1 && nChargedB > 0)) {
              _h_sigma[ipsi][ieta]->fill(_ecms);
            }
            if (ipsi > 0) break;
            // angular dists for J/psi only
            // production
            const double cProd = p.p3().z() / p.p3().mod();
            _h_angle[0][ieta]->fill(abs(cProd));
            _p[0][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cProd)));
            _p[2][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cProd)));
            // helicity angle
            const double cHel = helicityAngle(p);
            if (cHel > 1.) break;
            _h_angle[1][ieta]->fill(abs(cHel));
            _p[1][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cHel)));
            _p[2][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cHel)));
            break;
          }
        }
        if (matched) break;
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int ix = 0; ix < 2; ++ix) {
        scale(_h_sigma[ix], crossSection() / sumOfWeights() / femtobarn);
        normalize(_h_angle[ix]);
      }
      // extract the alpha parameters
      for (unsigned int ix = 0; ix < 3; ++ix) {
        for (unsigned int iy = 0; iy < 3; ++iy) {
          const double val = _p[ix][iy]->bin(1).mean(2);
          const double err = _p[ix][iy]->bin(1).stdErr(2);
          BinnedEstimatePtr<string> tmp;
          book(tmp, 3, 1 + ix, 1 + iy);
          const double alpha = 3. * val / (1 - val);
          const double error = 3. / sqr(1. - val) * err;
          tmp->bin(1).set(alpha, error);
        }
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    BinnedHistoPtr<string> _h_sigma[2][3];
    Histo1DPtr _h_angle[2][3];
    BinnedProfilePtr<string> _p[3][3];
    string _ecms = "10.6";
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BELLE_2004_I653673);

}