Rivet analyses

title: BELLE_2004_I653673

Double charmonium production in $e^+e^-$ collisions at $\sqrt{s}=10.6\,$GeV

Experiment: BELLE (KEKB)

Inspire ID: 653673

Status: VALIDATED NOHEPDATA

Authors: - Peter Richardson

References: - Phys.Rev.D 70 (2004) 071102

Beams: e+ e-

Beam energies: (5.3, 5.3)GeV

Run details: - e+e- > double charmonium

Double charmonium production at $\sqrt{s}=10.6\,$GeV. The cross sections and $\alpha$ parameters were taken from the tables in the paper and the corrected angular distributions read from the figures.

Source code:BELLE_2004_I653673.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief e+ e- -> double charmonium class BELLE_2004_I653673 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2004_I653673);


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

/// Book histograms and initialise projections before the run
void init() {
  // projections
  declare("FS", FinalState());
  declare("UFS",
          UnstableParticles(Cuts::pid == 443 || Cuts::pid == 100443 || Cuts::pid == 441
                            || Cuts::pid == 10441 || Cuts::pid == 100441));
  // histograms
  for (unsigned int ix = 0; ix < 3; ++ix) {
    for (unsigned int iy = 0; iy < 3; ++iy) {
      book(_p[ix][iy], "TMP/p_" + toString(ix + 1) + "_" + toString(iy + 1),
           refData<YODA::BinnedEstimate<string>>(3, 1 + ix, 1 + iy));
      if (ix == 2) continue;
      book(_h_sigma[ix][iy], 1 + ix, 1, 1 + 2 * iy);
      book(_h_angle[ix][iy], 4, 1 + ix, 1 + iy);
    }
  }
}

void findChildren(const Particle& p, map<long, int>& nRes, int& ncount, unsigned int& nCharged) {
  for (const Particle& child : p.children()) {
    if (child.children().empty()) {
      --nRes[child.pid()];
      --ncount;
      if (PID::isCharged(p.pid())) ++nCharged;
    }
    else {
      findChildren(child, nRes, ncount, nCharged);
    }
  }
}

double helicityAngle(const Particle& p) const {
  if (p.children().size() != 2) return 10.;
  if (p.children()[0].abspid() != PID::MUON || p.children()[0].pid() != -p.children()[1].pid()) {
    return 10.;
  }
  Vector3 axis = p.p3().unit();
  const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
  for (const Particle& child : p.children()) {
    if (child.pid() != PID::MUON) continue;
    return axis.dot(boost.transform(child.momentum()).p3().unit());
  }
  return 10;
}


/// Perform the per-event analysis
void analyze(const Event& event) {
  // final state particles
  const FinalState& fs = apply<FinalState>(event, "FS");
  map<long, int> nCount;
  int ntotal(0);
  for (const Particle& p : fs.particles()) {
    nCount[p.pid()] += 1;
    ++ntotal;
  }
  // loop over J/psi and psi(2S)
  const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
  bool matched = false;
  for (const Particle& p : ufs.particles(Cuts::pid == 443 || Cuts::pid == 100443)) {
    if (p.children().empty()) continue;
    map<long, int> nRes = nCount;
    int ncount = ntotal;
    unsigned int nCharged = 0;
    findChildren(p, nRes, ncount, nCharged);
    // eta_c, chi_c0, eta_c(2S)
    for (const Particle& p2 :
         ufs.particles(Cuts::pid == 441 || Cuts::pid == 10441 || Cuts::pid == 100441)) {
      map<long, int> nResB = nRes;
      int ncountB = ncount;
      unsigned int nChargedB = 0;
      findChildren(p2, nResB, ncountB, nChargedB);
      if (ncountB != 0) continue;
      matched = true;
      for (const auto& val : nResB) {
        if (val.second != 0) {
          matched = false;
          break;
        }
      }
      if (matched) {
        unsigned int ipsi = p.pid() / 100000;
        unsigned int ieta = p2.pid() / 10000;
        if (ieta > 1) ieta = 2;
        // fill the cross sections
        if ((ipsi == 0 && nChargedB > 2) || (ipsi == 1 && nChargedB > 0)) {
          _h_sigma[ipsi][ieta]->fill(_ecms);
        }
        if (ipsi > 0) break;
        // angular dists for J/psi only
        // production
        const double cProd = p.p3().z() / p.p3().mod();
        _h_angle[0][ieta]->fill(abs(cProd));
        _p[0][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cProd)));
        _p[2][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cProd)));
        // helicity angle
        const double cHel = helicityAngle(p);
        if (cHel > 1.) break;
        _h_angle[1][ieta]->fill(abs(cHel));
        _p[1][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cHel)));
        _p[2][ieta]->fill(_ecms, -1.25 * (1. - 3. * sqr(cHel)));
        break;
      }
    }
    if (matched) break;
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  for (unsigned int ix = 0; ix < 2; ++ix) {
    scale(_h_sigma[ix], crossSection() / sumOfWeights() / femtobarn);
    normalize(_h_angle[ix]);
  }
  // extract the alpha parameters
  for (unsigned int ix = 0; ix < 3; ++ix) {
    for (unsigned int iy = 0; iy < 3; ++iy) {
      const double val = _p[ix][iy]->bin(1).mean(2);
      const double err = _p[ix][iy]->bin(1).stdErr(2);
      BinnedEstimatePtr<string> tmp;
      book(tmp, 3, 1 + ix, 1 + iy);
      const double alpha = 3. * val / (1 - val);
      const double error = 3. / sqr(1. - val) * err;
      tmp->bin(1).set(alpha, error);
    }
  }
}

/// @}


/// @name Histograms
/// @{
BinnedHistoPtr<string> _h_sigma[2][3];
Histo1DPtr _h_angle[2][3];
BinnedProfilePtr<string> _p[3][3];
string _ecms = "10.6";
/// @}

};

RIVET_DECLARE_PLUGIN(BELLE_2004_I653673);

} ```