Belle II Software  release-08-01-10
EvtB0toKsKK.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <iomanip>
10 #include <cmath>
11 #include <string>
12 
13 #include "EvtGenBase/EvtCPUtil.hh"
14 #include "EvtGenBase/EvtPDL.hh"
15 #include "EvtGenBase/EvtTensor4C.hh"
16 
17 #include <generators/evtgen/EvtGenModelRegister.h>
18 #include "generators/evtgen/models/EvtB0toKsKK.h"
19 
20 namespace Belle2 {
28 
29  EvtB0toKsKK::~EvtB0toKsKK() {}
30 
31  std::string EvtB0toKsKK::getName()
32  {
33  return "B0toKsKK";
34  }
35 
36  EvtDecayBase* EvtB0toKsKK::clone()
37  {
38  return new EvtB0toKsKK;
39  }
40 
42  {
43  // Check number of arguments
44  checkNArg(32);
45 
46  // Check number of daughters
47  checkNDaug(3);
48 
49  // Check decay chain
50  if ((getParentId() != EvtPDL::getId("B0") &&
51  getParentId() != EvtPDL::getId("anti-B0")) ||
52  (getDaug(0) != EvtPDL::getId("K_S0")) ||
53  (getDaug(1) != EvtPDL::getId("K+") &&
54  getDaug(1) != EvtPDL::getId("K-")) ||
55  (getDaug(2) != EvtPDL::getId("K+") &&
56  getDaug(2) != EvtPDL::getId("K-"))) {
57  std::cout << "ERROR: Invalid decay" << std::endl;
58  std::cout << "USAGE: K_S0 K+ K-" << std::endl;
59  exit(1);
60  }
61 
62  a_f0ks_ =
63  (1.0 + getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
64  getArg(1) * sin(getArg(2) * M_PI / 180.0));
65  a_phiks_ =
66  (1.0 + getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
67  getArg(5) * sin(getArg(6) * M_PI / 180.0));
68  a_fxks_ =
69  (1.0 + getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
70  getArg(9) * sin(getArg(10) * M_PI / 180.0));
71  a_chic0ks_ =
72  (1.0 + getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
73  getArg(13) * sin(getArg(14) * M_PI / 180.0));
74  a_kpkmnr_ =
75  (1.0 + getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
76  getArg(17) * sin(getArg(18) * M_PI / 180.0));
77  a_kskpnr_ =
78  (1.0 + getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
79  getArg(22) * sin(getArg(23) * M_PI / 180.0));
80  a_kskmnr_ =
81  (1.0 + getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
82  getArg(27) * sin(getArg(28) * M_PI / 180.0));
83 
84  abar_f0ks_ =
85  (1.0 - getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
86  getArg(1) * sin(getArg(2) * M_PI / 180.0));
87  abar_phiks_ =
88  (1.0 - getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
89  getArg(5) * sin(getArg(6) * M_PI / 180.0));
90  abar_fxks_ =
91  (1.0 - getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
92  getArg(9) * sin(getArg(10) * M_PI / 180.0));
94  (1.0 - getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
95  getArg(13) * sin(getArg(14) * M_PI / 180.0));
96  abar_kpkmnr_ =
97  (1.0 - getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
98  getArg(17) * sin(getArg(18) * M_PI / 180.0));
99  abar_kskpnr_ =
100  (1.0 - getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
101  getArg(22) * sin(getArg(23) * M_PI / 180.0));
102  abar_kskmnr_ =
103  (1.0 - getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
104  getArg(27) * sin(getArg(28) * M_PI / 180.0));
105 
106  alpha_kpkmnr = getArg(21);
107  alpha_kskpnr = getArg(26);
108  alpha_kskmnr = getArg(31);
109 
110  std::cout << setiosflags(std::ios::left) << std::endl
111  << "B0 Channel " << std::setw(20) << "Relative amplitude" << std::setw(20) << "Relative phase" << std::endl
112  << "f0ks " << std::setw(20) << real(a_f0ks_) << std::setw(20) << imag(a_f0ks_) << std::endl
113  << "phiks " << std::setw(20) << real(a_phiks_) << std::setw(20) << imag(a_phiks_) << std::endl
114  << "fxks " << std::setw(20) << real(a_fxks_) << std::setw(20) << imag(a_fxks_) << std::endl
115  << "chic0ks " << std::setw(20) << real(a_chic0ks_) << std::setw(20) << imag(a_chic0ks_) << std::endl
116  << "kpkmnr " << std::setw(20) << real(a_kpkmnr_) << std::setw(20) << imag(a_kpkmnr_) << std::endl
117  << "kskpnr " << std::setw(20) << real(a_kskpnr_) << std::setw(20) << imag(a_kskpnr_) << std::endl
118  << "kskmnr " << std::setw(20) << real(a_kskmnr_) << std::setw(20) << imag(a_kskmnr_) << std::endl
119  << std::endl;
120 
121  std::cout << setiosflags(std::ios::left) << std::endl
122  << "B0B Channel " << std::setw(20) << "Relative amplitude" << std::setw(20) << "Relative phase" << std::endl
123  << "f0ks " << std::setw(20) << real(abar_f0ks_) << std::setw(20) << imag(abar_f0ks_) << std::endl
124  << "phiks " << std::setw(20) << real(abar_phiks_) << std::setw(20) << imag(abar_phiks_) << std::endl
125  << "fxks " << std::setw(20) << real(abar_fxks_) << std::setw(20) << imag(abar_fxks_) << std::endl
126  << "chic0ks " << std::setw(20) << real(abar_chic0ks_) << std::setw(20) << imag(abar_chic0ks_) << std::endl
127  << "kpkmnr " << std::setw(20) << real(abar_kpkmnr_) << std::setw(20) << imag(abar_kpkmnr_) << std::endl
128  << "kskpnr " << std::setw(20) << real(abar_kskpnr_) << std::setw(20) << imag(abar_kskpnr_) << std::endl
129  << "kskmnr " << std::setw(20) << real(abar_kskmnr_) << std::setw(20) << imag(abar_kskmnr_) << std::endl
130  << std::endl;
131 
132  // Open debugging file
133  debugfile_.open("debug.dat", std::ios::out);
134  }
135 
137  {
138  //setProbMax(50000.0);
139  //setProbMax(100000.0);
140  //setProbMax(200000.0);
141  //setProbMax(400000.0);
142  //setProbMax(600000.0);
143  setProbMax(1100000.0);
144  }
145 
146  void EvtB0toKsKK::decay(EvtParticle* p)
147  {
148  // Btag
149  static EvtId B0 = EvtPDL::getId("B0");
150  static EvtId B0B = EvtPDL::getId("anti-B0");
151 
152  double t;
153  EvtId other_b;
154 
155  //std::cout << EvtCPUtil::getInstance()->getMixingType() << std::endl;
156  //EvtCPUtil::getInstance()->setMixingType(0);
157  EvtCPUtil::getInstance()->OtherB(p, t, other_b, 0.5);
158  //EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.4);
159  //EvtCPUtil::getInstance()->OtherB(p,t,other_b);
160 
161  // Brec
162  p->initializePhaseSpace(getNDaug(), getDaugs());
163  EvtVector4R p4ks, p4kp, p4km;
164  if (p->getId() == B0) {
165  p4ks = p->getDaug(0)->getP4();
166  p4kp = p->getDaug(1)->getP4();
167  p4km = p->getDaug(2)->getP4();
170  /*p4ks = p->getDaug(0)->getP4();
171  p4kp = p->getDaug(2)->getP4();
172  p4km = p->getDaug(1)->getP4();*/
173  } else {
174  p4ks = p->getDaug(0)->getP4();
175  p4kp = p->getDaug(2)->getP4();
176  p4km = p->getDaug(1)->getP4();
178  /*p4ks = p->getDaug(0)->getP4();
179  p4kp = p->getDaug(1)->getP4();
180  p4km = p->getDaug(2)->getP4();*/
181  }
182 
183  /*std::cout << (p4ks + p4kp).mass() << " " << (p4ks + p4km).mass() << " "
184  << (p4kp + p4km).mass() << std::endl;
185  if( p->getId() == other_b )
186  std::cout << "same flavour" << std::endl;
187  std::cout << p->getId() << " --> " << p->getDaug(0)->getId() << " "
188  << p->getDaug(1)->getId() << " " << p->getDaug(2)->getId()
189  << std::endl;
190  std::cout << B0 << std::endl;
191  std::cout << B0B << std::endl;
192  std::cout << other_b << std::endl;*/
193  /*std::cout << p->getLifetime() << std::endl;
194  std::cout << p->get4Pos() << std::endl;
195  std::cout << p->getP4() << std::endl;
196  std::cout << p->getP4Lab() << std::endl;
197  EvtParticle *parent=p->getParent();
198  if (parent->getDaug(0)!=p){
199  std::cout << parent->getDaug(0)->getLifetime() << std::endl;
200  std::cout << parent->getDaug(0)->get4Pos() << std::endl;
201  std::cout << parent->getDaug(0)->getP4() << std::endl;
202  std::cout << parent->getDaug(0)->getP4Lab() << std::endl;
203  std::cout << parent->getDaug(0)->getP4Lab()+p->getP4Lab() << std::endl;
204  }
205  else{
206  std::cout << parent->getDaug(1)->getLifetime() << std::endl;
207  std::cout << parent->getDaug(1)->get4Pos() << std::endl;
208  std::cout << parent->getDaug(1)->getP4() << std::endl;
209  std::cout << parent->getDaug(1)->getP4Lab() << std::endl;
210  std::cout << parent->getDaug(1)->getP4Lab()+p->getP4Lab() << std::endl;
211  }
212  std::cout << parent->getP4() << std::endl;
213  std::cout << t << std::endl;*/
214 
215  // Relative amplides and phases with direct CP violation
216  a_f0ks_ =
217  (1.0 + getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
218  getArg(1) * sin(getArg(2) * M_PI / 180.0));
219  a_phiks_ =
220  (1.0 + getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
221  getArg(5) * sin(getArg(6) * M_PI / 180.0));
222  a_fxks_ =
223  (1.0 + getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
224  getArg(9) * sin(getArg(10) * M_PI / 180.0));
225  a_chic0ks_ =
226  (1.0 + getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
227  getArg(13) * sin(getArg(14) * M_PI / 180.0));
228  a_kpkmnr_ =
229  (1.0 + getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
230  getArg(17) * sin(getArg(18) * M_PI / 180.0));
231  a_kskpnr_ =
232  (1.0 + getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
233  getArg(22) * sin(getArg(23) * M_PI / 180.0));
234  a_kskmnr_ =
235  (1.0 + getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
236  getArg(27) * sin(getArg(28) * M_PI / 180.0));
237 
238  abar_f0ks_ =
239  (1.0 - getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
240  getArg(1) * sin(getArg(2) * M_PI / 180.0));
241  abar_phiks_ =
242  (1.0 - getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
243  getArg(5) * sin(getArg(6) * M_PI / 180.0));
244  abar_fxks_ =
245  (1.0 - getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
246  getArg(9) * sin(getArg(10) * M_PI / 180.0));
247  abar_chic0ks_ =
248  (1.0 - getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
249  getArg(13) * sin(getArg(14) * M_PI / 180.0));
250  abar_kpkmnr_ =
251  (1.0 - getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
252  getArg(17) * sin(getArg(18) * M_PI / 180.0));
253  abar_kskpnr_ =
254  (1.0 - getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
255  getArg(22) * sin(getArg(23) * M_PI / 180.0));
256  abar_kskmnr_ =
257  (1.0 - getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
258  getArg(27) * sin(getArg(28) * M_PI / 180.0));
259 
260  // Mixing-induced CP asymmetry
261  const double pCP_f0ks = getArg(4) * M_PI / 180.0;
262  const double pCP_phiks = getArg(8) * M_PI / 180.0;
263  const double pCP_fxks = getArg(12) * M_PI / 180.0;
264  const double pCP_chic0ks = getArg(16) * M_PI / 180.0;
265  const double pCP_kpkmnr = getArg(20) * M_PI / 180.0;
266  const double pCP_kskpnr = getArg(25) * M_PI / 180.0;
267  const double pCP_kskmnr = getArg(30) * M_PI / 180.0;
268 
269  if (other_b == B0) {
270  a_f0ks_ *=
271  EvtComplex(cos(+2.0 * pCP_f0ks), sin(+2.0 * pCP_f0ks));
272  a_phiks_ *=
273  EvtComplex(cos(+2.0 * pCP_phiks), sin(+2.0 * pCP_phiks));
274  a_fxks_ *=
275  EvtComplex(cos(+2.0 * pCP_fxks), sin(+2.0 * pCP_fxks));
276  a_chic0ks_ *=
277  EvtComplex(cos(+2.0 * pCP_chic0ks), sin(+2.0 * pCP_chic0ks));
278  a_kpkmnr_ *=
279  EvtComplex(cos(+2.0 * pCP_kpkmnr), sin(+2.0 * pCP_kpkmnr));
280  a_kskpnr_ *=
281  EvtComplex(cos(+2.0 * pCP_kskpnr), sin(+2.0 * pCP_kskpnr));
282  a_kskmnr_ *=
283  EvtComplex(cos(+2.0 * pCP_kskmnr), sin(+2.0 * pCP_kskmnr));
284  }
285  if (other_b == B0B) {
286  abar_f0ks_ *=
287  EvtComplex(cos(-2.0 * pCP_f0ks), sin(-2.0 * pCP_f0ks));
288  abar_phiks_ *=
289  EvtComplex(cos(-2.0 * pCP_phiks), sin(-2.0 * pCP_phiks));
290  abar_fxks_ *=
291  EvtComplex(cos(-2.0 * pCP_fxks), sin(-2.0 * pCP_fxks));
292  abar_chic0ks_ *=
293  EvtComplex(cos(-2.0 * pCP_chic0ks), sin(-2.0 * pCP_chic0ks));
294  abar_kpkmnr_ *=
295  EvtComplex(cos(-2.0 * pCP_kpkmnr), sin(-2.0 * pCP_kpkmnr));
296  abar_kskpnr_ *=
297  EvtComplex(cos(-2.0 * pCP_kskpnr), sin(-2.0 * pCP_kskpnr));
298  abar_kskmnr_ *=
299  EvtComplex(cos(-2.0 * pCP_kskmnr), sin(-2.0 * pCP_kskmnr));
300  }
301 
302  //Form Factors
303  EvtComplex Amp_f0ks = A_f0ks(p4ks, p4kp, p4km);
304  EvtComplex Amp_phiks = A_phiks(p4ks, p4kp, p4km);
305  EvtComplex Amp_fxks = A_fxks(p4ks, p4kp, p4km);
306  EvtComplex Amp_chic0ks = A_chic0ks(p4ks, p4kp, p4km);
307  EvtComplex Amp_kpkmnr = A_kknr(p4kp, p4km, alpha_kpkmnr);
308  EvtComplex Amp_kskpnr = A_kknr(p4ks, p4kp, alpha_kskpnr);
309  EvtComplex Amp_kskmnr = A_kknr(p4ks, p4km, alpha_kskmnr);
310 
311  EvtComplex Ampbar_f0ks = A_f0ks(p4ks, p4km, p4kp);
312  EvtComplex Ampbar_phiks = A_phiks(p4ks, p4km, p4kp);
313  EvtComplex Ampbar_fxks = A_fxks(p4ks, p4km, p4kp);
314  EvtComplex Ampbar_chic0ks = A_chic0ks(p4ks, p4km, p4kp);
315  EvtComplex Ampbar_kpkmnr = A_kknr(p4km, p4kp, alpha_kpkmnr);
316  EvtComplex Ampbar_kskpnr = A_kknr(p4ks, p4km, alpha_kskpnr);
317  EvtComplex Ampbar_kskmnr = A_kknr(p4ks, p4kp, alpha_kskmnr);
318 
319  const EvtComplex A_B0toKsKK =
320  (a_f0ks_ * Amp_f0ks) +
321  (a_phiks_ * Amp_phiks) +
322  (a_fxks_ * Amp_fxks) +
323  (a_chic0ks_ * Amp_chic0ks) +
324  (a_kpkmnr_ * Amp_kpkmnr) +
325  (a_kskpnr_ * Amp_kskpnr) +
326  (a_kskmnr_ * Amp_kskmnr);
327 
328  const EvtComplex Abar_B0toKsKK =
329  (abar_f0ks_ * Ampbar_f0ks) +
330  (abar_phiks_ * Ampbar_phiks) +
331  (abar_fxks_ * Ampbar_fxks) +
332  (abar_chic0ks_ * Ampbar_chic0ks) +
333  (abar_kpkmnr_ * Ampbar_kpkmnr) +
334  (abar_kskpnr_ * Ampbar_kskpnr) +
335  (abar_kskmnr_ * Ampbar_kskmnr);
336 
337  // CP asymmetry
338  const double dm = getArg(0);
339 
340  EvtComplex amp;
341  if (other_b == B0B) {
342  amp =
343  A_B0toKsKK * cos(dm * t / (2.0 * EvtConst::c)) +
344  EvtComplex(0.0, 1.0) * Abar_B0toKsKK * sin(dm * t / (2.0 * EvtConst::c));
345  }
346  if (other_b == B0) {
347  amp =
348  A_B0toKsKK *
349  EvtComplex(0.0, 1.0) * sin(dm * t / (2.0 * EvtConst::c)) +
350  Abar_B0toKsKK * cos(dm * t / (2.0 * EvtConst::c));
351  }
352 
353  if (abs2(amp) > 50000.0)
354  debugfile_ << abs2(amp) << std::endl;
355  vertex(amp);
356 
357  return;
358  }
359 
360  EvtVector4R EvtB0toKsKK::umu(const EvtVector4R& p4a, const EvtVector4R& p4b,
361  const EvtVector4R& p4c)
362  {
363  const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
364 
365  const EvtVector4R umu = (p4a + p4b + p4c) / sqrt(s);
366 
367  return umu;
368  }
369 
370  EvtVector4R EvtB0toKsKK::Smu(const EvtVector4R& p4a, const EvtVector4R& p4b,
371  const EvtVector4R& /*p4c*/)
372  {
373  const double ma = p4a.mass();
374  const double mb = p4b.mass();
375  const double mres = (p4a + p4b).mass();
376  const double ma2 = ma * ma;
377  const double mb2 = mb * mb;
378  const double mres2 = mres * mres;
379 
380  const double N1 =
381  mres /
382  (sqrt(mres2 - ((ma + mb) * (ma + mb))) * sqrt(mres2 - ((ma - mb) * (ma - mb))));
383 
384  const EvtVector4R Smu =
385  N1 * ((p4a - p4b) - ((p4a + p4b) * (ma2 - mb2) / mres2));
386 
387  return Smu;
388  }
389 
390  EvtVector4R EvtB0toKsKK::Lmu(const EvtVector4R& p4a, const EvtVector4R& p4b,
391  const EvtVector4R& p4c)
392  {
393  const double mc = p4c.mass();
394  const double mres = (p4a + p4b).mass();
395  const double mc2 = mc * mc;
396  const double mres2 = mres * mres;
397  const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
398 
399  const double N2 =
400  sqrt(s) /
401  (sqrt(s - ((mres + mc) * (mres + mc))) * sqrt(s - ((mres - mc) * (mres - mc))));
402 
403  const EvtVector4R Lmu =
404  N2 * ((p4c - (p4a + p4b)) - ((p4c + (p4a + p4b)) * (mc2 - mres2) / s));
405 
406  return Lmu;
407  }
408 
409  EvtTensor4C EvtB0toKsKK::gmunu_tilde(const EvtVector4R& p4a,
410  const EvtVector4R& p4b,
411  const EvtVector4R& p4c)
412  {
413  const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
414 
415  const EvtVector4R umu = (p4a + p4b + p4c) / sqrt(s);
416 
417  const EvtTensor4C gmunu_tilde =
418  EvtTensor4C::g() - EvtGenFunctions::directProd(umu, umu);
419 
420  return gmunu_tilde;
421  }
422 
423  EvtTensor4C EvtB0toKsKK::Tmunu(const EvtVector4R& p4a, const EvtVector4R& p4b,
424  const EvtVector4R& p4c)
425  {
426  const double mres = (p4a + p4b).mass();
427  const EvtVector4R wmu = (p4a + p4b) / mres;
428 
429  const EvtTensor4C Tmunu =
430  sqrt(3.0 / 2.0) *
431  (EvtGenFunctions::directProd(Smu(p4a, p4b, p4c), Smu(p4a, p4b, p4c)) +
432  (1.0 / 3.0 * (EvtTensor4C::g() - EvtGenFunctions::directProd(wmu, wmu))));
433 
434  return Tmunu;
435  }
436 
437  EvtTensor4C EvtB0toKsKK::Multiply(const EvtTensor4C& t1,
438  const EvtTensor4C& t2)
439  {
440  EvtTensor4C t;
441  for (unsigned int i = 0; i < 4; i++)
442  for (unsigned int j = 0; j < 4; j++) {
443  const EvtComplex element =
444  t1.get(i, 0) * t2.get(0, j) +
445  t1.get(i, 1) * t2.get(1, j) +
446  t1.get(i, 2) * t2.get(2, j) +
447  t1.get(i, 3) * t2.get(3, j);
448  t.set(i, j, element);
449  }
450  return t;
451  }
452 
453  EvtTensor4C EvtB0toKsKK::RaiseIndices(const EvtTensor4C& t)
454  {
455  const EvtTensor4C tmp = Multiply(t, EvtTensor4C::g());
456  const EvtTensor4C t_raised = Multiply(EvtTensor4C::g(), tmp);
457 
458  return t_raised;
459  }
460 
461  void EvtB0toKsKK::RaiseIndex(EvtVector4R& vector)
462  {
463  vector.set(1, -vector.get(1));
464  vector.set(2, -vector.get(2));
465  vector.set(3, -vector.get(3));
466  }
467 
468  EvtTensor4C EvtB0toKsKK::Mmunu(const EvtVector4R& p4a, const EvtVector4R& p4b,
469  const EvtVector4R& p4c)
470  {
471  const EvtTensor4C Mmunu =
472  sqrt(3.0 / 2.0) *
473  (EvtGenFunctions::directProd(Lmu(p4a, p4b, p4c), Lmu(p4a, p4b, p4c)) +
474  (1.0 / 3.0 * gmunu_tilde(p4a, p4b, p4c)));
475 
476  return Mmunu;
477  }
478 
479  double EvtB0toKsKK::BWBF(const double& q, const unsigned int& L)
480  {
481  //Meson radius
482  const double d = 1.0;
483 
484  const double z = q * q * d * d;
485 
486  double bwbf = 1.0;
487 
488  if (L == 1)
489  bwbf = sqrt(2.0 * z / (1.0 + z));
490  if (L == 2)
491  bwbf = sqrt(13.0 * z * z / (((z - 3.0) * (z - 3.0)) + (9.0 * z)));
492  if (L > 2) {
493  std::cout << "ERROR: BWBF not implemented for L>2" << std::endl;
494  exit(1);
495  }
496  return bwbf;
497  }
498 
499  double EvtB0toKsKK::BWBF(const double& q, const double& q0,
500  const unsigned int& L)
501  {
502  //Meson radius
503  const double d = 1.0;
504 
505  const double z = q * q * d * d;
506  const double z0 = q0 * q0 * d * d;
507 
508  double bwbf = 1.0;
509 
510  if (L == 1)
511  bwbf = sqrt((1.0 + z0) / (1.0 + z));
512  if (L == 2)
513  bwbf = sqrt((((z0 - 3.0) * (z0 - 3.0)) + (9.0 * z0)) / (((z - 3.0) * (z - 3.0)) + (9.0 * z)));
514  if (L > 2) {
515  std::cout << "ERROR: BWBF not implemented for L>2" << std::endl;
516  exit(1);
517  }
518  return bwbf;
519  }
520 
521  EvtComplex EvtB0toKsKK::BreitWigner(const double& m, const double& m0,
522  const double& Gamma0,
523  const double& q, const double& q0,
524  const unsigned int& L)
525  {
526  const double s = m * m;
527  const double s0 = m0 * m0;
528 
529  const double Gamma =
530  Gamma0 * m0 / m * pow(q / q0, (2 * L) + 1) *
531  BWBF(q, q0, L) * BWBF(q, q0, L);
532 
533  const EvtComplex BreitWigner = 1.0 / EvtComplex(s0 - s, -m0 * Gamma);
534 
535  return BreitWigner;
536  }
537 
538  EvtVector4R EvtB0toKsKK::Boost(const EvtVector4R& p4,
539  const EvtVector4R& boost)
540  {
541  const double bx = boost.get(1) / boost.get(0);
542  const double by = boost.get(2) / boost.get(0);
543  const double bz = boost.get(3) / boost.get(0);
544 
545  const double bx2 = bx * bx;
546  const double by2 = by * by;
547  const double bz2 = bz * bz;
548 
549  const double b2 = bx2 + by2 + bz2;
550  if (b2 == 0.0)
551  return p4;
552  assert(b2 < 1.0);
553 
554  const double gamma = 1.0 / sqrt(1 - b2);
555 
556  const double gb2 = (gamma - 1.0) / b2;
557 
558  const double gb2xy = gb2 * bx * by;
559  const double gb2xz = gb2 * bx * bz;
560  const double gb2yz = gb2 * by * bz;
561 
562  const double gbx = gamma * bx;
563  const double gby = gamma * by;
564  const double gbz = gamma * bz;
565 
566  const double e_b =
567  (gamma * p4.get(0)) - (gbx * p4.get(1)) - (gby * p4.get(2)) - (gbz * p4.get(3));
568  const double x_b =
569  (-gbx * p4.get(0)) + ((1.0 + (gb2 * bx2)) * p4.get(1)) +
570  (gb2xy * p4.get(2)) + (gb2xz * p4.get(3));
571  const double y_b =
572  (-gby * p4.get(0)) + (gb2xy * p4.get(1)) +
573  ((1.0 + (gb2 * by2)) * p4.get(2)) + (gb2yz * p4.get(3));
574  const double z_b =
575  (-gbz * p4.get(0)) + (gb2xz * p4.get(1)) +
576  (gb2yz * p4.get(2)) + ((1.0 + (gb2 * bz2)) * p4.get(3));
577 
578  const EvtVector4R p4_b(e_b, x_b, y_b, z_b);
579 
580  return p4_b;
581  }
582 
583  double EvtB0toKsKK::p(const double& mab, const double& M, const double& mc)
584  {
585  const double sab = mab * mab;
586  const double M_p_mc2 = (M + mc) * (M + mc);
587  const double M_m_mc2 = (M - mc) * (M - mc);
588 
589  const double p = sqrt((sab - M_p_mc2) * (sab - M_m_mc2)) / (2.0 * mab);
590 
591  return p;
592  }
593 
594  double EvtB0toKsKK::q(const double& mab, const double& ma, const double& mb)
595  {
596  const double mab2 = mab * mab;
597  const double ma_p_mb2 = (ma + mb) * (ma + mb);
598  const double ma_m_mb2 = (ma - mb) * (ma - mb);
599 
600  const double q = sqrt((mab2 - ma_p_mb2) * (mab2 - ma_m_mb2)) / (2.0 * mab);
601 
602  return q;
603  }
604 
605  EvtComplex EvtB0toKsKK::Flatte_k(const double& s, const double& m_h)
606  {
607  const double k2 = 1.0 - (4.0 * m_h * m_h / s);
608  EvtComplex k;
609  if (k2 < 0.0)
610  k = EvtComplex(0.0, sqrt(fabs(k2)));
611  else
612  k = sqrt(k2);
613 
614  return k;
615  }
616 
617  EvtComplex EvtB0toKsKK::Flatte(const double& m, const double& m0)
618  {
619  const double g_pipi = 0.165;
620  const double g_kk = 4.21 * g_pipi;
621 
622  static EvtId pip = EvtPDL::getId("pi+");
623  static EvtId kp = EvtPDL::getId("K+");
624 
625  const double s = m * m;
626  const double s0 = m0 * m0;
627 
628  const EvtComplex rho_pipi = 2.0 * Flatte_k(s, EvtPDL::getMeanMass(pip)) / m;
629  const EvtComplex rho_kk = 2.0 * Flatte_k(s, EvtPDL::getMeanMass(kp)) / m;
630 
631  const EvtComplex Gamma =
632  EvtComplex(0.0, 1.0) * ((g_pipi * rho_pipi) + (g_kk * rho_kk));
633 
634  const EvtComplex Flatte = 1.0 / (s0 - s - Gamma);
635 
636  return Flatte;
637  }
638 
639  EvtComplex EvtB0toKsKK::A_f0ks(const EvtVector4R& /*p4ks*/,
640  const EvtVector4R& p4kp,
641  const EvtVector4R& p4km)
642  {
643  static EvtId f0 = EvtPDL::getId("f_0");
644 
645  const double f0_m = EvtPDL::getMeanMass(f0);
646 
647  const EvtVector4R p4kpkm = p4kp + p4km;
648  const double mkpkm = p4kpkm.mass();
649 
650  //Angular distribution
651  const EvtComplex H_f0ks = 1.0;
652 
653  //Barrier factors
654  const EvtComplex D_f0ks = 1.0;
655 
656  //Line shape
657  const EvtComplex L_f0ks = Flatte(mkpkm, f0_m);
658 
659  //Amplitude
660  const EvtComplex A_f0ks = D_f0ks * H_f0ks * L_f0ks;
661 
662  return A_f0ks;
663  }
664 
666  double s13_min(const double& s12, const double& M,
667  const double& m1, const double& m2, const double& m3)
668  {
669  const double E2star = (s12 - (m2 * m2) + (m1 * m1)) / 2.0 / sqrt(s12);
670  const double E3star = (M * M - s12 - (m3 * m3)) / 2.0 / sqrt(s12);
671 
672  const double s23_min =
673  ((E2star + E3star) * (E2star + E3star)) -
674  ((sqrt((E2star * E2star) - (m1 * m1)) + sqrt((E3star * E3star) - (m3 * m3))) *
675  (sqrt((E2star * E2star) - (m1 * m1)) + sqrt((E3star * E3star) - (m3 * m3))));
676 
677  return s23_min;
678  }
679 
681  double s13_max(const double& s12, const double& M,
682  const double& m1, const double& m2, const double& m3)
683  {
684  const double E2star = (s12 - (m2 * m2) + (m1 * m1)) / 2.0 / sqrt(s12);
685  const double E3star = (M * M - s12 - (m3 * m3)) / 2.0 / sqrt(s12);
686 
687  const double s23_max =
688  ((E2star + E3star) * (E2star + E3star)) -
689  ((sqrt((E2star * E2star) - (m1 * m1)) - sqrt((E3star * E3star) - (m3 * m3))) *
690  (sqrt((E2star * E2star) - (m1 * m1)) - sqrt((E3star * E3star) - (m3 * m3))));
691 
692  return s23_max;
693  }
694 
695  EvtComplex EvtB0toKsKK::A_phiks(const EvtVector4R& p4ks,
696  const EvtVector4R& p4kp,
697  const EvtVector4R& p4km)
698  {
699  static EvtId phi = EvtPDL::getId("phi");
700 
701  const double phi_m = EvtPDL::getMeanMass(phi);
702  const double phi_w = EvtPDL::getWidth(phi);
703 
704  const EvtVector4R p4kpkm = p4kp + p4km;
705  const double mkpkm = p4kpkm.mass();
706 
707  //Angular distribution
708  const EvtVector4R S_mu = Smu(p4kp, p4km, p4ks);
709  const EvtVector4R L_mu = Lmu(p4kp, p4km, p4ks);
710 
711  const EvtTensor4C g_munu_tilde = gmunu_tilde(p4kp, p4km, p4ks);
712  const EvtTensor4C g__munu_tilde = RaiseIndices(g_munu_tilde);
713 
714  EvtComplex H_phiks = 0.0;
715  for (unsigned int i = 0; i < 4; i++)
716  for (unsigned int j = 0; j < 4; j++)
717  H_phiks += g__munu_tilde.get(i, j) * S_mu.get(i) * L_mu.get(j);
718 
719  //Barrier factors
720  const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
721 
722  const EvtComplex D_phiks = BWBF(p4kp_kpkm.d3mag(), 1);
723 
724  //Line shape
725  const double q0_phi = q(phi_m, p4kp.mass(), p4km.mass());
726  const EvtComplex L_phiks =
727  BreitWigner(mkpkm, phi_m, phi_w, p4kp_kpkm.d3mag(), q0_phi, 1);
728 
729  //Amplitude
730  const EvtComplex A_phiks = D_phiks * H_phiks * L_phiks;
731 
733  /*const EvtVector4R p4ks_kpkm = Boost(p4ks, p4kpkm);
734  const double p4ks_kpkm_alt = p(mkpkm, (p4ks+p4kp+p4km).mass(), p4ks.mass());
735  const double p4kp_kpkm_alt = q(mkpkm, p4kp.mass(), p4km.mass());
736  //std::cout << p4ks_kpkm.d3mag() << std::endl;
737  //std::cout << p4ks_kpkm_alt << std::endl;
738  //std::cout << p4kp_kpkm.d3mag() << std::endl;
739  //std::cout << p4kp_kpkm_alt << std::endl;
740  const double costheta =
741  p4ks_kpkm.dot(p4kp_kpkm)/p4ks_kpkm.d3mag()/p4kp_kpkm.d3mag();
742  //const double z = p4ks_kpkm.d3mag()/(p4ks+p4kp+p4km).mass();
743  //const double H_phiks_alt = -sqrt(1.0+(z*z))*costheta;
744  const double sp_min =
745  s13_min( mkpkm*mkpkm, (p4ks+p4kp+p4km).mass(),
746  p4kp.mass(), p4km.mass(), p4ks.mass() );
747  const double sp_max =
748  s13_max( mkpkm*mkpkm, (p4ks+p4kp+p4km).mass(),
749  p4kp.mass(), p4km.mass(), p4ks.mass() );
750  const double sp = (p4ks+p4kp).mass2();
751  const double costheta_alt = (sp_max + sp_min - (2.0*sp))/(sp_max - sp_min);
752  //std::cout << costheta << std::endl;
753  //std::cout << costheta_alt << std::endl;
754  //exit(1);
755  const double z_alt = p4ks_kpkm_alt/(p4ks+p4kp+p4km).mass();
756  const double H_phiks_alt = -sqrt(1.0+(z_alt*z_alt))*costheta_alt;
757  //std::cout << H_phiks << std::endl;
758  //std::cout << H_phiks_alt << std::endl;
759  if( real(H_phiks) != H_phiks_alt )
760  {
761  std::cout << H_phiks << std::endl;
762  std::cout << H_phiks_alt << std::endl;
763  }
764  const EvtComplex A_phiks = D_phiks*H_phiks_alt*L_phiks;*/
766 
767  return A_phiks;
768  }
769 
770  EvtComplex EvtB0toKsKK::A_fxks(const EvtVector4R& /*p4ks*/,
771  const EvtVector4R& p4kp,
772  const EvtVector4R& p4km)
773  {
774  static EvtId fx = EvtPDL::getId("f_0(1500)");
775 
776  const double fx_m = EvtPDL::getMeanMass(fx);
777  const double fx_w = EvtPDL::getWidth(fx);
778 
779  const EvtVector4R p4kpkm = p4kp + p4km;
780  const double mkpkm = p4kpkm.mass();
781 
782  //Angular distribution
783  const EvtComplex H_fxks = 1.0;
784 
785  //Barrier factors
786  const EvtComplex D_fxks = 1.0;
787 
788  //Line shape
789  const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
790  const double q0_fx = q(fx_m, p4kp.mass(), p4km.mass());
791  const EvtComplex L_fxks =
792  BreitWigner(mkpkm, fx_m, fx_w, p4kp_kpkm.d3mag(), q0_fx, 0);
793 
794  //Amplitude
795  const EvtComplex A_fxks = D_fxks * H_fxks * L_fxks;
796 
797  return A_fxks;
798  }
799 
800  EvtComplex EvtB0toKsKK::A_chic0ks(const EvtVector4R& /*p4ks*/,
801  const EvtVector4R& p4kp,
802  const EvtVector4R& p4km)
803  {
804  static EvtId chic0 = EvtPDL::getId("chi_c0");
805 
806  const double chic0_m = EvtPDL::getMeanMass(chic0);
807  const double chic0_w = EvtPDL::getWidth(chic0);
808 
809  const EvtVector4R p4kpkm = p4kp + p4km;
810  const double mkpkm = p4kpkm.mass();
811 
812  //Angular distribution
813  const EvtComplex H_chic0ks = 1.0;
814 
815  //Barrier factors
816  const EvtComplex D_chic0ks = 1.0;
817 
818  //Line shape
819  const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
820  const double q0_chic0 = q(chic0_m, p4kp.mass(), p4km.mass());
821  const EvtComplex L_chic0ks =
822  BreitWigner(mkpkm, chic0_m, chic0_w, p4kp_kpkm.d3mag(), q0_chic0, 0);
823 
824  //Amplitude
825  const EvtComplex A_chic0ks = D_chic0ks * H_chic0ks * L_chic0ks;
826 
827  return A_chic0ks;
828  }
829 
830  EvtComplex EvtB0toKsKK::A_kknr(const EvtVector4R& p4k1,
831  const EvtVector4R& p4k2,
832  const double& alpha_kk)
833  {
834  const EvtVector4R p4kk = p4k1 + p4k2;
835  const double m2kk = p4kk.mass2();
836 
837  const EvtComplex A_kknr = exp(-alpha_kk * m2kk);
838 
839  return A_kknr;
840  }
841 
843 } // Belle 2 Namespace
Register Decay model EvtB0toKsKK.
Definition: EvtB0toKsKK.h:21
EvtComplex abar_kskpnr_
Variable member abar_kskpnr_
Definition: EvtB0toKsKK.h:100
double alpha_kskmnr
Variable member alpha_kskmnr.
Definition: EvtB0toKsKK.h:105
EvtComplex a_chic0ks_
Variable member a_chic0ks_.
Definition: EvtB0toKsKK.h:90
EvtComplex a_fxks_
Variable member a_fxks_
Definition: EvtB0toKsKK.h:89
std::ofstream debugfile_
debuging stream
Definition: EvtB0toKsKK.h:107
EvtComplex a_f0ks_
<Variable names for form factors
Definition: EvtB0toKsKK.h:87
double alpha_kpkmnr
Variable member alpha_kpkmnr.
Definition: EvtB0toKsKK.h:103
EvtComplex abar_fxks_
Variable member abar_fxks_
Definition: EvtB0toKsKK.h:97
EvtComplex abar_f0ks_
Variable member abar_f0ks_
Definition: EvtB0toKsKK.h:95
double alpha_kskpnr
Variable member alpha_kskpnr.
Definition: EvtB0toKsKK.h:104
EvtComplex abar_chic0ks_
Variable member abar_chic0ks_.
Definition: EvtB0toKsKK.h:98
EvtComplex a_kskmnr_
Variable member a_kskmnr_.
Definition: EvtB0toKsKK.h:93
EvtComplex abar_kskmnr_
Variable member abar_kskmnr_
Definition: EvtB0toKsKK.h:101
EvtComplex a_kpkmnr_
Variable member a_kpkmnr_.
Definition: EvtB0toKsKK.h:91
EvtComplex abar_phiks_
Variable member abar_phiks_.
Definition: EvtB0toKsKK.h:96
EvtComplex a_kskpnr_
Variable member a_kskpnr_.
Definition: EvtB0toKsKK.h:92
EvtComplex a_phiks_
Variable member a_phiks_
Definition: EvtB0toKsKK.h:88
EvtComplex abar_kpkmnr_
Variable member abar_kpkmnr_
Definition: EvtB0toKsKK.h:99
void init()
Initialize standard stream objects
Definition: EvtB0toKsKK.cc:41
EvtVector4R Smu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector Smu.
Definition: EvtB0toKsKK.cc:370
EvtComplex A_f0ks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_f0ks is amplitude of f0.
Definition: EvtB0toKsKK.cc:639
EvtComplex A_phiks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_phiks is amplitude of phi.
Definition: EvtB0toKsKK.cc:695
EvtVector4R umu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector umu.
Definition: EvtB0toKsKK.cc:360
EvtComplex Flatte(const double &m, const double &m0)
Constant Flatte.
Definition: EvtB0toKsKK.cc:617
double s13_max(const double &s12, const double &M, const double &m1, const double &m2, const double &m3)
maximum s13
Definition: EvtB0toKsKK.cc:681
double q(const double &mab, const double &ma, const double &mb)
Constants q.
Definition: EvtB0toKsKK.cc:594
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
EvtTensor4C gmunu_tilde(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor gmunu
Definition: EvtB0toKsKK.cc:409
EvtComplex BreitWigner(const double &m, const double &m0, const double &Gamma0, const double &q, const double &q0, const unsigned int &L)
BreitWigner Shape.
Definition: EvtB0toKsKK.cc:521
EvtTensor4C Mmunu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor Mmunu.
Definition: EvtB0toKsKK.cc:468
EvtTensor4C RaiseIndices(const EvtTensor4C &t)
Function RaiseIndices
Definition: EvtB0toKsKK.cc:453
void RaiseIndex(EvtVector4R &vector)
Member function RaiseIndices.
Definition: EvtB0toKsKK.cc:461
EvtDecayBase * clone()
Clone the decay of B0toKsKK.
Definition: EvtB0toKsKK.cc:36
EvtVector4R Lmu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector Lmu.
Definition: EvtB0toKsKK.cc:390
EvtComplex A_fxks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_fxks is amplitude of fxks.
Definition: EvtB0toKsKK.cc:770
double p(const double &mab, const double &M, const double &mc)
Constants p
Definition: EvtB0toKsKK.cc:583
EvtComplex Flatte_k(const double &s, const double &m_h)
Constant Flatte_k.
Definition: EvtB0toKsKK.cc:605
EvtComplex A_kknr(const EvtVector4R &p4k1, const EvtVector4R &p4k2, const double &alpha_kk)
A_kknr is amplitude of kknr.
Definition: EvtB0toKsKK.cc:830
void initProbMax()
Initialize standard stream objects for probability function
Definition: EvtB0toKsKK.cc:136
double BWBF(const double &q, const unsigned int &L)
Meson radius
Definition: EvtB0toKsKK.cc:479
EvtTensor4C Tmunu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor Tmunu
Definition: EvtB0toKsKK.cc:423
double s13_min(const double &s12, const double &M, const double &m1, const double &m2, const double &m3)
minimum s13
Definition: EvtB0toKsKK.cc:666
EvtVector4R Boost(const EvtVector4R &p4, const EvtVector4R &boost)
Parameter for boost frame
Definition: EvtB0toKsKK.cc:538
std::string getName()
Get function Name
Definition: EvtB0toKsKK.cc:31
void decay(EvtParticle *p)
Member of particle in EvtGen.
Definition: EvtB0toKsKK.cc:146
B2_EVTGEN_REGISTER_MODEL(EvtB0toKsKK)
register the model in EvtGen
EvtTensor4C Multiply(const EvtTensor4C &t1, const EvtTensor4C &t2)
Function Tensor Multiply
Definition: EvtB0toKsKK.cc:437
EvtComplex A_chic0ks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_chic0ks is amplitude of chic0ks.
Definition: EvtB0toKsKK.cc:800
Abstract base class for different kinds of events.