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 Bbar0 -> D*+ semileptonic
 10  class BELLE_2017_I1512299 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2017_I1512299);
 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      declare(UnstableParticles(), "UFS");
 25
 26      // Book histograms
 27      book(_h_w     , "TMP/h_w"     , refData(1, 1, 1));
 28      book(_h_costhv, "TMP/h_costhv", refData(2, 1, 1));
 29      book(_h_costhl, "TMP/h_costhl", refData(3, 1, 1));
 30      book(_h_chi   , "TMP/h_chi"   , refData(4, 1, 1));
 31      book(_nB,"TMP/nB");
 32    }
 33
 34
 35    /// Perform the per-event analysis
 36    bool analyzeDecay(Particle mother, vector<int> ids) {
 37      // There is no point in looking for decays with less particles than to be analysed
 38      if (mother.children().size() == ids.size()) {
 39        bool decayfound = true;
 40        for (int id : ids) {
 41          if (!contains(mother, id)) decayfound = false;
 42        }
 43        return decayfound;
 44      }
 45      return false;
 46    }
 47
 48    bool contains(Particle& mother, int id) {
 49      return any(mother.children(), HasPID(id));
 50    }
 51
 52
 53    double recoilW(const Particle& mother) {
 54      FourMomentum lepton, neutrino, meson, q;
 55      for(const Particle& c : mother.children()) {
 56        if (c.isNeutrino()) neutrino=c.mom();
 57        if (c.isLepton() && !c.isNeutrino()) lepton =c.mom();
 58        if (c.isHadron()) meson=c.mom();
 59      }
 60      q = lepton + neutrino; //no hadron before
 61      double mb2= mother.mom()*mother.mom();
 62      double mD2 = meson*meson;
 63      return (mb2 + mD2 - q*q )/ (2. * sqrt(mb2) * sqrt(mD2) );
 64    }
 65
 66    /// Perform the per-event analysis
 67    void analyze(const Event& event) {
 68      FourMomentum pl, pnu, pB, pD, pDs, ppi;
 69      // Iterate of B0bar mesons
 70      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==-511)) {
 71        if(p.children().size()==1 && p.children()[0].abspid()==511) continue;
 72        _nB->fill();
 73        pB = p.momentum();
 74        // Find semileptonic decays
 75        if (analyzeDecay(p, {PID::DSTARPLUS,-12,11}) || analyzeDecay(p, {PID::DSTARPLUS,-14,13}) ) {
 76          _h_w->fill(recoilW(p));
 77          // Get the necessary momenta for the angles
 78          bool foundDdecay=false;
 79          for (const Particle & c : p.children()) {
 80            if ( (c.pid() == PID::DSTARPLUS)  && (analyzeDecay(c, {PID::PIPLUS, PID::D0}) || analyzeDecay(c, {PID::PI0, PID::DPLUS})) ) {
 81              foundDdecay=true;
 82              pDs = c.momentum();
 83              for (const Particle & dc : c.children()) {
 84                if (dc.hasCharm()) pD = dc.momentum();
 85                else ppi = dc.momentum();
 86              }
 87            }
 88            if (c.pid() ==  11 || c.pid() ==  13) pl  = c.momentum();
 89            if (c.pid() == -12 || c.pid() == -14) pnu = c.momentum();
 90          }
 91          // This is the angle analysis
 92          if (foundDdecay) {
 93
 94            // First boost all relevant momenta into the B-rest frame
 95            const LorentzTransform B_boost = LorentzTransform::mkFrameTransformFromBeta(pB.betaVec());
 96            // Momenta in B rest frame
 97            FourMomentum lv_brest_Dstar = B_boost.transform(pDs);
 98            FourMomentum lv_brest_w     = B_boost.transform(pB - pDs);
 99            FourMomentum lv_brest_D     = B_boost.transform(pD);
100            FourMomentum lv_brest_lep   = B_boost.transform(pl);
101            
102            const LorentzTransform Ds_boost = LorentzTransform::mkFrameTransformFromBeta(lv_brest_Dstar.betaVec());
103            FourMomentum lv_Dstarrest_D     = Ds_boost.transform(lv_brest_D);
104            const LorentzTransform W_boost  = LorentzTransform::mkFrameTransformFromBeta(lv_brest_w.betaVec());
105            FourMomentum lv_wrest_lep       = W_boost.transform(lv_brest_lep);
106
107            double cos_thetaV = cos(lv_brest_Dstar.p3().angle(lv_Dstarrest_D.p3()));
108            _h_costhv->fill(cos_thetaV);
109
110            double cos_thetaL = cos(lv_brest_w.p3().angle(lv_wrest_lep.p3()));
111            _h_costhl->fill(cos_thetaL);
112
113            Vector3 LTrans = lv_wrest_lep.p3()   - cos_thetaL*lv_wrest_lep.p3().perp()*lv_brest_w.p3().unit();
114            Vector3 VTrans = lv_Dstarrest_D.p3() - cos_thetaV*lv_Dstarrest_D.p3().perp()*lv_brest_Dstar.p3().unit();
115            float chi = atan2(LTrans.cross(VTrans).dot(lv_brest_w.p3().unit()), LTrans.dot(VTrans));
116            if(chi < 0) chi += TWOPI;
117
118            _h_chi->fill(chi);
119          }
120        }
121      }
122    }
123
124    /// Normalise histograms etc., after the run
125    void finalize() {
126
127      double GAMMA_B0 = 4.32e-13 *1e15; // Normalise histos to partial width ifetime converted to GeV (x10^-15)
128      // factor 1/2 as e+mu in hists 
129      scale(_h_w,      0.5*GAMMA_B0/ *_nB);
130      scale(_h_costhv, 0.5*GAMMA_B0/ *_nB);
131      scale(_h_costhl, 0.5*GAMMA_B0/ *_nB);
132      scale(_h_chi,    0.5*GAMMA_B0/ *_nB);
133      Estimate1DPtr tmp;
134      book(tmp,1,1,1);
135      barchart(_h_w,tmp);
136      book(tmp,2,1,1);
137      barchart(_h_costhv,tmp);
138      book(tmp,3,1,1);
139      barchart(_h_costhl,tmp);
140      book(tmp,4,1,1);
141      barchart(_h_chi,tmp);
142    }
143
144    /// @}
145
146
147    /// @name Histograms
148    /// @{
149    Histo1DPtr _h_w;
150    Histo1DPtr _h_costhv;
151    Histo1DPtr _h_costhl;
152    Histo1DPtr _h_chi;
153    CounterPtr _nB;
154    /// @}
155
156
157  };
158
159
160  RIVET_DECLARE_PLUGIN(BELLE_2017_I1512299);
161
162
163}