1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Bs0 -> Ds* Ds*
  9  class BELLE_2013_I1124584 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2013_I1124584);
 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      UnstableParticles ufs = UnstableParticles(Cuts::abspid==531);
 23      declare(ufs, "UFS");
 24      // histograms
 25      for(unsigned int ix=0;ix<2;++ix)
 26	book(_h[ix],1,1,1+ix);
 27    }
 28
 29
 30    /// Perform the per-event analysis
 31    void analyze(const Event& event) {
 32      Particles BS0 = apply<UnstableParticles>(event, "UFS").particles();
 33      for(const Particle & p : BS0) {
 34	if(p.children().size()!=2) continue;
 35	if(p.children()[0].pid()!=-p.children()[1].pid()) continue;
 36	if(p.children()[0].abspid()!=433) continue;
 37	Particle Dp = p.children()[0];
 38	Particle Dm = p.children()[1];
 39	if     (p.pid()>0 && Dp.pid()<0) swap(Dp,Dm);
 40	else if(p.pid()<0 && Dp.pid()>0) swap(Dp,Dm);
 41	// boost to rest frame
 42	LorentzTransform boostB = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 43	FourMomentum pB =  boostB.transform(p.momentum());
 44	if(Dp.children().size()==2) {
 45	  Particle gamma;
 46	  bool found = true;
 47	  if(Dp.children()[0].pid()==PID::GAMMA &&
 48	     Dp.children()[1].abspid()==431)
 49	    gamma = Dp.children()[0];
 50	  else if (Dp.children()[1].pid()==PID::GAMMA &&
 51		   Dp.children()[0].abspid()==431)
 52	    gamma = Dp.children()[1];
 53	  else
 54	    found = false;
 55	  if( found) {
 56	    FourMomentum pD     = boostB.transform(Dp.momentum());
 57	    FourMomentum pgamma = boostB.transform(gamma.momentum());
 58	    LorentzTransform boostD = LorentzTransform::mkFrameTransformFromBeta(pD.betaVec());
 59	    Vector3 axisB = boostD.transform(pB).p3().unit();
 60	    Vector3 axisG = boostD.transform(pgamma).p3().unit();
 61	    _h[0]->fill(axisB.dot(axisG));
 62	  }
 63	}
 64	if(Dm.children().size()==2) {
 65	  Particle gamma;
 66	  bool found = true;
 67	  if(Dm.children()[0].pid()==PID::GAMMA &&
 68	     Dm.children()[1].abspid()==431)
 69	    gamma = Dm.children()[0];
 70	  else if (Dm.children()[1].pid()==PID::GAMMA &&
 71		   Dm.children()[0].abspid()==431)
 72	    gamma = Dm.children()[1];
 73	  else
 74	    found = false;
 75	  if( found) {
 76	    FourMomentum pD     = boostB.transform(Dm.momentum());
 77	    FourMomentum pgamma = boostB.transform(gamma.momentum());
 78	    LorentzTransform boostD = LorentzTransform::mkFrameTransformFromBeta(pD.betaVec());
 79	    Vector3 axisB = boostD.transform(pB).p3().unit();
 80	    Vector3 axisG = boostD.transform(pgamma).p3().unit();
 81	    _h[1]->fill(axisB.dot(axisG));
 82	  }
 83	}
 84      }
 85    }
 86
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      for(unsigned int ix=0;ix<2;++ix)
 91	normalize(_h[ix],1.,false);
 92    }
 93
 94    /// @}
 95
 96
 97    /// @name Histograms
 98    /// @{
 99    Histo1DPtr _h[2];
100    /// @}
101
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(BELLE_2013_I1124584);
107
108}