Rivet analyses

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

namespace Rivet {


  /// @brief Charm meson spectra in bottom decays
  class BABAR_2006_I719111 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2006_I719111);


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

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

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

      // Book histograms
      book(_h_Bm_D0      , 1, 1, 1);
      book(_h_Bm_Dbar0   , 1, 1, 2);
      book(_h_Bm_Dp      , 2, 1, 1);
      book(_h_Bm_Dm      , 2, 1, 2);
      book(_h_Bm_Dsp     , 3, 1, 1);
      book(_h_Bm_Dsm     , 3, 1, 2);
      book(_h_Bm_lam     , 4, 1, 1);
      book(_h_Bm_lbar    , 4, 1, 2);
      book(_h_Bbar0_D0   , 5, 1, 1);
      book(_h_Bbar0_Dbar0, 5, 1, 2);
      book(_h_Bbar0_Dp   , 6, 1, 1);
      book(_h_Bbar0_Dm   , 6, 1, 2);
      book(_h_Bbar0_Dsp  , 7, 1, 1);
      book(_h_Bbar0_Dsm  , 7, 1, 2);
      book(_h_Bbar0_lam  , 8, 1, 1);
      book(_h_Bbar0_lbar , 8, 1, 2);

      book(_c_Bm   ,"/TMP/Bm");
      book(_c_Bbar0,"/TMP/B0");
    }

    void findDecayProducts(Particle p,
               Particles & D0, Particles & Dbar0,
               Particles & Dp, Particles & Dm,
               Particles & Dsp, Particles & Dsm,
               Particles & lam, Particles & lbar) {
      for(const Particle & child : p.children()) {
    if(child.pid()==421)
      D0.push_back(child);
    else if(child.pid()==-421)
      Dbar0.push_back(child);
    else if(child.pid()==411)
      Dp.push_back(child);
    else if(child.pid()==-411)
      Dm.push_back(child);
    else if(child.pid()==431)
      Dsp.push_back(child);
    else if(child.pid()==-431)
      Dsm.push_back(child);
    else if(child.pid()==4122)
      lam.push_back(child);
    else if(child.pid()==-4122)
      lbar.push_back(child);
    else if(!child.children().empty())
      findDecayProducts(child,D0,Dbar0,Dp,Dm, Dsp, Dsm, lam, lbar);
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      for(const Particle & p : ufs.particles(Cuts::abspid==511 || Cuts::abspid==521)) {
        if(p.abspid()==511)
          _c_Bbar0->fill();
        else if(p.abspid()==521)
          _c_Bm   ->fill();
        Particles D0, Dbar0, Dp, Dm, Dsp, Dsm, lam, lbar;
        findDecayProducts(p,D0,Dbar0,Dp,Dm, Dsp, Dsm, lam, lbar);
        LorentzTransform boost =
          LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
        if(p.pid()>0) {
          swap(D0, Dbar0);
          swap(Dp, Dm);
          swap(Dsp, Dsm);
          swap(lam, lbar);
        }
        for(const Particle & child : D0) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_D0   ->fill(pChild);
          else
            _h_Bbar0_D0->fill(pChild);
        }
        for(const Particle & child : Dbar0) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_Dbar0   ->fill(pChild);
          else
            _h_Bbar0_Dbar0->fill(pChild);
        }
        for(const Particle & child : Dp) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_Dp   ->fill(pChild);
          else
            _h_Bbar0_Dp->fill(pChild);
        }
        for(const Particle & child : Dm) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_Dm   ->fill(pChild);
          else
            _h_Bbar0_Dm->fill(pChild);
        }
        for(const Particle & child : Dsp) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_Dsp   ->fill(pChild);
          else
            _h_Bbar0_Dsp->fill(pChild);
        }
        for(const Particle & child : Dsm) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_Dsm   ->fill(pChild);
          else
            _h_Bbar0_Dsm->fill(pChild);
        }
        for(const Particle & child : lam) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_lam   ->fill(pChild);
          else
            _h_Bbar0_lam->fill(pChild);
        }
        for(const Particle & child : lbar) {
          double pChild = boost.transform(child.momentum()).p3().mod();
          if(p.abspid()==521)
            _h_Bm_lbar   ->fill(pChild);
          else
            _h_Bbar0_lbar->fill(pChild);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h_Bm_D0      ,100./ _c_Bm->val());
      scale(_h_Bm_Dbar0   ,100./ _c_Bm->val());
      scale(_h_Bm_Dp      ,100./ _c_Bm->val());
      scale(_h_Bm_Dm      ,100./ _c_Bm->val());
      scale(_h_Bm_Dsp     ,100./ _c_Bm->val());
      scale(_h_Bm_Dsm     ,100./ _c_Bm->val());
      scale(_h_Bm_lam     ,100./ _c_Bm->val());
      scale(_h_Bm_lbar    ,100./ _c_Bm->val());
      scale(_h_Bbar0_D0   ,100./ _c_Bbar0->val());
      scale(_h_Bbar0_Dbar0,100./ _c_Bbar0->val());
      scale(_h_Bbar0_Dp   ,100./ _c_Bbar0->val());
      scale(_h_Bbar0_Dm   ,100./ _c_Bbar0->val());
      scale(_h_Bbar0_Dsp  ,100./ _c_Bbar0->val());
      scale(_h_Bbar0_Dsm  ,100./ _c_Bbar0->val());
      scale(_h_Bbar0_lam  ,100./ _c_Bbar0->val());
      scale(_h_Bbar0_lbar ,100./ _c_Bbar0->val());
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_Bm_D0,_h_Bm_Dbar0,_h_Bm_Dp,_h_Bm_Dm,_h_Bm_Dsp,_h_Bm_Dsm,
      _h_Bm_lam,_h_Bm_lbar,_h_Bbar0_D0;
    Histo1DPtr _h_Bbar0_Dbar0,_h_Bbar0_Dp,_h_Bbar0_Dm,
      _h_Bbar0_Dsp,_h_Bbar0_Dsm,_h_Bbar0_lam,_h_Bbar0_lbar;

    CounterPtr _c_Bm ,_c_Bbar0;
    /// @}

  };


  RIVET_DECLARE_PLUGIN(BABAR_2006_I719111);


}