1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief D_s+ -> K+K- mu+ nu_mu
 10  class BESIII_2023_I2674768 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2674768);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      UnstableParticles ufs = UnstableParticles(Cuts::pid==431);
 25      declare(ufs, "UFS");
 26      DecayedParticles DS(ufs);
 27      DS.addStable(PID::PI0);
 28      DS.addStable(PID::K0S);
 29      DS.addStable(PID::ETA);
 30      DS.addStable(PID::ETAPRIME);
 31      declare(DS, "DS");
 32      // Book histograms
 33      for (unsigned int ix=0; ix<5; ++ix) {
 34        book(_h[ix], 1, 1, 1+ix);
 35      }
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41      DecayedParticles DS = apply<DecayedParticles>(event, "DS");
 42      // loop over particles
 43      for (unsigned int ix=0; ix<DS.decaying().size(); ++ix) {
 44        if (!DS.modeMatches(ix, 4, mode)) continue;
 45        const Particle& Kp = DS.decayProducts()[ix].at( 321)[0];
 46        const Particle& Km = DS.decayProducts()[ix].at(-321)[0];
 47        const Particle& ep = DS.decayProducts()[ix].at( -13)[0];
 48        const Particle& nue= DS.decayProducts()[ix].at(  14)[0];
 49        FourMomentum pPhi = Kp.mom()+Km.mom();
 50        _h[0]->fill(pPhi.mass());
 51        FourMomentum qq = DS.decaying()[ix].mom()-pPhi;
 52        _h[1]->fill(qq.mass2());
 53        // boost momenta to DS rest frame
 54        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(DS.decaying()[ix].mom().betaVec());
 55        FourMomentum pPHI = boost.transform(pPhi);
 56        Matrix3 ptoz(-pPHI.p3().unit(), Vector3(0,0,1));
 57        boost.preMult(ptoz);
 58        // the momenta in frane to W along z
 59        FourMomentum pD  = boost.transform(DS.decaying()[ix].mom());
 60        FourMomentum pKp = boost.transform(Kp .mom());
 61        FourMomentum pKm = boost.transform(Km .mom());
 62        FourMomentum pe  = boost.transform(ep .mom());
 63        FourMomentum pnu = boost.transform(nue.mom());
 64        pPhi = pKp+pKm;
 65        qq = pD-pPhi;
 66        LorentzTransform boostK = LorentzTransform::mkFrameTransformFromBeta(pPhi.betaVec());
 67        Vector3 axisK = boostK.transform(pKp).p3().unit();
 68        _h[2]->fill(axisK.dot(pPhi.p3().unit()));
 69        LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(  qq.betaVec());
 70        Vector3 axisE = boostW.transform(pe).p3().unit();
 71        _h[3]->fill(axisE.dot(qq.p3().unit()));
 72        axisK.setZ(0.);
 73        axisE.setZ(0.);
 74        const double chi = atan2(axisE.cross(axisK).dot(qq.p3().unit()), axisE.dot(axisK));
 75        _h[4]->fill(chi);
 76      }
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82      normalize(_h);
 83    }
 84
 85    /// @}
 86
 87
 88    /// @name Histograms
 89    /// @{
 90    Histo1DPtr _h[5];
 91    const map<PdgId,unsigned int> mode = { { 321,1}, {-321,1}, {-13,1}, { 14,1}};
 92    /// @}
 93
 94
 95  };
 96
 97
 98  RIVET_DECLARE_PLUGIN(BESIII_2023_I2674768);
 99
100}