TEvtGen/EvtGenModels/EvtBTo3pi.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: EvtBTo3pi.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 B -->-- 3pions decay generator
  22 C Its is composed of the following subroutines:
  23 C
  24 C [*] EVTHowToUse
  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] EVT3pions
  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 EVT3pions 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
  43 C                 Generation of the kinematics of the 3 pions
  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     -->-- 3pions
  49 C                                   and of the B0_bar -->-- 3pions
  50 C                 corrected for the generation mechanism
  51 C                 The notations used are the ones of the paper of
  52 C                 A. Snyder and H. Quinn [Phys.Rev.D48 (1993) 2139]
  53 C [3] BreitWigner
  54 C                 compute the Breit-Wigner of the contributing rho s
  55 C                 taking into account the cosine term linked to the
  56 C                 zero-helicity of the rho s. There is three forms of
  57 C                 Breit-Wigners available. The first one is the simple
  58 C                 non-relativistic form, while the second and third
  59 C                 ones are more involved relativistic expressions which
  60 C                 in addition incorporate the contributions of three
  61 C                 rho resonances [rho(770:1450:1700)]. The parameters
  62 C                 used are the ones resulting from the ALEPH analysis
  63 C                 which might be found in CERN-PPE:97013 (submitted to
  64 C                 Zeitschrift für Physik C). The two parametrizations
  65 C                 of the relativistic Breit-Wigners are the ones of
  66 C                 Kuhn-SantaMaria [default] and of Gounaris-Sakurai.
  67 C                 The default setting is the non-relativistic Breit-
  68 C                 Wigner form.
  69 C [4] Set_constants
  70 C                 Set the constants values, from the pion mass to the
  71 C                 penguin contributions. It is called by EVT3pions
  72 C
  73 C And other routines which do not deserve comment here.
  74 C===================================================================
  75 c
  76 c      call EvtHowToUse
  77 c      Stop
  78 c      End
  79
  80       subroutine EvtHowToUse
  81
  82       Implicit none
  83       Real*8  alphaCP
  84       Integer iset,number,i,j,jset,N_gener,N_asked
  85       Real*8  p_pi_plus(4),p_pi_minus(4),p_gamma_1(4),p_gamma_2(4)
  86       Real*8  q_pi_plus(4),q_pi_minus(4),q_gamma_1(4),q_gamma_2(4)
  87       Real*8  CPweight
  88       Real*8  Real_B0,Imag_B0,Real_B0bar,Imag_B0bar
  89       Real*8  t3pions,ttag,t,Weight,Weight_max
  90       Real*8  AB0,AB0bar,Ainter,xd
  91       Real*8  m_rhop_2,m_rhom_2
  92       Real*4  Evtranf,Tag
  93       Real*8  Kin_moy,Kin_moy2,Kin,Lambda
  94
  95       alphaCP      = 0.4
  96       Lambda     = Dsin(2.*alphaCP)
  97       Kin_moy    =0.
  98       Kin_moy2   =0.
  99       N_gener    =0
 100       N_asked    =100000
 101
 102       weight_max = 1.0
 103
 104 c xd is not needed by the generator but by the time-dependant
 105 C sector of the general BaBar Generator.
 106       xd         = 0.71
 107
 108 c run : Simulation of the Anders Ryd Generator
 109
 110       Do number=1,N_asked ! weighted events as generated here
 111
 112       If(number.eq.1) then
 113       iset=10000          ! 10^4 events are used to normalize amplitudes
 114       Else
 115       iset=0              ! iset must be reset to zero after the first call
 116       End If
 117
 118 c Here is the call to EVT3pions !!!!!!!!!!!!!!!!!!!!!!!!!
 119 c communication of data is done by argument only <<<<<<<<
 120        call EVT3pions(
 121      + alphaCP,iset,
 122      + p_pi_plus,p_pi_minus,p_gamma_1,p_gamma_2,
 123      + Real_B0,Imag_B0,Real_B0bar,Imag_B0bar)
 124
 125 C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 126 c what follows is just to show how the ouput of EVT3pions
 127 C should be used at the time dependant level.
 128 c Decay time of the B->-3 pions
 129  1    Continue
 130       t3pions=-Alog(Evtranf())
 131       If(t3pions.gt.10.) Go to 1
 132 c Decay time of the tagging B
 133  2    Continue
 134       ttag   =-Alog(Evtranf())
 135       If(ttag   .gt.10.) Go to 2
 136 c time-difference between the 3pi decay and the tagging decay
 137       t=t3pions-ttag
 138 c select the Tag
 139       Tag   =evtranf()
 140 C get the relevant quantities
 141       AB0   = Real_B0   **2 + Imag_B0   **2
 142       AB0bar= Real_B0bar**2 + Imag_B0bar**2
 143       Ainter= Real_B0*Imag_B0bar - Imag_B0*Real_B0bar
 144 c generate acording to the tag
 145       If(Tag.gt.0.5) Then
 146
 147 c a B0bar tag => at t3pions=ttag the decay is one from a B0
 148       Weight= AB0   *Dcos(xd*t/2.D+00)**2
 149      +      + AB0bar*Dsin(xd*t/2.D+00)**2
 150      +      - Ainter*Dsin(xd*t       )
 151
 152            Else
 153
 154 c a B0    tag => at t3pions=ttag the decay is one from a B0bar
 155       Weight= AB0bar*Dcos(xd*t/2.D+00)**2
 156      +      + AB0   *Dsin(xd*t/2.D+00)**2
 157      +      + Ainter*Dsin(xd*t       )
 158
 159       End If
 160 c       Print*,'weight',weight
 161       If(Weight.Gt.evtranf()) Then
 162 c       Print*,'unweighted events'
 163 c----------------------------------------------------------------------------
 164 c unweighted event production
 165 c----------------------------------------------------------------------------
 166       N_gener=N_gener+1
 167 C here is just a Dalitz plot and a few prints
 168
 169       m_rhop_2=(p_pi_plus (4)+p_gamma_1(4)+p_gamma_2(4))**2
 170       m_rhom_2=(p_pi_minus(4)+p_gamma_1(4)+p_gamma_2(4))**2
 171       do j=1,3
 172       m_rhop_2=m_rhop_2-(p_pi_plus (j)+p_gamma_1(j)+p_gamma_2(j))**2
 173       m_rhom_2=m_rhom_2-(p_pi_minus(j)+p_gamma_1(j)+p_gamma_2(j))**2
 174       end do
 175
 176 c here is a check that weight_max is one
 177
 178       If(Weight.gt.Weight_max) Then
 179       Weight_max=Weight
 180       Print*,' overweighted event found at weight = ',Weight_max
 181       End If
 182
 183 c here is an example of use of the iset<0 option.
 184 c One does not want to generate a new event but just to
 185 c compute the matrix element of the CP conjugate event:
 186 c the example provided below uses the Kin variable, but
 187 c this kind of facility will be very usefull for other
 188 c types of applications. For example, to deal with events
 189 c in a likelihood analysis, one needs to compute matrix
 190 c elements, and it would be nice to be able to do that from
 191 c the BaBar generator itself, to be sure that the Physics
 192 c is the same.
 193
 194       q_pi_plus (4)=   p_pi_minus(4)
 195       q_pi_minus(4)=   p_pi_plus (4)
 196       q_gamma_1 (4)=   p_gamma_1 (4)
 197       q_gamma_2 (4)=   p_gamma_2 (4)
 198
 199       do i=1,3
 200       q_pi_plus (i)= - p_pi_minus(i)
 201       q_pi_minus(i)= - p_pi_plus (i)
 202       q_gamma_1 (i)= - p_gamma_1 (i)
 203       q_gamma_2 (i)= - p_gamma_2 (i)
 204       end do
 205
 206 c Here is the call to EVT3pions !!!!!!!!!!!!!!!!!!!!!!!!!
 207 c with the iset<0 option
 208
 209        jset=-1
 210
 211        call EVT3pions(
 212      + alphaCP,jset,
 213      + q_pi_plus,q_pi_minus,q_gamma_1,q_gamma_2,
 214      + Real_B0,Imag_B0,Real_B0bar,Imag_B0bar)
 215 c that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 216 C get the relevant quantities
 217       AB0   = Real_B0   **2 + Imag_B0   **2
 218       AB0bar= Real_B0bar**2 + Imag_B0bar**2
 219       Ainter= Real_B0*Imag_B0bar - Imag_B0*Real_B0bar
 220
 221       If(Tag.lt.0.5) Then
 222
 223 c a B0bar tag => at t3pions=ttag the decay is one from a B0
 224       CPWeight= AB0   *Dcos(xd*t/2.D+00)**2
 225      +      + AB0bar*Dsin(xd*t/2.D+00)**2
 226      +      - Ainter*Dsin(xd*t       )
 227
 228            Else
 229
 230 c a B0    tag => at t3pions=ttag the decay is one from a B0bar
 231       CPWeight= AB0bar*Dcos(xd*t/2.D+00)**2
 232      +      + AB0   *Dsin(xd*t/2.D+00)**2
 233      +      + Ainter*Dsin(xd*t       )
 234
 235       End If
 236
 237       Kin = (Weight-CPweight)/(Weight+CPweight)/Lambda
 238       Kin_moy  = Kin_moy  + Kin
 239       Kin_moy2 = Kin_moy2 + Kin**2
 240
 241 c----------------------------------------------------------------------------
 242       End If
 243
 244 c end of the loop over events
 245       End Do
 246
 247       Print*,'number of unity-weight events generated : ',N_gener
 248       Print*,'number of trials                        : ',N_asked
 249       Print*,' measurement of sin(2 alpha)            = ',
 250      +      Kin_moy/Kin_moy2
 251       Print*,' true value                             = ',Lambda
 252
 253       End
 254 C===================================================================
 255       subroutine Evt3pions(
 256      + alpha_input,iset,
 257      + p_pi_plus,p_pi_minus,p_gamma_1,p_gamma_2,
 258      + Real_B0,Imag_B0,Real_B0bar,Imag_B0bar)
 259 c-----------------------------------------------------------------
 260 c       ----------------------------------------------------------
 261 c       --- This is the routine to be called by the Main generator
 262 c****************************************************************************
 263 c       --- Inputs are  :
 264 c
 265 c       ---               alpha     : the value of the alpha angle (radians)
 266 c             Note that the amplitude which are returned incorporate the
 267 c             CP violating phase present in the mixing.
 268 c
 269 c       ---               iset      : to initialize the generator
 270 c                                     iset should be set to the number
 271 c                                     of events to be used for the pre-run
 272 c                                     devoted to the computation of the
 273 c                                     normalization factor designed to help
 274 c                                     the subsequent time-dependent generation
 275 c                                     of events.
 276 c                                     After the first call iset should be
 277 c                                     set to 0.
 278 c****************************************************************************                                         set to 0.
 279 c       --- Outputs are :
 280 c
 281 c       ---               p_pi_plus : the four momentum of the pi+
 282 c       ---               p_pi_minus: the four momentum of the pi-
 283 c       ---               p_gamma_1 : the four momentum of the first  photon
 284 c       ---               p_gamma_2 : the four momentum of the second photon
 285 c             Note that : the energy is stored in the fourth component
 286 c                         the values are the ones of the B rest frame
 287 c                         a random rotation has been applied
 288 c
 289 c       ---               Real_B0   : The real      part of the amplitude of
 290 c                                    the B0    ->- 3 pions decay
 291 c       ---               Imag_B0   : The imaginary part of the amplitude of
 292 c                                    the B0    ->- 3 pions decay
 293 c       ---               Real_B0bar: The real      part of the amplitude of
 294 c                                    the B0bar ->- 3 pions decay
 295 c       ---               Imag_B0bar: The imaginary part of the amplitude of
 296 c                                    the B0bar ->- 3 pions decay
 297 c****************************************************************************
 298 c-----------------------------------------------------------------
 299       Implicit none
 300 #include "EvtGenModels/EvtBTo3pi.inc"
 301       Real*8  alpha_input
 302       Real*8  beta_input
 303       Integer iset
 304       Real*8  p_pi_plus(4),p_pi_minus(4),p_gamma_1(4),p_gamma_2(4)
 305       Real*8  Real_B0,Imag_B0,Real_B0bar,Imag_B0bar
 306
 307 c Working quantities
 308       Integer i,number
 309       Real*8  p1(5),p2(5),p3(5),Gamma1(5),Gamma2(5)
 310       Real*8 factor,factor_max,AB0,AB0bar,Ainter,R1,R2
 311       data factor_max/1.D+00/
 312       Integer ierr
 313       ierr =0
 314 c-------------------------------------------------------------------
 315       If(iset.eq.0) Then
 316 c-------------------------------------------------------------------
 317 c     this is the normal mode of operation
 318 c     First, generate the kinematics
 319
 320       p1(5)=M_pip**2
 321       p2(5)=M_pi0**2
 322       p3(5)=M_pim**2
 323
 324  10   continue
 325
 326
 327       call Evtfirst_step(p1,p2,p3)
 328
 329 c     Then, compute the amplitudes
 330
 331
 332       Call EvtCompute(p1,p2,p3,
 333      +    Real_B0,Imag_B0,Real_B0bar,Imag_B0bar,iset,ierr)
 334       if(ierr.ne.0)  Go To 10
 335
 336 c-------------------------------------------------------------------
 337       ElseIf(iset.lt.0) Then
 338 c-------------------------------------------------------------------
 339 c     This is an user mode of operation where the kinematics is
 340 c     provided by the user who only wants the corresponding amplitudes
 341 c     to be computed
 342
 343       Do i=1,4
 344       p1(i)= p_pi_plus (i)
 345       p2(i)= p_gamma_1 (i) + p_gamma_2 (i)
 346       p3(i)= p_pi_minus(i)
 347       End Do
 348       p1(5)= M_pip**2
 349       p2(5)= M_pi0**2
 350       p3(5)= M_pim**2
 351
 352         Call EvtCompute(p1,p2,p3,
 353      +          Real_B0,Imag_B0,Real_B0bar,Imag_B0bar,iset,ierr)
 354
 355       if(ierr.ne.0) Then
 356         Print*,'the provided kinematics are not physical'
 357         Print*,'ierr=',ierr
 358         Print*,'the program will stop'
 359         Stop
 360       endif
 361 c-------------------------------------------------------------------
 362       ElseIf(iset.gt.0) Then
 363 c-------------------------------------------------------------------
 364 c     This is the pre-run mode of operation where initializations are
 365 c     performed.
 366
 367       factor_max= 0
 368
 369 C Ghm : beta is not needed for this generation - put a default value
 370       beta_input = 0.362
 371 C
 372       call Evtset_constants(alpha_input, beta_input)
 373       p1(5)=M_pip**2
 374       p2(5)=M_pi0**2
 375       p3(5)=M_pim**2
 376
 377 c     pre-run
 378       Do number=1,iset
 379
 380  20   continue
 381
 382       call Evtfirst_step(p1,p2,p3)
 383
 384       Call EvtCompute(p1,p2,p3,Real_B0,Imag_B0,Real_B0bar,
 385      +  Imag_B0bar,iset,ierr)
 386          if(ierr.ne.0) Go To 20
 387       AB0   = Real_B0   **2 + Imag_B0   **2
 388       AB0bar= Real_B0bar**2 + Imag_B0bar**2
 389       Ainter= Real_B0*Imag_B0bar - Imag_B0*Real_B0bar
 390       R1    =  (AB0-AB0bar) /(AB0+AB0bar)
 391       R2    =(2.D+00*Ainter)/(AB0+AB0bar)
 392       factor=(1.D+00+Dsqrt(R1**2+R2**2))*(AB0+AB0bar)/2.D+00
 393       If(factor.gt.factor_max) factor_max=factor
 394
 395       End Do
 396 c     end of the pre-run
 397 c     Normalization factor
 398       factor_max=1.D+00/Dsqrt(factor_max)
 399 c       Print*,'factor_max',factor_max,p1,p2,p3
 400 c-------------------------------------------------------------------
 401       End If
 402 c-------------------------------------------------------------------
 403
 404       Real_B0   =Real_B0    * factor_max
 405       Imag_B0   =Imag_B0    * factor_max
 406       Real_B0bar=Real_B0bar * factor_max
 407       Imag_B0bar=Imag_B0bar * factor_max
 408
 409       if(iset.lt.0) return
 410
 411 c     P1,p2,p3 ---> random rotation in B rest frame
 412
 413       Call EvtRotation(p1,1)
 414       Call EvtRotation(p2,0)
 415       Call EvtRotation(p3,0)
 416
 417 C     Desintegrate the pi_0
 418
 419       Call EvtGammaGamma(p2,Gamma1,Gamma2)
 420
 421 C     Feed the output four vectors
 422
 423       Do i=1,4
 424
 425       p_pi_plus (i)=p1    (i)
 426       p_pi_minus(i)=p3    (i)
 427       p_gamma_1 (i)=Gamma1(i)
 428       p_gamma_2 (i)=Gamma2(i)
 429
 430       End Do
 431
 432       Return
 433
 434       End
 435
 436 c===================================================================
 437       subroutine Evtfirst_step(P1,P2,P3)
 438 c-----------------------------------------------------------------
 439 c       ----------------------------------------------------------
 440 c       --- This routine generates the 5-vectors P1,P2,P3
 441 c       --- Associated respectively with the Pi+,Pi0 and Pi-
 442 c       ---             P1(1) = Px
 443 c       ---             P1(2) = Py
 444 c       ---             P1(3) = Pz
 445 c       ---             P1(4) = E
 446 c       ---             P1(5) = M**2
 447 c       ----------------------------------------------------------
 448 c       ---     Input Four Vectors
 449 C       ---     Particle [1] is the pi+
 450 C       ---     Particle [2] is the pi0
 451 C       ---     Particle [3] is the pi-
 452 c       ----------------------------------------------------------
 453
 454 c       ----------------------------------------------------------
 455 c       --- commons
 456 c       ----------------------------------------------------------
 457 #include "EvtGenModels/EvtBTo3pi.inc"
 458 c       ----------------------------------------------------------
 459 c       --- Used Functions
 460 c       ----------------------------------------------------------
 461
 462         real evtranf
 463
 464 c       ----------------------------------------------------------
 465 c       --- Variables in Argument
 466 c       ----------------------------------------------------------
 467
 468          real*8 P1(5),P2(5),P3(5)
 469
 470 c       ----------------------------------------------------------
 471 c       --- Local Variables
 472 c       ----------------------------------------------------------
 473
 474         real*8 m12,min_m12, max_m12
 475         real*8 m13,min_m13, max_m13
 476         real*8 m23,min_m23, max_m23
 477         Real*8 cost13,cost12,cost23
 478
 479         real*8 p1mom,p2mom,p3mom
 480
 481         real*8 x, y, z, mass
 482         integer i
 483 C There is two ways of generating the events:
 484 C The first one used a pole-compensation method to generate the
 485 C events efficiently taking into account the poles due to the
 486 C Breit-Wigners of the rho s. It is activated by setting
 487 C Phase_Space to .false.
 488 C The second one generates events according to phase space. It is
 489 C inneficient but allows the exploration of the full Dalitz plot
 490 C in an uniform way. It was found to be usefull fopr some peculiar
 491 C applications. It is activated by setting
 492 C Phase_space to .true.
 493 C Note that in that case, the generation is no longer correct.
 494
 495         Logical Phase_space
 496         data Phase_space/.false./
 497
 498 c       initialize to avoid warning on Linux.
 499         mass=0.0
 500
 501 c       ----------------------------------------------------------
 502 c       --- Computation
 503 c       ----------------------------------------------------------
 504         max_m12 = M_B**2
 505         min_m12 = P1(5) + P2(5)
 506
 507         max_m13 = M_B**2
 508         min_m13 = P1(5) + P3(5)
 509
 510         max_m23 = M_B**2
 511         min_m23 = P2(5) + P3(5)
 512
 513
 514 100     Continue
 515
 516 c       ----------------------------------------------------------
 517 c       --- Generation of the Mass of the Rho(+,-,0)
 518 c       ----------------------------------------------------------
 519
 520         If(.not.Phase_space) Then
 521         y = evtranf()*PI - PI/2.
 522         x = Dtan(y)
 523         mass = x*Gam_rho/2. +Mass_rho
 524         End If
 525
 526 c       ----------------------------------------------------------
 527 c       --- z is the Flag needed to choose between the generation
 528 c       --- of a Rho+, a Rho-,or a Rho0
 529 c       ----------------------------------------------------------
 530
 531         z = 3.* evtranf()
 532
 533         if(z.lt.1.) Then
 534
 535         If(Phase_space) Then
 536                 m12 = evtranf()*(max_m12-min_m12)+min_m12
 537         Else
 538                 m12 = mass**2
 539         End If
 540                 m13 = evtranf()*(max_m13-min_m13)+min_m13
 541                 m23 = MB2 - m12 - m13
 542         else if (z.lt.2.) Then
 543
 544         If(Phase_space) Then
 545                 m13 = evtranf()*(max_m13-min_m13)+min_m13
 546         Else
 547                 m13 =  mass**2
 548         End If
 549                 m12 = evtranf()*(max_m12-min_m12)+min_m12
 550                 m23 = MB2 - m12 - m13
 551         else
 552
 553         If(Phase_space) Then
 554                 m23 = evtranf()*(max_m23-min_m23)+min_m23
 555         Else
 556                 m23 = mass**2
 557         End If
 558                 m12 = evtranf()*(max_m12-min_m12)+min_m12
 559                 m13 = MB2 - m12 - m23
 560         endif
 561
 562 c       ----------------------------------------------------------
 563 c       --- Check that the physics is OK :
 564 c       --- Are the invariant Masses in allowed ranges ?
 565 c       ----------------------------------------------------------
 566
 567         If(m23.lt.min_m23.or.m23.gt.max_m23) Go to 100
 568         If(m13.lt.min_m13.or.m13.gt.max_m13) Go to 100
 569         If(m12.lt.min_m12.or.m12.gt.max_m12) Go to 100
 570
 571 c       ----------------------------------------------------------
 572 c       --- Are the Cosines of the angles between particles
 573 c       --- Between -1 and +1 ?
 574 c       ----------------------------------------------------------
 575
 576         P1(4)=(M_B**2+P1(5)-m23)/(2.*M_B)
 577         P2(4)=(M_B**2+P2(5)-m13)/(2.*M_B)
 578         P3(4)=(M_B**2+P3(5)-m12)/(2.*M_B)
 579         if (p1(4)**2-P1(5).le.0.0) go to 100
 580         if (p2(4)**2-P2(5).le.0.0) go to 100
 581         if (p3(4)**2-P3(5).le.0.0) go to 100
 582         p1mom=Dsqrt(p1(4)**2-P1(5))
 583         p2mom=Dsqrt(p2(4)**2-P2(5))
 584         p3mom=Dsqrt(p3(4)**2-P3(5))
 585         cost13=(2.*p1(4)*p3(4)+P1(5)+p3(5)-m13)/(2.*p1mom*p3mom)
 586         cost12=(2.*p1(4)*p2(4)+P1(5)+p2(5)-m12)/(2.*p1mom*p2mom)
 587         cost23=(2.*p2(4)*p3(4)+P2(5)+p3(5)-m23)/(2.*p2mom*p3mom)
 588         If(Dabs(cost13).gt.1.) Go to 100
 589         If(Dabs(cost12).gt.1.) Go to 100
 590         If(Dabs(cost23).gt.1.) Go to 100
 591
 592 c       ----------------------------------------------------------
 593 c       --- Filling the 5-vectors P1,P2,P3
 594 c       ----------------------------------------------------------
 595
 596         P3(1) = 0
 597         P3(2) = 0
 598         p3(3) = p3mom
 599
 600         P1(3) = p1mom*cost13
 601         P1(1) = p1mom*Dsqrt(1.D+00-cost13**2)
 602         p1(2) = 0.
 603
 604         Do i=1,3
 605         P2(i)=-p1(i)-p3(i)
 606         End do
 607
 608         END
 609
 610
 611 c======================================================================
 612       Subroutine EvtCompute(p1,p2,p3,
 613      +           Real_B0,Imag_B0,Real_B0bar,Imag_B0bar,iset,ierr)
 614 c-----------------------------------------------------------------------
 615         IMPLICIT None
 616 #include "EvtGenModels/EvtBTo3pi.inc"
 617
 618
 619         Real*8 m12, m13, m23, W12, W13, W23, Wtot
 620         Real*4 evt_gmas
 621         Complex*16 MatBp,MatBm
 622         Real*8  Real_B0,Imag_B0,Real_B0bar,Imag_B0bar
 623         Integer Iset,ierr
 624         Real*8     p1(5),p2(5),p3(5)
 625
 626         COMPLEX*16 Mat_rhop,Mat_rhom,Mat_rho0
 627         Complex*16 BreitWigner,Mat_1,Mat_2,Mat_3
 628 c       ----------------------------------------------------------------
 629 C       ---     Account for the pole compensation
 630 c       ----------------------------------------------------------------
 631         ierr = 0
 632
 633         if(evt_gmas(p1,p2).lt.0.or.evt_gmas(p1,p3)
 634      +          .lt.0.or.evt_gmas(p2,p3).lt.0) Then
 635                 ierr=1
 636                 return
 637         endif
 638
 639         m12 = sqrt(evt_gmas(p1,p2))
 640         m13 = sqrt(evt_gmas(p1,p3))
 641         m23 = sqrt(evt_gmas(p2,p3))
 642
 643         W12 = (1./m12)*1./((Mass_rho - m12)**2+(Gam_rho/2.)**2)
 644         W13 = (1./m13)*1./((Mass_rho - m13)**2+(Gam_rho/2.)**2)
 645         W23 = (1./m23)*1./((Mass_rho - m23)**2+(Gam_rho/2.)**2)
 646
 647         If(iset.ge.0) Then
 648         Wtot = 1./Dsqrt(W12 + W13 + W23)
 649         else
 650         Wtot = 1.
 651         endif
 652
 653 c       ----------------------------------------------------------------
 654 C       ---     Compute Breit-Wigners
 655 c       ----------------------------------------------------------------
 656
 657         Mat_rhop=BreitWigner(p1,p2,p3,ierr)
 658         Mat_rhom=BreitWigner(p2,p3,p1,ierr)
 659         Mat_rho0=BreitWigner(p1,p3,p2,ierr)
 660
 661
 662 c       ----------------------------------------------------------------
 663 C       ---     Build up the amplitudes
 664 c       ----------------------------------------------------------------
 665
 666         Mat_1 =  Mat_s3*Mat_rhop
 667         Mat_2 =  Mat_s4*Mat_rhom
 668         Mat_3 =  Mat_s5*Mat_rho0*0.5D+00
 669
 670         MatBp    = (Mat_1+Mat_2+Mat_3) * Wtot
 671
 672         Mat_1 =  Nat_s3*Mat_rhom
 673         Mat_2 =  Nat_s4*Mat_rhop
 674         Mat_3 =  Nat_s5*Mat_rho0*0.5D+00
 675
 676         MatBm    = (Mat_1+Mat_2+Mat_3) * Wtot
 677
 678 c       Pick up the Real and Imaginary parts
 679
 680 c       changed Real to DBLE and Imag DIMAG to make it complie
 681 c       with Absoft and g77 (ryd)
 682
 683         Real_B0    = DBLE(MatBp)
 684         Imag_B0    = DIMAG(MatBp)
 685
 686         Real_B0bar = DBLE(MatBm)
 687         Imag_B0bar = DIMAG(MatBm)
 688
 689         Return
 690         End
 691
 692 c======================================================================
 693       Function BreitWigner(p1,p2,p3,ierr)
 694 c----------------------------------------------------------------------
 695         IMPLICIT None
 696 #include "EvtGenModels/EvtBTo3pi.inc"
 697         Complex *16 BreitWigner,EvtCRhoF_W
 698         Logical Aleph
 699         Integer ierr
 700         Data Aleph /.true./
 701
 702 c       ---------------------------------------------------------------
 703 c       ---     Input Four Vectors
 704 c       ---------------------------------------------------------------
 705         Real*8     p1(5),p2(5),p3(5)
 706
 707 c       ---------------------------------------------------------------
 708 C       ---     intermediate variables
 709 c       ---------------------------------------------------------------
 710         Real*8 E12_2,m12_2,beta,gamma,argu,m13_2,costet,coscms,m12
 711         Real*8 Factor,Num_real,Num_imag
 712         Integer i
 713         Real *8 p1z,p1zcms12,e1cms12,p1cms12
 714
 715 c       ---------------------------------------------------------------
 716 C       ---     Boost factor
 717 c       ---------------------------------------------------------------
 718
 719         BreitWigner=Dcmplx(0,0)
 720
 721
 722         ierr  = 0
 723                 E12_2=(p1(4)+p2(4))**2
 724         m12_2=E12_2
 725         Do i=1,3
 726                 m12_2=m12_2-(p1(i)+p2(i))**2
 727         End Do
 728         Argu = 1.D+00 - m12_2 / E12_2
 729
 730         If(argu.gt.0.) Then
 731                 beta = Dsqrt(Argu)
 732         Else
 733                 Print *,'Abnormal beta ! Argu  = ',Argu
 734
 735                 Argu = 0.
 736                 Beta = 0.
 737         End If
 738
 739         If(m12_2.gt.0.)Then
 740                 m12     = Dsqrt(m12_2)
 741         Else
 742                 Print *,'Abnormal m12  ! m12_2 = ',m12_2
 743                 Print*,'p1 = ',p1
 744                 Print*,'p2 = ',p2
 745                 Print*,'p3 = ',p3
 746                 Stop
 747         End if
 748
 749         gamma=Dsqrt(E12_2/m12_2)
 750 c       ---------------------------------------------------------------
 751 C       ---     get the cosine in the B CMS
 752 c       ---------------------------------------------------------------
 753
 754         m13_2=(p1(4)+p3(4))**2
 755         Do i=1,3
 756                 m13_2=m13_2-(p1(i)+p3(i))**2
 757         End Do
 758         if(m13_2.lt.0)            Go To 50
 759         if((p1(4)**2-p1(5)).lt.0) Go To 50
 760         if((p3(4)**2-p3(5)).lt.0) Go To 50
 761         costet= (2.D+00*p1(4)*p3(4)-m13_2+p1(5)+p3(5))
 762      >      /
 763      >        (2.D+00*Dsqrt( (p1(4)**2-p1(5)) * (p3(4)**2-p3(5)) ))
 764
 765 c       ---------------------------------------------------------------
 766 C       ---     get the costet in the 1-2 CMS
 767 c       ---------------------------------------------------------------
 768
 769         p1z=dsqrt(P1(4)**2-p1(5))*costet
 770         p1zcms12=gamma*(p1z+beta*P1(4))
 771         e1cms12 =gamma*(p1(4)+beta*p1z)
 772         p1cms12 =Dsqrt(e1cms12**2-p1(5))
 773         coscms=p1zcms12/p1cms12
 774
 775 c       ---------------------------------------------------------------
 776 C       ---     Build the Breit Wigner
 777 c       ---------------------------------------------------------------
 778
 779         If(aleph) then
 780                 BreitWigner= coscms * EvtCRhoF_W( m12_2 )
 781         else
 782          Factor  = 2.D+00* ( (Mass_rho-m12)**2+(0.5D+00*Gam_rho)**2 )
 783          Factor  = coscms * Gam_rho / Factor
 784          Num_real=  (Mass_rho-m12)                    * Factor
 785          Num_imag=                    0.5D+00*Gam_rho * Factor
 786          BreitWigner=Dcmplx(Num_real,Num_imag)
 787         End if
 788
 789         Return
 790  50     continue
 791         ierr = 2
 792         Return
 793         End
 794
 795 c======================================================================
 796       Subroutine EvtSet_Constants(balpha, bbeta)
 797 c----------------------------------------------------------------------
 798       IMPLICIT None
 799 #include "EvtGenModels/EvtBTo3pi.inc"
 800
 801 c       -------------------------------------------------------------
 802 c       ---     Basic matrix elements
 803 c       -------------------------------------------------------------
 804
 805         Complex *16   Mat_Tp0,Mat_Tm0,Mat_T0p,Mat_T0m,Mat_Tpm,Mat_Tmp
 806      >             ,Mat_P1,Mat_P0
 807      >             ,Nat_Tp0,Nat_Tm0,Nat_T0p,Nat_T0m,Nat_Tpm,Nat_Tmp
 808      >             ,Nat_P1,Nat_P0, Mat_Ppm, Mat_Pmp
 809      >             ,StrongExp
 810       Real*8        balpha, bbeta
 811       Real*8        calpha,salpha
 812       Real*8        cbeta,sbeta
 813       Real*8        StrongPhase
 814       Real*8  Rho_gen,Eta_gen,c_factor
 815 c       -------------------------------------------------------------
 816       pi       = 3.141592653 D+00
 817       alphaCP  = balpha
 818       calpha   = Dcos(alphaCP)
 819       salpha   = Dsin(alphaCP)
 820       betaCP   = bbeta
 821       cbeta    = Dcos(betaCP)
 822       sbeta    = Dsin(betaCP)
 823       M_B      = 5.2794D+00
 824       M_pip    = 0.13957D0
 825       M_pim    = M_pip
 826       M_pi0    = 0.134976D0
 827       M_Kp     = 0.49368
 828       Mass_rho = 0.770D0
 829       Gam_rho  = 0.150D0
 830       Mass_Kstarp = 0.8916
 831       Gam_Kstarp  = 0.0498
 832       Mass_Kstar0 = 0.8961
 833       Gam_Kstar0  = 0.0505
 834
 835       MB2         = M_B**2 + M_pip**2 + M_pim**2 + M_pi0**2
 836       MC2         = M_B**2 + M_Kp **2 + M_pim**2 + M_pi0**2
 837
 838       StrongPhase = 0.
 839
 840 c       -------------------------------------------------------------
 841       StrongExp = Dcmplx(Dcos(StrongPhase),Dsin(StrongPhase))
 842
 843       rho_gen = 0.D+00
 844 c       if(balpha.ne.0) Then
 845 c      eta_gen = cos(balpha)/sin(balpha)
 846 c       else
 847 c      eta_gen =0.
 848 c       endif
 849         eta_gen = 0.5D+00
 850       c_factor=
 851      >Sqrt((rho_gen**2+eta_gen**2)/((1.-rho_gen)**2+eta_gen**2))
 852
 853       Mat_Tp0 = Dcmplx(calpha, -salpha) *1.09D+00
 854       Mat_Tm0 = Dcmplx(calpha, -salpha) *1.09D+00
 855       Mat_T0p = Dcmplx(calpha, -salpha) *0.66D+00
 856       Mat_T0m = Dcmplx(calpha, -salpha) *0.66D+00
 857       Mat_Tpm = Dcmplx(calpha, -salpha) *1.00D+00
 858       Mat_Tmp = Dcmplx(calpha, -salpha) *0.47D+00
 859
 860       Mat_Ppm = -Dcmplx(Dcos(-0.5D0),Dsin(-0.5D0))*0.2D0
 861       Mat_Pmp =  Dcmplx(Dcos(2.D0),Dsin(2.D0))*0.15D0
 862
 863       Mat_P1  = 0.5D0*(Mat_Ppm-Mat_Pmp)
 864       Mat_P0  = 0.5D0*(Mat_Ppm+Mat_Pmp)
 865       Nat_Tp0 = Dcmplx(calpha, salpha) *1.09D+00
 866       Nat_Tm0 = Dcmplx(calpha, salpha) *1.09D+00
 867       Nat_T0p = Dcmplx(calpha, salpha) *0.66D+00
 868       Nat_T0m = Dcmplx(calpha, salpha) *0.66D+00
 869       Nat_Tpm = Dcmplx(calpha, salpha) *1.00D+00
 870       Nat_Tmp = Dcmplx(calpha, salpha) *0.47D+00
 871       Nat_P1  = Mat_P1
 872       Nat_P0  = Mat_P0
 873
 874       Mat_Tpm = StrongExp * Mat_Tpm
 875       Nat_Tpm = StrongExp * Nat_Tpm
 876
 877       Mat_S1  = Mat_Tp0 + 2.D+00 * Mat_P1
 878       Mat_S2  = Mat_T0p - 2.D+00 * Mat_P1
 879       Mat_S3  = Mat_Tpm + Mat_P1 + Mat_P0
 880       Mat_S4  = Mat_Tmp - Mat_P1 + Mat_P0
 881       Mat_S5  = - Mat_Tpm - Mat_Tmp + Mat_Tp0 + Mat_T0p
 882      >        - 2.D+00 * Mat_P0
 883
 884       Nat_S1  = Nat_Tp0 + 2.D+00 * Nat_P1
 885       Nat_S2  = Nat_T0p - 2.D+00 * Nat_P1
 886       Nat_S3  = Nat_Tpm + Nat_P1 + Nat_P0
 887       Nat_S4  = Nat_Tmp - Nat_P1 + Nat_P0
 888       Nat_S5  = - Nat_Tpm - Nat_Tmp + Nat_Tp0 + Nat_T0p
 889      >        - 2.D+00 * Nat_P0
 890
 891
 892 c=====================================================================
 893 c=====================================================================
 894
 895
 896 C Gautier 14-Jan-98 :
 897 C  set matrix elements according to Jerome s calculations
 898 C  (so called Reference set)
 899 c B0    -->-- K*+ pi- Amplitudes (Trees + Penguins)
 900       MatKstarp = Dcmplx(calpha, -salpha) * Dcmplx( 0.220,0.)
 901      >          + Dcmplx(cbeta,   sbeta)  * Dcmplx(-1.200,0.)
 902 c B0    -->-- K*0 pi0 Amplitudes (Trees + Penguins)
 903       MatKstar0 = Dcmplx(calpha, -salpha) * Dcmplx( 0.015,0.)
 904      >          + Dcmplx(cbeta,   sbeta)  * Dcmplx( 0.850,0.)
 905 c B0    -->-- K+ rho- Amplitudes (Trees + Penguins)
 906       MatKrho   = Dcmplx(calpha, -salpha) * Dcmplx( 0.130,0.)
 907      >          + Dcmplx(cbeta,   sbeta)  * Dcmplx( 0.160,0.)
 908 c=====================================================================
 909 c B0bar -->-- K*+ pi- Amplitudes (Trees + Penguins)
 910       NatKstarp = Dcmplx(0.,0.)
 911 c B0bar -->-- K*0 pi0 Amplitudes (Trees + Penguins)
 912       NatKstar0 = Dcmplx(0.,0.)
 913 c B0bar -->-- K+ rho- Amplitudes (Trees + Penguins)
 914       NatKrho   = Dcmplx(0.,0.)
 915 c=====================================================================
 916 c=====================================================================
 917 c       Print*,'Mat_Tp0',Mat_Tp0
 918 c       Print*,'Mat_T0p',Mat_T0p
 919 c       Print*,'Mat_Tpm',Mat_Tpm
 920 c       Print*,'Mat_Tmp',Mat_Tmp
 921 c       Print*,'Mat_P1',Mat_P1
 922 c       Print*,'Mat_P0',Mat_P0
 923
 924       Return
 925       End
 926
 927
 928
 929 C=====================================================================
 930         subroutine EvtGammaGamma(P2,PG1Boost,PG2Boost)
 931 c---------------------------------------------------------------------
 932 c       --------------------------------------------------------------
 933 c       ---     The aim of this routine is to generate two photons
 934 c       ---     for the decay of the pi0
 935 c       --------------------------------------------------------------
 936         implicit none
 937 c --- commons
 938
 939 #include "EvtGenModels/EvtBTo3pi.inc"
 940
 941 c --- Used Functions
 942
 943         Real evtranf
 944
 945 c --- Variables in Argument
 946
 947         Real*8 PG1Boost(5), PG2Boost(5)
 948         Real*8 P2(5)
 949
 950 c --- Local Variables
 951
 952         Real*8 EGammaCmsPi0
 953         Real*8 PGamma1(5), PGamma2(5)
 954         Real*8 sinThetaRot, cosThetaRot, PhiRot
 955         Real*8 Beta(4)
 956         Integer i
 957
 958         Do i=1, 3
 959           Beta(i) = P2(i)/p2(4)
 960         Enddo
 961
 962         EGammaCmsPi0 = Dsqrt(P2(5))/2.
 963
 964
 965         cosThetaRot = evtranf()*2.D+00-1.D+00
 966         sinThetaRot = Dsqrt(1.D+00 - cosThetaRot**2)
 967         PhiRot = evtranf()*2.D+00*pi
 968
 969         PGamma1(1) = EGammaCmsPi0 * sinThetaRot * Dcos(PhiRot)
 970         PGamma1(2) = EGammaCmsPi0 * sinThetaRot * Dsin(PhiRot)
 971         PGamma1(3) = EGammaCmsPi0 * cosThetaRot
 972         PGamma1(4) = EGammaCmsPi0
 973         PGamma1(5) = 0.
 974
 975         PGamma2(1) = - PGamma1(1)
 976         PGamma2(2) = - PGamma1(2)
 977         PGamma2(3) = - PGamma1(3)
 978         PGamma2(4) =   PGamma1(4)
 979         PGamma2(5) = 0.
 980
 981         Call EvtREFREF(BETA,PGamma1,PG1Boost,1)
 982         Call EvtREFREF(BETA,PGamma2,PG2Boost,1)
 983
 984         end
 985
 986 C===================================================================
 987        SUBROUTINE EvtREFREF(BETA,PIN,POUT,S)
 988 c-------------------------------------------------------------------
 989 c       ------------------------------------------------------------
 990 c       ---      GIVEN A FINAL STATE , TO BOOST IT IN ITS C.M.S TAKE
 991 c       ---      BETA  =-SIGMA(P(K))/SIGMA(E(K))
 992 c       ---      S     =+1.
 993 c       ---      TO GO BACK TO THE INITIAL FRAME : SET
 994 c       ---      S     =-1.
 995 c       ------------------------------------------------------------
 996
 997         implicit None
 998         Real*8 Beta,Pin,Pout
 999         Integer s,i
