Rivet analyses

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

namespace Rivet {


  /// @brief Xi_c' spectrum
  class BABAR_2007_I722622 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2007_I722622);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // Initialise and register projections
      declare(UnstableParticles(), "UFS");
      // Book histograms
      book(_h_p_0, 3, 1, 2);
      book(_h_p_p, 3, 1, 1);
      book(_h_ctheta, 4, 1, 1);
      book(_b_p, 1, 1, 1);
      book(_b_0, 1, 1, 2);
      book(_r_p, 2, 1, 1);
      book(_r_0, 2, 1, 2);
      book(_ups, "/TMP/ups");
    }

    void findChildren(Particle parent, unsigned int& nStable, Particles& Xi, unsigned int& nPi) {
      for (const Particle& p : parent.children()) {
        if (p.abspid() == PID::PIPLUS) {
          ++nPi;
          ++nStable;
        }
        else if (p.abspid() == PID::XIMINUS) {
          Xi.push_back(p);
          ++nStable;
        }
        else if (!p.children().empty()) {
          findChildren(p, nStable, Xi, nPi);
        }
        else {
          ++nStable;
        }
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      static const int id0 = 4312, idp = 4322;
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      bool ups = !ufs.particles(Cuts::pid == 300553).empty();
      if (ups) _ups->fill();
      for (const Particle& p : ufs.particles(Cuts::abspid == idp or Cuts::abspid == id0)) {
        int idXic = 4232;
        if (p.abspid() == idp) {
          _h_p_p->fill(p.momentum().p3().mod());
        }
        else {
          _h_p_0->fill(p.momentum().p3().mod());
          idXic = 4132;
        }
        // first find the Xi_c in the decay
        if (p.children().size() != 2) continue;
        int sign = p.pid() / p.abspid();
        Particle Xi_c;
        if (p.children()[0].pid() == sign * idXic && p.children()[1].pid() == 22) {
          Xi_c = p.children()[0];
        }
        else if (p.children()[1].pid() == sign * idXic && p.children()[0].pid() == 22) {
          Xi_c = p.children()[1];
        }
        else
          continue;
        // and the children of the Xi_c
        Particles Xi;
        unsigned int nStable(0), nPi(0);
        findChildren(Xi_c, nStable, Xi, nPi);
        if (Xi.size() != 1) continue;
        if (!(p.abspid() == idp && nPi == 2 && nStable == 3)
            && !(p.abspid() == id0 && nPi == 1 && nStable == 2))
          continue;
        if (p.abspid() == idp) {
          if (ups)
            _b_p->fill();
          else
            _r_p->fill("10.58"s);
        }
        else {
          if (ups)
            _b_0->fill();
          else
            _r_0->fill("10.58"s);
          idXic = 4132;
        }
        // boost to Xi'_c rest frame
        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
        FourMomentum pXic = boost1.transform(Xi_c.momentum());
        FourMomentum pXi = boost1.transform(Xi[0].momentum());
        // then to Xi_c rest frame
        LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pXic.betaVec());
        Vector3 axis = pXic.p3().unit();
        FourMomentum pp = boost2.transform(pXi);
        // calculate angle
        double cTheta = pp.p3().unit().dot(axis);
        _h_ctheta->fill(cTheta);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h_p_0);
      normalize(_h_p_p);
      normalize(_h_ctheta);
      if (_ups->effNumEntries() != 0) {
        scale(_b_p, 0.5 / *_ups);
        scale(_b_0, 0.5 / *_ups);
      }
      scale(_r_p, crossSection() / sumOfWeights() / femtobarn);
      scale(_r_0, crossSection() / sumOfWeights() / femtobarn);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_p_0, _h_p_p, _h_ctheta;
    CounterPtr _b_p, _b_0;
    BinnedHistoPtr<string> _r_p, _r_0;
    CounterPtr _ups;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BABAR_2007_I722622);

}