evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtbsToLLLLHyperCPAmp.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
21 #include "EvtGenModels/EvtbsToLLLLHyperCPAmp.hh"
22 
23 #include "EvtGenBase/EvtAmp.hh"
24 #include "EvtGenBase/EvtComplex.hh"
25 #include "EvtGenBase/EvtDiracSpinor.hh"
26 #include "EvtGenBase/EvtGenKine.hh"
27 #include "EvtGenBase/EvtId.hh"
28 #include "EvtGenBase/EvtIdSet.hh"
29 #include "EvtGenBase/EvtPDL.hh"
30 #include "EvtGenBase/EvtParticle.hh"
31 #include "EvtGenBase/EvtPatches.hh"
32 #include "EvtGenBase/EvtReport.hh"
33 #include "EvtGenBase/EvtScalarParticle.hh"
34 #include "EvtGenBase/EvtTensor4C.hh"
35 #include "EvtGenBase/EvtVector4C.hh"
36 
37 #include <cstdlib>
38 
39 // input: *parent - the pointer to the parent particle (B-meson, the
40 // object of the EvtParticle class);
41 // mS - the mass of the scalar sgoldstino "S" (GeV);
42 // mP - the mass of the pseudoscalar sgoldstino "P" (GeV);
43 // gammaS - the decay width of the scalar sgoldstino "S" (GeV);
44 // gammaP - the decay width of the pseudoscalar sgoldstino "P" (GeV);
45 // mLiiLR -
46 // Fc - coupling constant (GeV^2);
47 // mDijLL(RR) - parameters for \bar Bq-decays
48 // mDjiLL(RR) - parameters for Bq-decays (i <-> j!)
49 // d==1, s==2, b==3
50 //
51 void EvtbsToLLLLHyperCPAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp,
52  double mS, double mP, double gammaS,
53  double gammaP, double mLiiLR, double Fc,
54  double mD23LL, double mD23RR, double mD32LL,
55  double mD32RR, double mD13LL, double mD13RR,
56  double mD31LL, double mD31RR )
57 {
58  // FILE *mytest;
59 
60  int il1 = 0, il2 = 1, il3 = 2,
61  il4 = 3; // leptons are the first, second, thirds
62  // and fourth daughter particles
63 
64  EvtComplex unit1( 1.0, 0.0 ); // real unit
65  EvtComplex uniti( 0.0, 1.0 ); // imaginary unit
66 
67  parent->mass(); // B - meson mass, GeV
68  double fb = 0.0; // leptonic decay constant of B-meson, GeV
69 
70  double Cl = 0.0; // LPL and LSL - vertexes
71  if ( Fc != 0.0 ) {
72  Cl = mLiiLR * mLiiLR / ( sqrt( 2 ) * Fc );
73  }
74  if ( Cl == 0.0 ) {
75  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
76  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcAmp(...)"
77  << "\n Error in the Cl setting!"
78  << "\n Cl = " << Cl << "\n mLiiLR = " << mLiiLR
79  << "\n Fc = " << Fc << std::endl;
80  ::abort();
81  }
82 
83  EvtComplex MS = unit1 * mS -
84  uniti * gammaS /
85  2.0; // complex mass of the scalar sgoldstino
86  EvtComplex MP = unit1 * mP -
87  uniti * gammaP /
88  2.0; // complex mass of the pseudoscalar sgoldstino
89 
90  //
91  // Setting of the different Bq-mesons tipes
92  //
93  EvtId idparent = parent->getId(); // Bq-meson Id
94  EvtId IdMu1, IdMu2, IdMu3, IdMu4;
95 
96  double CB = 0.0;
97 
98  if ( idparent == EvtPDL::getId( std::string( "B_s0" ) ) ) {
99  fb = 0.24; // leptonic decay constant
100  CB = mD32LL * mD32LL + mD32RR * mD32RR;
101  }
102 
103  if ( idparent == EvtPDL::getId( std::string( "anti-B_s0" ) ) ) {
104  fb = 0.24; // leptonic decay constant
105  CB = mD23LL * mD23LL + mD23RR * mD23RR;
106  }
107 
108  if ( idparent == EvtPDL::getId( std::string( "B0" ) ) ) {
109  fb = 0.20; // leptonic decay constant
110  CB = mD31LL * mD31LL + mD31RR * mD31RR;
111  }
112 
113  if ( idparent == EvtPDL::getId( std::string( "anti-B0" ) ) ) {
114  fb = 0.20; // leptonic decay constant
115  CB = mD13LL * mD13LL + mD13RR * mD13RR;
116  }
117 
118  if ( CB == 0.0 ) {
119  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
120  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcAmp(...)"
121  << "\n Error in the CB setting!"
122  << "\n CB = " << CB << "\n mD32LL = " << mD32LL
123  << "\n mD32RR = " << mD32RR << "\n mD23LL = " << mD23LL
124  << "\n mD23RR = " << mD23RR << "\n mD31LL = " << mD31LL
125  << "\n mD31RR = " << mD31RR << "\n mD13LL = " << mD13LL
126  << "\n mD13RR = " << mD13RR << "\n idparent = " << idparent
127  << std::endl;
128  ::abort();
129  }
130 
131  //
132  // Setting the leptonic kinematical properties
133  //
134 
135  // to find charges of ell^+ and ell^- in the B-meson daughters
136  int charge1 = ( EvtPDL::chg3( parent->getDaug( il1 )->getId() ) ) / 3;
137  int charge2 = ( EvtPDL::chg3( parent->getDaug( il2 )->getId() ) ) / 3;
138  int charge3 = ( EvtPDL::chg3( parent->getDaug( il3 )->getId() ) ) / 3;
139  int charge4 = ( EvtPDL::chg3( parent->getDaug( il4 )->getId() ) ) / 3;
140  if ( ( abs( charge1 ) != 1 ) || ( abs( charge2 ) != 1 ) ||
141  ( abs( charge3 ) != 1 ) || ( abs( charge4 ) != 1 ) ||
142  ( charge1 + charge2 + charge3 + charge4 != 0 ) ) {
143  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
144  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcAmp(...)"
145  << "\n Error in the leptonic charge definition!"
146  << "\n charge1 =" << charge1
147  << "\n charge2 =" << charge2
148  << "\n charge3 =" << charge3
149  << "\n charge4 =" << charge4
150  << "\n number of daughters =" << parent->getNDaug() << std::endl;
151  ::abort();
152  }
153 
154  EvtParticle* lep1Plus = 0;
155  EvtParticle* lep1Minus = 0;
156  EvtParticle* lep2Plus = 0;
157  EvtParticle* lep2Minus = 0;
158 
159  EvtVector4R p; // B-meson momentum in the B-rest frame
160 
161  EvtVector4R q; // first transition 4-momentum in the B-rest frame
162  EvtVector4R k; // second transition 4-momentum in the B-rest frame
163  double q2; // Mandelstam variable s=q^2
164  double k2; // Mandelstam variable t=k^2
165 
166  EvtVector4R qsecond; // first transition 4-momentum in the B-rest frame
167  EvtVector4R ksecond; // second transition 4-momentum in the B-rest frame
168  double q2second; // Mandelstam variable s=q^2
169  double k2second; // Mandelstam variable t=k^2
170 
171  EvtVector4R k_1; // 4-momentum of ell^+ in the B-rest frame
172  EvtVector4R k_2; // 4-momentum of ell^- in the B-rest frame
173  EvtVector4R k_3; // 4-momentum of ell^+ in the B-rest frame
174  EvtVector4R k_4; // 4-momentum of ell^- in the B-rest frame
175 
176  k_1.set( 0.0, 0.0, 0.0, 0.0 );
177  k_2.set( 0.0, 0.0, 0.0, 0.0 );
178  k_3.set( 0.0, 0.0, 0.0, 0.0 );
179  k_4.set( 0.0, 0.0, 0.0, 0.0 );
180 
181  if ( ( charge1 + charge2 == 0 ) && ( charge3 + charge4 == 0 ) ) {
182  // positive charged lepton 1
183  lep1Plus = ( charge1 > charge2 ) ? parent->getDaug( il1 )
184  : parent->getDaug( il2 );
185  // negative charged lepton 1
186  lep1Minus = ( charge1 < charge2 ) ? parent->getDaug( il1 )
187  : parent->getDaug( il2 );
188  if ( charge1 > charge2 ) {
189  k_1 = parent->getDaug( il1 )->getP4();
190  k_2 = parent->getDaug( il2 )->getP4();
191  IdMu1 = parent->getDaug( il1 )->getId();
192  IdMu2 = parent->getDaug( il2 )->getId();
193  } else {
194  k_1 = parent->getDaug( il2 )->getP4();
195  k_2 = parent->getDaug( il1 )->getP4();
196  IdMu1 = parent->getDaug( il2 )->getId();
197  IdMu2 = parent->getDaug( il1 )->getId();
198  }
199  // positive charged lepton 2
200  lep2Plus = ( charge3 > charge4 ) ? parent->getDaug( il3 )
201  : parent->getDaug( il4 );
202  // negative charged lepton 2
203  lep2Minus = ( charge3 < charge4 ) ? parent->getDaug( il3 )
204  : parent->getDaug( il4 );
205  if ( charge3 > charge4 ) {
206  k_3 = parent->getDaug( il3 )->getP4();
207  k_4 = parent->getDaug( il4 )->getP4();
208  IdMu3 = parent->getDaug( il3 )->getId();
209  IdMu4 = parent->getDaug( il4 )->getId();
210  } else {
211  k_3 = parent->getDaug( il4 )->getP4();
212  k_4 = parent->getDaug( il3 )->getP4();
213  IdMu3 = parent->getDaug( il4 )->getId();
214  IdMu4 = parent->getDaug( il3 )->getId();
215  }
216  }
217  if ( ( charge1 + charge3 == 0 ) && ( charge2 + charge4 == 0 ) ) {
218  // positive charged lepton 1
219  lep1Plus = ( charge1 > charge3 ) ? parent->getDaug( il1 )
220  : parent->getDaug( il3 );
221  // negative charged lepton 1
222  lep1Minus = ( charge1 < charge3 ) ? parent->getDaug( il1 )
223  : parent->getDaug( il3 );
224  if ( charge1 > charge3 ) {
225  k_1 = parent->getDaug( il1 )->getP4();
226  k_2 = parent->getDaug( il3 )->getP4();
227  IdMu1 = parent->getDaug( il1 )->getId();
228  IdMu2 = parent->getDaug( il3 )->getId();
229  } else {
230  k_1 = parent->getDaug( il3 )->getP4();
231  k_2 = parent->getDaug( il1 )->getP4();
232  IdMu1 = parent->getDaug( il3 )->getId();
233  IdMu2 = parent->getDaug( il1 )->getId();
234  }
235  // positive charged lepton 2
236  lep2Plus = ( charge2 > charge4 ) ? parent->getDaug( il2 )
237  : parent->getDaug( il4 );
238  // negative charged lepton 2
239  lep2Minus = ( charge2 < charge4 ) ? parent->getDaug( il2 )
240  : parent->getDaug( il4 );
241  if ( charge2 > charge4 ) {
242  k_3 = parent->getDaug( il2 )->getP4();
243  k_4 = parent->getDaug( il4 )->getP4();
244  IdMu3 = parent->getDaug( il2 )->getId();
245  IdMu4 = parent->getDaug( il4 )->getId();
246  } else {
247  k_3 = parent->getDaug( il4 )->getP4();
248  k_4 = parent->getDaug( il2 )->getP4();
249  IdMu3 = parent->getDaug( il4 )->getId();
250  IdMu4 = parent->getDaug( il2 )->getId();
251  }
252  }
253 
254  p = parent->getP4Restframe();
255 
256  //
257  // The calculation of the FIRST part of the amplitude
258  //
259 
260  q = k_1 + k_2;
261  k = k_3 + k_4;
262  q2 = q.mass2(); // Mandelstam variable s=q^2
263  k2 = k.mass2(); // Mandelstam variable t=k^2
264 
265  //
266  // The calculation of the SECOND part of the amplitude
267  //
268 
269  qsecond = k_1 + k_4;
270  ksecond = k_3 + k_2;
271  q2second = qsecond.mass2(); // Mandelstam variable s=q^2
272  k2second = ksecond.mass2(); // Mandelstam variable t=k^2
273 
274  // B - and barB - mesons descriptors
275  static EvtIdSet bmesons( "anti-B0", "anti-B_s0" );
276  static EvtIdSet bbarmesons( "B0", "B_s0" );
277 
278  EvtId parentID = parent->getId();
279 
280  if ( bmesons.contains( parentID ) ) {
281  // The amplitude for the decay barB -> ell^+ ell^- ell^+ ell^- or
282  // b \bar q -> ell^+ ell^- ell^+ ell^-
283 
284  int i1, i2, i3, i4; // leptonic spin structures counters
285  int leptonicspin[4]; // array for the saving of the leptonic spin configuration
286 
287  // Tables for correspondings
288  // l^+(k_1) && lep1Plus && k_1 && i1
289  // l^-(k_2) && lep1Minus && k_2 && i2
290  // l^+(k_3) && lep2Plus && k_3 && i3
291  // l^-(k_4) && lep2Minus && k_4 && i4
292 
293  for ( i2 = 0; i2 < 2; i2++ ) {
294  leptonicspin[0] = i2;
295  for ( i1 = 0; i1 < 2; i1++ ) {
296  leptonicspin[1] = i1;
297  for ( i4 = 0; i4 < 2; i4++ ) {
298  leptonicspin[2] = i4;
299  for ( i3 = 0; i3 < 2; i3++ ) {
300  leptonicspin[3] = i3;
301 
302  EvtComplex SL2L1, PL4L3;
303  EvtComplex SL2L1second, PL4L3second;
304 
305  SL2L1 = EvtLeptonSCurrent( lep1Minus->spParent( i2 ),
306  lep1Plus->spParent( i1 ) );
307  PL4L3 = EvtLeptonPCurrent( lep2Minus->spParent( i4 ),
308  lep2Plus->spParent( i3 ) );
309 
310  SL2L1second = EvtLeptonSCurrent(
311  lep2Minus->spParent( i2 ), lep1Plus->spParent( i1 ) );
312  PL4L3second = EvtLeptonPCurrent(
313  lep1Minus->spParent( i4 ), lep2Plus->spParent( i3 ) );
314 
315  amp.vertex(
316  leptonicspin,
317  Cl * Cl * CB * fb *
318  ( SL2L1 * PL4L3 * ( q2 - k2 ) /
319  ( ( q2 - MS * MS ) * ( k2 - MP * MP ) ) -
320  SL2L1second * PL4L3second *
321  ( q2second - k2second ) /
322  ( ( q2second - MS * MS ) *
323  ( k2second - MP * MP ) ) ) /
324  ( 4.0 * Fc * Fc ) );
325  }
326  }
327  }
328  }
329 
330  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n The function EvtbsToLLLLHyperCPAmp::CalcAmp(...) passed with arguments:"
331  // << "\n ============================================================================"
332  // << "\n Input parameters:"
333  // << "\n mS = " << mS
334  // << "\n mP = " << mP
335  // << "\n gammaS = " << gammaS
336  // << "\n gammaP = " << gammaP
337  // << "\n mLiiLR = " << mLiiLR
338  // << "\n Fc = " << Fc
339  // << "\n mD23LL = " << mD23LL
340  // << "\n mD23RR = " << mD23RR
341  // << "\n mD32LL = " << mD32LL
342  // << "\n mD32RR = " << mD32RR
343  // << "\n mD13LL = " << mD13LL
344  // << "\n mD13RR = " << mD13RR
345  // << "\n mD31LL = " << mD31LL
346  // << "\n mD31RR = " << mD31RR
347  // << "\n ============================================================================"
348  // << "\n Kinematics:"
349  // << "\n k_1 = " << k_1
350  // << "\n m_ell_1 =" << parent->getDaug(il1)->mass()
351  // << "\n k_2 = " << k_2
352  // << "\n m_ell_2 =" << parent->getDaug(il2)->mass()
353  // << "\n k_3 = " << k_3
354  // << "\n m_ell_3 =" << parent->getDaug(il3)->mass()
355  // << "\n k_4 = " << k_4
356  // << "\n m_ell_4 =" << parent->getDaug(il4)->mass()
357  // << "\n p = " << p
358  // << "\n q = " << q
359  // << "\n k = " << k
360  // << "\n qsecond = " << qsecond
361  // << "\n ksecond = " << ksecond
362  // << "\n ============================================================================"
363  // << "\n Form-factors"
364  // << "\n fb = " << fb
365  // << "\n ============================================================================"
366  // << "\n Coefficients:"
367  // << "\n Cl = " << Cl
368  // << "\n CB = " << CB
369  // << "\n ============================================================================"
370  // << "\n Partical Properties:"
371  // << "\n IdB = " << idparent << " == " << EvtPDL::getId(std::string("anti-B_s0"))
372  // << "\n IdMu1 = " << IdMu1 << " == " << EvtPDL::getId(std::string("mu+"))
373  // << "\n IdMu2 = " << IdMu2 << " == " << EvtPDL::getId(std::string("mu-"))
374  // << "\n IdMu3 = " << IdMu3 << " == " << EvtPDL::getId(std::string("mu+"))
375  // << "\n IdMu4 = " << IdMu4 << " == " << EvtPDL::getId(std::string("mu-"))
376  // << "\n\n\n\n"
377  // << std::endl;
378 
379  } else {
380  if ( bbarmesons.contains( parentID ) ) {
381  // The amplitude for the decay B -> ell^+ ell^- ell^+ ell^- or
382  // q bar b -> ell^+ ell^- ell^+ ell^-
383 
384  int i1, i2, i3, i4; // leptonic spin structures counters
385  int leptonicspin[4]; // array for the saving of the leptonic spin configuration
386 
387  // Tables for correspondings
388  // l^+(k_1) && lep1Plus && k_1 && i1
389  // l^-(k_2) && lep1Minus && k_2 && i2
390  // l^+(k_3) && lep2Plus && k_3 && i3
391  // l^-(k_4) && lep2Minus && k_4 && i4
392 
393  for ( i2 = 1; i2 < 0; i2-- ) {
394  leptonicspin[0] = i2;
395  for ( i1 = 1; i1 < 0; i1-- ) {
396  leptonicspin[1] = i1;
397  for ( i4 = 1; i4 < 0; i4-- ) {
398  leptonicspin[2] = i4;
399  for ( i3 = 1; i3 < 0; i3-- ) {
400  leptonicspin[3] = i3;
401 
402  EvtComplex SL2L1, PL4L3;
403  EvtComplex SL2L1second, PL4L3second;
404 
405  SL2L1 = EvtLeptonSCurrent( lep1Minus->spParent( i2 ),
406  lep1Plus->spParent( i1 ) );
407  PL4L3 = EvtLeptonPCurrent( lep2Minus->spParent( i4 ),
408  lep2Plus->spParent( i3 ) );
409 
410  SL2L1second =
411  EvtLeptonSCurrent( lep2Minus->spParent( i2 ),
412  lep1Plus->spParent( i1 ) );
413  PL4L3second =
414  EvtLeptonPCurrent( lep1Minus->spParent( i4 ),
415  lep2Plus->spParent( i3 ) );
416 
417  amp.vertex( leptonicspin,
418  Cl * Cl * CB * fb *
419  ( SL2L1 * PL4L3 * ( q2 - k2 ) /
420  ( ( q2 - MS * MS ) *
421  ( k2 - MP * MP ) ) -
422  SL2L1second * PL4L3second *
423  ( q2second - k2second ) /
424  ( ( q2second - MS * MS ) *
425  ( k2second - MP * MP ) ) ) /
426  ( 4.0 * Fc * Fc ) );
427  }
428  }
429  }
430  }
431 
432  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n The function EvtbsToLLLLHyperCPAmp::CalcAmp(...) passed with arguments:"
433  // << "\n ============================================================================"
434  // << "\n Input parameters:"
435  // << "\n mS = " << mS
436  // << "\n mP = " << mP
437  // << "\n gammaS = " << gammaS
438  // << "\n gammaP = " << gammaP
439  // << "\n mLiiLR = " << mLiiLR
440  // << "\n Fc = " << Fc
441  // << "\n mD23LL = " << mD23LL
442  // << "\n mD23RR = " << mD23RR
443  // << "\n mD32LL = " << mD32LL
444  // << "\n mD32RR = " << mD32RR
445  // << "\n mD13LL = " << mD13LL
446  // << "\n mD13RR = " << mD13RR
447  // << "\n mD31LL = " << mD31LL
448  // << "\n mD31RR = " << mD31RR
449  // << "\n ============================================================================"
450  // << "\n Kinematics:"
451  // << "\n k_1 = " << k_1
452  // << "\n m_ell_1 =" << parent->getDaug(il1)->mass()
453  // << "\n k_2 = " << k_2
454  // << "\n m_ell_2 =" << parent->getDaug(il2)->mass()
455  // << "\n k_3 = " << k_3
456  // << "\n m_ell_3 =" << parent->getDaug(il3)->mass()
457  // << "\n k_4 = " << k_4
458  // << "\n m_ell_4 =" << parent->getDaug(il4)->mass()
459  // << "\n p = " << p
460  // << "\n q = " << q
461  // << "\n k = " << k
462  // << "\n qsecond = " << qsecond
463  // << "\n ksecond = " << ksecond
464  // << "\n ============================================================================"
465  // << "\n Form-factors"
466  // << "\n fb = " << fb
467  // << "\n ============================================================================"
468  // << "\n Coefficients:"
469  // << "\n Cl = " << Cl
470  // << "\n CB = " << CB
471  // << "\n ============================================================================"
472  // << "\n Partical Properties:"
473  // << "\n IdB = " << idparent << " == " << EvtPDL::getId(std::string("anti-B_s0"))
474  // << "\n IdMu1 = " << IdMu1 << " == " << EvtPDL::getId(std::string("mu+"))
475  // << "\n IdMu2 = " << IdMu2 << " == " << EvtPDL::getId(std::string("mu-"))
476  // << "\n IdMu3 = " << IdMu3 << " == " << EvtPDL::getId(std::string("mu+"))
477  // << "\n IdMu4 = " << IdMu4 << " == " << EvtPDL::getId(std::string("mu-"))
478  // << "\n\n\n\n"
479  // << std::endl;
480 
481  } else {
482  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
483  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcAmp(...)"
484  << "\n Wrong Bq-meson number" << std::endl;
485  ::abort();
486  }
487  }
488 }
489 
490 //
491 // The decays Bq -> ell^+ ell^- ell^+ ell^- maximum probability calculation
492 //
493 double EvtbsToLLLLHyperCPAmp::CalcMaxProb(
494  EvtId parnum, EvtId l1num, EvtId /*l2num*/, EvtId /*l3num*/,
495  EvtId /*l4num*/, double mS, double mP, double gammaS, double gammaP,
496  double mLiiLR, double Fc, double mD23LL, double mD23RR, double mD32LL,
497  double mD32RR, double mD13LL, double mD13RR, double mD31LL, double mD31RR )
498 {
499  if ( Fc == 0.0 ) {
500  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
501  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcMaxProb"
502  << "\n Error in the Fc setting!"
503  << "\n Fc = " << Fc << "\n mD32LL = " << mD32LL
504  << "\n mD32RR = " << mD32RR << "\n mD23LL = " << mD23LL
505  << "\n mD23RR = " << mD23RR << "\n mD31LL = " << mD31LL
506  << "\n mD31RR = " << mD31RR << "\n mD13LL = " << mD13LL
507  << "\n mD13RR = " << mD13RR << "\n parnum = " << parnum
508  << std::endl;
509  ::abort();
510  }
511 
512  double Cl = 0.0; // LPL and LSL - vertexes
513  if ( Fc != 0.0 ) {
514  Cl = mLiiLR * mLiiLR / ( sqrt( 2 ) * Fc );
515  }
516  if ( Cl == 0.0 ) {
517  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
518  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcMaxProb"
519  << "\n Error in the Cl setting!"
520  << "\n Cl = " << Cl << "\n mLiiLR = " << mLiiLR
521  << "\n Fc = " << Fc << std::endl;
522  ::abort();
523  }
524 
525  //
526  // Setting of the different Bq-mesons tipes
527  //
528 
529  double fb = 0.0;
530  double CB = 0.0;
531 
532  if ( parnum == EvtPDL::getId( std::string( "B_s0" ) ) ) {
533  fb = 0.24; // leptonic decay constant
534  CB = mD32LL * mD32LL + mD32RR * mD32RR;
535  }
536 
537  if ( parnum == EvtPDL::getId( std::string( "anti-B_s0" ) ) ) {
538  fb = 0.24; // leptonic decay constant
539  CB = mD23LL * mD23LL + mD23RR * mD23RR;
540  }
541 
542  if ( parnum == EvtPDL::getId( std::string( "B0" ) ) ) {
543  fb = 0.20; // leptonic decay constant
544  CB = mD31LL * mD31LL + mD31RR * mD31RR;
545  }
546 
547  if ( parnum == EvtPDL::getId( std::string( "anti-B0" ) ) ) {
548  fb = 0.20; // leptonic decay constant
549  CB = mD13LL * mD13LL + mD13RR * mD13RR;
550  }
551 
552  if ( CB == 0.0 ) {
553  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
554  << "\n\n The function EvtbsToLLLLHyperCPAmp::CalcMaxProb"
555  << "\n Error in the CB setting!"
556  << "\n CB = " << CB << "\n mD32LL = " << mD32LL
557  << "\n mD32RR = " << mD32RR << "\n mD23LL = " << mD23LL
558  << "\n mD23RR = " << mD23RR << "\n mD31LL = " << mD31LL
559  << "\n mD31RR = " << mD31RR << "\n mD13LL = " << mD13LL
560  << "\n mD13RR = " << mD13RR << "\n parnum = " << parnum
561  << std::endl;
562  ::abort();
563  }
564 
565  double M1 = EvtPDL::getMeanMass( parnum ); // B - meson mass
566  EvtPDL::getMeanMass( l1num ); // leptonic mass
567 
568  // We find the maximum amplitude probability
569  double maxfoundprob = Cl * Cl * CB * fb *
570  fabs( mS * mS + mP * mP + M1 * M1 ) * 10000000.0 /
571  ( 4.0 * Fc * Fc * mS * gammaS * mP * gammaP );
572 
573  if ( maxfoundprob <= 0.0 ) {
574  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
575  << "\n\n In the function EvtbsToLLLLHyperCPAmp::CalcMaxProb"
576  << "\n maxfoundprob = " << maxfoundprob << " < 0 or =0!"
577  << "\n mS = " << mS << "\n mP = " << mP
578  << "\n gammaS = " << gammaS << "\n gammaP = " << gammaP
579  << "\n mLiiLR = " << mLiiLR << "\n Fc = " << Fc
580  << "\n mD32LL = " << mD32LL << "\n mD32RR = " << mD32RR
581  << "\n mD23LL = " << mD23LL << "\n mD23RR = " << mD23RR
582  << "\n mD31LL = " << mD31LL << "\n mD31RR = " << mD31RR
583  << "\n mD13LL = " << mD13LL << "\n mD13RR = " << mD13RR
584  << "\n parnum = " << parnum << std::endl;
585  ::abort();
586  }
587 
588  EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
589  << "\n maxfoundprob (...) = " << maxfoundprob << std::endl;
590 
591  maxfoundprob *= 1.01;
592 
593  return maxfoundprob;
594 }
595 
596 // Triangular function
597 double EvtbsToLLLLHyperCPAmp::lambda( double a, double b, double c )
598 {
599  double l;
600 
601  l = pow( a, 2.0 ) + pow( b, 2.0 ) + pow( c, 2.0 ) - 2.0 * a * b -
602  2.0 * a * c - 2.0 * b * c;
603 
604  return l;
605 }
EvtbsToLLLLHyperCPAmp::CalcMaxProb
double CalcMaxProb(EvtId parnum, EvtId l1num, EvtId l2num, EvtId l3num, EvtId l4num, double mS, double mP, double gammaS, double gammaP, double mLiiLR, double Fc, double mD23LL, double mD23RR, double mD32LL, double mD32RR, double mD13LL, double mD13RR, double mD31LL, double mD31RR)
Definition: EvtbsToLLLLHyperCPAmp.cpp:493
EvtbsToLLLLHyperCPAmp::lambda
double lambda(double a, double b, double c)
Definition: EvtbsToLLLLHyperCPAmp.cpp:597
EvtLeptonPCurrent
EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:260
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.cpp:314
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cpp:130
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:167
EvtbsToLLLLHyperCPAmp::CalcAmp
void CalcAmp(EvtParticle *parent, EvtAmp &amp, double mS, double mP, double gammaS, double gammaP, double mLiiLR, double Fc, double mD23LL, double mD23RR, double mD32LL, double mD32RR, double mD13LL, double mD13RR, double mD31LL, double mD31RR)
Definition: EvtbsToLLLLHyperCPAmp.cpp:51
EvtPDL::chg3
static int chg3(EvtId i)
Definition: EvtPDL.cpp:372
EvtParticle::spParent
virtual EvtDiracSpinor spParent(int) const
Definition: EvtParticle.cpp:665
EvtId
Definition: EvtId.hh:27
EvtParticle::getNDaug
size_t getNDaug() const
Definition: EvtParticle.cpp:160
EvtAmp::vertex
void vertex(const EvtComplex &amp)
Definition: EvtAmp.cpp:454
EvtVector4R::mass2
double mass2() const
Definition: EvtVector4R.hh:100
EvtAmp
Definition: EvtAmp.hh:30
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:201
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cpp:150
EvtParticle::mass
double mass() const
Definition: EvtParticle.cpp:165
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtIdSet::contains
int contains(const EvtId id)
Definition: EvtIdSet.cpp:422
EvtLeptonSCurrent
EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:248
EvtParticle::getP4Restframe
EvtVector4R getP4Restframe() const
Definition: EvtParticle.cpp:793
Generated by   doxygen 1.8.15