Rivet analyses

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

namespace Rivet {


  /// @brief D*(sJ) production
  class BABAR_2006_I714447 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2006_I714447);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare(UnstableParticles(), "UFS");
      // histograms
      // dists
      book(_s_2317, 2, 1, 1);
      book(_s_2460_1, 3, 1, 1);
      book(_s_2460_2, 4, 1, 1);
      book(_hel, 5, 1, 1);
      // rates
      book(_r_2317, 1, 1, 1);
      book(_r_2460_1, 1, 1, 2);
      book(_r_2460_2, 1, 1, 3);
      book(_r_2460_3, 1, 1, 4);
      book(_r_2460_4, 1, 1, 5);
      book(_r_2536, 1, 1, 6);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      for (const Particle& Ds :
           apply<UnstableParticles>(event, "UFS")
               .particles(Cuts::pid == 10433 or Cuts::pid == 20433 or Cuts::pid == 10431)) {
        // cut of 3.2 GeV of the momentum
        double mom = Ds.momentum().p3().mod();
        if (mom <= 3.2) continue;
        if (Ds.abspid() == 10431) {
          _s_2317->fill(mom);
        }
        else if (Ds.abspid() == 20433) {
          _s_2460_1->fill(mom);
          _s_2460_2->fill(mom);
        }
        int sign = Ds.pid() / Ds.abspid();
        Particle DsStar;
        bool doHelicity = false;
        // identify the decay mode
        if (Ds.children().size() == 2) {
          if (Ds.children()[0].pid() == sign * 431 && Ds.children()[1].pid() == 111) {
            if (Ds.abspid() == 10431) _r_2317->fill("10.6"s);
          }
          else if (Ds.children()[1].pid() == sign * 431 && Ds.children()[0].pid() == 111) {
            if (Ds.abspid() == 10431) _r_2317->fill("10.6"s);
          }
          else if (Ds.children()[0].pid() == sign * 433 && Ds.children()[1].pid() == 111) {
            if (Ds.abspid() == 20433) {
              _r_2460_3->fill("10.6"s);
              DsStar = Ds.children()[0];
              doHelicity = true;
            }
          }
          else if (Ds.children()[1].pid() == sign * 433 && Ds.children()[0].pid() == 111) {
            if (Ds.abspid() == 20433) {
              _r_2460_3->fill("10.6"s);
              DsStar = Ds.children()[1];
              doHelicity = true;
            }
          }
          else if (Ds.children()[0].pid() == sign * 431 && Ds.children()[1].pid() == 22) {
            if (Ds.abspid() == 20433) _r_2460_1->fill("10.6"s);
          }
          else if (Ds.children()[1].pid() == sign * 431 && Ds.children()[0].pid() == 22) {
            if (Ds.abspid() == 20433) _r_2460_1->fill("10.6"s);
          }
          if (doHelicity && DsStar.children().size() == 2) {
            Particle gamma;
            doHelicity = false;
            for (const Particle& p : DsStar.children()) {
              if (p.pid() == 22) {
                doHelicity = true;
                gamma = p;
              }
            }
            if (doHelicity) {
              LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(Ds.momentum().betaVec());
              FourMomentum pDsStar = boost.transform(DsStar.momentum());
              FourMomentum pGamma = boost.transform(gamma.momentum());
              LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDsStar.betaVec());
              FourMomentum pGamma2 = boost2.transform(pGamma);
              double cTheta = pGamma2.p3().unit().dot(pDsStar.p3().unit());
              _hel->fill(cTheta);
            }
          }
        }
        Particles D, pip, pi0, gamma;
        unsigned int nstable(0);
        findDecayProducts(Ds, D, pip, pi0, gamma, nstable);
        if (nstable == 3) {
          if (D.size() == 1 && pip.size() == 2) {
            if (Ds.abspid() == 20433) {
              _r_2460_4->fill("10.6"s);
            }
            else if (Ds.abspid() == 10433) {
              _r_2536->fill("10.6"s);
            }
          }
          else if (D.size() == 1 && pi0.size() == 1 && gamma.size() == 1) {
            if (Ds.abspid() == 20433) _r_2460_2->fill("10.6"s);
          }
        }
      }
    }

    void findDecayProducts(Particle parent,
                           Particles& Ds,
                           Particles& pip,
                           Particles& pi0,
                           Particles& gamma,
                           unsigned int& nstable) {
      for (const Particle& p : parent.children()) {
        if (p.pid() == 111) {
          pi0.push_back(p);
          ++nstable;
        }
        else if (p.pid() == 22) {
          gamma.push_back(p);
          ++nstable;
        }
        else if (p.abspid() == 211) {
          pip.push_back(p);
          ++nstable;
        }
        else if (p.abspid() == 431) {
          Ds.push_back(p);
          ++nstable;
        }
        else if (p.children().empty()) {
          ++nstable;
        }
        else
          findDecayProducts(p, Ds, pip, pi0, gamma, nstable);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_s_2317, 1., false);
      normalize(_s_2460_1, 1., false);
      normalize(_s_2460_2, 1., false);
      normalize(_hel, 1., false);
      scale(_r_2317, crossSection() / femtobarn / sumOfWeights());
      scale(_r_2460_1, crossSection() / femtobarn / sumOfWeights());
      scale(_r_2460_2, crossSection() / femtobarn / sumOfWeights());
      scale(_r_2460_3, crossSection() / femtobarn / sumOfWeights());
      scale(_r_2460_4, crossSection() / femtobarn / sumOfWeights());
      scale(_r_2536, crossSection() / femtobarn / sumOfWeights());
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _s_2317, _s_2460_1, _s_2460_2, _hel;
    BinnedHistoPtr<string> _r_2317, _r_2460_1, _r_2460_2, _r_2460_3, _r_2460_4, _r_2536;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BABAR_2006_I714447);

}