Rivet analyses

title: BELLE_2002_I563840

Prompt Charmonium production at 10.6 GeV

Experiment: BELLE (KEKB)

Inspire ID: 563840

Status: VALIDATED

Authors: - Peter Richardson

References: - Phys.Rev.Lett. 88 (2002) 052001

Beams: e+ e-

Beam energies: (5.3, 5.3)GeV

Run details: - e+e- to hadrons, only the continuuum contributes not Upslion(4S) decays

Measurement of the spectra for prompt $J/\psi$ and $\psi(2S)$ production at 10.6 GeV by BELLE.

Source code:BELLE_2002_I563840.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief charmonium production at 10.6 GeV class BELLE_2002_I563840 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2002_I563840);


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

/// Book histograms and initialise projections before the run
void init() {
  // projections
  declare(UnstableParticles(), "UFS");
  // book the histograms
  // spectra
  book(_h_Jpsi, 3, 1, 1);
  book(_h_feed, 3, 1, 2);
  book(_h_Psi2, 3, 2, 1);
  // cross sections
  book(_h_sig_JPsi_all, 1, 1, 1);
  book(_h_sig_Jpsi_high, 1, 2, 1);
  book(_h_sig_Jpsi_feed, 1, 2, 2);
  book(_h_sig_Psi2_high, 1, 2, 3);
  // angular distributions
  const vector<double> bins = {2., 2.6, 3.4, 4.9};
  book(_h_cThetaStar, bins);
  book(_h_cThetaH, bins);
  for (size_t ix = 1; ix < _h_cThetaH->numBins() + 1; ++ix) {
    if (ix <= 2) {
      const string suff = to_string(ix + 1);
      book(_h_cThetaStar->bin(ix), "/TMP/cThetaStar_" + suff, 5, -1.0, 1.0);
      book(_h_cThetaH->bin(ix), "/TMP/cThetaH_" + suff, 5, -1.0, 1.0);
    }
    else {
      book(_h_cThetaStar->bin(ix), 4, 1, 2);
      book(_h_cThetaH->bin(ix), 4, 2, 2);
    }
  }
  book(_h_cS_low, 4, 1, 1);
  book(_h_cS_high, "/TMP/cS_high", 5, -1.0, 1.0);
  book(_h_cH_low, 4, 2, 1);
  book(_h_cH_high, "/TMP/cH_high", 5, -1.0, 1.0);
}

void findLeptons(const Particle& mother, unsigned int& nstable, Particles& lp, Particles& lm) {
  for (const Particle& p : mother.children()) {
    int id = p.pid();
    if (id == 11 || id == 13) {
      lm.push_back(p);
      ++nstable;
    }
    else if (id == -11 || id == -13) {
      lp.push_back(p);
      ++nstable;
    }
    else if (!p.children().empty()) {
      findLeptons(p, nstable, lp, lm);
    }
    else {
      ++nstable;
    }
  }
}

/// Perform the per-event analysis
void analyze(const Event& event) {
  for (const Particle& p :
       apply<UnstableParticles>("UFS", event).particles(Cuts::pid == 443 or Cuts::pid == 100443)) {
    LorentzTransform boost = cmsTransform(beams());
    // check if prompt (i.e. not from B decay)
    if (p.fromBottom()) continue;
    bool feedDown = false;
    FourMomentum mom = boost.transform(p.momentum());
    double pStar = mom.p3().mod();
    if (p.pid() == 443) {
      Particle parent = p;
      while (!parent.parents().empty()) {
        parent = parent.parents()[0];
        if (p.pid() == parent.pid()) continue;
        if ((parent.abspid() % 1000) / 10 == 44) {
          feedDown = true;
          break;
        }
      }
      _h_Jpsi->fill(pStar);
      _h_sig_JPsi_all->fill("10.6"s);
      if (pStar > 2.) {
        _h_sig_Jpsi_high->fill("10.6"s);
        if (feedDown) {
          _h_feed->fill(pStar);
          _h_sig_Jpsi_feed->fill("10.6"s);
        }
        double cThetaS = cos(mom.p3().polarAngle());
        _h_cThetaStar->fill(pStar, cThetaS);
        if (pStar < 3.4)
          _h_cS_low->fill(cThetaS);
        else
          _h_cS_high->fill(cThetaS);
        // leptons from J/psi decay
        unsigned int nstable = 0;
        Particles lp, lm;
        findLeptons(p, nstable, lp, lm);
        if (nstable == 2 && lp.size() == 1 && lm.size() == 1) {
          FourMomentum pl = boost.transform(lp[0].momentum());
          LorentzTransform b2 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
          pl = b2.transform(pl);
          double cThetaH = pl.p3().unit().dot(p.p3().unit());
          _h_cThetaH->fill(pStar, cThetaH);
          if (pStar < 3.4)
            _h_cH_low->fill(cThetaH);
          else
            _h_cH_high->fill(cThetaH);
        }
      }
    }
    else {
      _h_Psi2->fill(pStar);
      if (pStar > 2.) _h_sig_Psi2_high->fill("10.6"s);
    }
  }
}

