Rivet analyses

Cross sections for e+e → b → (Ds+, D0)X for $\sqrt{s}$ between 10.63 and 11.02 GeV

Experiment: BELLE (KEKB)

Inspire ID: 2660525

Status: VALIDATED NOHEPDATA

Authors: - Peter Richardson

References: - arxiv:2305.10098 [hep-ex]

Beams: * *

Beam energies: ANY

Run details: - Any process producing Upsilon(4,5S) for decays and e+e–> hadrons for the cross sections

The cross sections for e+e → b → (Ds+, D0)X are measured for $\sqrt{s}$ between 10.63 and 11.02 GeV. The cross sections for e+e → b → BX and e+e → b → Bs0s0X are then extracted. In addition the Ds+, D0 spectra in Υ(4, 5S) decays are measured. There are two modes, one DECAY (the default), for the spectra in the decays and a second SIGMA for the cross sections.

Source code:BELLE_2023_I2660525.cc

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

namespace Rivet {


  /// @brief D_s D0 production in e+e-
  class BELLE_2023_I2660525 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2023_I2660525);


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

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

      // Initialise and register projections
      declare(UnstableParticles(), "UFS");
      _mode = getOption("MODE") == "SIGMA";
      // Upsilon decays
      if (_mode==0) {
        for (size_t ix=0; ix<2; ++ix) {
          book(_c[ix],"TMP/n_Ups_"+toString(ix+4));
          book(_h_x [0][ix], 3, 1, ix+3);
          book(_h_x [1][ix], 3, 1, ix+1);
          for (size_t iy=0; iy<2; ++iy) {
            book(_h_br[ix][iy], 1, 1, 2*ix+iy+1);
          }
        }
        book(_h_fs, 1, 1, 5);
      }
      else {
        for (size_t ix=0; ix<4; ++ix) {
          book(_sigma[ix], 2, 1, 1+ix);
        }
        for (const string& en : _sigma[0].binning().edges<0>()) {
          double eval = stod(en)*GeV;
          if (isCompatibleWithSqrtS(eval)) {
            _sqs = en; break;
          }
        }
        raiseBeamErrorIf(_sqs.empty());
      }
    }

    void findDecayProducts(const Particle& p, bool& hasBs, Particles& d0, Particles& ds) const {
      for (const Particle& child : p.children()) {
        hasBs |= child.abspid()==531;
        if (child.abspid()==431) {
          ds.push_back(child);
        }
        else if (child.abspid()==421) {
          d0.push_back(child);
        }
        else if (!child.children().empty()) {
          findDecayProducts(child,hasBs,d0,ds);
        }
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Upsilon(4,5S) decay
      if (_mode==0) {
        const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
        Particles upsilons = ufs.particles(Cuts::pid==300553 || Cuts::pid==400553 || Cuts::pid==9000553);
        for (const Particle& ups : upsilons) {
          unsigned int iups = ups.pid()==300553 ? 0 : 1;
          _c[iups]->fill();
          LorentzTransform cms_boost;
          if (ups.p3().mod() > 1*MeV)
            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.mom().betaVec());
          Particles ds,d0;
          bool hasBs=false;
          findDecayProducts(ups, hasBs, d0, ds);
          if(iups==1&&hasBs) _h_fs->fill();
          for (const Particle & p : ds) {
            FourMomentum p2 = cms_boost.transform(p.mom());
            double x = p2.p3().mod()/sqrt(0.25*ups.mass2()-p2.mass2());
            _h_br[iups][0]->fill();
            _h_x [iups][0]->fill(x  );

          }
          for (const Particle& p : d0) {
            FourMomentum p2 = cms_boost.transform(p.mom());
            double x = p2.p3().mod()/sqrt(0.25*ups.mass2()-p2.mass2());
            _h_br[iups][1]->fill();
            _h_x [iups][1]->fill(x);
          }
        }
      }
      else {
        const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
        Particles B  = ufs.particles(Cuts::abspid==511 || Cuts::abspid==521);
        Particles Bs = ufs.particles(Cuts::abspid==531);
        // require B mesons
        if (B.empty() && Bs.empty()) vetoEvent;
        if (!B.empty()) _sigma[3]->fill(_sqs);
        Particles Ds = ufs.particles(Cuts::abspid==431);
        for (unsigned int ix=0;ix<Ds.size();++ix) {
          _sigma[0]->fill(_sqs);
        }
        Particles D0 = ufs.particles(Cuts::abspid==421);
        for (unsigned int ix=0;ix<D0.size();++ix) {
          _sigma[1]->fill(_sqs);
        }
        if (!Bs.empty() && !Ds.empty()) {
          for (const Particle& p : Bs) {
            if (p.children().size()==1) continue;
            Particles ds,d0;
            bool hasBs=false;
            findDecayProducts(p, hasBs, d0, ds);
            for (unsigned int ix=0;ix<ds.size();++ix) _sigma[2]->fill(_sqs);
          }
        }
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      if (_mode==0) {
        for (size_t ix=0; ix<2; ++ix) {
          if (_c[ix]->numEntries()==0) continue;
          scale(_h_x[ix], 1.0/ *_c[ix]);
          // brs in % and at 4S /2 as br is for B decays
          if (ix==0) scale(_h_br[ix], 50  / *_c[ix]);
          else       scale(_h_br[ix], 100./ *_c[ix]);
          if (ix==1) scale(_h_fs,100./ *_c[1]);
        }
      }
      else {
        scale(_sigma, crossSection()/ sumOfWeights() /picobarn);
      }
    }
    /// @}


    /// @name Histograms
    /// @{
    size_t _mode;
    Histo1DPtr _h_x[2][2];
    CounterPtr _h_br[2][2],_h_fs,_c[2];
    BinnedHistoPtr<string> _sigma[4];
    string _sqs = "";
    /// @}


  };


  RIVET_DECLARE_PLUGIN(BELLE_2023_I2660525);

}