EvtGen 2.2.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
Loading...
Searching...
No Matches
EvtBTo3hCP.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/EvtBTo3hCP.hh"
22
23#include "EvtGenBase/EvtCPUtil.hh"
24#include "EvtGenBase/EvtComplex.hh"
25#include "EvtGenBase/EvtGenKine.hh"
26#include "EvtGenBase/EvtId.hh"
27#include "EvtGenBase/EvtPDL.hh"
28#include "EvtGenBase/EvtParticle.hh"
29#include "EvtGenBase/EvtRandom.hh"
30#include "EvtGenBase/EvtReport.hh"
31#include "EvtGenBase/EvtVector4R.hh"
32
33#include <cmath>
34#include <stdlib.h>
35#include <string>
36
37#define square( x ) ( ( x ) * ( x ) )
38
39/*
40 * MK - 25/Aug/2016
41 * The code bellow is not necessarilly well writen as it is 1-to-1 rewrite
42 * from FORTRAN code (which is spread over 4 different source codes in not
43 * fully obvious way). Once it is rewriten and giving correct results, I
44 * will think how to give it more proper structure (mainly by checking what
45 * is duplicated and how to simplify it).
46 */
47
48void EvtBTo3hCP::setConstants( double balpha, double bbeta )
49{
50 m_alphaCP = balpha;
51 double calpha = cos( m_alphaCP );
52 double salpha = sin( m_alphaCP );
53 m_betaCP = bbeta;
54 double cbeta = cos( m_betaCP );
55 double sbeta = sin( m_betaCP );
56
57 m_MA2 = square( m_M_B ) + square( m_M_pip ) + square( m_M_pi0 ) +
58 square( m_M_pi0 );
59 m_MB2 = square( m_M_B ) + square( m_M_pip ) + square( m_M_pim ) +
60 square( m_M_pi0 );
61 m_MC2 = square( m_M_B ) + square( m_M_Kp ) + square( m_M_pim ) +
62 square( m_M_pi0 );
63
64 double StrongPhase = 0;
65 EvtComplex StrongExp( cos( StrongPhase ), sin( StrongPhase ) );
66
67 EvtComplex Mat_Tp0( calpha, -salpha );
68 Mat_Tp0 *= 1.09;
69 EvtComplex Mat_Tm0( calpha, -salpha );
70 Mat_Tm0 *= 1.09;
71 EvtComplex Mat_T0p( calpha, -salpha );
72 Mat_T0p *= 0.66;
73 EvtComplex Mat_T0m( calpha, -salpha );
74 Mat_T0m *= 0.66;
75 EvtComplex Mat_Tpm( calpha, -salpha );
76 Mat_Tpm *= 1.00;
77 EvtComplex Mat_Tmp( calpha, -salpha );
78 Mat_Tmp *= 0.47;
79
80 EvtComplex Mat_Ppm( cos( -0.5 ), sin( -0.5 ) );
81 Mat_Ppm *= -0.2;
82 EvtComplex Mat_Pmp( cos( 2. ), sin( 2. ) );
83 Mat_Pmp *= 0.15;
84
85 EvtComplex Mat_P1 = 0.5 * ( Mat_Ppm - Mat_Pmp );
86 EvtComplex Mat_P0 = 0.5 * ( Mat_Ppm + Mat_Pmp );
87 EvtComplex Nat_Tp0( calpha, salpha );
88 Nat_Tp0 *= 1.09;
89 EvtComplex Nat_Tm0( calpha, salpha );
90 Nat_Tm0 *= 1.09;
91 EvtComplex Nat_T0p( calpha, salpha );
92 Nat_T0p *= 0.66;
93 EvtComplex Nat_T0m( calpha, salpha );
94 Nat_T0m *= 0.66;
95 EvtComplex Nat_Tpm( calpha, salpha );
96 Nat_Tpm *= 1.00;
97 EvtComplex Nat_Tmp( calpha, salpha );
98 Nat_Tmp *= 0.47;
99 EvtComplex Nat_P1 = Mat_P1;
100 EvtComplex Nat_P0 = Mat_P0;
101
102 Mat_Tpm = StrongExp * Mat_Tpm;
103 Nat_Tpm = StrongExp * Nat_Tpm;
104
105 m_Mat_S1 = Mat_Tp0 + 2. * Mat_P1;
106 m_Mat_S2 = Mat_T0p - 2. * Mat_P1;
107 m_Mat_S3 = Mat_Tpm + Mat_P1 + Mat_P0;
108 m_Mat_S4 = Mat_Tmp - Mat_P1 + Mat_P0;
109 m_Mat_S5 = -Mat_Tpm - Mat_Tmp + Mat_Tp0 + Mat_T0p - 2. * Mat_P0;
110
111 m_Nat_S1 = Nat_Tp0 + 2. * Nat_P1;
112 m_Nat_S2 = Nat_T0p - 2. * Nat_P1;
113 m_Nat_S3 = Nat_Tpm + Nat_P1 + Nat_P0;
114 m_Nat_S4 = Nat_Tmp - Nat_P1 + Nat_P0;
115 m_Nat_S5 = -Nat_Tpm - Nat_Tmp + Nat_Tp0 + Nat_T0p - 2. * Nat_P0;
116
117 // B0 -->-- K*+ pi- Amplitudes (Trees + Penguins)
118 m_MatKstarp = EvtComplex( calpha, -salpha ) * EvtComplex( 0.220, 0. ) +
119 EvtComplex( cbeta, sbeta ) * EvtComplex( -1.200, 0. );
120 // B0 -->-- K*0 pi0 Amplitudes (Trees + Penguins)
121 m_MatKstar0 = EvtComplex( calpha, -salpha ) * EvtComplex( 0.015, 0. ) +
122 EvtComplex( cbeta, sbeta ) * EvtComplex( 0.850, 0. );
123 // B0 -->-- K+ rho- Amplitudes (Trees + Penguins)
124 m_MatKrho = EvtComplex( calpha, -salpha ) * EvtComplex( 0.130, 0. ) +
125 EvtComplex( cbeta, sbeta ) * EvtComplex( 0.160, 0. );
126 // B0bar -->-- K*+ pi- Amplitudes (Trees + Penguins)
127 m_NatKstarp = EvtComplex( 0., 0. );
128 // B0bar -->-- K*0 pi0 Amplitudes (Trees + Penguins)
129 m_NatKstar0 = EvtComplex( 0., 0. );
130 // B0bar -->-- K+ rho- Amplitudes (Trees + Penguins)
131 m_NatKrho = EvtComplex( 0., 0. );
132}
133
134void EvtBTo3hCP::Evt3pi( double alpha, int iset, EvtVector4R& p_pi_plus,
135 EvtVector4R& p_pi_minus, EvtVector4R& p_gamma_1,
136 EvtVector4R& p_gamma_2, double& Real_B0,
137 double& Imag_B0, double& Real_B0bar, double& Imag_B0bar )
138{
139 EvtVector4R p_p2;
140 double AB0, AB0bar, Ainter, R1, R2;
141 double factor;
142 int ierr = 0;
143
144 // Ghm : beta is not needed for this generation - put a default value
145 setConstants( alpha, 0.362 );
146
147 if ( iset == 0 ) {
148 p_pi_plus.set( m_M_pip, 0, 0, 0 );
149 p_p2.set( m_M_pi0, 0, 0, 0 );
150 p_pi_minus.set( m_M_pim, 0, 0, 0 );
151
152 do {
153 firstStep( p_pi_plus, p_p2, p_pi_minus, 1 );
154 ierr = compute3pi( p_pi_plus, p_p2, p_pi_minus, Real_B0, Imag_B0,
155 Real_B0bar, Imag_B0bar, iset );
156 } while ( ierr != 0 );
157 } else if ( iset < 0 ) {
158 p_p2 = p_gamma_1 + p_gamma_2;
159 ierr = compute3pi( p_pi_plus, p_p2, p_pi_minus, Real_B0, Imag_B0,
160 Real_B0bar, Imag_B0bar, iset );
161 if ( ierr != 0 ) {
162 std::cout << "Provided kinematics is not physical\n";
163 std::cout << "Program will stop\n";
164 exit( 1 );
165 }
166 } else // iset > 0
167 {
168 m_factor_max = 0;
169
170 int endLoop = iset;
171 for ( int i = 0; i < endLoop; ++i ) {
172 p_pi_plus.set( m_M_pip, 0, 0, 0 );
173 p_p2.set( m_M_pi0, 0, 0, 0 );
174 p_pi_minus.set( m_M_pim, 0, 0, 0 );
175
176 firstStep( p_pi_plus, p_p2, p_pi_minus, 1 );
177 ierr = compute3pi( p_pi_plus, p_p2, p_pi_minus, Real_B0, Imag_B0,
178 Real_B0bar, Imag_B0bar, iset );
179 if ( ierr != 0 ) {
180 continue;
181 }
182 AB0 = square( Real_B0 ) + square( Imag_B0 );
183 AB0bar = square( Real_B0bar ) + square( Imag_B0bar );
184 Ainter = Real_B0 * Imag_B0bar - Imag_B0 * Real_B0bar;
185 R1 = ( AB0 - AB0bar ) / ( AB0 + AB0bar );
186 R2 = ( 2.0 * Ainter ) / ( AB0 + AB0bar );
187 factor = ( 1.0 + sqrt( square( R1 ) + square( R2 ) ) ) *
188 ( AB0 + AB0bar ) / 2.0;
189 if ( factor > m_factor_max )
190 m_factor_max = factor;
191 }
192 m_factor_max = 1.0 / std::sqrt( m_factor_max );
193 }
194
195 Real_B0 *= m_factor_max;
196 Imag_B0 *= m_factor_max;
197 Real_B0bar *= m_factor_max;
198 Imag_B0bar *= m_factor_max;
199
200 if ( iset < 0 ) {
201 return;
202 }
203
204 rotation( p_pi_plus, 1 );
205 rotation( p_p2, 0 );
206 rotation( p_pi_minus, 0 );
207
208 gammaGamma( p_p2, p_gamma_1, p_gamma_2 );
209}
210
211void EvtBTo3hCP::Evt3piMPP( double alpha, int iset, EvtVector4R& p_p1,
212 EvtVector4R& p_p2, EvtVector4R& p_p3,
213 double& Real_B0, double& Imag_B0,
214 double& Real_B0bar, double& Imag_B0bar )
215{
216 double ABp, ABm;
217 int ierr = 0;
218
219 // Ghm : beta is not needed for this generation - put a default value
220 setConstants( alpha, 0.362 );
221
222 if ( iset == 0 ) {
223 p_p1.set( m_M_pim, 0, 0, 0 );
224 p_p2.set( m_M_pip, 0, 0, 0 );
225 p_p3.set( m_M_pip, 0, 0, 0 );
226
227 do {
228 firstStep( p_p1, p_p2, p_p3, 2 );
229 ierr = compute3piMPP( p_p1, p_p2, p_p3, Real_B0, Imag_B0,
230 Real_B0bar, Imag_B0bar, iset );
231 } while ( ierr != 0 );
232 } else if ( iset < 0 ) {
233 ierr = compute3piMPP( p_p1, p_p2, p_p3, Real_B0, Imag_B0, Real_B0bar,
234 Imag_B0bar, iset );
235 if ( ierr != 0 ) {
236 std::cout << "Provided kinematics is not physical\n";
237 std::cout << "Program will stop\n";
238 exit( 1 );
239 }
240 } else // iset > 0
241 {
242 m_factor_max = 0;
243
244 int endLoop = iset;
245 for ( int i = 0; i < endLoop; ++i ) {
246 p_p1.set( m_M_pim, 0, 0, 0 );
247 p_p2.set( m_M_pip, 0, 0, 0 );
248 p_p3.set( m_M_pip, 0, 0, 0 );
249
250 firstStep( p_p1, p_p2, p_p3, 2 );
251 ierr = compute3piMPP( p_p1, p_p2, p_p3, Real_B0, Imag_B0,
252 Real_B0bar, Imag_B0bar, iset );
253 if ( ierr != 0 ) {
254 continue;
255 }
256 ABp = square( Real_B0 ) + square( Imag_B0 );
257 ABm = square( Real_B0bar ) + square( Imag_B0bar );
258 if ( ABp > m_factor_max )
259 m_factor_max = ABp;
260 if ( ABm > m_factor_max )
261 m_factor_max = ABm;
262 }
263 m_factor_max = 1.0 / std::sqrt( m_factor_max );
264 }
265
266 Real_B0 *= m_factor_max;
267 Imag_B0 *= m_factor_max;
268 Real_B0bar *= m_factor_max;
269 Imag_B0bar *= m_factor_max;
270
271 if ( iset < 0 ) {
272 return;
273 }
274
275 rotation( p_p1, 1 );
276 rotation( p_p2, 0 );
277 rotation( p_p3, 0 );
278}
279
280void EvtBTo3hCP::Evt3piP00( double alpha, int iset, EvtVector4R& p_p1,
281 EvtVector4R& p_p1_gamma1, EvtVector4R& p_p1_gamma2,
282 EvtVector4R& p_p2_gamma1, EvtVector4R& p_p2_gamma2,
283 double& Real_B0, double& Imag_B0,
284 double& Real_B0bar, double& Imag_B0bar )
285{
286 double ABp, ABm;
287 EvtVector4R p_p2, p_p3;
288 int ierr = 0;
289
290 // Ghm : beta is not needed for this generation - put a default value
291 setConstants( alpha, 0.362 );
292
293 if ( iset == 0 ) {
294 p_p1.set( m_M_pip, 0, 0, 0 );
295 p_p2.set( m_M_pi0, 0, 0, 0 );
296 p_p3.set( m_M_pi0, 0, 0, 0 );
297
298 do {
299 firstStep( p_p1, p_p2, p_p3, 3 );
300 ierr = compute3piP00( p_p1, p_p2, p_p3, Real_B0, Imag_B0,
301 Real_B0bar, Imag_B0bar, iset );
302 } while ( ierr != 0 );
303 } else if ( iset < 0 ) {
304 p_p2 = p_p1_gamma1 + p_p1_gamma2;
305 p_p3 = p_p2_gamma1 + p_p2_gamma2;
306 ierr = compute3piP00( p_p1, p_p2, p_p3, Real_B0, Imag_B0, Real_B0bar,
307 Imag_B0bar, iset );
308 if ( ierr != 0 ) {
309 std::cout << "Provided kinematics is not physical\n";
310 std::cout << "Program will stop\n";
311 exit( 1 );
312 }
313 } else // iset > 0
314 {
315 m_factor_max = 0;
316
317 int endLoop = iset;
318 for ( int i = 0; i < endLoop; ++i ) {
319 p_p1.set( m_M_pip, 0, 0, 0 );
320 p_p2.set( m_M_pi0, 0, 0, 0 );
321 p_p3.set( m_M_pi0, 0, 0, 0 );
322
323 firstStep( p_p1, p_p2, p_p3, 3 );
324 ierr = compute3piP00( p_p1, p_p2, p_p3, Real_B0, Imag_B0,
325 Real_B0bar, Imag_B0bar, iset );
326 if ( ierr != 0 ) {
327 continue;
328 }
329 ABp = square( Real_B0 ) + square( Imag_B0 );
330 ABm = square( Real_B0bar ) + square( Imag_B0bar );
331 if ( ABp > m_factor_max )
332 m_factor_max = ABp;
333 if ( ABm > m_factor_max )
334 m_factor_max = ABm;
335 }
336 m_factor_max = 1.0 / std::sqrt( m_factor_max );
337 }
338
339 Real_B0 *= m_factor_max;
340 Imag_B0 *= m_factor_max;
341 Real_B0bar *= m_factor_max;
342 Imag_B0bar *= m_factor_max;
343
344 if ( iset < 0 ) {
345 return;
346 }
347
348 rotation( p_p1, 1 );
349 rotation( p_p2, 0 );
350 rotation( p_p3, 0 );
351
352 gammaGamma( p_p2, p_p1_gamma1, p_p1_gamma2 );
353 gammaGamma( p_p3, p_p2_gamma1, p_p2_gamma2 );
354}
355
356void EvtBTo3hCP::EvtKpipi( double alpha, double beta, int iset,
357 EvtVector4R& p_K_plus, EvtVector4R& p_pi_minus,
358 EvtVector4R& p_gamma_1, EvtVector4R& p_gamma_2,
359 double& Real_B0, double& Imag_B0, double& Real_B0bar,
360 double& Imag_B0bar )
361{
362 EvtVector4R p_p3;
363 double ABp, ABm;
364 int ierr = 0;
365
366 setConstants( alpha, beta );
367
368 if ( iset == 0 ) {
369 p_K_plus.set( m_M_Kp, 0, 0, 0 );
370 p_pi_minus.set( m_M_pim, 0, 0, 0 );
371 p_p3.set( m_M_pi0, 0, 0, 0 );
372
373 do {
374 firstStep( p_K_plus, p_pi_minus, p_p3, 0 );
375 ierr = computeKpipi( p_K_plus, p_pi_minus, p_p3, Real_B0, Imag_B0,
376 Real_B0bar, Imag_B0bar, iset );
377 } while ( ierr != 0 );
378 } else if ( iset < 0 ) {
379 p_p3 = p_gamma_1 + p_gamma_2;
380 ierr = computeKpipi( p_K_plus, p_pi_minus, p_p3, Real_B0, Imag_B0,
381 Real_B0bar, Imag_B0bar, iset );
382 if ( ierr != 0 ) {
383 std::cout << "Provided kinematics is not physical\n";
384 std::cout << "Program will stop\n";
385 exit( 1 );
386 }
387 } else // iset > 0
388 {
389 m_factor_max = 0;
390
391 int endLoop = iset;
392 for ( int i = 0; i < endLoop; ++i ) {
393 p_K_plus.set( m_M_Kp, 0, 0, 0 );
394 p_pi_minus.set( m_M_pim, 0, 0, 0 );
395 p_p3.set( m_M_pi0, 0, 0, 0 );
396 firstStep( p_K_plus, p_pi_minus, p_p3, 0 );
397 ierr = computeKpipi( p_K_plus, p_pi_minus, p_p3, Real_B0, Imag_B0,
398 Real_B0bar, Imag_B0bar, iset );
399 if ( ierr != 0 ) {
400 continue;
401 }
402 ABp = square( Real_B0 ) + square( Imag_B0 );
403 ABm = square( Real_B0bar ) + square( Imag_B0bar );
404
405 if ( ABp > m_factor_max ) {
406 m_factor_max = ABp;
407 }
408 if ( ABm > m_factor_max ) {
409 m_factor_max = ABm;
410 }
411 }
412 m_factor_max = 1.0 / std::sqrt( m_factor_max );
413 }
414
415 Real_B0 *= m_factor_max;
416 Imag_B0 *= m_factor_max;
417 Real_B0bar *= m_factor_max;
418 Imag_B0bar *= m_factor_max;
419
420 if ( iset < 0 ) {
421 return;
422 }
423
424 rotation( p_K_plus, 1 );
425 rotation( p_pi_minus, 0 );
426 rotation( p_p3, 0 );
427
428 gammaGamma( p_p3, p_gamma_1, p_gamma_2 );
429}
430
431void EvtBTo3hCP::firstStep( EvtVector4R& p1, EvtVector4R& p2, EvtVector4R& p3,
432 int mode )
433{
434 const double m1sq = p1.mass2();
435 const double m2sq = p2.mass2();
436 const double m3sq = p3.mass2();
437 double min_m12, min_m13, min_m23;
438
439 double max_m12 = square( m_M_B );
440 double max_m13 = square( m_M_B );
441 double max_m23 = square( m_M_B );
442
443 if ( mode == 0 ) {
444 min_m12 = m1sq + m2sq + 2 * sqrt( m1sq * m2sq );
445 min_m13 = m1sq + m3sq + 2 * sqrt( m1sq * m3sq );
446 min_m23 = m2sq + m3sq + 2 * sqrt( m2sq * m3sq );
447 } else {
448 min_m12 = m1sq + m2sq;
449 min_m13 = m1sq + m3sq;
450 min_m23 = m2sq + m3sq;
451 }
452
453 bool eventOK;
454 double m13, m12, m23;
455 double E1;
456 double E2;
457 double E3;
458 double p1mom;
459 double p2mom;
460 double p3mom;
461 double cost13;
462 double cost12;
463 double cost23;
464 eventOK = false;
465
466 do {
467 switch ( mode ) {
468 case 0:
469 generateSqMasses_Kpipi( m12, m13, m23, m_MC2, m1sq, m2sq, m3sq );
470 break;
471 case 1:
472 generateSqMasses_3pi( m12, m13, m23, m_MB2, m1sq, m2sq, m3sq );
473 break;
474 case 2:
475 generateSqMasses_3piMPP( m12, m13, m23, m_MB2, m1sq, m2sq, m3sq );
476 break;
477 case 3:
478 generateSqMasses_3piP00( m12, m13, m23, m_MA2, m1sq, m2sq, m3sq );
479 break;
480 default:
481 break;
482 }
483 // Check whether event is physical
484 if ( ( m23 < min_m23 ) || ( m23 > max_m23 ) )
485 continue;
486 if ( ( m13 < min_m13 ) || ( m13 > max_m13 ) )
487 continue;
488 if ( ( m12 < min_m12 ) || ( m12 > max_m12 ) )
489 continue;
490
491 // Now check the cosines of the angles
492 E1 = ( square( m_M_B ) + m1sq - m23 ) / ( 2. * m_M_B );
493 E2 = ( square( m_M_B ) + m2sq - m13 ) / ( 2. * m_M_B );
494 E3 = ( square( m_M_B ) + m3sq - m12 ) / ( 2. * m_M_B );
495 p1mom = square( E1 ) - m1sq;
496 p2mom = square( E2 ) - m2sq;
497 p3mom = square( E3 ) - m3sq;
498 if ( p1mom < 0 || p2mom < 0 || p3mom < 0 ) {
499 // std::cout<<"Momenta magnitude negative\n";
500 continue;
501 }
502 p1mom = sqrt( p1mom );
503 p2mom = sqrt( p2mom );
504 p3mom = sqrt( p3mom );
505
506 cost13 = ( 2. * E1 * E3 + m1sq + m3sq - m13 ) / ( 2. * p1mom * p3mom );
507 cost12 = ( 2. * E1 * E2 + m1sq + m2sq - m12 ) / ( 2. * p1mom * p2mom );
508 cost23 = ( 2. * E2 * E3 + m2sq + m3sq - m23 ) / ( 2. * p2mom * p3mom );
509 if ( cost13 < -1. || cost13 > 1. || cost12 < -1. || cost12 > 1. ||
510 cost23 < -1. || cost23 > 1. ) {
511 continue;
512 }
513 eventOK = true;
514 } while ( eventOK == false );
515
516 // Now is time to fill 4-vectors
517 p3.set( E3, 0, 0, p3mom );
518 p1.set( E1, p1mom * sqrt( 1 - square( cost13 ) ), 0, p1mom * cost13 );
519 p2.set( 0, E2 );
520 for ( int i = 1; i < 4; ++i ) {
521 p2.set( i, -p1.get( i ) - p3.get( i ) );
522 }
523 if ( p1.get( 0 ) < p1.d3mag() ) {
524 std::cout << "Unphysical p1 generated: " << p1 << std::endl;
525 }
526 if ( p2.get( 0 ) < p2.d3mag() ) {
527 std::cout << "Unphysical p2 generated: " << p2 << std::endl;
528 }
529 if ( p3.get( 0 ) < p3.d3mag() ) {
530 std::cout << "Unphysical p3 generated: " << p3 << std::endl;
531 }
532 double testMB2 = m_MB2;
533 switch ( mode ) {
534 case 0:
535 testMB2 = m_MC2;
536 break;
537 case 1:
538 case 2:
539 testMB2 = m_MB2;
540 break;
541 case 3:
542 testMB2 = m_MA2;
543 break;
544 }
545
546 if ( fabs( m12 + m13 + m23 - testMB2 ) > 1e-4 ) {
547 std::cout << "Unphysical event generated: " << m12 << " " << m13 << " "
548 << m23 << std::endl;
549 }
550}
551
552void EvtBTo3hCP::generateSqMasses_Kpipi( double& m12, double& m13, double& m23,
553 double MB2, double m1sq, double m2sq,
554 double m3sq )
555{
556 /*
557 C There is two ways of generating the events:
558 C The first one used a pole-compensation method to generate the
559 C events efficiently taking into account the poles due to the
560 C Breit-Wigners of the rho s. It is activated by setting
561 C Phase_Space to .false.
562 C The second one generates events according to phase space. It is
563 C inneficient but allows the exploration of the full Dalitz plot
564 C in an uniform way. It was found to be usefull fopr some peculiar
565 C applications. It is activated by setting
566 C Phase_space to .true.
567 C Note that in that case, the generation is no longer correct.
568 */
569 static const bool phaseSpace = false;
570
571 double max_m12 = square( m_M_B );
572 double min_m12 = m1sq + m2sq + 2 * sqrt( m1sq * m2sq );
573
574 double max_m13 = square( m_M_B );
575 double min_m13 = m1sq + m3sq + 2 * sqrt( m1sq * m3sq );
576
577 double max_m23 = square( m_M_B );
578 double min_m23 = m2sq + m3sq + 2 * sqrt( m2sq * m3sq );
579
580 double z = 3. * EvtRandom::Flat();
581 if ( z < 1. ) // K*+
582 {
583 if ( phaseSpace ) {
584 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
585 } else {
586 double y = EvtRandom::Flat() * m_pi - m_pi / 2;
587 double x = std::tan( y );
588 double mass = x * m_Gam_Kstarp / 2. + m_Mass_Kstarp;
589 m13 = square( mass );
590 }
591 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
592 m23 = MB2 - m12 - m13;
593 } else if ( z < 2. ) // K*0
594 {
595 if ( phaseSpace ) {
596 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
597 } else {
598 double y = EvtRandom::Flat() * m_pi - m_pi / 2;
599 double x = std::tan( y );
600 double mass = x * m_Gam_Kstar0 / 2. + m_Mass_Kstar0;
601 m12 = square( mass );
602 }
603 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
604 m23 = MB2 - m12 - m13;
605 } else // rho-
606 {
607 if ( phaseSpace ) {
608 m23 = EvtRandom::Flat() * ( max_m23 - min_m23 ) + min_m23;
609 } else {
610 double y = EvtRandom::Flat() * m_pi - m_pi / 2;
611 double x = std::tan( y );
612 double mass = x * m_Gam_rho / 2. + m_Mass_rho;
613 m23 = square( mass );
614 }
615 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
616 m12 = MB2 - m23 - m13;
617 }
618}
619
620void EvtBTo3hCP::generateSqMasses_3pi( double& m12, double& m13, double& m23,
621 double MB2, double m1sq, double m2sq,
622 double m3sq )
623{
624 /*
625 C There is two ways of generating the events:
626 C The first one used a pole-compensation method to generate the
627 C events efficiently taking into account the poles due to the
628 C Breit-Wigners of the rho s. It is activated by setting
629 C Phase_Space to .false.
630 C The second one generates events according to phase space. It is
631 C inneficient but allows the exploration of the full Dalitz plot
632 C in an uniform way. It was found to be usefull fopr some peculiar
633 C applications. It is activated by setting
634 C Phase_space to .true.
635 C Note that in that case, the generation is no longer correct.
636 */
637 static const bool phaseSpace = false;
638
639 double max_m12 = square( m_M_B );
640 double min_m12 = m1sq + m2sq;
641
642 double max_m13 = square( m_M_B );
643 double min_m13 = m1sq + m3sq;
644
645 double max_m23 = square( m_M_B );
646 double min_m23 = m2sq + m3sq;
647 double mass = 0;
648
649 if ( !phaseSpace ) {
650 double y = EvtRandom::Flat() * m_pi - m_pi / 2;
651 double x = std::tan( y );
652 mass = x * m_Gam_rho / 2. + m_Mass_rho;
653 }
654
655 double z = 3. * EvtRandom::Flat();
656 if ( z < 1. ) {
657 if ( phaseSpace ) {
658 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
659 } else {
660 m12 = square( mass );
661 }
662 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
663 m23 = MB2 - m12 - m13;
664 } else if ( z < 2. ) {
665 if ( phaseSpace ) {
666 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
667 } else {
668 m13 = square( mass );
669 }
670 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
671 m23 = MB2 - m12 - m13;
672 } else {
673 if ( phaseSpace ) {
674 m23 = EvtRandom::Flat() * ( max_m23 - min_m23 ) + min_m23;
675 } else {
676 m23 = square( mass );
677 }
678 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
679 m13 = MB2 - m12 - m23;
680 }
681}
682
683void EvtBTo3hCP::generateSqMasses_3piMPP( double& m12, double& m13, double& m23,
684 double MB2, double m1sq, double m2sq,
685 double m3sq )
686{
687 /*
688 C There is two ways of generating the events:
689 C The first one used a pole-compensation method to generate the
690 C events efficiently taking into account the poles due to the
691 C Breit-Wigners of the rho s. It is activated by setting
692 C Phase_Space to .false.
693 C The second one generates events according to phase space. It is
694 C inneficient but allows the exploration of the full Dalitz plot
695 C in an uniform way. It was found to be usefull fopr some peculiar
696 C applications. It is activated by setting
697 C Phase_space to .true.
698 C Note that in that case, the generation is no longer correct.
699 */
700 static const bool phaseSpace = false;
701
702 double max_m12 = square( m_M_B );
703 double min_m12 = m1sq + m2sq;
704
705 double max_m13 = square( m_M_B );
706 double min_m13 = m1sq + m3sq;
707
708 double mass = 0;
709
710 if ( !phaseSpace ) {
711 double y = EvtRandom::Flat() * m_pi - m_pi / 2;
712 double x = std::tan( y );
713 mass = x * m_Gam_rho / 2. + m_Mass_rho;
714 }
715
716 double z = EvtRandom::Flat();
717 if ( z < 0.5 ) {
718 if ( phaseSpace ) {
719 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
720 } else {
721 m12 = square( mass );
722 }
723 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
724 m23 = MB2 - m12 - m13;
725 } else {
726 if ( phaseSpace ) {
727 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
728 } else {
729 m13 = square( mass );
730 }
731 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
732 m23 = MB2 - m12 - m13;
733 }
734}
735
736void EvtBTo3hCP::generateSqMasses_3piP00( double& m12, double& m13, double& m23,
737 double MB2, double m1sq, double m2sq,
738 double m3sq )
739{
740 /*
741 C There is two ways of generating the events:
742 C The first one used a pole-compensation method to generate the
743 C events efficiently taking into account the poles due to the
744 C Breit-Wigners of the rho s. It is activated by setting
745 C Phase_Space to .false.
746 C The second one generates events according to phase space. It is
747 C inneficient but allows the exploration of the full Dalitz plot
748 C in an uniform way. It was found to be usefull fopr some peculiar
749 C applications. It is activated by setting
750 C Phase_space to .true.
751 C Note that in that case, the generation is no longer correct.
752 */
753 static const bool phaseSpace = false;
754
755 double max_m12 = square( m_M_B );
756 double min_m12 = m1sq + m2sq;
757
758 double max_m13 = square( m_M_B );
759 double min_m13 = m1sq + m3sq;
760
761 double mass = 0;
762
763 if ( !phaseSpace ) {
764 double y = EvtRandom::Flat() * m_pi - m_pi / 2;
765 double x = std::tan( y );
766 mass = x * m_Gam_rho / 2. + m_Mass_rho;
767 }
768
769 double z = EvtRandom::Flat();
770 if ( z < 0.5 ) {
771 if ( phaseSpace ) {
772 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
773 } else {
774 m12 = square( mass );
775 }
776 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
777 m23 = MB2 - m12 - m13;
778 } else {
779 if ( phaseSpace ) {
780 m13 = EvtRandom::Flat() * ( max_m13 - min_m13 ) + min_m13;
781 } else {
782 m13 = square( mass );
783 }
784 m12 = EvtRandom::Flat() * ( max_m12 - min_m12 ) + min_m12;
785 m23 = MB2 - m12 - m13;
786 }
787}
788int EvtBTo3hCP::compute3pi( EvtVector4R& p1, EvtVector4R& p2, EvtVector4R& p3,
789 double& real_B0, double& imag_B0,
790 double& real_B0bar, double& imag_B0bar, int iset )
791{
792 int ierr = 0;
793
794 double m12 = ( p1 + p2 ).mass();
795 double m13 = ( p1 + p3 ).mass();
796 double m23 = ( p2 + p3 ).mass();
797
798 double W12 =
799 1. / ( ( square( m_Mass_rho - m12 ) + square( m_Gam_rho / 2. ) ) * m12 );
800 double W13 =
801 1. / ( ( square( m_Mass_rho - m13 ) + square( m_Gam_rho / 2. ) ) * m13 );
802 double W23 =
803 1. / ( ( square( m_Mass_rho - m23 ) + square( m_Gam_rho / 2. ) ) * m23 );
804
805 double Wtot = 1.;
806 if ( iset >= 0 ) {
807 Wtot = 1. / sqrt( W12 + W13 + W23 );
808 }
809
810 EvtComplex Mat_rhop = BreitWigner( p1, p2, p3, ierr );
811 EvtComplex Mat_rhom = BreitWigner( p2, p3, p1, ierr );
812 EvtComplex Mat_rho0 = BreitWigner( p1, p3, p2, ierr );
813
814 EvtComplex Mat_1 = m_Mat_S3 * Mat_rhop;
815 EvtComplex Mat_2 = m_Mat_S4 * Mat_rhom;
816 EvtComplex Mat_3 = m_Mat_S5 * Mat_rho0 * 0.5;
817
818 EvtComplex MatBp = ( Mat_1 + Mat_2 + Mat_3 ) * Wtot;
819
820 Mat_1 = m_Nat_S3 * Mat_rhom;
821 Mat_2 = m_Nat_S4 * Mat_rhop;
822 Mat_3 = m_Nat_S5 * Mat_rho0 * 0.5;
823
824 EvtComplex MatBm = ( Mat_1 + Mat_2 + Mat_3 ) * Wtot;
825
826 real_B0 = real( MatBp );
827 imag_B0 = imag( MatBp );
828
829 real_B0bar = real( MatBm );
830 imag_B0bar = imag( MatBm );
831
832 return ierr;
833}
834
835int EvtBTo3hCP::compute3piMPP( EvtVector4R& p1, EvtVector4R& p2,
836 EvtVector4R& p3, double& real_B0, double& imag_B0,
837 double& real_B0bar, double& imag_B0bar, int iset )
838{
839 int ierr = 0;
840 const double ASHQ = sqrt( 2. );
841 double m12 = ( p1 + p2 ).mass();
842 double m13 = ( p1 + p3 ).mass();
843
844 double W12 =
845 1. / ( ( square( m_Mass_rho - m12 ) + square( m_Gam_rho / 2. ) ) * m12 );
846 double W13 =
847 1. / ( ( square( m_Mass_rho - m13 ) + square( m_Gam_rho / 2. ) ) * m13 );
848
849 double Wtot = 1.;
850 if ( iset >= 0 ) {
851 Wtot = 1. / sqrt( W12 + W13 );
852 }
853
854 EvtComplex Mat_rhop = BreitWigner( p1, p2, p3, ierr ) +
855 BreitWigner( p1, p3, p2, ierr );
856
857 EvtComplex MatBp = m_Mat_S2 * Mat_rhop * Wtot * ASHQ;
858 EvtComplex MatBm = m_Nat_S2 * Mat_rhop * Wtot * ASHQ;
859
860 real_B0 = real( MatBp );
861 imag_B0 = imag( MatBp );
862
863 real_B0bar = real( MatBm );
864 imag_B0bar = imag( MatBm );
865
866 return ierr;
867}
868
869int EvtBTo3hCP::compute3piP00( EvtVector4R& p1, EvtVector4R& p2,
870 EvtVector4R& p3, double& real_B0, double& imag_B0,
871 double& real_B0bar, double& imag_B0bar, int iset )
872{
873 int ierr = 0;
874 const double ASHQ = sqrt( 2. );
875 double m12 = ( p1 + p2 ).mass();
876 double m13 = ( p1 + p3 ).mass();
877
878 double W12 =
879 1. / ( ( square( m_Mass_rho - m12 ) + square( m_Gam_rho / 2. ) ) * m12 );
880 double W13 =
881 1. / ( ( square( m_Mass_rho - m13 ) + square( m_Gam_rho / 2. ) ) * m13 );
882
883 double Wtot = 1.;
884 if ( iset >= 0 ) {
885 Wtot = 1. / sqrt( W12 + W13 );
886 }
887
888 EvtComplex Mat_rhop = BreitWigner( p1, p2, p3, ierr ) +
889 BreitWigner( p1, p3, p2, ierr );
890
891 EvtComplex MatBp = m_Mat_S1 * Mat_rhop * Wtot * ASHQ;
892 EvtComplex MatBm = m_Nat_S1 * Mat_rhop * Wtot * ASHQ;
893
894 real_B0 = real( MatBp );
895 imag_B0 = imag( MatBp );
896
897 real_B0bar = real( MatBm );
898 imag_B0bar = imag( MatBm );
899
900 return ierr;
901}
902
903int EvtBTo3hCP::computeKpipi( EvtVector4R& p1, EvtVector4R& p2, EvtVector4R& p3,
904 double& real_B0, double& imag_B0,
905 double& real_B0bar, double& imag_B0bar, int iset )
906{
907 int ierr = 0;
908
909 double m12 = ( p1 + p2 ).mass();
910 double m13 = ( p1 + p3 ).mass();
911 double m23 = ( p2 + p3 ).mass();
912
913 double W12 =
914 1. / ( ( square( m_Mass_Kstar0 - m12 ) + square( m_Gam_Kstar0 / 2. ) ) *
915 m12 );
916 double W13 =
917 1. / ( ( square( m_Mass_Kstarp - m13 ) + square( m_Gam_Kstarp / 2. ) ) *
918 m13 );
919 double W23 =
920 1. / ( ( square( m_Mass_rho - m23 ) + square( m_Gam_rho / 2. ) ) * m23 );
921
922 double Wtot = 1.;
923 if ( iset >= 0 ) {
924 Wtot = 1. / sqrt( W12 + W13 + W23 );
925 }
926
927 EvtComplex BW13 = BreitWigner( p1, p3, p2, ierr, m_Mass_Kstarp, m_Gam_Kstarp );
928 if ( ierr != 0 )
929 return ierr;
930 EvtComplex BW12 = BreitWigner( p1, p2, p3, ierr, m_Mass_Kstar0, m_Gam_Kstar0 );
931 if ( ierr != 0 )
932 return ierr;
933 /*
934 If the rho is to be treated on the same footing as K* ==> use the line below
935 EvtComplex BW23=BreitWigner(p2, p3, p1, ierr, Mass_Rho, Gam_Rho);
936 */
937 EvtComplex BW23 = BreitWigner( p2, p3, p1, ierr );
938 if ( ierr != 0 )
939 return ierr;
940
941 // Build up amplitudes
942 EvtComplex MatB0 = m_MatKstarp * BW13 + m_MatKstar0 * BW12 + m_MatKrho * BW23;
943 EvtComplex MatB0bar = m_NatKstarp * BW13 + m_NatKstar0 * BW12 +
944 m_NatKrho * BW23;
945
946 real_B0 = real( MatB0 ) * Wtot;
947 imag_B0 = imag( MatB0 ) * Wtot;
948 real_B0bar = real( MatB0bar ) * Wtot;
949 imag_B0bar = imag( MatB0bar ) * Wtot;
950
951 return ierr;
952}
953
954void EvtBTo3hCP::rotation( EvtVector4R& p, int newRot )
955{
956 if ( newRot ) {
957 double phi2 = EvtRandom::Flat() * 2. * m_pi;
958 double phi3 = EvtRandom::Flat() * 2. * m_pi;
959
960 double c1 = 2. * EvtRandom::Flat() - 1.;
961 double c2 = cos( phi2 );
962 double c3 = cos( phi3 );
963
964 double s1 = sqrt( 1. - square( c1 ) );
965 double s2 = sin( phi2 );
966 double s3 = sin( phi3 );
967
968 m_rotMatrix[0][0] = c1;
969 m_rotMatrix[0][1] = s1 * c3;
970 m_rotMatrix[0][2] = s1 * s3;
971 m_rotMatrix[1][0] = -s1 * c2;
972 m_rotMatrix[1][1] = c1 * c2 * c3 - s2 * s3;
973 m_rotMatrix[1][2] = c1 * c2 * s3 + s2 * c3;
974 m_rotMatrix[2][0] = s1 * s2;
975 m_rotMatrix[2][1] = -c1 * s2 * c3 - c2 * s3;
976 m_rotMatrix[2][2] = -c1 * s2 * s3 + c2 * c3;
977 }
978
979 double mom[3];
980 for ( int i = 1; i < 4; ++i ) {
981 mom[i - 1] = p.get( i );
982 p.set( i, 0 );
983 }
984 for ( int i = 0; i < 3; ++i ) {
985 for ( int j = 0; j < 3; ++j ) {
986 p.set( i + 1, p.get( i + 1 ) + m_rotMatrix[i][j] * mom[j] );
987 }
988 }
989}
990
991void EvtBTo3hCP::gammaGamma( EvtVector4R& p, EvtVector4R& pgamma1,
992 EvtVector4R& pgamma2 )
993{
994 double EGammaCmsPi0 = sqrt( p.mass2() ) / 2.;
995
996 double cosThetaRot = EvtRandom::Flat() * 2. - 1.;
997 double sinThetaRot = sqrt( 1. - square( cosThetaRot ) );
998 double PhiRot = EvtRandom::Flat() * 2. * m_pi;
999
1000 pgamma1.set( 1, EGammaCmsPi0 * sinThetaRot * cos( PhiRot ) );
1001 pgamma1.set( 2, EGammaCmsPi0 * sinThetaRot * sin( PhiRot ) );
1002 pgamma1.set( 3, EGammaCmsPi0 * cosThetaRot );
1003 pgamma1.set( 0, EGammaCmsPi0 );
1004
1005 for ( int i = 1; i < 4; ++i ) {
1006 pgamma2.set( i, -pgamma1.get( i ) );
1007 }
1008 pgamma2.set( 0, pgamma1.get( 0 ) );
1009
1010 pgamma1.applyBoostTo( p );
1011 pgamma2.applyBoostTo( p );
1012}
1013
1014EvtComplex EvtBTo3hCP::BreitWigner( EvtVector4R& p1, EvtVector4R& p2,
1015 EvtVector4R& p3, int& ierr, double Mass,
1016 double Width )
1017{
1018 bool pipiMode = true;
1019
1020 if ( Mass > 1e-5 ) {
1021 pipiMode = false;
1022 }
1023
1024 bool relatBW = true;
1025 bool aleph = true;
1026 EvtComplex result( 0, 0 );
1027 ierr = 0;
1028
1029 double m12 = ( p1 + p2 ).mass();
1030 double e12 = ( p1 + p2 ).get( 0 );
1031 double argu = 1. - square( m12 ) / square( e12 );
1032 double beta = 0;
1033 if ( argu > 0 ) {
1034 beta = sqrt( argu );
1035 } else {
1036 std::cout << "Abnormal beta ! Argu = " << argu << std::endl;
1037 argu = 0;
1038 }
1039 double gamma = e12 / m12;
1040
1041 double m13sq = ( p1 + p3 ).mass2();
1042 double costet = ( 2. * p1.get( 0 ) * p3.get( 0 ) - m13sq + p1.mass2() +
1043 p3.mass2() ) /
1044 ( 2. * p1.d3mag() * p3.d3mag() );
1045
1046 double p1z = p1.d3mag() * costet;
1047 double p1zcms12 = gamma * ( p1z + beta * p1.get( 0 ) );
1048 double e1cms12 = gamma * ( p1.get( 0 ) + beta * p1z );
1049 double p1cms12 = sqrt( square( e1cms12 ) - p1.mass2() );
1050 double coscms = p1zcms12 / p1cms12;
1051
1052 if ( pipiMode ) {
1053 if ( aleph ) {
1054 double m12_2 = square( m12 );
1055 result = coscms * EvtCRhoF_W( m12_2 );
1056 } else {
1057 double factor = 2 * ( square( m_Mass_rho - m12 ) +
1058 square( 0.5 * m_Gam_rho ) );
1059 factor = coscms * m_Gam_rho / factor;
1060 double numReal = ( m_Mass_rho - m12 ) * factor;
1061 double numImg = 0.5 * m_Gam_rho * factor;
1062 result = EvtComplex( numReal, numImg );
1063 }
1064 } else {
1065 if ( relatBW ) {
1066 double Am2Min = p1.mass2() + p2.mass2() + 2 * p1.mass() * p2.mass();
1067 result = coscms *
1068 EvtRBW( square( m12 ), square( Mass ), Width, Am2Min );
1069 } else {
1070 double factor = 2 * ( square( Mass - m12 ) + square( 0.5 * Width ) );
1071 factor = coscms * Width / factor;
1072 double numReal = ( Mass - m12 ) * factor;
1073 double numImg = 0.5 * Width * factor;
1074 result = EvtComplex( numReal, numImg );
1075 }
1076 }
1077
1078 return result;
1079}
1080
1081EvtComplex EvtBTo3hCP::EvtCRhoF_W( double s )
1082{
1083 const bool kuhn_santa = true; // type of Breit-Wigner formula
1084 // double lambda = 1.0;
1085
1086 double AmRho, GamRho, AmRhoP, GamRhoP, beta, AmRhoPP, GamRhoPP, gamma;
1087 if ( kuhn_santa ) {
1088 //...rho(770)
1089 AmRho = 0.7734;
1090 GamRho = 0.1477;
1091 //...rho(1450)
1092 AmRhoP = 1.465;
1093 GamRhoP = 0.696;
1094 beta = -0.229;
1095 //...rho(1700)
1096 AmRhoPP = 1.760;
1097 GamRhoPP = 0.215;
1098 gamma = 0.075;
1099 } else {
1100 //...rho(770)
1101 AmRho = 0.7757;
1102 GamRho = 0.1508;
1103 //...rho(1450)
1104 AmRhoP = 1.448;
1105 GamRhoP = 0.503;
1106 beta = -0.161;
1107 //...rho(1700)
1108 AmRhoPP = 1.757;
1109 GamRhoPP = 0.237;
1110 gamma = 0.076;
1111 }
1112
1113 EvtComplex result( 0, 0 );
1114
1115 if ( kuhn_santa ) {
1116 result = ( EvtcBW_KS( s, square( AmRho ), GamRho ) +
1117 EvtcBW_KS( s, square( AmRhoP ), GamRhoP ) *
1118 ( beta ) +
1119 EvtcBW_KS( s, square( AmRhoPP ), GamRhoPP ) *
1120 ( gamma ) ) /
1121 ( 1. + beta + gamma );
1122 } else {
1123 result = ( EvtcBW_GS( s, square( AmRho ), GamRho ) +
1124 EvtcBW_GS( s, square( AmRhoP ), GamRhoP ) * ( beta ) +
1125 EvtcBW_GS( s, square( AmRhoPP ), GamRhoPP ) * ( gamma ) ) /
1126 ( 1. + beta + gamma );
1127 }
1128 return result;
1129}
1130
1131EvtComplex EvtBTo3hCP::EvtRBW( double s, double Am2, double Gam, double Am2Min )
1132{
1133 EvtComplex result( 0, 0 );
1134
1135 if ( s < Am2Min ) {
1136 return result;
1137 }
1138
1139 double tmp = ( ( s - Am2Min ) / ( Am2 - Am2Min ) );
1140 double G = Gam * ( Am2 / s ) * sqrt( square( tmp ) * tmp );
1141 double D = square( Am2 - s ) + s * square( G );
1142 double X = Am2 * ( Am2 - s );
1143 double Y = Am2 * sqrt( s ) * G;
1144
1145 result = EvtComplex( X / D, Y / D );
1146 return result;
1147}
1148
1149EvtComplex EvtBTo3hCP::EvtcBW_KS( double s, double Am2, double Gam )
1150{
1151 EvtComplex result( 0, 0 );
1152 const double AmPi2 = square( 0.13956995 );
1153 return EvtRBW( s, Am2, Gam, 4. * AmPi2 );
1154}
1155
1156EvtComplex EvtBTo3hCP::EvtcBW_GS( double s, double Am2, double Gam )
1157{
1158 EvtComplex result( 0, 0 );
1159 const double AmPi2 = square( 0.13956995 );
1160
1161 if ( s < 4. * AmPi2 ) {
1162 return result;
1163 }
1164
1165 double tmp = ( ( s - 4. * AmPi2 ) / ( Am2 - 4. * AmPi2 ) );
1166
1167 double G = Gam * ( Am2 / s ) * sqrt( square( tmp ) * tmp );
1168 double z1 = Am2 - s + Evtfs( s, Am2, Gam );
1169 double z2 = sqrt( s ) * G;
1170 double z3 = Am2 + d( Am2 ) * Gam * sqrt( Am2 );
1171
1172 double X = z3 * z1;
1173 double Y = z3 * z2;
1174 double N = square( z1 ) + square( z2 );
1175
1176 result = EvtComplex( X / N, Y / N );
1177
1178 return result;
1179}
1180
1181double EvtBTo3hCP::d( double AmRho2 )
1182{
1183 const double AmPi = 0.13956995;
1184 const double AmPi2 = square( AmPi );
1185 double AmRho = sqrt( AmRho2 );
1186 double k_AmRho2 = k( AmRho2 );
1187 double result = 3. / m_pi * AmPi2 / square( k_AmRho2 ) *
1188 log( ( AmRho + 2. * k_AmRho2 ) / ( 2. * AmPi ) ) +
1189 AmRho / ( 2. * m_pi * k_AmRho2 ) -
1190 AmPi2 * AmRho / ( m_pi * ( square( k_AmRho2 ) * k_AmRho2 ) );
1191 return result;
1192}
1193
1194double EvtBTo3hCP::k( double s )
1195{
1196 const double AmPi2 = square( 0.13956995 );
1197 return 0.5 * sqrt( s - 4. * AmPi2 );
1198}
1199
1200double EvtBTo3hCP::Evtfs( double s, double AmRho2, double GamRho )
1201{
1202 double k_s = k( s );
1203 double k_Am2 = k( AmRho2 );
1204
1205 return GamRho * AmRho2 / ( square( k_Am2 ) * k_Am2 ) *
1206 ( square( k_s ) * ( h( s ) - h( AmRho2 ) ) +
1207 ( AmRho2 - s ) * square( k_Am2 ) * dh_ds( AmRho2 ) );
1208}
1209
1210double EvtBTo3hCP::h( double s )
1211{
1212 const double AmPi = 0.13956995;
1213 double sqrts = sqrt( s );
1214 double k_s = k( s );
1215 return 2. / m_pi * ( k_s / sqrts ) *
1216 log( ( sqrts + 2. * k_s ) / ( 2. * AmPi ) );
1217}
1218
1219double EvtBTo3hCP::dh_ds( double s )
1220{
1221 return h( s ) * ( 1. / ( 8. * square( k( s ) ) ) - 1. / ( 2 * s ) ) +
1222 1. / ( 2. * m_pi * s );
1223}
square
#define square(x)
Definition EvtBTo3hCP.cpp:37
imag
double imag(const EvtComplex &c)
Definition EvtComplex.hh:237
real
double real(const EvtComplex &c)
Definition EvtComplex.hh:232
EvtBTo3hCP::k
double k(double s)
Definition EvtBTo3hCP.cpp:1194
EvtBTo3hCP::m_Nat_S3
EvtComplex m_Nat_S3
Definition EvtBTo3hCP.hh:95
EvtBTo3hCP::m_betaCP
double m_betaCP
Definition EvtBTo3hCP.hh:98
EvtBTo3hCP::compute3pi
int compute3pi(EvtVector4R &p1, EvtVector4R &p2, EvtVector4R &p3, double &real_B0, double &imag_B0, double &real_B0bar, double &imag_B0bar, int set)
Definition EvtBTo3hCP.cpp:788
EvtBTo3hCP::m_Mat_S1
EvtComplex m_Mat_S1
Definition EvtBTo3hCP.hh:94
EvtBTo3hCP::EvtcBW_GS
EvtComplex EvtcBW_GS(double s, double Am2, double Gam)
Definition EvtBTo3hCP.cpp:1156
EvtBTo3hCP::m_Nat_S4
EvtComplex m_Nat_S4
Definition EvtBTo3hCP.hh:95
EvtBTo3hCP::EvtCRhoF_W
EvtComplex EvtCRhoF_W(double s)
Definition EvtBTo3hCP.cpp:1081
EvtBTo3hCP::m_rotMatrix
double m_rotMatrix[3][3]
Definition EvtBTo3hCP.hh:115
EvtBTo3hCP::m_Gam_rho
double m_Gam_rho
Definition EvtBTo3hCP.hh:104
EvtBTo3hCP::compute3piP00
int compute3piP00(EvtVector4R &p1, EvtVector4R &p2, EvtVector4R &p3, double &real_B0, double &imag_B0, double &real_B0bar, double &imag_B0bar, int set)
Definition EvtBTo3hCP.cpp:869
EvtBTo3hCP::m_Mat_S3
EvtComplex m_Mat_S3
Definition EvtBTo3hCP.hh:94
EvtBTo3hCP::m_NatKrho
EvtComplex m_NatKrho
Definition EvtBTo3hCP.hh:96
EvtBTo3hCP::d
double d(double AmRho2)
Definition EvtBTo3hCP.cpp:1181
EvtBTo3hCP::Evt3piP00
void Evt3piP00(double alpha, int iset, EvtVector4R &p_p1, EvtVector4R &p_p1_gamma1, EvtVector4R &p_p1_gamma2, EvtVector4R &p_p2_gamma1, EvtVector4R &p_p2_gamma2, double &Real_B0, double &Imag_B0, double &Real_B0bar, double &Imag_B0bar)
Definition EvtBTo3hCP.cpp:280
EvtBTo3hCP::gammaGamma
void gammaGamma(EvtVector4R &p, EvtVector4R &pgamma1, EvtVector4R &pgamma2)
Definition EvtBTo3hCP.cpp:991
EvtBTo3hCP::generateSqMasses_Kpipi
void generateSqMasses_Kpipi(double &m12, double &m13, double &m23, double MB2, double m1sq, double m2sq, double m3sq)
Definition EvtBTo3hCP.cpp:552
EvtBTo3hCP::m_MatKstarp
EvtComplex m_MatKstarp
Definition EvtBTo3hCP.hh:95
EvtBTo3hCP::m_Mat_S5
EvtComplex m_Mat_S5
Definition EvtBTo3hCP.hh:94
EvtBTo3hCP::m_Mat_S4
EvtComplex m_Mat_S4
Definition EvtBTo3hCP.hh:94
EvtBTo3hCP::h
double h(double s)
Definition EvtBTo3hCP.cpp:1210
EvtBTo3hCP::generateSqMasses_3piMPP
void generateSqMasses_3piMPP(double &m12, double &m13, double &m23, double MB2, double m1sq, double m2sq, double m3sq)
Definition EvtBTo3hCP.cpp:683
EvtBTo3hCP::m_Gam_Kstarp
double m_Gam_Kstarp
Definition EvtBTo3hCP.hh:112
EvtBTo3hCP::Evtfs
double Evtfs(double s, double AmRho2, double GamRho)
Definition EvtBTo3hCP.cpp:1200
EvtBTo3hCP::m_factor_max
double m_factor_max
Definition EvtBTo3hCP.hh:116
EvtBTo3hCP::m_alphaCP
double m_alphaCP
Definition EvtBTo3hCP.hh:97
EvtBTo3hCP::m_NatKstarp
EvtComplex m_NatKstarp
Definition EvtBTo3hCP.hh:96
EvtBTo3hCP::m_pi
double m_pi
Definition EvtBTo3hCP.hh:102
EvtBTo3hCP::BreitWigner
EvtComplex BreitWigner(EvtVector4R &p1, EvtVector4R &p2, EvtVector4R &p3, int &ierr, double Mass=0, double Width=0)
Definition EvtBTo3hCP.cpp:1014
EvtBTo3hCP::m_Mass_rho
double m_Mass_rho
Definition EvtBTo3hCP.hh:103
EvtBTo3hCP::m_Nat_S2
EvtComplex m_Nat_S2
Definition EvtBTo3hCP.hh:95
EvtBTo3hCP::m_Mass_Kstar0
double m_Mass_Kstar0
Definition EvtBTo3hCP.hh:111
EvtBTo3hCP::firstStep
void firstStep(EvtVector4R &p1, EvtVector4R &p2, EvtVector4R &p3, int mode)
Definition EvtBTo3hCP.cpp:431
EvtBTo3hCP::m_Nat_S5
EvtComplex m_Nat_S5
Definition EvtBTo3hCP.hh:95
EvtBTo3hCP::EvtKpipi
void EvtKpipi(double alpha, double beta, int iset, EvtVector4R &p_K_plus, EvtVector4R &p_pi_minus, EvtVector4R &p_gamma_1, EvtVector4R &p_gamma_2, double &Real_B0, double &Imag_B0, double &Real_B0bar, double &Imag_B0bar)
Definition EvtBTo3hCP.cpp:356
EvtBTo3hCP::rotation
void rotation(EvtVector4R &p, int newRot)
Definition EvtBTo3hCP.cpp:954
EvtBTo3hCP::computeKpipi
int computeKpipi(EvtVector4R &p1, EvtVector4R &p2, EvtVector4R &p3, double &real_B0, double &imag_B0, double &real_B0bar, double &imag_B0bar, int set)
Definition EvtBTo3hCP.cpp:903
EvtBTo3hCP::generateSqMasses_3pi
void generateSqMasses_3pi(double &m12, double &m13, double &m23, double MB2, double m1sq, double m2sq, double m3sq)
Definition EvtBTo3hCP.cpp:620
EvtBTo3hCP::m_MC2
double m_MC2
Definition EvtBTo3hCP.hh:101
EvtBTo3hCP::m_M_B
double m_M_B
Definition EvtBTo3hCP.hh:105
EvtBTo3hCP::m_M_pi0
double m_M_pi0
Definition EvtBTo3hCP.hh:108
EvtBTo3hCP::m_M_pip
double m_M_pip
Definition EvtBTo3hCP.hh:106
EvtBTo3hCP::EvtRBW
EvtComplex EvtRBW(double s, double Am2, double Gam, double Am2Min)
Definition EvtBTo3hCP.cpp:1131
EvtBTo3hCP::generateSqMasses_3piP00
void generateSqMasses_3piP00(double &m12, double &m13, double &m23, double MB2, double m1sq, double m2sq, double m3sq)
Definition EvtBTo3hCP.cpp:736
EvtBTo3hCP::m_Mass_Kstarp
double m_Mass_Kstarp
Definition EvtBTo3hCP.hh:110
EvtBTo3hCP::m_MA2
double m_MA2
Definition EvtBTo3hCP.hh:99
EvtBTo3hCP::setConstants
void setConstants(double balpha, double bbeta)
Definition EvtBTo3hCP.cpp:48
EvtBTo3hCP::EvtcBW_KS
EvtComplex EvtcBW_KS(double s, double Am2, double Gam)
Definition EvtBTo3hCP.cpp:1149
EvtBTo3hCP::m_MB2
double m_MB2
Definition EvtBTo3hCP.hh:100
EvtBTo3hCP::compute3piMPP
int compute3piMPP(EvtVector4R &p1, EvtVector4R &p2, EvtVector4R &p3, double &real_B0, double &imag_B0, double &real_B0bar, double &imag_B0bar, int set)
Definition EvtBTo3hCP.cpp:835
EvtBTo3hCP::m_MatKrho
EvtComplex m_MatKrho
Definition EvtBTo3hCP.hh:96
EvtBTo3hCP::m_Gam_Kstar0
double m_Gam_Kstar0
Definition EvtBTo3hCP.hh:113
EvtBTo3hCP::m_M_pim
double m_M_pim
Definition EvtBTo3hCP.hh:107
EvtBTo3hCP::m_M_Kp
double m_M_Kp
Definition EvtBTo3hCP.hh:109
EvtBTo3hCP::m_Mat_S2
EvtComplex m_Mat_S2
Definition EvtBTo3hCP.hh:94
EvtBTo3hCP::m_NatKstar0
EvtComplex m_NatKstar0
Definition EvtBTo3hCP.hh:96
EvtBTo3hCP::dh_ds
double dh_ds(double s)
Definition EvtBTo3hCP.cpp:1219
EvtBTo3hCP::Evt3piMPP
void Evt3piMPP(double alpha, int iset, EvtVector4R &p_p1, EvtVector4R &p_p2, EvtVector4R &p_p3, double &Real_B0, double &Imag_B0, double &Real_B0bar, double &Imag_B0bar)
Definition EvtBTo3hCP.cpp:211
EvtBTo3hCP::Evt3pi
void Evt3pi(double alpha, int iset, EvtVector4R &p_K_plus, EvtVector4R &p_pi_minus, EvtVector4R &p_gamma_1, EvtVector4R &p_gamma_2, double &Real_B0, double &Imag_B0, double &Real_B0bar, double &Imag_B0bar)
Definition EvtBTo3hCP.cpp:134
EvtBTo3hCP::m_MatKstar0
EvtComplex m_MatKstar0
Definition EvtBTo3hCP.hh:95
EvtBTo3hCP::m_Nat_S1
EvtComplex m_Nat_S1
Definition EvtBTo3hCP.hh:94
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cpp:95
EvtVector4R::mass
double mass() const
Definition EvtVector4R.cpp:48
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:163
EvtVector4R::d3mag
double d3mag() const
Definition EvtVector4R.cpp:192
EvtVector4R::mass2
double mass2() const
Definition EvtVector4R.hh:100
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:168
EvtVector4R::applyBoostTo
void applyBoostTo(const EvtVector4R &p4, bool inverse=false)
Definition EvtVector4R.cpp:111
Generated by doxygen 1.16.1