Rivet analyses

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

namespace Rivet {


  /// @brief Single and di-hadron spectra
  class BELLE_2020_I1777678 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2020_I1777678);


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

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      declare(ChargedFinalState(Cuts::abspid==211 or Cuts::abspid==321 or Cuts::abspid==2212),"CFS");
      // projections
      FinalState fs;
      declare(fs,"FS");
      declare(Thrust(fs),"Thrust");
      // single particle hists
      vector<int> pdg={211,321,2212};
      for (size_t ix=0; ix<3; ++ix) {
        book(_s_all   [pdg[ix]], 1, ix+1, 1);
        book(_s_strong[pdg[ix]], 1, ix+1, 2);
      }
      // dihadron histograms
      const vector<double> edges{0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,
                                 0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00};
      size_t i0=1;
      for (size_t defn=0; defn<3; ++defn) {
        for (size_t hemi=0; hemi<3; ++hemi) {
          for (size_t ip=0; ip<6; ++ip) {
            ++i0;
            size_t ymax=16;
            if(i0==7 || i0==19)      ymax=15;
            else if(i0>=8 && i0<=12) ymax=14;
            else if(i0==13)          ymax=13;
            else if(i0==26)          ymax=10;
            else if(i0==27||i0==30)  ymax= 9;
            else if(i0==28)          ymax= 7;
            else if(i0==29)          ymax= 8;
            else if(i0==31||i0==44)  ymax= 6;
            else if(i0==45||i0==48)  ymax= 5;
            else if(i0==46||i0==47)  ymax= 4;
            else if(i0==49        )  ymax= 3;
            book(_d_all[ip][defn][hemi], edges);
            book(_d_strong[ip][defn][hemi], edges);
            for (size_t iy=1; iy < _d_all[ip][defn][hemi]->numBins()+1; ++iy) {
              if (iy <= ymax) {
                book(_d_all[ip][defn][hemi]->bin(iy), i0, 1, iy);
                book(_d_strong[ip][defn][hemi]->bin(iy), i0, 2, iy);
              }
              else {
                _d_all[ip][defn][hemi]->maskBin(iy);
                _d_strong[ip][defn][hemi]->maskBin(iy);
              }

            }
          }
        }
      }
    }

    bool isWeak(const Particle & p) {
      bool weak = false;
      if(p.parents().empty()) return weak;
      Particle parent = p.parents()[0];
      while (!parent.parents().empty()) {
        if(parent.abspid()==411  || parent.abspid()==421  || parent.abspid()==431  ||
           parent.abspid()==4122 || parent.abspid()==4232 || parent.abspid()==4132 ||
           parent.abspid()==4332) {
          weak=true;
          break;
        }
        parent = parent.parents()[0];
      }
      return weak;
    }

    void fillHistos(int ip,bool strong,bool same,bool opp,
                    const Particle & p1, const Particle & p2) {
      for (size_t def=0; def<3; ++def) {
        double z1 = 0., z2 = 0.;
        if(def==0) {
          z1 = 2.*p1.momentum().t()/sqrtS();
          z2 = 2.*p2.momentum().t()/sqrtS();
        }
        else if(def==1) {
          z1 = 2.*p1.momentum().t()/sqrtS();
          z2 = (p1.momentum()*p2.momentum())/p1.momentum().t()/sqrtS();
        }
        else if(def==2) {
          double p1p2 = p1.momentum()*p2.momentum();
          double p1q = p1.momentum().t()*sqrtS();
          double p2q = p2.momentum().t()*sqrtS();
          z1 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p2q-p2.mass2()*p1q/p1p2);
          z2 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p1q-p1.mass2()*p2q/p1p2);
        }
        _d_all[ip][def][0]->fill(z1,z2,0.5);
        if (strong) _d_strong[ip][def][0]->fill(z1,z2,0.5);
        if (same) {
          _d_all[ip][def][1]->fill(z1,z2,0.5);
          if (strong) _d_strong[ip][def][1]->fill(z1,z2,0.5);
        }
        if (opp) {
          _d_all[ip][def][2]->fill(z1,z2,0.5);
          if (strong) _d_strong[ip][def][2]->fill(z1,z2,0.5);
        }
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // apply projection
      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
      // fill single particle histos
      for (const Particle& p : cfs.particles()) {
        const double z = 2.*p.momentum().t()/sqrtS();
        _s_all[p.abspid()]->fill(z);
        if (!isWeak(p))  _s_strong[p.abspid()]->fill(z);
      }
      // get thrust
      const Thrust thrust = apply<Thrust>(event,"Thrust");
      ThreeVector axis = thrust.thrustAxis();
      Particles piK = cfs.particles(Cuts::abspid==PID::KPLUS or Cuts::abspid==PID::PIPLUS);
      for (size_t ix=0; ix<piK.size(); ++ix) {
        double dot1 = axis.dot(piK[ix].momentum().p3());
        bool weak1 = isWeak(piK[ix]);
        for (size_t iy=0; iy<piK.size(); ++iy) {
          if (ix==iy) continue;
          double dot2 = axis.dot(piK[iy].momentum().p3());
          bool weak2 = isWeak(piK[iy]);
          bool strong = !weak1 && !weak2;
          bool same = thrust.thrust()>0.8 && dot1*dot2>0.;
          bool opp  = thrust.thrust()>0.8 && dot1*dot2<0.;
          unsigned int ip=0;
          if (piK[ix].pid()==PID::PIPLUS) {
            if (piK[iy].pid()==PID::PIPLUS)       ip=1;
            else if (piK[iy].pid()==PID::PIMINUS) ip=0;
            else if (piK[iy].pid()==PID::KPLUS  ) ip=3;
            else if (piK[iy].pid()==PID::KMINUS)  ip=2;
          }
          else if (piK[ix].pid()==PID::PIMINUS) {
            if (piK[iy].pid()==PID::PIPLUS)       ip=0;
            else if (piK[iy].pid()==PID::PIMINUS) ip=1;
            else if (piK[iy].pid()==PID::KPLUS  ) ip=2;
            else if (piK[iy].pid()==PID::KMINUS ) ip=3;
          }
          else if(piK[ix].pid()==PID::KPLUS) {
            if(piK[iy].pid()==PID::PIPLUS)       ip=3;
            else if(piK[iy].pid()==PID::PIMINUS) ip=2;
            else if(piK[iy].pid()==PID::KPLUS)   ip=5;
            else if(piK[iy].pid()==PID::KMINUS)  ip=4;
          }
          else if(piK[ix].pid()==PID::KMINUS) {
            if(piK[iy].pid()==PID::PIPLUS)       ip=2;
            else if(piK[iy].pid()==PID::PIMINUS) ip=3;
            else if(piK[iy].pid()==PID::KPLUS)   ip=4;
            else if(piK[iy].pid()==PID::KMINUS)  ip=5;
          }
          fillHistos(ip,strong,same,opp,piK[ix],piK[iy]);
        }
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {

      const double sf = crossSection()/femtobarn/sumOfWeights();
      scale(_s_all, sf);
      scale(_s_strong, sf);
      for (size_t ix=0; ix<6; ++ix) {
        for (size_t iy=0; iy<3; ++iy) {
          scale(_d_all[ix][iy], sf);
          divByGroupWidth(_d_all[ix][iy]);
          scale(_d_strong[ix][iy], sf);
          divByGroupWidth(_d_strong[ix][iy]);
        }
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    map<int,Histo1DPtr> _s_all,_s_strong;
    Histo1DGroupPtr _d_all[6][3][3], _d_strong[6][3][3];
    /// @}


  };


  RIVET_DECLARE_PLUGIN(BELLE_2020_I1777678);

}