TEvtGen/EvtGenModels/EvtBToKpipi.F

   1 C--------------------------------------------------------------------------
   2 C
   3 C Environment:
   4 C      This software is part of the EvtGen package developed jointly
   5 C      for the BaBar and CLEO collaborations.  If you use all or part
   6 C      of it, please give an appropriate acknowledgement.
   7 C
   8 C Copyright Information: See EvtGen/COPYRIGHT
   9 C      Copyright (C) 1998      Caltech, UCSB
  10 C
  11 C Module: EvtBToKpipi.F
  12 C
  13 C Description:
  14 C
  15 C Modification history:
  16 C
  17 C    DJL/RYD     August 11, 1998         Module created
  18 C
  19 C------------------------------------------------------------------------
  20 C===================================================================
  21 C This package is providing a neutral B -->-- K pi pi decay generator
  22 C Its is composed of the following subroutines:
  23 C
  24 C [*] HowToUse
  25 C                 This is an How To Use routine where one may find the
  26 C                 implementation of the time dependance: That is to
  27 C                 say that it shows how the output of the routine is
  28 C                 supposed to be used in the mind of the authors.
  29 C
  30 C===================================================================
  31 C [0] EVTKpipi
  32 C                 The routine to be called. Note that on the first call
  33 C                 some initialization will be made, in particular the
  34 C                 computation of a normalization factor designed to help
  35 C                 the subsequent time-dependent generation of events.
  36 C                 The normalisation done inside EVTKpipi is such that
  37 C                 at the level of the time implementation, the maximum
  38 C                 time-dependant weight is unity : nothing is to be
  39 C                 computed to generate unity-weight events. The exact
  40 C                 meaning of the above is made explicit in the HowToUse
  41 C                 routine.
  42 C [1] first_step_Kpipi
  43 C                 Generation of the kinematics of the 3 prongs
  44 C                 It uses the function evtranf which is a random number
  45 C                 generator providing an uniform distribution
  46 C                 of Real*4 random number between 0 and 1
  47 C [2] compute
  48 C                 return the amplitudes of the B0 -->-- K+pi-pi0
  49 C                 corrected for the generation mechanism.
  50 c                 Note that this is a Tagging Mode. The CP conjugate
  51 c                 mode (B0bar -->-- K-pi+pi0) is treated at the same time.
  52 C [3] BreitWigner
  53 C                 compute the Breit-Wigner of the contributing K* and rho
  54 C                 taking into account the cosine term linked to the
  55 C                 zero-helicity of the spin-1 resonances. There is two
  56 c                 forms of Breit-Wigners available. The first one is the
  57 c                 simple non-relativistic form, while the second is the
  58 C                 relativistic expressions.
  59 C                 The default setting is the relativistic one.
  60 C [4] Set_constants
  61 C                 Set the constants values, from the pion mass to the
  62 C                 penguin contributions. It is called by EVTKpipi
  63 C
  64 C And other routines which do not deserve comment here.
  65 C===================================================================
  66 c      Implicit none
  67 C      Real Hmemor
  68 C      Common/Pawc/Hmemor(1000000)
  69 C      Call Hlimit(1000000)
  70 C      Call EvtHowToUse_Kpipi
  71 C      Stop
  72 C      End
  73
  74
  75       subroutine EvtHowToUse_Kpipi
  76
  77       Implicit none
  78       Real*8  alphaCP, betaCP
  79       Integer iset,number,j,N_gener,N_asked
  80
  81       Real*8  p_K_plus(4),p_pi_minus(4)
  82       Real*8  p_gamma_1(4),p_gamma_2(4)
  83
  84       Real*8  Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar
  85       Real*8  Weight,Weight_max
  86       Real*8  m_Kstarp,m_Kstar0,m_rhom
  87       Real*8  Wrong
  88       Real*4  Evtranf,Tag
  89
  90       alphaCP      = 1.35
  91       betaCP       = 0.362
  92       N_gener    = 0
  93       N_asked    = 100000
  94
  95       weight_max = 1.0
  96
  97 c run : Simulation of the Anders Ryd Generator
  98
  99       Do number=1,N_asked ! weighted events as generated here
 100
 101       If(number.eq.1) then
 102       iset=10000          ! 10^4 events are used to normalize amplitudes
 103       Else
 104       iset=0              ! iset must be reset to zero after the first call
 105       End If
 106
 107 c Here is the call to EVTKpipi  !!!!!!!!!!!!!!!!
 108 c communication of data is done by argument only <<<<<<<<
 109        call EVTKpipi(
 110      + alphaCP,betaCP,iset,
 111      + p_K_plus,p_pi_minus,
 112      + p_gamma_1,p_gamma_2,
 113      + Real_B0,Imag_B0,real_B0bar,Imag_B0bar)
 114 C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 115 c select the Tag
 116       Tag   =evtranf()
 117 c generate acording to the tag
 118       If(Tag.gt.0.5) Then
 119
 120 c a B0bar tag =>  the decay is one from a B0
 121       Weight= Real_B0      **2 + Imag_B0      **2
 122
 123       Else
 124
 125 c a B0    tag =>  the decay is one from a B0bar
 126       Weight= Real_B0bar   **2 + Imag_B0bar   **2
 127
 128       End If
 129
 130       If(Weight.Gt.evtranf()) Then
 131 c----------------------------------------------------------------------------
 132 c unweighted event production
 133 c----------------------------------------------------------------------------
 134       N_gener=N_gener+1
 135 C here is just a Dalitz plot and a few prints
 136
 137       m_Kstarp=(p_K_plus   (4)+p_gamma_1(4)+p_gamma_2(4))**2
 138       m_rhom  =(p_pi_minus (4)+p_gamma_1(4)+p_gamma_2(4))**2
 139       m_Kstar0=(p_K_plus   (4)+p_pi_minus (4))           **2
 140       do j=1,3
 141       m_Kstarp=m_Kstarp -(p_K_plus   (j)+p_gamma_1(j)+p_gamma_2(j))**2
 142       m_rhom  =m_rhom   -(p_pi_minus (j)+p_gamma_1(j)+p_gamma_2(j))**2
 143       m_Kstar0=m_Kstar0 -(p_K_plus   (j)+p_pi_minus (j))           **2
 144       end do
 145
 146 c here is the Dalitz plot when assuming the pion mass for the Kaon
 147 c the energy is redefined
 148       Wrong =0.
 149       do j=1,3
 150       Wrong = Wrong + p_K_plus(j)**2
 151       end do
 152       Wrong = Dsqrt(Wrong+0.139**2)
 153
 154       m_Kstarp=(Wrong         +p_gamma_1(4)+p_gamma_2(4))**2
 155       m_Kstar0=(Wrong         +p_pi_minus (4))           **2
 156       do j=1,3
 157       m_Kstarp=m_Kstarp -(p_K_plus   (j)+p_gamma_1(j)+p_gamma_2(j))**2
 158       m_Kstar0=m_Kstar0 -(p_K_plus   (j)+p_pi_minus (j))           **2
 159       end do
 160
 161 c here is a check that weight_max is one
 162
 163       If(Weight.gt.Weight_max) Then
 164       Weight_max=Weight
 165       Print*,' overweighted event found at weight = ',Weight_max
 166       End If
 167
 168 c----------------------------------------------------------------------------
 169       End If
 170
 171 c end of the loop over events
 172       End Do
 173
 174       Print*,'number of unity-weight events generated : ',N_gener
 175       Print*,'number of trials                        : ',N_asked
 176
 177       End
 178 C===================================================================
 179       subroutine EvtKpipi(
 180      + alpha_input,beta_input,iset,
 181      + p_K_plus,p_pi_minus,
 182      + p_gamma_1,p_gamma_2,
 183      + Real_B0,Imag_B0,Real_B0bar,Imag_B0bar)
 184 c-----------------------------------------------------------------
 185 c       ----------------------------------------------------------
 186 c       --- This is the routine to be called by the Main generator
 187 c           to get the decay of B0    -->-- K+ pi- pi0
 188 c           The decay proceeeds through three channels:
 189 c           a) B0 -->-- K*+ pi-  ; K*+    -->-- K+ pi0
 190 c           b)          K*0 pi0  ; K*0bar -->-- K+ pi-
 191 c           c)          K-  rho+ ; rho+   -->-- pi+ pi0
 192 c           .) The K0 rho0 channel is not implemented since it does
 193 c              not lead to Kpipi final state, but it is interesting
 194 c              in itself.
 195 c       It provides at the same time the CP conjugate decay
 196 c                               B0bar -->-- K- pi+ pi0
 197 c****************************************************************************
 198 c       --- Outputs are :
 199 c
 200 c       ---               p_K_plus      : the four momentum of the K+
 201 c       ---               p_pi_minus    : ........................ pi-
 202 c       ---               p_gamma_1     : ........................ first  gamma of the pi0
 203 c       ---               p_gamma_2     : ........................ second gamma of the pi0
 204 c
 205 c             Note that : the energy is stored in the fourth component
 206 c                         the values are the ones of the B rest frame
 207 c                         a random rotation has been applied
 208 c
 209 c       ---               Real_B0      : The real      part of the amplitude of the decay
 210 c                                                                 B0    -->-- K+ pi- pi0
 211 c       ---               Imag_B0      : ... imaginary ..................................
 212 c            similarly
 213 c       ---               Real_B0bar   : The real      part of the amplitude of the decay
 214 c                                                                 B0bar -->-- K- pi+ pi0
 215 c       ---               Imag_B0bar   : ... imaginary ..................................
 216
 217 c****************************************************************************
 218 c-----------------------------------------------------------------
 219       Implicit none
 220 #include "EvtGenModels/EvtBTo3pi.inc"
 221       Real*8  alpha_input, beta_input
 222       Integer iset
 223       Real*8  p_K_plus(4),p_pi_minus(4)
 224       Real*8  p_gamma_1(4),p_gamma_2(4)
 225       Real*8  Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar
 226
 227 c Working quantities
 228       Integer i,number
 229       Real*8  p1(5),p2(5),p3(5)
 230       Real*8  Gamma1(5),Gamma2(5)
 231       Real*8  factor_max,ABp,ABm
 232       Integer ierr
 233         data factor_max/1.D+00/
 234         ierr =0
 235 c-------------------------------------------------------------------
 236       If(iset.eq.0) Then
 237 c-------------------------------------------------------------------
 238 c     this is the normal mode of operation
 239 c     First, generate the kinematics
 240
 241       p1(5)= M_Kp **2
 242       p2(5)= M_pim**2
 243       p3(5)= M_pi0**2
 244
 245  10   continue
 246       call Evtfirst_step_Kpipi(p1,p2,p3)
 247
 248 c     Then, compute the amplitudes
 249
 250       Call EvtCompute_Kpipi(p1,p2,p3,
 251      +    Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
 252         if(ierr.ne.0 ) Go To 10
 253 c----------------nedit EvtBto---------------------------------------------------
 254       ElseIf(iset.lt.0) Then
 255 c-------------------------------------------------------------------
 256 c     This is an user mode of operation where the kinematics is
 257 c     provided by the user who only wants the corresponding amplitudes
 258 c     to be computed
 259
 260       Do i=1,4
 261       p1(i)= p_K_plus  (i)
 262       p2(i)= p_pi_minus(i)
 263       p3(i)= p_gamma_1 (i) + p_gamma_2 (i)
 264       End Do
 265       p1(5)= M_Kp **2
 266       p2(5)= M_pim**2
 267       p3(5)= M_pi0**2
 268
 269       Call EvtCompute_Kpipi(p1,p2,p3,
 270      +    Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
 271
 272        if(ierr.ne.0) Then
 273         Print*,'the provided kinematics are not physical'
 274         Print*,'ierr=',ierr
 275         Print*,'the program will stop'
 276         Stop
 277        endif
 278 c-------------------------------------------------------------------
 279       ElseIf(iset.gt.0) Then
 280 c-------------------------------------------------------------------
 281 c     This is the pre-run mode of operation where initializations are
 282 c     performed.
 283
 284       factor_max= 0
 285       call Evtset_constants(alpha_input, beta_input)
 286       p1(5)= M_Kp **2
 287       p2(5)= M_pim**2
 288       p3(5)= M_pi0**2
 289
 290 c     pre-run
 291       Do number=1,iset
 292
 293   20  continue
 294       call Evtfirst_step_Kpipi(p1,p2,p3)
 295
 296       Call EvtCompute_Kpipi(p1,p2,p3,
 297      +    Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
 298         if(ierr.ne.0) Go To 20
 299       ABp   = Real_B0      **2 + Imag_B0      **2
 300       ABm   = Real_B0bar   **2 + Imag_B0Bar   **2
 301
 302       If(ABp.gt.factor_max) factor_max=ABp
 303       If(ABm.gt.factor_max) factor_max=ABm
 304
 305       End Do
 306 c     end of the pre-run
 307
 308       factor_max=1.D+00/Dsqrt(factor_max)
 309
 310 c-------------------------------------------------------------------
 311       End If
 312 c-------------------------------------------------------------------
 313
 314       Real_B0   =Real_B0    * factor_max
 315       Imag_B0   =Imag_B0    * factor_max
 316       Real_B0bar=Real_B0bar * factor_max
 317       Imag_B0Bar=Imag_B0Bar * factor_max
 318
 319       if(iset.lt.0) return
 320 c     P1,p2,p3 ---> random rotation in B rest frame
 321
 322       Call EvtRotation(p1,1)
 323       Call EvtRotation(p2,0)
 324       Call EvtRotation(p3,0)
 325
 326 C     Desintegrate the pi_0 s
 327
 328       Call EvtGammaGamma(p3,Gamma1,Gamma2)
 329
 330 C     Feed the output four vectors
 331
 332       Do i=1,4
 333
 334       p_K_plus  (i)=p1    (i)
 335       p_pi_minus(i)=p2    (i)
 336       p_gamma_1 (i)=Gamma1(i)
 337       p_gamma_2 (i)=Gamma2(i)
 338       End Do
 339
 340       Return
 341
 342       End
 343
 344 c===================================================================
 345       subroutine Evtfirst_step_Kpipi(P1,P2,P3)
 346 c-----------------------------------------------------------------
 347 c       ----------------------------------------------------------
 348 c       --- This routine generates the 5-vectors P1,P2,P3
 349 c       --- Associated respectively with the Pi+ and two Pi0 s
 350 c       ---             P1(1) = Px
 351 c       ---             P1(2) = Py
 352 c       ---             P1(3) = Pz
 353 c       ---             P1(4) = E
 354 c       ---             P1(5) = M**2
 355 c       ----------------------------------------------------------
 356 c       ---     Input Four Vectors
 357 C       ---     Particle [1] is the K+
 358 C       ---     Particle [2] is the pi-
 359 C       ---     Particle [3] is the pi0
 360 c       ----------------------------------------------------------
 361
 362 c       ----------------------------------------------------------
 363 c       --- commons
 364 c       ----------------------------------------------------------
 365 #include "EvtGenModels/EvtBTo3pi.inc"
 366 c       ----------------------------------------------------------
 367 c       --- Used Functions
 368 c       ----------------------------------------------------------
 369
 370         real evtranf
 371
 372 c       ----------------------------------------------------------
 373 c       --- Variables in Argument
 374 c       ----------------------------------------------------------
 375
 376          real*8 P1(5),P2(5),P3(5)
 377
 378 c       ----------------------------------------------------------
 379 c       --- Local Variables
 380 c       ----------------------------------------------------------
 381
 382         real*8 m12,min_m12, max_m12
 383         real*8 m13,min_m13, max_m13
 384         real*8 m23,min_m23, max_m23
 385         Real*8 cost13,cost12,cost23
 386
 387         real*8 p1mom,p2mom,p3mom
 388
 389         real*8 x, y, z, mass
 390         integer i
 391
 392         Logical Phase_space
 393         data Phase_space/.false./
 394
 395 c       ----------------------------------------------------------
 396 c       --- Computation
 397 c       ----------------------------------------------------------
 398         max_m12 = M_B**2
 399         min_m12 = P1(5) + P2(5) + 2.*Dsqrt(p1(5)*p2(5))
 400
 401         max_m13 = M_B**2
 402         min_m13 = P1(5) + P3(5) + 2.*Dsqrt(p1(5)*p3(5))
 403
 404         max_m23 = M_B**2
 405         min_m23 = P2(5) + P3(5) + 2.*Dsqrt(p2(5)*p3(5))
 406
 407 100     Continue
 408
 409 c       ----------------------------------------------------------
 410 c       --- Generation of the Mass of the Rho or Kstar
 411 c       ----------------------------------------------------------
 412
 413
 414 c       ----------------------------------------------------------
 415 c       --- z is the Flag needed to choose between the generation
 416 c       --- of a K*+, K*0 or rho- resonance
 417 c       ----------------------------------------------------------
 418
 419         z = 3.*evtranf()
 420
 421         MC2 = M_B**2
 422
 423         If(z.lt.1.) Then
 424 c K*+ production
 425                 If(Phase_space) Then
 426                 m13 = evtranf()*(max_m13-min_m13)+min_m13
 427                 Else
 428                 y = evtranf()*PI - PI/2.
 429                 x = Dtan(y)
 430                 mass = x*Gam_Kstarp/2. +Mass_Kstarp
 431                 m13 = mass**2
 432                 End If
 433
 434         m12 = evtranf()*(max_m12-min_m12)+min_m12
 435         m23 = MC2 - m12 - m13
 436
 437         ElseIf(z.lt.2.) Then
 438 c K*0 production
 439                 If(Phase_space) Then
 440                 m12 = evtranf()*(max_m12-min_m12)+min_m12
 441                 Else
 442                 y = evtranf()*PI - PI/2.
 443                 x = Dtan(y)
 444                 mass = x*Gam_Kstar0/2. +Mass_Kstar0
 445                 m12 =  mass**2
 446                 End If
 447
 448         m13 = evtranf()*(max_m13-min_m13)+min_m13
 449         m23 = MC2 - m12 - m13
 450
 451         Else
 452 c rho- production
 453                 If(Phase_space) Then
 454                 m23 = evtranf()*(max_m23-min_m23)+min_m23
 455                 Else
 456                 y = evtranf()*PI - PI/2.
 457                 x = Dtan(y)
 458                 mass = x*Gam_rho/2. +Mass_rho
 459                 m23 =  mass**2
 460                 End If
 461
 462         m13 = evtranf()*(max_m13-min_m13)+min_m13
 463         m12 = MC2 - m23 - m13
 464
 465         Endif
 466
 467 c       ----------------------------------------------------------
 468 c       --- Check that the physics is OK :
 469 c       --- Are the invariant Masses in allowed ranges ?
 470 c       ----------------------------------------------------------
 471
 472         If(m23.lt.min_m23.or.m23.gt.max_m23) Go to 100
 473         If(m13.lt.min_m13.or.m13.gt.max_m13) Go to 100
 474         If(m12.lt.min_m12.or.m12.gt.max_m12) Go to 100
 475
 476 c       ----------------------------------------------------------
 477 c       --- Are the Cosines of the angles between particles
 478 c       --- Between -1 and +1 ?
 479 c       ----------------------------------------------------------
 480
 481         P1(4)=(M_B**2+P1(5)-m23)/(2.*M_B)
 482         P2(4)=(M_B**2+P2(5)-m13)/(2.*M_B)
 483         P3(4)=(M_B**2+P3(5)-m12)/(2.*M_B)
 484
 485         p1mom=p1(4)**2-P1(5)
 486         p2mom=p2(4)**2-P2(5)
 487         p3mom=p3(4)**2-P3(5)
 488         If(p1mom.lt.0) Go to 100
 489         If(p2mom.lt.0) Go to 100
 490         If(p3mom.lt.0) Go to 100
 491         p1mom=Dsqrt(p1mom)
 492         p2mom=Dsqrt(p2mom)
 493         p3mom=Dsqrt(p3mom)
 494
 495         cost13=(2.*p1(4)*p3(4)+P1(5)+p3(5)-m13)/(2.*p1mom*p3mom)
 496         cost12=(2.*p1(4)*p2(4)+P1(5)+p2(5)-m12)/(2.*p1mom*p2mom)
 497         cost23=(2.*p2(4)*p3(4)+P2(5)+p3(5)-m23)/(2.*p2mom*p3mom)
 498         If(Dabs(cost13).gt.1.) Go to 100
 499         If(Dabs(cost12).gt.1.) Go to 100
 500         If(Dabs(cost23).gt.1.) Go to 100
 501
 502 c       ----------------------------------------------------------
 503 c       --- Filling the 5-vectors P1,P2,P3
 504 c       ----------------------------------------------------------
 505
 506         P3(1) = 0
 507         P3(2) = 0
 508         p3(3) = p3mom
 509
 510         P1(3) = p1mom*cost13
 511         P1(1) = p1mom*Dsqrt(1.D+00-cost13**2)
 512         p1(2) = 0.
 513
 514         Do i=1,3
 515         P2(i)=-p1(i)-p3(i)
 516         End do
 517
 518         END
 519
 520
 521 c======================================================================
 522       Subroutine EvtCompute_Kpipi(p1,p2,p3,
 523      +           Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
 524 c-----------------------------------------------------------------------
 525         IMPLICIT None
 526 #include "EvtGenModels/EvtBTo3pi.inc"
 527
 528
 529         Real*8 m12, m13, m23, W12, W13, W23, Wtot
 530         Real*4 evt_gmas
 531         Complex*16 MatB0,MatB0bar,BW12,BW13,BW23
 532         Real*8  Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar
 533         Real*8     p1(5),p2(5),p3(5)
 534
 535
 536         Integer ierr,iset
 537         Complex*16 BrightWagner,BreitWigner
 538
 539         ierr = 0
 540 c       ----------------------------------------------------------------
 541 C       ---     Account for the pole compensation
 542 c       ----------------------------------------------------------------
 543
 544         m12 = evt_gmas(p1,p2)
 545         m13 = evt_gmas(p1,p3)
 546         m23 = evt_gmas(p2,p3)
 547
 548         if(m12.lt.0. .or. m13.lt.0. .or. m23.lt.0.) Then
 549                 ierr=1
 550                 Print*,'ierr = ',ierr
 551                 return
 552         endif
 553
 554         m12 = sqrt(m12)
 555         m13 = sqrt(m13)
 556         m23 = sqrt(m23)
 557
 558         W12 = 1. / (((Mass_Kstar0 - m12)**2+(Gam_Kstar0/2.)**2)*m12)
 559         W13 = 1. / (((Mass_Kstarp - m13)**2+(Gam_Kstarp/2.)**2)*m13)
 560         W23 = 1. / (((Mass_rho    - m23)**2+(Gam_rho   /2.)**2)*m23)
 561
 562         if(iset.ge.0) Then
 563         Wtot = 1.D+00/Dsqrt(W12 + W13 + W23)
 564             else
 565         Wtot =1.
 566         Endif
 567 c       ----------------------------------------------------------------
 568 C       ---     Compute Breit-Wigners
 569 c       ----------------------------------------------------------------
 570
 571         BW13=BrightWagner(p1,p3,p2,Mass_Kstarp,Gam_Kstarp,ierr)
 572         If(ierr.ne.0) Return
 573         BW12=BrightWagner(p1,p2,p3,Mass_Kstar0,Gam_Kstar0,ierr)
 574         If(ierr.ne.0) Return
 575 c If the rho is to be treated on the same footing as K* ==> use the line below
 576 c        BW23=BrightWagner(p2,p3,p1,Mass_Rho   ,Gam_Rho   ,ierr)
 577         BW23=BreitWigner(p2,p3,p1,ierr)
 578         If(ierr.ne.0) Return
 579
 580 c       ----------------------------------------------------------------
 581 c        -              and
 582 C       ---     Build up the amplitudes
 583 c       ----------------------------------------------------------------
 584
 585 c Here come the relative amplitudes of the three decay channels
 586 c First, one computes the B0     decay B0    ->- K+ pi- pi0
 587         MatB0    = MatKstarp * BW13
 588      >           + MatKstar0 * BW12
 589      >           + MatKrho   * BW23
 590 c Second, one computes the B0bar decay B0bar ->- K- pi+ pi0
 591         MatB0bar = NatKstarp * BW13
 592      >           + NatKstar0 * BW12
 593      >           + NatKrho   * BW23
 594
 595
 596 c       Pick up the Real and Imaginary parts
 597
 598         Real_B0    = dreal(MatB0   )*Wtot
 599         Imag_B0    = Imag(MatB0   )*Wtot
 600
 601         Real_B0bar = dreal(MatB0bar)*Wtot
 602         Imag_B0Bar = Imag(MatB0bar)*Wtot
 603
 604         Return
 605         End
 606 c======================================================================
 607       Function BrightWagner(p1,p2,p3,Mass,Width,ierr)
 608 c----------------------------------------------------------------------
 609         IMPLICIT None
 610 #include "EvtGenModels/EvtBTo3pi.inc"
 611         Complex *16 BrightWagner,EvtRBW
 612         Integer ierr
 613         Logical RelatBW
 614         Data    RelatBW/.true./
 615         Real*8  Mass,Width,Am2Min
 616
 617 c       ---------------------------------------------------------------
 618 c       ---     Input Four Vectors
 619 c       ---------------------------------------------------------------
 620         Real*8     p1(5),p2(5),p3(5)
 621
 622 c       ---------------------------------------------------------------
 623 C       ---     intermediate variables
 624 c       ---------------------------------------------------------------
 625         Real*8 E12_2,m12_2,beta,gamma,argu,m13_2,costet,coscms,m12
 626         Real*8 Factor,Num_real,Num_imag
 627         Integer i
 628         Real *8 p1z,p1zcms12,e1cms12,p1cms12
 629
 630 c       ---------------------------------------------------------------
 631 C       ---     Boost factor
 632 c       ---------------------------------------------------------------
 633
 634         BrightWagner=Dcmplx(0,0)
 635
 636         ierr  = 0
 637                 E12_2=(p1(4)+p2(4))**2
 638         m12_2=E12_2
 639         Do i=1,3
 640                 m12_2=m12_2-(p1(i)+p2(i))**2
 641         End Do
 642         Argu = 1.D+00 - m12_2 / E12_2
 643
 644         If(argu.gt.0.) Then
 645                 beta = Dsqrt(Argu)
 646         Else
 647                 Print *,'Abnormal beta ! Argu  = ',Argu
 648
 649                 Argu = 0.
 650                 Beta = 0.
 651         End If
 652
 653         If(m12_2.gt.0.)Then
 654                 m12     = Dsqrt(m12_2)
 655         Else
 656                 Print *,'Abnormal m12  ! m12_2 = ',m12_2
 657                 Print*,'p1 = ',p1
 658                 Print*,'p2 = ',p2
 659                 Print*,'p3 = ',p3
 660                 Stop
 661         End if
 662
 663         gamma=Dsqrt(E12_2/m12_2)
 664 c       ---------------------------------------------------------------
 665 C       ---     get the cosine in the B CMS
 666 c       ---------------------------------------------------------------
 667
 668         m13_2=(p1(4)+p3(4))**2
 669         Do i=1,3
 670                 m13_2=m13_2-(p1(i)+p3(i))**2
 671         End Do
 672         if(m13_2.lt.0)            Go To 50
 673         if((p1(4)**2-p1(5)).lt.0) Go To 50
 674         if((p3(4)**2-p3(5)).lt.0) Go To 50
 675         costet= (2.D+00*p1(4)*p3(4)-m13_2+p1(5)+p3(5))
 676      >      /
 677      >        (2.D+00*Dsqrt( (p1(4)**2-p1(5)) * (p3(4)**2-p3(5)) ))
 678
 679 c       ---------------------------------------------------------------
 680 C       ---     get the costet in the 1-2 CMS
 681 c       ---------------------------------------------------------------
 682
 683         p1z=dsqrt(P1(4)**2-p1(5))*costet
 684         p1zcms12=gamma*(p1z+beta*P1(4))
 685         e1cms12 =gamma*(p1(4)+beta*p1z)
 686         p1cms12 =Dsqrt(e1cms12**2-p1(5))
 687         coscms=p1zcms12/p1cms12
 688
 689 c       ---------------------------------------------------------------
 690 C       ---     Build the Breit Wigner
 691 c       ---------------------------------------------------------------
 692
 693         If(RelatBW) then
 694                 Am2Min       = p1(5) + p2(5) + 2.*Dsqrt( p1(5)*p2(5) )
 695                 BrightWagner = coscms * EvtRBW(m12_2,Mass**2,Width,Am2Min)
 696         Else
 697          Factor  = 2.D+00* ( (Mass-m12)**2+(0.5D+00*Width)**2 )
 698          Factor  = coscms * Width / Factor
 699          Num_real=  (Mass-m12)                    * Factor
 700          Num_imag=                  0.5D+00*Width * Factor
 701          BrightWagner=Dcmplx(Num_real,Num_imag)
 702         End if
 703
 704         Return
 705  50     continue
 706         ierr = 2
 707         Return
 708         End
 709
 710
 711       FUNCTION EvtRBW(s,Am2,Gam,Am2Min)
 712
 713       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 714       COMPLEX*16 EvtRBW
 715
 716       EvtRBW = 0.
 717       IF (s.le.Am2Min) RETURN
 718
 719       G  =  Gam* (Am2/s) * ((s-Am2Min)/(Am2-Am2Min))**1.5
 720       D  =  (Am2-s)**2 + s*G**2
 721       X  =  Am2*(Am2-s)
 722       Y  =  Am2*SQRT(s)*G
 723
 724       EvtRBW = DCMPLX(X/D,Y/D)
 725
 726       RETURN
 727       END