1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief B+ -> psi(3770) K+
 9  class BELLE_2004_I623994 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2004_I623994);
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==521);
23      declare(ufs, "UFS");
24      // histos
25      book(_h,1,1,1);
26    }
27
28
29    /// Perform the per-event analysis
30    void analyze(const Event& event) {
31      UnstableParticles ufs = apply<UnstableParticles>(event,"UFS");
32      for (const Particle& B : ufs.particles()) {
33      	if (B.children().size()!=2) continue;
34      	int sign = B.pid()/B.abspid();
35        // boost to B rest frame
36        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(B.mom().betaVec());
37        Particle resonance;
38        if (B.children()[0].pid()==30443 && B.children()[1].pid()==sign*321) {
39          resonance=B.children()[0];
40        }
41        else if (B.children()[1].pid()==30443 && B.children()[0].pid()==sign*321) {
42          resonance=B.children()[1];
43        }
44        else {
45          continue;
46        }
47        if (resonance.children().size()!=2) continue;
48        FourMomentum pRes = boost1.transform(resonance.mom());
49        LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pRes.betaVec());
50        Vector3 axis1 = pRes.p3().unit();
51        Particle child;
52        if (resonance.children()[0].pid()!=-resonance.children()[1].pid()) continue;
53        if (resonance.children()[0].abspid()!=411 && resonance.children()[0].abspid()!=421) continue;
54        child = resonance.children()[0].pid()>0 ? resonance.children()[0] : resonance.children()[1];
55        Vector3 axis2 = boost2.transform(boost1.transform(child.mom())).p3().unit();
56      	_h->fill(axis1.dot(axis2));
57      }
58    }
59
60
61    /// Normalise histograms etc., after the run
62    void finalize() {
63      normalize(_h, 1.0, false);
64    }
65
66    /// @}
67
68
69    /// @name Histograms
70    /// @{
71    Histo1DPtr _h;
72    /// @}
73
74
75  };
76
77
78  RIVET_DECLARE_PLUGIN(BELLE_2004_I623994);
79
80}