Rivet analyses

title: BELLE_2019_I1752523

Azimuthal asymmetries in bact-to-back $\pi^\pm$ and $(\pi^0,\eta,\pi^\pm)$ production at $\sqrt{s}=10.58$GeV

Experiment: BELLE (KEKB)

Inspire ID: 1752523

Status: VALIDATED NOHEPDATA

Authors: - Peter Richardson

References: - Phys.Rev.D 100 (2019) 9, 092008

Beams: e+ e-

Beam energies: (5.3, 5.3)GeV

Run details: - e+e- to hadrons

Azimuthal asymmetries in bact-to-back $\pi^\pm$ and $(\pi^0,\eta,\pi^\pm)$ production at $\sqrt{s}=10.58$GeV

Source code:BELLE_2019_I1752523.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/Thrust.hh"

include "Rivet/Projections/UnstableParticles.hh"

include "Rivet/Tools/Random.hh"

namespace Rivet {

/// @brief azimuthal asymmetries in pi+/- (pi0, eta, pi+-) production class BELLE_2019_I1752523 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2019_I1752523);


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

/// Book histograms and initialise projections before the run
void init() {
  // projections
  const FinalState fs;
  declare(fs, "FS");
  declare(UnstableParticles(Cuts::pid == 111 || Cuts::pid == 221), "UFS");
  declare(Thrust(fs), "Thrust");
  declare(Beam(), "Beams");
  // histograms
  unsigned int nbin = 20;
  // charged
  for (unsigned int ix = 0; ix < 3; ++ix) {
    // charged in pT bins
    for (unsigned int iy = 0; iy < 4; ++iy) {
      book(_h_charged_pt[ix][iy], "TMP/h_charged_pt_" + toString(ix + 1) + "_" + toString(iy + 1), nbin,
           0., M_PI);
      for (unsigned int iz = 0; iz < 4; ++iz) {
        book(_b_charged_pt[ix][iy][iz],
             "TMP/b_charged_pt_" + toString(ix) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
      }
    }
    // charged in z bins
    for (unsigned int iy = 0; iy < 6; ++iy) {
      book(_h_charged_z[ix][iy], "TMP/h_charged_z_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0.,
           M_PI);
      if (iy == 5) continue;
      for (unsigned int iz = 0; iz < 5; ++iz) {
        book(_b_charged_z[ix][iy][iz],
             "TMP/b_charged_z_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1),
             nbin, 0., M_PI);
      }
    }
  }
  for (unsigned int ix = 0; ix < 2; ++ix) {
    // pi0
    for (unsigned int iy = 0; iy < 4; ++iy) {
      book(_h_pi0_pt[ix][iy], "TMP/h_pi0_pt_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0.,
           M_PI);
      book(_h_eta_pt[ix][iy], "TMP/h_eta_pt_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0.,
           M_PI);
      book(_h_pi0B_pt[ix][iy], "TMP/h_pi0B_pt_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0.,
           M_PI);
      for (unsigned int iz = 0; iz < 4; ++iz) {
        book(_b_pi0_pt[ix][iy][iz],
             "TMP/b_pi0_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
        book(_b_eta_pt[ix][iy][iz],
             "TMP/b_eta_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
        book(_b_pi0B_pt[ix][iy][iz],
             "TMP/b_pi0B_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
        book(_b_charged_z_pt[ix][iy][iz],
             "TMP/b_charged_z_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1),
             nbin, 0., M_PI);
        book(_b_pi0_z_pt[ix][iy][iz],
             "TMP/b_pi0_z_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
        if (iy > 2) continue;
        book(_b_eta_z_pt[ix][iy][iz],
             "TMP/b_eta_z_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
        book(_b_pi0B_z_pt[ix][iy][iz],
             "TMP/b_pi0B_z_pt_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1),
             nbin, 0., M_PI);
      }
    }
    // pi0 in z bins
    for (unsigned int iy = 0; iy < 6; ++iy) {
      book(_h_pi0_z[ix][iy], "TMP/h_pi0_z_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0., M_PI);
      if (iy == 5) continue;
      for (unsigned int iz = 0; iz < 5; ++iz) {
        book(_b_pi0_z[ix][iy][iz],
             "TMP/b_pi0_z_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
      }
    }
    // eta in z bins
    for (unsigned int iy = 0; iy < 5; ++iy) {
      book(_h_eta_z[ix][iy], "TMP/h_eta_z_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0., M_PI);
      book(_h_pi0B_z[ix][iy], "TMP/h_pi0B_z_" + toString(ix + 1) + "_" + toString(iy + 1), nbin, 0.,
           M_PI);
      if (iy > 2) continue;
      for (unsigned int iz = 0; iz < 3; ++iz) {
        book(_b_eta_z[ix][iy][iz],
             "TMP/b_eta_z_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
        book(_b_pi0B_z[ix][iy][iz],
             "TMP/b_pi0B_z_" + toString(ix + 1) + "_" + toString(iy + 1) + "_" + toString(iz + 1), nbin,
             0., M_PI);
      }
    }
  }
}

unsigned int iBin_z0(double z) {
  if (z < .3)
    return 0;
  else if (z < .4)
    return 1;
  else if (z < .5)
    return 2;
  else if (z < .6)
    return 3;
  else if (z < .7)
    return 4;
  else
    return 5;
}

unsigned int iBin_z1(double z) {
  if (z < .2)
    return 0;
  else if (z < .3)
    return 1;
  else if (z < .5)
    return 2;
  else if (z < .7)
    return 3;
  else
    return 4;
}

unsigned int iBin_pT1(double pT) {
  if (pT < .15)
    return 0;
  else if (pT < .3)
    return 1;
  else if (pT < .5)
    return 2;
  else
    return 3;
}

/// Perform the per-event analysis
void analyze(const Event& event) {
  static const double c03 = cos(0.3);
  // get the axis, direction of incoming electron
  const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
  Vector3 axis1;
  if (beams.first.pid() > 0) {
    axis1 = beams.second.mom().p3().unit();
  }
  else {
    axis1 = beams.first.mom().p3().unit();
  }
  // apply thrust cuts  T > 0.8
  Thrust thrust = apply<Thrust>(event, "Thrust");
  if (thrust.thrust() <= 0.8) vetoEvent;
  ThreeVector t_z = thrust.thrustAxis();
  ThreeVector t_x = (axis1 - t_z.dot(axis1) * t_z).unit();
  ThreeVector t_y = t_z.cross(t_x);
  double cThetaT = t_z.dot(axis1);
  static double cTheta[2] = {cos(1.34), cos(2.03)};
  if (cThetaT < cTheta[1] || cThetaT > cTheta[0]) vetoEvent;
  Vector3 cross = axis1.cross(t_z);
  Particles charged = apply<FinalState>(event, "FS").particles(Cuts::abspid == PID::PIPLUS);
  Particles neutral = apply<UnstableParticles>(event, "UFS").particles();
  // first charged particle
  for (unsigned int ix = 0; ix < charged.size(); ++ix) {
    // z and angle cut
    const double x1 = 2. * charged[ix].mom().t() / sqrtS();
    double dot1 = t_z.dot(charged[ix].p3().unit());
    if (abs(dot1) < c03 || x1 < 0.1) continue;
    // second charged particle
    for (unsigned int iy = ix + 1; iy < charged.size(); ++iy) {
      const double x2 = 2. * charged[iy].mom().t() / sqrtS();
      double dot2 = t_z.dot(charged[iy].p3().unit());
      if (abs(dot2) < c03 || x2 < 0.1 || dot1 * dot2 > 0) continue;
      // randomly order the particles
      Particle p1 = charged[ix], p2 = charged[iy];
      double z1(x1), z2(x2);
      if (rand01() < 0.5) {
        swap(p1, p2);
        swap(z1, z2);
      }
      const double pT1 = sqrt(sqr(dot(p1.mom().p3(), thrust.thrustMajorAxis()))
                              + sqr(dot(p1.mom().p3(), thrust.thrustMinorAxis())));
      const double pT2 = sqrt(sqr(dot(p2.mom().p3(), thrust.thrustMajorAxis()))
                              + sqr(dot(p2.mom().p3(), thrust.thrustMinorAxis())));
      if (pT1 > 3. || pT2 > 3.) continue;
      double phi12 = atan2(p1.p3().dot(t_y), p1.p3().dot(t_x))
          + atan2(p2.p3().dot(t_y), p2.p3().dot(t_x));
      if (phi12 > M_PI) phi12 -= 2 * M_PI;
      if (phi12 < -M_PI) phi12 += 2 * M_PI;
      if (phi12 < 0.) phi12 = -phi12;
      if (z1 > .2 && z2 > .2) {
        unsigned int ibin = iBin_z0(z1);
        unsigned int ipT1 = iBin_pT1(pT1);
        unsigned int ipT2 = iBin_pT1(pT2);
        unsigned int ibin1 = iBin_z1(z1) - 1;
        if (p1.pid() * p2.pid() > 0) {
          _h_charged_z[0][ibin]->fill(phi12);
          _h_charged_pt[0][ipT1]->fill(phi12);
          _b_charged_pt[0][ipT1][ipT2]->fill(phi12);
          _h_pi0_z[1][ibin]->fill(phi12);
          _h_pi0_pt[1][ipT1]->fill(phi12);
          _b_pi0_pt[1][ipT1][ipT2]->fill(phi12);
          _b_pi0_z_pt[1][ibin1][ipT1]->fill(phi12);
          if (z1 > .3 && z2 > .3) {
            _h_eta_z[1][ibin - 1]->fill(phi12);
            _h_pi0B_z[1][ibin - 1]->fill(phi12);
            _h_eta_pt[1][ipT1]->fill(phi12);
            _h_pi0B_pt[1][ipT1]->fill(phi12);
            _b_eta_pt[1][ipT1][ipT2]->fill(phi12);
            _b_pi0B_pt[1][ipT1][ipT2]->fill(phi12);
            _b_eta_z_pt[1][ibin1 - 1][ipT1]->fill(phi12);
            _b_pi0B_z_pt[1][ibin1 - 1][ipT1]->fill(phi12);
          }
        }
        else {
          _h_charged_z[1][ibin]->fill(phi12);
          _h_charged_pt[1][ipT1]->fill(phi12);
          _b_charged_pt[1][ipT1][ipT2]->fill(phi12);
          _b_charged_z_pt[0][ibin1][ipT1]->fill(phi12);
        }
        _h_charged_z[2][ibin]->fill(phi12);
        _h_charged_pt[2][ipT1]->fill(phi12);
        _b_charged_pt[2][ipT1][ipT2]->fill(phi12);
        _b_charged_z_pt[1][ibin1][ipT1]->fill(phi12);
      }
      unsigned int ibin1 = iBin_z1(z1), ibin2 = iBin_z1(z2);
      if (p1.pid() * p2.pid() > 0) {
        _b_charged_z[0][ibin1][ibin2]->fill(phi12);
        _b_pi0_z[1][ibin1][ibin2]->fill(phi12);
        if (z1 > .3 && z2 > .3) {
          _b_eta_z[1][ibin1 - 2][ibin2 - 2]->fill(phi12);
          _b_pi0B_z[1][ibin1 - 2][ibin2 - 2]->fill(phi12);
        }
      }
      else {
        _b_charged_z[1][ibin1][ibin2]->fill(phi12);
      }
      _b_charged_z[2][ibin1][ibin2]->fill(phi12);
    }
    // pi0 or eta
    for (unsigned int iy = 0; iy < neutral.size(); ++iy) {
      const double x2 = 2. * neutral[iy].mom().t() / sqrtS();
      double dot2 = t_z.dot(neutral[iy].p3().unit());
      if (abs(dot2) < c03 || x2 < 0.2 || dot1 * dot2 > 0) continue;
      // neutral particle first
      Particle p2 = charged[ix], p1 = neutral[iy];
      double z1(x2), z2(x1);
      const double pT1 = sqrt(sqr(dot(p1.mom().p3(), thrust.thrustMajorAxis()))
                              + sqr(dot(p1.mom().p3(), thrust.thrustMinorAxis())));
      const double pT2 = sqrt(sqr(dot(p2.mom().p3(), thrust.thrustMajorAxis()))
                              + sqr(dot(p2.mom().p3(), thrust.thrustMinorAxis())));
      if (pT1 > 3. || pT2 > 3.) continue;
      double phi12 = atan2(p1.p3().dot(t_y), p1.p3().dot(t_x))
          + atan2(p2.p3().dot(t_y), p2.p3().dot(t_x));
      if (phi12 > M_PI) phi12 -= 2 * M_PI;
      if (phi12 < -M_PI) phi12 += 2 * M_PI;
      if (phi12 < 0.) phi12 = -phi12;
      if (z1 > .2) {
        unsigned int ibin = iBin_z0(z1);
        unsigned int ipT1 = iBin_pT1(pT1);
        unsigned int ipT2 = iBin_pT1(pT2);
        unsigned int ibin1 = iBin_z1(z1) - 1;
        if (p1.pid() == 111) {
          _h_pi0_z[0][ibin]->fill(phi12);
          _h_pi0_pt[0][ipT1]->fill(phi12);
          _b_pi0_pt[0][ipT1][ipT2]->fill(phi12);
          _b_pi0_z_pt[0][ibin1][ipT1]->fill(phi12);
          if (z1 > .3 && z2 > .3) {
            _h_pi0B_z[0][ibin - 1]->fill(phi12);
            _h_pi0B_pt[0][ipT1]->fill(phi12);
            _b_pi0B_pt[0][ipT1][ipT2]->fill(phi12);
            _b_pi0B_z_pt[0][ibin1 - 1][ipT1]->fill(phi12);
          }
        }
        else if (z1 > .3 && z2 > .3) {
          _h_eta_z[0][ibin - 1]->fill(phi12);
          _h_eta_pt[0][ipT1]->fill(phi12);
          _b_eta_pt[0][ipT1][ipT2]->fill(phi12);
          _b_eta_z_pt[0][ibin1 - 1][ipT1]->fill(phi12);
        }
      }
      unsigned int ibin1 = iBin_z1(z1), ibin2 = iBin_z1(z2);
      if (p1.pid() == 111) {
        _b_pi0_z[0][ibin1][ibin2]->fill(phi12);
        if (z1 > .3 && z2 > .3) _b_pi0B_z[0][ibin1 - 2][ibin2 - 2]->fill(phi12);
      }
      else if (z1 > .3 && z2 > .3) {
        _b_eta_z[0][ibin1 - 2][ibin2 - 2]->fill(phi12);
      }
    }
  }
}

pair<double, double> calcAsymmetry(Estimate1DPtr hist, double fact = 1.) {
  double sum1(0.), sum2(0.);
  for (const auto& bin : hist->bins()) {
    double Oi = bin.val();
    if (Oi == 0. || std::isnan(Oi)) continue;
    double ai = 1.;
    double bi = (sin(fact * bin.xMax()) - sin(fact * bin.xMin())) / (bin.xMax() - bin.xMin()) / fact;
    double Ei = bin.errAvg();
    sum1 += sqr(bi / Ei);
    sum2 += bi / sqr(Ei) * (Oi - ai);
  }
  if (sum1 == 0.) return make_pair(0., 0.);
  return make_pair(sum2 / sum1 * 1e2, sqrt(1. / sum1) * 1e2);
}

/// Normalise histograms etc., after the run
void finalize() {
  // charged in pT bins
  Estimate1DPtr h_charged_UL, h_charged_UC;
  book(h_charged_UL, 1, 1, 1);
  book(h_charged_UC, 1, 1, 2);
  Estimate1DPtr h_pi0_UL, h_eta_UL, h_pi0B_UL;
  book(h_pi0_UL, 6, 1, 1);
  book(h_eta_UL, 7, 1, 1);
  book(h_pi0B_UL, 8, 1, 1);
  for (unsigned int ix = 0; ix < 4; ++ix) {
    // first the 1d dists
    for (unsigned int ii = 0; ii < 3; ++ii) {
      normalize(_h_charged_pt[ii][ix]);
      if (ii == 2) continue;
      normalize(_h_pi0_pt[ii][ix]);
      normalize(_h_eta_pt[ii][ix]);
      normalize(_h_pi0B_pt[ii][ix]);
    }
    Estimate1DPtr htemp;
    // UL
    // charged
    book(htemp, "TMP/R_charged_pt_UL_" + toString(ix), _h_charged_pt[1][ix]->xEdges());
    divide(_h_charged_pt[1][ix], _h_charged_pt[0][ix], htemp);
    pair<double, double> asym = calcAsymmetry(htemp);
    h_charged_UL->bin(ix + 1).setVal(asym.first);
    h_charged_UL->bin(ix + 1).setErr(asym.second);
    // UC
    book(htemp, "TMP/R_charged_pt_UC_" + toString(ix), _h_charged_pt[1][ix]->xEdges());
    divide(_h_charged_pt[1][ix], _h_charged_pt[2][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_charged_UC->bin(ix + 1).setVal(asym.first);
    h_charged_UC->bin(ix + 1).setErr(asym.second);
    // pi0
    book(htemp, "TMP/R_pi0_pt_UL_" + toString(ix), _h_pi0_pt[0][ix]->xEdges());
    divide(_h_pi0_pt[0][ix], _h_pi0_pt[1][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_pi0_UL->bin(ix + 1).setVal(asym.first);
    h_pi0_UL->bin(ix + 1).setErr(asym.second);
    // eta
    book(htemp, "TMP/R_eta_pt_UL_" + toString(ix), _h_eta_pt[0][ix]->xEdges());
    divide(_h_eta_pt[0][ix], _h_eta_pt[1][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_eta_UL->bin(ix + 1).setVal(asym.first);
    h_eta_UL->bin(ix + 1).setErr(asym.second);
    // pi0 z>.3
    book(htemp, "TMP/R_pi0B_pt_UL_" + toString(ix), _h_pi0B_pt[0][ix]->xEdges());
    divide(_h_pi0B_pt[0][ix], _h_pi0B_pt[1][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_pi0B_UL->bin(ix + 1).setVal(asym.first);
    h_pi0B_UL->bin(ix + 1).setErr(asym.second);
    // then 2nd dists
    Estimate1DPtr h_charged_UL2, h_charged_UC2;
    book(h_charged_UL2, 1, 2 + ix, 1);
    book(h_charged_UC2, 1, 2 + ix, 2);
    Estimate1DPtr h_pi0_UL2, h_eta_UL2, h_pi0B_UL2;
    book(h_pi0_UL2, 6, 2 + ix, 1);
    book(h_eta_UL2, 7, 2 + ix, 1);
    book(h_pi0B_UL2, 8, 2 + ix, 1);
    for (unsigned int iy = 0; iy < 4; ++iy) {
      for (unsigned int ii = 0; ii < 3; ++ii) {
        normalize(_b_charged_pt[ii][ix][iy]);
        if (ii == 2) continue;
        normalize(_b_pi0_pt[ii][ix][iy]);
        normalize(_b_eta_pt[ii][ix][iy]);
        normalize(_b_pi0B_pt[ii][ix][iy]);
      }
      // UL
      book(htemp, "TMP/R_charged_pt_UL_" + toString(ix) + "_" + toString(iy),
           _b_charged_pt[1][ix][iy]->xEdges());
      divide(_b_charged_pt[1][ix][iy], _b_charged_pt[0][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_charged_UL2->bin(iy + 1).setVal(asym.first);
      h_charged_UL2->bin(iy + 1).setErr(asym.second);
      // UC
      book(htemp, "TMP/R_charged_pt_UC_" + toString(ix) + "_" + toString(iy),
           _b_charged_pt[1][ix][iy]->xEdges());
      divide(_b_charged_pt[1][ix][iy], _b_charged_pt[2][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_charged_UC2->bin(iy + 1).setVal(asym.first);
      h_charged_UC2->bin(iy + 1).setErr(asym.second);
      // pi0
      book(htemp, "TMP/R_pi0_pt_UL_" + toString(ix) + "_" + toString(iy), _b_pi0_pt[0][ix][iy]->xEdges());
      divide(_b_pi0_pt[0][ix][iy], _b_pi0_pt[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_pi0_UL2->bin(iy + 1).setVal(asym.first);
      h_pi0_UL2->bin(iy + 1).setErr(asym.second);
      // eta
      book(htemp, "TMP/R_eta_pt_UL_" + toString(ix) + "_" + toString(iy), _b_eta_pt[0][ix][iy]->xEdges());
      divide(_b_eta_pt[0][ix][iy], _b_eta_pt[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_eta_UL2->bin(iy + 1).setVal(asym.first);
      h_eta_UL2->bin(iy + 1).setErr(asym.second);
      // pi0 z>0.3
      book(htemp, "TMP/R_pi0B_pt_UL_" + toString(ix) + "_" + toString(iy),
           _b_pi0B_pt[0][ix][iy]->xEdges());
      divide(_b_pi0B_pt[0][ix][iy], _b_pi0B_pt[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_pi0B_UL2->bin(iy + 1).setVal(asym.first);
      h_pi0B_UL2->bin(iy + 1).setErr(asym.second);
    }
  }
  // charged and pi0 in z bins
  book(h_charged_UL, 2, 1, 1);
  book(h_charged_UC, 2, 1, 2);
  book(h_pi0_UL, 3, 1, 1);
  book(h_eta_UL, 4, 1, 1);
  book(h_pi0B_UL, 5, 1, 1);
  for (unsigned int ix = 0; ix < 6; ++ix) {
    // first the 1d dists
    for (unsigned int ii = 0; ii < 3; ++ii) {
      normalize(_h_charged_z[ii][ix]);
      if (ii == 2) continue;
      normalize(_h_pi0_z[ii][ix]);
      if (ix == 5) continue;
      normalize(_h_eta_z[ii][ix]);
      normalize(_h_pi0B_z[ii][ix]);
    }
    Estimate1DPtr htemp;
    // UL
    // charged
    book(htemp, "TMP/R_charged_z_UL_" + toString(ix), _h_charged_z[1][ix]->xEdges());
    divide(_h_charged_z[1][ix], _h_charged_z[0][ix], htemp);
    pair<double, double> asym = calcAsymmetry(htemp);
    h_charged_UL->bin(ix + 1).setVal(asym.first);
    h_charged_UL->bin(ix + 1).setErr(asym.second);
    // UC
    book(htemp, "TMP/R_charged_z_UC_" + toString(ix), _h_charged_z[1][ix]->xEdges());
    divide(_h_charged_z[1][ix], _h_charged_z[2][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_charged_UC->bin(ix + 1).setVal(asym.first);
    h_charged_UC->bin(ix + 1).setErr(asym.second);
    // pi0
    book(htemp, "TMP/R_pi0_z_UL_" + toString(ix), _h_pi0_z[0][ix]->xEdges());
    divide(_h_pi0_z[0][ix], _h_pi0_z[1][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_pi0_UL->bin(ix + 1).setVal(asym.first);
    h_pi0_UL->bin(ix + 1).setErr(asym.second);
    if (ix == 5) continue;
    // eta
    book(htemp, "TMP/R_eta_z_UL_" + toString(ix), _h_eta_z[0][ix]->xEdges());
    divide(_h_eta_z[0][ix], _h_eta_z[1][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_eta_UL->bin(ix + 1).setVal(asym.first);
    h_eta_UL->bin(ix + 1).setErr(asym.second);
    // pi0 z>.3
    book(htemp, "TMP/R_pi0B_z_UL_" + toString(ix), _h_pi0B_z[0][ix]->xEdges());
    divide(_h_pi0B_z[0][ix], _h_pi0B_z[1][ix], htemp);
    asym = calcAsymmetry(htemp);
    h_pi0B_UL->bin(ix + 1).setVal(asym.first);
    h_pi0B_UL->bin(ix + 1).setErr(asym.second);
    // then 2nd dists
    Estimate1DPtr h_charged_UL2, h_charged_UC2;
    book(h_charged_UL2, 2, 2 + ix, 1);
    book(h_charged_UC2, 2, 2 + ix, 2);
    Estimate1DPtr h_pi0_UL2, h_eta_UL2, h_pi0B_UL2;
    book(h_pi0_UL2, 3, 2 + ix, 1);
    if (ix < 3) {
      book(h_eta_UL2, 4, 2 + ix, 1);
      book(h_pi0B_UL2, 5, 2 + ix, 1);
    }
    for (unsigned int iy = 0; iy < 5; ++iy) {
      for (unsigned int ii = 0; ii < 3; ++ii) {
        normalize(_b_charged_z[ii][ix][iy]);
        if (ii == 2) continue;
        normalize(_b_pi0_z[ii][ix][iy]);
        if (ix > 2 || iy > 2) continue;
        normalize(_b_eta_z[ii][ix][iy]);
        normalize(_b_pi0B_z[ii][ix][iy]);
      }
      // UL
      book(htemp, "TMP/R_charged_z_UL_" + toString(ix) + "_" + toString(iy),
           _b_charged_z[1][ix][iy]->xEdges());
      divide(_b_charged_z[1][ix][iy], _b_charged_z[0][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_charged_UL2->bin(iy + 1).setVal(asym.first);
      h_charged_UL2->bin(iy + 1).setErr(asym.second);
      // UC
      book(htemp, "TMP/R_charged_z_UC_" + toString(ix) + "_" + toString(iy),
           _b_charged_z[1][ix][iy]->xEdges());
      divide(_b_charged_z[1][ix][iy], _b_charged_z[2][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_charged_UC2->bin(iy + 1).setVal(asym.first);
      h_charged_UC2->bin(iy + 1).setErr(asym.second);
      // pi0
      book(htemp, "TMP/R_pi0_z_UL_" + toString(ix) + "_" + toString(iy), _b_pi0_z[0][ix][iy]->xEdges());
      divide(_b_pi0_z[0][ix][iy], _b_pi0_z[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_pi0_UL2->bin(iy + 1).setVal(asym.first);
      h_pi0_UL2->bin(iy + 1).setErr(asym.second);
      if (ix > 2 || iy > 2) continue;
      // eta
      book(htemp, "TMP/R_eta_z_UL_" + toString(ix) + "_" + toString(iy), _b_eta_z[0][ix][iy]->xEdges());
      divide(_b_eta_z[0][ix][iy], _b_eta_z[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_eta_UL2->bin(iy + 1).setVal(asym.first);
      h_eta_UL2->bin(iy + 1).setErr(asym.second);
      // pi0 z>0.3
      book(htemp, "TMP/R_pi0B_z_UL_" + toString(ix) + "_" + toString(iy), _b_pi0B_z[0][ix][iy]->xEdges());
      divide(_b_pi0B_z[0][ix][iy], _b_pi0B_z[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_pi0B_UL2->bin(iy + 1).setVal(asym.first);
      h_pi0B_UL2->bin(iy + 1).setErr(asym.second);
    }
  }
  // finally z and pT
  // only 2d dists
  for (unsigned int ix = 0; ix < 4; ++ix) {
    Estimate1DPtr h_charged_UC2;
    book(h_charged_UC2, 9, 1 + ix, 1);
    Estimate1DPtr h_pi0_UL2, h_eta_UL2, h_pi0B_UL2;
    book(h_pi0_UL2, 10, 1 + ix, 1);
    if (ix < 3) {
      book(h_eta_UL2, 11, 1 + ix, 1);
      book(h_pi0B_UL2, 12, 1 + ix, 1);
    }
    for (unsigned int iy = 0; iy < 4; ++iy) {
      for (unsigned int ii = 0; ii < 2; ++ii) {
        normalize(_b_charged_z_pt[ii][ix][iy]);
        normalize(_b_pi0_z_pt[ii][ix][iy]);
        if (ix > 2) continue;
        normalize(_b_eta_z_pt[ii][ix][iy]);
        normalize(_b_pi0B_z_pt[ii][ix][iy]);
      }
      pair<double, double> asym;
      Estimate1DPtr htemp;
      // UC (missing last point for first histo)
      if (ix != 0 || iy != 3) {
        book(htemp, "TMP/R_charged_z_pt_UC_" + toString(ix) + "_" + toString(iy),
             _b_charged_z_pt[0][ix][iy]->xEdges());
        divide(_b_charged_z_pt[0][ix][iy], _b_charged_z_pt[1][ix][iy], htemp);
        asym = calcAsymmetry(htemp);
        h_charged_UC2->bin(iy + 1).setVal(asym.first);
        h_charged_UC2->bin(iy + 1).setErr(asym.second);
        // pi0
        book(htemp, "TMP/R_pi0_z_pt_UL_" + toString(ix) + "_" + toString(iy),
             _b_pi0_z_pt[0][ix][iy]->xEdges());
        divide(_b_pi0_z_pt[0][ix][iy], _b_pi0_z_pt[1][ix][iy], htemp);
        asym = calcAsymmetry(htemp);
        h_pi0_UL2->bin(iy + 1).setVal(asym.first);
        h_pi0_UL2->bin(iy + 1).setErr(asym.second);
      }
      if (ix > 2) continue;
      // eta
      book(htemp, "TMP/R_eta_z_pt_UL_" + toString(ix) + "_" + toString(iy),
           _b_eta_z_pt[0][ix][iy]->xEdges());
      divide(_b_eta_z_pt[0][ix][iy], _b_eta_z_pt[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_eta_UL2->bin(iy + 1).setVal(asym.first);
      h_eta_UL2->bin(iy + 1).setErr(asym.second);
      // pi0 z>0.3
      book(htemp, "TMP/R_pi0B_z_pt_UL_" + toString(ix) + "_" + toString(iy),
           _b_pi0B_z_pt[0][ix][iy]->xEdges());
      divide(_b_pi0B_z_pt[0][ix][iy], _b_pi0B_z_pt[1][ix][iy], htemp);
      asym = calcAsymmetry(htemp);
      h_pi0B_UL2->bin(iy + 1).setVal(asym.first);
      h_pi0B_UL2->bin(iy + 1).setErr(asym.second);
    }
  }
}

/// @}


/// @name Histograms
/// @{
Histo1DPtr _h_charged_pt[3][4], _h_charged_z[3][6];
Histo1DPtr _b_charged_pt[3][4][4], _b_charged_z[3][5][5], _b_charged_z_pt[2][4][4];
Histo1DPtr _h_pi0_z[2][6], _h_eta_z[2][5], _h_pi0B_z[2][5], _h_pi0_pt[2][4], _h_eta_pt[2][4],
    _h_pi0B_pt[2][4];
Histo1DPtr _b_pi0_z[2][5][5], _b_eta_z[2][3][3], _b_pi0B_z[2][3][3], _b_pi0_pt[2][4][4];
Histo1DPtr _b_eta_pt[2][4][4], _b_pi0B_pt[2][4][4], _b_pi0_z_pt[2][4][4], _b_pi0B_z_pt[2][3][4],
    _b_eta_z_pt[2][3][4];
/// @}

};

RIVET_DECLARE_PLUGIN(BELLE_2019_I1752523);

} ```