Rivet analyses

// -*- 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);

}