1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief B -> D* semileptonic
 10  class BELLE_2023_I2624324 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2023_I2624324);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(UnstableParticles(Cuts::abspid==511 || Cuts::abspid==521), "UFS");
 24      // histograms
 25      for (unsigned int ix=0; ix<4; ++ix) {
 26        book(_h_aver[ix], 2, 1+ix, 1);
 27        for (unsigned int iy=0; iy<4; ++iy) {
 28          book(_h[ix][iy], 1, 1+ix, 1+iy);
 29        }
 30      }
 31    }
 32
 33    /// Perform the per-event analysis
 34    bool analyzeDecay(const Particle& mother, const vector<int>& ids) {
 35      // There is no point in looking for decays with less particles than to be analysed
 36      if (mother.children().size() == ids.size()) {
 37        bool decayfound = true;
 38        for (int id : ids) {
 39          if (!contains(mother, id)) decayfound = false;
 40        }
 41        return decayfound;
 42      }
 43      return false;
 44    }
 45
 46    bool contains(const Particle& mother, int id) {
 47      return any(mother.children(), HasPID(id));
 48    }
 49
 50    double recoilW(const Particle& mother) {
 51      FourMomentum lepton, neutrino, meson, q;
 52      for(const Particle& c : mother.children()) {
 53        if (c.isNeutrino()) neutrino=c.mom();
 54        if (c.isLepton() && !c.isNeutrino()) lepton =c.mom();
 55        if (c.isHadron()) meson=c.mom();
 56      }
 57      q = lepton + neutrino; //no hadron before
 58      double mb2= mother.mom()*mother.mom();
 59      double mD2 = meson*meson;
 60      return (mb2 + mD2 - q*q )/ (2. * sqrt(mb2) * sqrt(mD2) );
 61    }
 62
 63    /// Perform the per-event analysis
 64    void analyze(const Event& event) {
 65      FourMomentum pl, pnu, pB, pD, pDs, ppi;
 66      // Iterate of B mesons
 67      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 68        pB = p.mom();
 69        // Find semileptonic decays
 70        int sign = p.pid()/p.abspid();
 71        int iDStar = sign*(p.abspid()==511 ? -413 : -423);
 72        int iloc=-1;
 73        if (analyzeDecay(p, {iDStar,12*sign,-11*sign})) iloc = 0;
 74        else if(analyzeDecay(p, {iDStar,14*sign,-13*sign})) iloc = 1;
 75        else continue;
 76        if (p.abspid()==521) iloc+=2;
 77        double w = recoilW(p);
 78        _h[0][iloc]->fill(w);
 79        _h_aver[0] ->fill(w);
 80        // Get the necessary momenta for the angles
 81        bool foundDdecay=false;
 82        for (const Particle& c : p.children()) {
 83          if (c.abspid()==413 || c.abspid()==423) {
 84            if ((c.pid() == -413  && (analyzeDecay(c, {PID::PIMINUS, PID::D0BAR}) || analyzeDecay(c, {PID::PI0, PID::DMINUS})) ) ||
 85                (c.pid() ==  413  && (analyzeDecay(c, {PID::PIPLUS , PID::D0   }) || analyzeDecay(c, {PID::PI0, PID::DPLUS })) ) ||
 86                (c.pid() == -423  && analyzeDecay(c, {PID::PI0, PID::D0BAR }) ) ||
 87                (c.pid() ==  423  && analyzeDecay(c, {PID::PI0, PID::D0    }) )) {
 88              foundDdecay=true;
 89              pDs = c.mom();
 90              for (const Particle & dc : c.children()) {
 91                if (dc.hasCharm()) pD = dc.mom();
 92                else ppi = dc.mom();
 93              }
 94            }
 95          }
 96          else if (c.abspid() ==  11 || c.abspid() ==  13) pl  = c.mom();
 97          else if (c.abspid() ==  12 || c.abspid() ==  14) pnu = c.mom();
 98        }
 99        // This is the angle analysis
100        if (!foundDdecay) continue;
101        // First boost all relevant momenta into the B-rest frame
102        const LorentzTransform B_boost = LorentzTransform::mkFrameTransformFromBeta(pB.betaVec());
103        // Momenta in B rest frame:
104        FourMomentum lv_brest_Dstar = B_boost.transform(pDs);
105        FourMomentum lv_brest_w     = B_boost.transform(pB - pDs);
106        FourMomentum lv_brest_D     = B_boost.transform(pD);
107        FourMomentum lv_brest_lep   = B_boost.transform(pl);
108
109        const LorentzTransform Ds_boost = LorentzTransform::mkFrameTransformFromBeta(lv_brest_Dstar.betaVec());
110        FourMomentum lv_Dstarrest_D     = Ds_boost.transform(lv_brest_D);
111        const LorentzTransform W_boost  = LorentzTransform::mkFrameTransformFromBeta(lv_brest_w.betaVec());
112        FourMomentum lv_wrest_lep       = W_boost.transform(lv_brest_lep);
113
114        double cos_thetaV = cos(lv_brest_Dstar.p3().angle(lv_Dstarrest_D.p3()));
115        _h[2][iloc]->fill(cos_thetaV);
116        _h_aver[2] ->fill(cos_thetaV);
117
118        double cos_thetaL = cos(lv_brest_w.p3().angle(lv_wrest_lep.p3()));
119        _h[1][iloc]->fill(cos_thetaL);
120        _h_aver[1] ->fill(cos_thetaL);
121
122        Vector3 LTrans = lv_wrest_lep.p3()   - cos_thetaL*lv_wrest_lep.p3().perp()*lv_brest_w.p3().unit();
123        Vector3 VTrans = lv_Dstarrest_D.p3() - cos_thetaV*lv_Dstarrest_D.p3().perp()*lv_brest_Dstar.p3().unit();
124        float chi = atan2(LTrans.cross(VTrans).dot(lv_brest_w.p3().unit()), LTrans.dot(VTrans));
125        if (chi<0.) chi+=2.*M_PI;
126      	_h[3][iloc]->fill(chi);
127      	_h_aver[3] ->fill(chi);
128      }
129    }
130
131
132    /// Normalise histograms etc., after the run
133    void finalize() {
134      normalize(_h     , 1.0);
135      normalize(_h_aver, 1.0);
136    }
137
138    /// @}
139
140
141    /// @name Histograms
142    /// @{
143    Histo1DPtr _h[4][4];
144    Histo1DPtr _h_aver[4];
145    /// @}
146
147
148  };
149
150
151  RIVET_DECLARE_PLUGIN(BELLE_2023_I2624324);
152
153}