1000         dimension BETA(4),PIN(5),POUT(5)
1001         Real*8 GAMMA,BETPIN,BETA2
1002
1003
1004         BETA2 =BETA(1)**2+BETA(2)**2+BETA(3)**2
1005         IF(BETA2.GT.1.D0) Print*,'Beta**2 .gt.1!!!'
1006         IF(BETA2.GT..99999D0) BETA2=0.99999
1007         BETA(4)=DSQRT(BETA2)
1008         GAMMA=1./DSQRT(1.D0-BETA2)
1009
1010         BetPin = 0.
1011         Do i=1,3
1012                 BetPin = Beta(i)*Pin(i) + BetPin
1013         Enddo
1014         BetPin = s*BetPin
1015
1016         DO 1 I=1,3
1017         POUT(I)=PIN(I)+S*BETA(I)*GAMMA*
1018      >                  (PIN(4)+GAMMA/(GAMMA+1.)*BETPIN)
1019   1     CONTINUE
1020         POUT(4) = GAMMA*(PIN(4)+BETPIN)
1021         POUT(5) = PIN(5)
1022
1023         RETURN
1024         END
1025
1026 C******************************************************************
1027         function evt_gmas(P1,P2)
1028 C******************************************************************
1029 c       -----------------------------------------------------------
1030 c       ---     compute the invariant mass between particle 1 and 2
1031 c       -----------------------------------------------------------
1032         implicit none
1033         real*8 p1(5),P2(5)
1034         real*4 evt_gmas
1035         integer i
1036
1037         evt_gmas = (P1(4)+P2(4))**2
1038
1039         do i=1,3
1040                 evt_gmas = evt_gmas - (P1(i) + P2(i))**2
1041         enddo
1042
1043         end
1044
1045 C===================================================================
1046         Subroutine EvtRotation(p,new)
1047 c------------------------------------------------------------------
1048         Implicit None
1049 #include "EvtGenModels/EvtBTo3pi.inc"
1050
1051         Real*8 p(5),pstor(3)
1052         Real*8 phi2,phi3
1053         Real*8 c1,c2,c3,s1,s2,s3
1054         Real*8 MatRot(3,3)
1055         save MatRot
1056         Real evtranf
1057         Integer new,i,j
1058
1059         If(new.ne.0) Then
1060 c       ---------------------------------------------------------
1061 C       --- generate a random rotation
1062 c       ---------------------------------------------------------
1063
1064                 phi2=evtranf()*2.*pi
1065                 phi3=evtranf()*2.*pi
1066
1067                 c1=2.D+00*evtranf()-1.D+00
1068                 c2=Dcos(phi2)
1069                 c3=Dcos(phi3)
1070
1071                 s1=Dsqrt(1.D+00-c1**2)
1072                 s2=Dsin(phi2)
1073                 s3=Dsin(phi3)
1074
1075 c       ---------------------------------------------------------
1076 C       ---     Compute the overall rotation matrix
1077 c       ---------------------------------------------------------
1078                 MatRot(1,1)=c1
1079                 MatRot(1,2)=s1*c3
1080                 MatRot(1,3)=s1*s3
1081                 MatRot(2,1)=-s1*c2
1082                 MatRot(2,2)=c1*c2*c3-s2*s3
1083                 MatRot(2,3)=c1*c2*s3+s2*c3
1084                 MatRot(3,1)=s1*s2
1085                 MatRot(3,2)=-c1*s2*c3-c2*s3
1086                 MatRot(3,3)=-c1*s2*s3+c2*c3
1087         End If
1088
1089 c       ---------------------------------------------------------
1090 C       ---     store the input 3-momentum
1091 c       ---------------------------------------------------------
1092         Do i=1,3
1093                 pstor(i)=p(i)
1094                 p(i) = 0.
1095         End Do
1096 c       ---------------------------------------------------------
1097 C       ---     Apply the rotation
1098 c       ---------------------------------------------------------
1099         Do i=1,3
1100           Do j=1,3
1101                 p(i)=p(i)+MatRot(i,j)*pstor(j)
1102           End do
1103         End do
1104
1105         end
1106
1107
1108
1109 c======================================================================
1110       FUNCTION EvtCRhoF_W( s )
1111 C -----------------------------------------------------------------------
1112 C
1113 C Author : Andreas Hoecker [Aleph Collaboration, BaBar]
1114 C
1115 C tau --> pi pi0 nu
1116 C
1117
1118       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1119       COMPLEX*16 EvtCRhoF_W, EvtcBW_KS, EvtcBW_GS
1120
1121       PARAMETER (AmPi2=0.13956995**2)
1122
1123       LOGICAL KUHN_SANTA    / .TRUE. / !...type of Breit-Wigner formula
1124       double precision lambda
1125          COMMON /LAM/ lambda
1126 C
1127 C store parameters
1128 C
1129       IF (KUHN_SANTA) THEN
1130 C...rho(770)
1131         AmRho     =  0.7734
1132         GamRho    =  0.1477
1133 C...rho(1450)
1134         AmRhoP    =  1.465
1135         GamRhoP   =  0.696
1136         beta      =  -0.229
1137 C...rho(1700)
1138         AmRhoPP   =  1.760
1139         GamRhoPP  =  0.215
1140         gamma     =  0.075
1141 C...exponent
1142         lambda    =  1.0
1143       ELSE
1144 C...rho(770)
1145         AmRho     =  0.7757
1146         GamRho    =  0.1508
1147 C...rho(1450)
1148         AmRhoP    =  1.448
1149         GamRhoP   =  0.503
1150         beta      =  -0.161
1151 C...rho(1700)
1152         AmRhoPP   =  1.757
1153         GamRhoPP  =  0.237
1154         gamma     =  0.076
1155 C...exponent
1156         lambda    =  1.0
1157       ENDIF
1158 C
1159 C init
1160 C
1161       EvtCRhoF_W = DCMPLX(0.,0.)
1162
1163       IF (KUHN_SANTA) THEN
1164 C
1165 C Kuehn-Santamaria model
1166 C
1167         EvtCRhoF_W =
1168      >       ( EvtcBW_KS(s,AmRho**2,GamRho)         +     !...BW-rho( 770)
1169      >       EvtcBW_KS(s,AmRhoP**2,GamRhoP)*(beta)  +     !...BW-rho(1450)
1170      >       EvtcBW_KS(s,AmRhoPP**2,GamRhoPP)*(gamma) ) / !...BW-rho(1700)
1171      >       (1.+beta+gamma)
1172       ELSE
1173 C
1174 C Gounaris-Sakurai model
1175 C
1176         EvtCRhoF_W =
1177      >       ( EvtcBW_GS(s,AmRho**2,GamRho)         +
1178      >       EvtcBW_GS(s,AmRhoP**2,GamRhoP)*(beta)  +
1179      >       EvtcBW_GS(s,AmRhoPP**2,GamRhoPP)*(gamma) ) /
1180      >       (1.+beta+gamma)
1181       ENDIF
1182
1183       RETURN
1184       END
1185 C
1186 C * ================================================================== *
1187 C
1188 C
1189 C * ------------------------------------------------------------- *
1190 C
1191       FUNCTION EvtcBW_KS( s,Am2,Gam )
1192
1193       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1194       COMPLEX*16 EvtcBW_KS
1195       PARAMETER (AmPi2=0.13956995**2)
1196
1197       DOUBLE PRECISION lambda
1198       COMMON /LAM/ lambda
1199
1200       EvtcBW_KS = 0.
1201       IF (s.le.4.*AmPi2) RETURN
1202
1203       G  =  Gam* (Am2/s)**lambda * ((s-4.*AmPi2)/(Am2-4.*AmPi2))**1.5
1204       D  =  (Am2-s)**2 + s*G**2
1205       X  =  Am2*(Am2-s)
1206       Y  =  Am2*SQRT(s)*G
1207
1208       EvtcBW_KS = DCMPLX(X/D,Y/D)
1209
1210       RETURN
1211       END
1212 C
1213 C * ------------------------------------------------------------- *
1214 C
1215       FUNCTION EvtcBW_GS( s,Am2,Gam )
1216
1217       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1218       COMPLEX*16 EvtcBW_GS
1219       PARAMETER (AmPi2=0.13956995**2)
1220
1221       DOUBLE PRECISION Evtfs,d,N
1222
1223       DOUBLE PRECISION lambda
1224       COMMON /LAM/ lambda
1225
1226       EvtcBW_GS = 0.
1227       IF (s.le.4.*AmPi2) RETURN
1228
1229       G  =  Gam* (Am2/s)**lambda * ((s-4.*AmPi2)/(Am2-4.*AmPi2))**1.5
1230       z1 =  Am2 - s + Evtfs(s,Am2,Gam)
1231       z2 =  SQRT(s)*G
1232       z3 =  Am2 + d(Am2)*Gam*SQRT(Am2)
1233
1234       X  =  z3 * z1
1235       Y  =  z3 * z2
1236       N  =  z1**2 + z2**2
1237
1238       EvtcBW_GS = DCMPLX(X/N,Y/N)
1239
1240       RETURN
1241       END
1242
1243 C---------------------------------------------------------------------
1244
1245       FUNCTION Evtfs(s,AmRho2,GamRho)
1246       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1247       DOUBLE PRECISION k,k_s,k_Am2
1248       k_s   = k(s)
1249       k_Am2 = k(AmRho2)
1250
1251       Evtfs = GamRho * AmRho2 / k_Am2**3 *
1252      >     (
1253      >          k_s**2 * (h(s)-h(AmRho2)) +
1254      >          (AmRho2-s) * k_Am2**2 * dh_ds(AmRho2)
1255      >     )
1256       RETURN
1257       END
1258 C --------------
1259       FUNCTION k(s)
1260       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1261       DOUBLE PRECISION k
1262       PARAMETER (AmPi2=0.13956995**2)
1263       k = 0.5 * SQRT(s - 4.*AmPi2)
1264       RETURN
1265       END
1266 C --------------
1267       FUNCTION h(s)
1268       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1269       DOUBLE PRECISION k,k_s
1270       PARAMETER (pi=3.141593)
1271       PARAMETER (AmPi=0.13957)
1272       SQRTs = SQRT(s)
1273       k_s = k(s)
1274       h = 2./pi * (k_s/SQRTs) * LOG((SQRTs+2.*k_s)/(2.*AmPi))
1275         RETURN
1276       END
1277 C --------------
1278       FUNCTION dh_ds(s)
1279       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1280       DOUBLE PRECISION k
1281       PARAMETER (pi=3.141593)
1282       dh_ds = h(s) * (1./(8.*k(s)**2) - 1./(2.*s)) + 1./(2.*pi*s)
1283         RETURN
1284       END
1285 C --------------
1286       FUNCTION d(AmRho2)
1287       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1288       DOUBLE PRECISION k,k_AmRho2
1289       PARAMETER (pi=3.141593)
1290       PARAMETER (AmPi=0.13956995)
1291       PARAMETER (AmPi2=AmPi**2)
1292       AmRho    = SQRT(AmRho2)
1293       k_AmRho2 = k(AmRho2)
1294       d = 3./pi * AmPi2/k_AmRho2**2 *
1295      >         LOG((AmRho+2.*k_AmRho2)/(2.*AmPi)) +
1296      >    AmRho/(2.*pi*k_AmRho2) - AmPi2*AmRho/(pi*k_AmRho2**3)
1297          RETURN
1298       END
1299
1300
1301