pair<double, pair<double, double>> calcAlpha(Histo1DPtr hist) {
  if (hist->numEntries() == 0.) return make_pair(0., make_pair(0., 0.));
  double sum1(0.), sum2(0.), sum3(0.), sum4(0.), sum5(0.);
  for (const auto& bin : hist->bins()) {
    double Oi = bin.sumW();
    if (Oi == 0.) continue;
    double a = 1.5 * (bin.xMax() - bin.xMin());
    double b = 0.5 * (pow(bin.xMax(), 3) - pow(bin.xMin(), 3));
    double Ei = bin.errW();
    sum1 += a * Oi / sqr(Ei);
    sum2 += b * Oi / sqr(Ei);
    sum3 += sqr(a) / sqr(Ei);
    sum4 += sqr(b) / sqr(Ei);
    sum5 += a * b / sqr(Ei);
  }
  // calculate alpha
  double alpha = (-3 * sum1 + 9 * sum2 + sum3 - 3 * sum5) / (sum1 - 3 * sum2 + 3 * sum4 - sum5);
  // and error
  double cc = -pow((sum3 + 9 * sum4 - 6 * sum5), 3);
  double bb = -2 * sqr(sum3 + 9 * sum4 - 6 * sum5) * (sum1 - 3 * sum2 + 3 * sum4 - sum5);
  double aa = sqr(sum1 - 3 * sum2 + 3 * sum4 - sum5)
      * (-sum3 - 9 * sum4 + sqr(sum1 - 3 * sum2 + 3 * sum4 - sum5) + 6 * sum5);
  double dis = sqr(bb) - 4. * aa * cc;
  if (dis > 0.) {
    dis = sqrt(dis);
    return make_pair(alpha, make_pair(0.5 * (-bb + dis) / aa, -0.5 * (-bb - dis) / aa));
  }
  else {
    return make_pair(alpha, make_pair(0., 0.));
  }
}

/// Normalise histograms etc., after the run
void finalize() {
  // spectra
  normalize(_h_Jpsi, 1., false);
  normalize(_h_feed, 1., false);
  normalize(_h_Psi2, 1., false);
  // cross sections
  double fact = 1. / sumOfWeights() * crossSection() / picobarn;
  scale(_h_sig_JPsi_all, fact);
  scale(_h_sig_Jpsi_high, fact);
  scale(_h_sig_Jpsi_feed, fact);
  scale(_h_sig_Psi2_high, fact);
  // angular dists and parameters from them
  vector<double> bins = {2., 2.6, 3.4, 4.9};
  Estimate1DPtr _h_A;
  book(_h_A, 2, 1, 1);
  Estimate1DPtr _h_alpha;
  book(_h_alpha, 2, 1, 2);
  for (size_t ix = 1; ix < _h_cThetaH->numBins() + 1; ++ix) {
    normalize(_h_cThetaStar->bin(ix));
    pair<double, pair<double, double>> alpha = calcAlpha(_h_cThetaStar->bin(ix));
    _h_A->bin(ix).set(alpha.first, alpha.second);
    normalize(_h_cThetaH->bin(ix));
    alpha = calcAlpha(_h_cThetaH->bin(ix));
    _h_alpha->bin(ix).set(alpha.first, alpha.second);
  }
  normalize(_h_cS_low);
  pair<double, pair<double, double>> alpha = calcAlpha(_h_cS_low);
  book(_h_A, 2, 2, 1);
  _h_A->bin(1).set(alpha.first, alpha.second);
  normalize(_h_cS_high);
  alpha = calcAlpha(_h_cS_high);
  book(_h_A, 2, 3, 1);
  _h_A->bin(1).set(alpha.first, alpha.second);
  normalize(_h_cH_low);
  alpha = calcAlpha(_h_cH_low);
  book(_h_alpha, 2, 2, 2);
  _h_alpha->bin(1).set(alpha.first, alpha.second);
  normalize(_h_cH_high);
  alpha = calcAlpha(_h_cH_high);
  book(_h_alpha, 2, 3, 2);
  _h_alpha->bin(1).set(alpha.first, alpha.second);
}

/// @}


/// @name Histograms
/// @{
Histo1DPtr _h_Jpsi, _h_Psi2, _h_feed;
BinnedHistoPtr<string> _h_sig_JPsi_all, _h_sig_Jpsi_high, _h_sig_Jpsi_feed, _h_sig_Psi2_high;
Histo1DPtr _h_cS_low, _h_cS_high, _h_cH_low, _h_cH_high;
Histo1DGroupPtr _h_cThetaStar, _h_cThetaH;


/// @}

};

RIVET_DECLARE_PLUGIN(BELLE_2002_I563840);

} ```