Rivet analyses

e+e → Dπ+π for $\sqrt{s}=4.0854$ to 4.6 GeV

Experiment: BESIII (BEPC)

Inspire ID: 1725786

Status: VALIDATED NOHEPDATA SINGLEWEIGHT

Authors: - Peter Richardson

References: - Phys.Rev.D 100 (2019) 3, 032005

Beams: e+ e-

Beam energies: (2.0, 2.0); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.2, 2.2); (2.2, 2.2); (2.2, 2.2); (2.2, 2.2); (2.2, 2.2); (2.3, 2.3); (2.3, 2.3); (2.3, 2.3)GeV

Run details: - e+ e- to hadrons

Analysis of the resonant contributions to e+e → Dπ+π for $\sqrt{s}=4.0854$ to 4.6 GeV, where D = D+D, D00. The cross sections for the production of π+πψ(3770) and D1(242000+c.c. are measured as a function of energy and the mass distrinutions are implemented for three energies.

Source code:BESIII_2019_I1725786.cc

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

namespace Rivet {


  /// @brief e+ e- > D D pi pi
  class BESIII_2019_I1725786 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2019_I1725786);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare(FinalState(), "FS");
      // histograms
      for (double eVal : allowedEnergies()) {

        const string en = toString(round(10 * eVal / MeV));
        if (isCompatibleWithSqrtS(eVal)) _sqs = en;

        if (en == "44156"s) {
          book(_h["Dpi0"], 1, 1, 1);
          book(_h["Dpi1"], 1, 1, 2);
          book(_h[en + "DD"], 2, 1, 3);
          book(_h[en + "mass0"], 3, 2, 1);
          book(_h[en + "mass1"], 3, 2, 2);
          book(_h[en + "mass2"], 3, 2, 3);
        }
        else if (en == "42580"s) {
          book(_h[en + "DD"], 2, 1, 1);
        }
        else if (en == "43583"s) {
          book(_h[en + "DD"], 2, 1, 2);
          book(_h[en + "mass0"], 3, 1, 1);
          book(_h[en + "mass1"], 3, 1, 2);
          book(_h[en + "mass2"], 3, 1, 3);
        }
        else if (en == "45995"s) {
          book(_h[en + "mass0"], 3, 3, 1);
          book(_h[en + "mass1"], 3, 3, 2);
          book(_h[en + "mass2"], 3, 3, 3);
        }
      }
      raiseBeamErrorIf(_sqs.empty());

      size_t ih = 0;
      for (size_t ix = 0; ix < 5; ++ix) {
        if (ix == 1) continue;
        book(_sigma[ix], 4 + ix, 1, 1);
        if (ix != 0 && ix != 2) continue;
        for (const string& en : _sigma[ix].binning().edges<0>()) {
          const double eval = stod(en) * GeV;
          if (isCompatibleWithSqrtS(eval)) {
            _edge[ih] = en;
            break;
          }
        }
        ++ih;
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      Particles fs = apply<FinalState>(event, "FS").particles();
      Particles DD, other;
      for (const Particle& p : fs) {
        Particle parent = p;
        while (!parent.parents().empty()) {
          parent = parent.parents()[0];
          if (parent.abspid() == PID::D0 || parent.abspid() == PID::DPLUS) break;
        }
        if (parent.abspid() != PID::D0 && parent.abspid() != PID::DPLUS) {
          other.push_back(p);
          continue;
        }
        bool found = false;
        for (const auto& D : DD) {
          // D already in list
          if (fuzzyEquals(D.mom(), parent.mom())) {
            found = true;
            break;
          }
        }
        if (!found) DD += parent;
      }
      //  d dbar + 2 other particles
      if (DD.size() != 2 || other.size() != 2) vetoEvent;
      // other particles pi+ pi-
      if (!(other[0].pid() == -other[1].pid() && other[0].abspid() == PID::PIPLUS)) vetoEvent;
      // d dbar pair
      if (DD[0].pid() != -DD[1].pid()) vetoEvent;
      if (other[0].pid() < 0) swap(other[0], other[1]);
      if (DD[0].pid() < 0) swap(DD[0], DD[1]);
      // fill the mass plots if needed
      if (_sqs == "44156"s && DD[0].pid() == 421) {
        _h["Dpi0"]->fill((DD[0].mom() + other[0].mom()).mass());
        _h["Dpi1"]->fill((DD[1].mom() + other[1].mom()).mass());
      }
      if (_sqs == "44156"s || _sqs == "42580"s || _sqs == "43583"s) {
        _h[_sqs + "DD"]->fill((DD[0].mom() + DD[1].mom()).mass());
      }
      const bool DstarP = (DD[0].pid() == 421 && DD[0].parents()[0].pid() == 413);
      const bool DstarM = (DD[1].pid() == -421 && DD[1].parents()[0].pid() == -413);
      const bool Dstar = DstarP || DstarM;
      if (_sqs == "44156"s || _sqs == "43583"s || _sqs == "45995"s) {
        if (DD[0].pid() == 421 && !Dstar) { // D0 pi+ pi- mass
          _h[_sqs + "mass0"]->fill((DD[0].mom() + other[0].mom() + other[1].mom()).mass());
        }
        if (DstarP) { // D*+ pi- mass
          _h[_sqs + "mass1"]->fill((DD[0].mom() + other[0].mom() + other[1].mom()).mass());
        }
        if (DD[0].pid() == 411) { // D+ pi+ pi- mass
          _h[_sqs + "mass2"]->fill((DD[0].mom() + other[0].mom() + other[1].mom()).mass());
        }
      }
      // now for the counters for the cross sections
      if (DD[0].parents()[0].pid() == 30443 && DD[1].parents()[0].pid() == 30443) {
        _sigma[0]->fill(_edge[0]);
      }
      else if (DD[0].pid() == 421
               && (DD[0].parents()[0].pid() == 10423 || DD[1].parents()[0].pid() == -10423)) {
        _sigma[2]->fill(_edge[1]);
      }
      else if (DD[0].pid() == 421
               && ((DstarP && DD[0].parents()[0].parents()[0].pid() == 10423)
                   || (DstarM && DD[1].parents()[0].parents()[0].pid() == -10423))) {
        _sigma[3]->fill(_edge[1]);
      }
      else if (DD[0].pid() == 411
               && (DD[0].parents()[0].pid() == 10413 || DD[1].parents()[0].pid() == -10413)) {
        _sigma[4]->fill(_edge[1]);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // dists
      normalize(_h, 1.0, false);
      // cross sections
      const vector<double> brs = {0.93, 1.0, 1.0, 0.677, 1.0};
      const double fact = crossSection() / sumOfWeights() / picobarn;
      ;
      for (size_t ih = 0; ih < 5; ++ih) {
        if (ih == 1) continue;
        scale(_sigma[ih], fact / brs[ih]);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    map<string, Histo1DPtr> _h;
    BinnedHistoPtr<string> _sigma[5];
    string _edge[2], _sqs = "";
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BESIII_2019_I1725786);

}