Rivet analyses

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

namespace Rivet {


  /// @brief tau polarization in B -> D* tau nu_tau
  class BELLE_2018_I1621272 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2018_I1621272);


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

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      declare(UnstableParticles(), "UFS");

      // Book histograms
      book(_h_pi, "/TMP/PI", 20, -1., 1.);
      book(_h_rho, "/TMP/RHO", 20, -1., 1.);
    }

    void findChildren(const Particle& p,
                      int& sign,
                      unsigned int& nprod,
                      Particles& Dstar,
                      Particles& tau,
                      Particles& nu) {
      for (const Particle& child : p.children()) {
        if (child.pid() == -sign * 413 || child.pid() == -sign * 423) {
          ++nprod;
          Dstar.push_back(child);
        }
        else if (child.pid() == -sign * 15) {
          ++nprod;
          tau.push_back(child);
        }
        else if (child.pid() == sign * 16) {
          ++nprod;
          nu.push_back(child);
        }
        else if (child.pid() == 22)
          continue;
        else if (child.children().empty() || child.pid() == 111 || child.pid() == 221 || child.pid() == 331) {
          ++nprod;
        }
        else {
          findChildren(child, sign, nprod, Dstar, tau, nu);
        }
      }
    }

    void findTau(const Particle& p,
                 int& sign,
                 unsigned int& nprod,
                 Particles& piP,
                 Particles& pi0,
                 Particles& nu) {
      for (const Particle& child : p.children()) {
        if (child.pid() == 111) {
          ++nprod;
          pi0.push_back(child);
        }
        else if (child.pid() == sign * 211) {
          ++nprod;
          piP.push_back(child);
        }
        else if (child.pid() == -sign * 16) {
          ++nprod;
          nu.push_back(child);
        }
        else if (child.pid() == 22)
          continue;
        else if (child.children().empty() || child.pid() == 221 || child.pid() == 331) {
          ++nprod;
        }
        else {
          findTau(child, sign, nprod, piP, pi0, nu);
        }
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Loop over B0 mesons
      for (const Particle& p : apply<UnstableParticles>(event, "UFS")
                                   .particles(Cuts::abspid == PID::B0 or Cuts::abspid == PID::BPLUS)) {
        // find the B decay
        int sign = p.pid() / p.abspid();
        unsigned int nprod = 0;
        Particles Dstar, tau, nu;
        findChildren(p, sign, nprod, Dstar, tau, nu);
        if (nprod != 3 || Dstar.size() != 1 || tau.size() != 1 || nu.size() != 1) continue;
        // check decay
        if (p.pid() == PID::B0) {
          if (Dstar[0].pid() != -sign * 413) vetoEvent;
        }
        else if (p.pid() == PID::BPLUS) {
          if (Dstar[0].pid() != -sign * 423) vetoEvent;
        }
        // find the tau decay
        nprod = 0;
        nu.clear();
        Particles piP, pi0;
        findTau(tau[0], sign, nprod, piP, pi0, nu);
        if (nu.size() != 1) continue;
        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
        FourMomentum ptau = boost1.transform(tau[0].momentum());
        LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(ptau.betaVec());
        // pion mode
        if (nprod == 2 && piP.size() == 1) {
          FourMomentum pPi = boost2.transform(boost1.transform(piP[0].momentum()));
          double cTheta = pPi.p3().unit().dot(ptau.p3().unit());
          _h_pi->fill(cTheta);
        }
        // rho mode
        else if (nprod == 3 && piP.size() == 1 && pi0.size() == 1) {
          FourMomentum pRho = boost2.transform(boost1.transform(piP[0].momentum() + pi0[0].momentum()));
          double cTheta = pRho.p3().unit().dot(ptau.p3().unit());
          _h_rho->fill(cTheta);
        }
      }
    }

    pair<double, double> calcAlpha(Histo1DPtr hist) {
      if (hist->numEntries() == 0.) return make_pair(0., 0.);
      double sum1(0.), sum2(0.);
      for (auto& bin : hist->bins()) {
        double Oi = bin.sumW();
        if (Oi == 0.) continue;
        double ai = 0.5 * (bin.xMax() - bin.xMin());
        double bi = 0.5 * ai * (bin.xMax() + bin.xMin());
        double Ei = bin.errW();
        sum1 += sqr(bi / Ei);
        sum2 += bi / sqr(Ei) * (Oi - ai);
      }
      return make_pair(sum2 / sum1, sqrt(1. / sum1));
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h_pi);
      normalize(_h_rho);
      // the polarization
      Estimate0DPtr _h_alpha;
      book(_h_alpha, 1, 1, 2);
      pair<double, double> alpha_pi = calcAlpha(_h_pi);
      pair<double, double> alpha_rho = calcAlpha(_h_rho);
      // 0.45 factor for rho
      alpha_rho.first /= 0.46;
      alpha_rho.second /= 0.46;
      pair<double, double> alpha;
      alpha.first = (alpha_pi.first * sqr(alpha_rho.second) + alpha_rho.first * sqr(alpha_pi.second))
          / (sqr(alpha_pi.second) + sqr(alpha_rho.second));
      alpha.second = alpha_pi.second * alpha_rho.second / sqrt(sqr(alpha_pi.second) + sqr(alpha_rho.second));
      _h_alpha->set(alpha.first, alpha.second);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_pi, _h_rho;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BELLE_2018_I1621272);

}