Rivet analyses referenceBELLE_2016_I1504055Angular Analysis of $B\to K^{*}\ell^-\ell^-$Experiment: BELLE (KEKB) Inspire ID: 1504055 Status: VALIDATED NOHEPDATA Authors:
Beam energies: ANY Run details:
Measurement of angular coefficients in $B\to K^{*}\ell^-\ell^-$, the code implements these by taking appropriate moments Source code: BELLE_2016_I1504055.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Projections/DecayedParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief B -> K* l+l-
10 class BELLE_2016_I1504055 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2016_I1504055);
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 UnstableParticles ufs = UnstableParticles(Cuts::abspid==511 or
24 Cuts::abspid==521);
25 declare(ufs, "UFS");
26 DecayedParticles BB(ufs);
27 BB.addStable( 443);
28 BB.addStable(100443);
29 BB.addStable( 313);
30 BB.addStable( 323);
31 BB.addStable(-313);
32 BB.addStable(-323);
33 declare(BB, "BB");
34 for(unsigned int ix=0;ix<2;++ix) {
35 for(unsigned int iy=0;iy<6;++iy) {
36 book(_p_P[ix][iy],"TMP/p_P_"+toString(ix)+"_"+toString(iy),refData(1,1+ix,1+iy));
37 if(iy>1) continue;
38 book(_p_Q[ix][iy],"TMP/p_Q_"+toString(ix)+"_"+toString(iy),refData(2,1+ix,1+iy));
39 }
40 }
41 book(_FL,"TMP/FL");
42 book(_norm,"TMP/norm");
43 }
44
45
46 /// Perform the per-event analysis
47 void analyze(const Event& event) {
48 static const map<PdgId,unsigned int> & mode1 = { { 323,1},{ 13,1}, {-13,1}};
49 static const map<PdgId,unsigned int> & mode1CC = { {-323,1},{ 13,1}, {-13,1}};
50 static const map<PdgId,unsigned int> & mode2 = { { 313,1},{ 13,1}, {-13,1}};
51 static const map<PdgId,unsigned int> & mode2CC = { {-313,1},{ 13,1}, {-13,1}};
52 static const map<PdgId,unsigned int> & mode3 = { { 323,1},{ 11,1}, {-11,1}};
53 static const map<PdgId,unsigned int> & mode3CC = { {-323,1},{ 11,1}, {-11,1}};
54 static const map<PdgId,unsigned int> & mode4 = { { 313,1},{ 11,1}, {-11,1}};
55 static const map<PdgId,unsigned int> & mode4CC = { {-313,1},{ 11,1}, {-11,1}};
56 DecayedParticles BB = apply<DecayedParticles>(event, "BB");
57 // loop over particles
58 for(unsigned int ix=0;ix<BB.decaying().size();++ix) {
59 int imode=0;
60 if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode1)) ||
61 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode1CC))) imode=0;
62 else if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode2)) ||
63 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode2CC))) imode=1;
64 else if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode3)) ||
65 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode3CC))) imode=2;
66 else if ((BB.decaying()[ix].pid()>0 && BB.modeMatches(ix,3,mode4)) ||
67 (BB.decaying()[ix].pid()<0 && BB.modeMatches(ix,3,mode4CC))) imode=3;
68 else continue;
69 int il = imode<2 ? 13 : 11;
70 int sign = BB.decaying()[ix].pid()>0 ? 1 : -1;
71 const Particle & lp = BB.decayProducts()[ix].at(-sign*il)[0];
72 const Particle & lm = BB.decayProducts()[ix].at( sign*il)[0];
73 double qq = (lp.momentum()+lm.momentum()).mass2();
74 int iK = BB.decaying()[ix].abspid()==521 ? 323 : 313;
75 iK *= BB.decaying()[ix].pid()/BB.decaying()[ix].abspid();
76 const Particle & Kstar = BB.decayProducts()[ix].at( iK)[0];
77 if(Kstar.children().size()!=2) continue;
78 Particle KK;
79 if(Kstar.abspid()==313) {
80 if(Kstar.children()[0].abspid()==321 &&
81 Kstar.children()[1].abspid()==211)
82 KK = Kstar.children()[0];
83 else if(Kstar.children()[1].abspid()==321 &&
84 Kstar.children()[0].abspid()==211)
85 KK = Kstar.children()[1];
86 else continue;
87 }
88 else {
89 if(Kstar.children()[0].abspid()==311 &&
90 Kstar.children()[1].abspid()==211)
91 KK = Kstar.children()[0];
92 else if(Kstar.children()[1].abspid()==311 &&
93 Kstar.children()[0].abspid()==211)
94 KK = Kstar.children()[1];
95 else if(Kstar.children()[0].abspid()==310 &&
96 Kstar.children()[1].abspid()==211)
97 KK = Kstar.children()[0];
98 else if(Kstar.children()[1].abspid()==310 &&
99 Kstar.children()[0].abspid()==211)
100 KK = Kstar.children()[1];
101 else if(Kstar.children()[0].abspid()==321 &&
102 Kstar.children()[1].abspid()==111 )
103 KK = Kstar.children()[0];
104 else if(Kstar.children()[1].abspid()==321 &&
105 Kstar.children()[0].abspid()==111 )
106 KK = Kstar.children()[1];
107 else continue;
108 if(KK.abspid()==311) {
109 if(KK.children().size()==1 && KK.children()[0].pid()==310)
110 KK = KK.children()[0];
111 else
112 continue;
113 }
114 }
115 // first boost to bottom frame
116 const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(BB.decaying()[ix].momentum().betaVec());
117 FourMomentum plp = boost.transform(lp .momentum());
118 FourMomentum plm = boost.transform(lm .momentum());
119 FourMomentum pKstar = boost.transform(Kstar.momentum());
120 FourMomentum pK = boost.transform(KK .momentum());
121 FourMomentum pB = boost.transform(BB.decaying()[ix].momentum());
122 // lepton stuff
123 const LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta((plp+plm).betaVec());
124 plp = boost2.transform(plp);
125 Vector3 axis1 = boost .transform(pB ).p3().unit();
126 double cThetaL = plp.p3().unit().dot(axis1);
127 Vector3 Trans1 = plp.p3() - cThetaL*plp.p3().mod()*axis1;
128 // kaon stuff
129 const LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pKstar.betaVec());
130 pK = boost3.transform(pK);
131 Vector3 axis2 = boost .transform(pB ).p3().unit();
132 double cThetaK = pK.p3().unit().dot(axis2);
133 double FL = .5*(5.*sqr(cThetaK)-1.);
134 Vector3 Trans2 = pK.p3() - cThetaK*pK.p3().mod()*axis2;
135 double phi = atan2(Trans1.cross(Trans2).dot(axis2),Trans1.dot(Trans2));
136 double sThetaL = sqrt(1.-sqr(cThetaL));
137 double sThetaK = sqrt(1.-sqr(cThetaK));
138 double S4 = 12.5*cThetaL*sThetaL*cThetaK*sThetaK*cos(phi);
139 double S5 = 5.*cThetaK*sThetaK*sThetaL*sin(phi);
140 _FL->fill(FL);
141 _norm->fill();
142 for(unsigned int ix=0;ix<2;++ix) {
143 _p_P[ix][0]->fill(qq,S4);
144 _p_P[ix][3]->fill(qq,S5);
145 if(il==11) {
146 _p_P[ix][1]->fill(qq,S4);
147 _p_P[ix][4]->fill(qq,S5);
148 _p_Q[ix][0]->fill(qq,-S4);
149 _p_Q[ix][1]->fill(qq,-S5);
150 }
151 else {
152 _p_P[ix][2]->fill(qq,S4);
153 _p_P[ix][5]->fill(qq,S5);
154 _p_Q[ix][0]->fill(qq,S4);
155 _p_Q[ix][1]->fill(qq,S5);
156 }
157 }
158 }
159 }
160
161
162 /// Normalise histograms etc., after the run
163 void finalize() {
164 Estimate0D R = *_FL/ *_norm;
165 double fl = R.val();
166 double fact = 1./sqrt(fl*(1.-fl));
167 for(unsigned int ix=0;ix<2;++ix) {
168 for(unsigned int iy=0;iy<6;++iy) {
169 Estimate1DPtr tmp;
170 book(tmp,1,1+ix,1+iy);
171 barchart(_p_P[ix][iy],tmp);
172 scale(tmp,fact);
173 if(iy>1) continue;
174 book(tmp,2,1+ix,1+iy);
175 barchart(_p_Q[ix][iy],tmp);
176 scale(tmp,fact);
177 }
178 }
179 }
180
181 /// @}
182
183
184 /// @name Histograms
185 /// @{
186 Profile1DPtr _p_P[2][6],_p_Q[2][2];
187 CounterPtr _FL,_norm;
188 /// @}
189
190 };
191
192
193 RIVET_DECLARE_PLUGIN(BELLE_2016_I1504055);
194
195}
|
