1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> Xs l+l-
  9  class BELLE_2016_I1283183 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2016_I1283183);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(Cuts::abspid==521 || Cuts::abspid==511), "UFS");
 23      // book the profile hist
 24      for(unsigned int ix=0;ix<4;++ix)
 25        book(_p[ix],1,1,ix+1);
 26    }
 27
 28    void findDecayProducts(bool & charm, const Particle& mother,
 29                           unsigned int& nK0, unsigned int& nKp,
 30			   unsigned int& nKm, Particles & lp, Particles & lm) {
 31      for (const Particle & p : mother.children()) {
 32        int id = p.pid();
 33	if(PID::isCharmHadron(p.pid())) charm = true;
 34	else if ( id == PID::POSITRON || id == PID::ANTIMUON) lp.push_back(p);
 35	else if ( id == PID::ELECTRON || id == PID::MUON    ) lm.push_back(p);
 36        else if ( id == PID::KPLUS )    ++nKp;
 37        else if (id == PID::KMINUS )    ++nKm;
 38        else if (id == PID::K0S)        ++nK0;
 39        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
 40	  continue;
 41        }
 42        else if ( !p.children().empty() ) {
 43          findDecayProducts(charm, p, nK0, nKp, nKm,lp,lm);
 44        }
 45      }
 46    }
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      if(_edges.empty()) {
 51        for(unsigned int ix=0;ix<4;++ix)
 52          _edges.push_back(_p[ix]->xEdges()[0]);
 53      }
 54      // Loop over bottoms
 55      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
 56       	// remove mixing entries etc
 57	if(bottom.children()[0].abspid()==bottom.abspid()) continue;
 58	bool charm = false;
 59	Particles lp,lm;
 60	unsigned int nK0(0),nKp(0),nKm(0);
 61	findDecayProducts(charm,bottom, nK0, nKp, nKm,lp,lm);
 62	if(charm) continue;
 63        unsigned int nk = nKp-nKm+nK0;
 64	if( nk % 2 == 0) continue;
 65	if (lp.size()!=1 || lm.size()!=1 || lp[0].pid()!=-lm[0].pid()) continue;
 66	if(bottom.pid()>0) swap(lp,lm);
 67	double q2 = (lp[0].momentum()+lm[0].momentum()).mass2();
 68	// veto region valid for muons but not electrons
 69	if(lm[0].pid()==PID::ELECTRON && ( (q2>7.3 && q2<10.5) || (q2>11.8 && q2<12.5) )) continue;
 70	// first boost to bottom frame
 71	const LorentzTransform boost  = LorentzTransform::mkFrameTransformFromBeta(bottom.momentum().betaVec());
 72	FourMomentum plp = boost.transform(lp[0] .momentum());
 73	FourMomentum plm = boost.transform(lm[0] .momentum());
 74	FourMomentum pB  = boost.transform(bottom.momentum());
 75	const LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta((plp+plm).betaVec());
 76	plp = boost2.transform(plp);
 77	pB  = boost .transform(pB );
 78	double cTheta = plp.p3().unit().dot(pB.p3().unit());
 79        if(q2>0.2&&q2<4.3)
 80          _p[0]->fill(_edges[0], cTheta>0 ? 1 : -1);
 81        else if(q2>4.3 && q2<8.1)
 82          _p[1]->fill(_edges[1], cTheta>0 ? 1 : -1);
 83        else if(q2>10.2 && q2<12.5)
 84          _p[2]->fill(_edges[2], cTheta>0 ? 1 : -1);
 85        else if(q2>14.3 && q2<25.0)
 86          _p[3]->fill(_edges[3], cTheta>0 ? 1 : -1);
 87      }
 88    }
 89
 90
 91    /// Normalise histograms etc., after the run
 92    void finalize() {
 93    }
 94
 95    /// @}
 96
 97
 98    /// @name Histograms
 99    /// @{
100    BinnedProfilePtr<string> _p[4];
101    vector<string> _edges;
102    /// @}
103
104
105  };
106
107
108  RIVET_DECLARE_PLUGIN(BELLE_2016_I1283183);
109
110}