TEvtGen/EvtGenModels/EvtBTo3piP00.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: EvtBTo3piP00.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 [*] 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] 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_P00
  43 C                 Generation of the kinematics of the 3 pions
  44 C                 It uses the function ranf 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_P00
  77 c      Stop
  78 c      End
  79
  80       subroutine EvtHowToUse_P00
  81
  82       Implicit none
  83       Real*8  alpha
  84       Integer iset,number,j,N_gener,N_asked
  85
  86       Real*8  p_pi_plus(4)
  87       Real*8  p_gamma_1(4),p_gamma_2(4),p_gamma_3(4),p_gamma_4(4)
  88
  89       Real*8  Real_Bp,Imag_Bp,Real_Bm,Imag_Bm
  90       Real*8  Weight,Weight_max
  91       Real*8  ABp,ABm
  92       Real*8  m_rho12,m_rho13
  93       Real*8  Evtranf,Tag
  94
  95       alpha      = 0.4
  96       N_gener    = 0
  97       N_asked    = 100000
  98
  99       weight_max = 1.0
 100
 101 c run : Simulation of the Anders Ryd Generator
 102
 103       Do number=1,N_asked ! weighted events as generated here
 104
 105       If(number.eq.1) then
 106       iset=10000          ! 10^4 events are used to normalize amplitudes
 107       Else
 108       iset=0              ! iset must be reset to zero after the first call
 109       End If
 110
 111 c Here is the call to EVT3pions_P00  !!!!!!!!!!!!!!!!
 112 c communication of data is done by argument only <<<<<<<<
 113        call EVT3pionsP00(
 114      + alpha,iset,
 115      + p_pi_plus,
 116      + p_gamma_1,p_gamma_2,p_gamma_3,p_gamma_4,
 117      + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm)
 118
 119 C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 120 c select the Tag
 121       Tag   =evtranf()
 122 C get the relevant quantities
 123       ABp   = Real_Bp   **2 + Imag_Bp   **2
 124       ABm   = Real_Bp   **2 + Imag_Bm   **2
 125 c generate acording to the tag
 126       If(Tag.gt.0.5) Then
 127
 128 c a Bm tag =>  the decay is one from a Bp
 129       Weight= ABp
 130
 131       Else
 132
 133 c a Bp tag =>  the decay is one from a Bm
 134       Weight= ABm
 135
 136       End If
 137       If(Weight.Gt.evtranf()) Then
 138 c----------------------------------------------------------------------------
 139 c unweighted event production
 140 c----------------------------------------------------------------------------
 141       N_gener=N_gener+1
 142 C here is just a Dalitz plot and a few prints
 143
 144       m_rho12=(p_pi_plus (4)+p_gamma_1(4)+p_gamma_2(4))**2
 145       m_rho13=(p_pi_plus (4)+p_gamma_3(4)+p_gamma_4(4))**2
 146       do j=1,3
 147       m_rho12=m_rho12-(p_pi_plus (j)+p_gamma_1(j)+p_gamma_2(j))**2
 148       m_rho13=m_rho13-(p_pi_plus (j)+p_gamma_3(j)+p_gamma_4(j))**2
 149       end do
 150
 151 c here is a check that weight_max is one
 152
 153       If(Weight.gt.Weight_max) Then
 154       Weight_max=Weight
 155       Print*,' overweighted event found at weight = ',Weight_max
 156       End If
 157
 158 c----------------------------------------------------------------------------
 159       End If
 160
 161 c end of the loop over events
 162       End Do
 163
 164       Print*,'number of unity-weight events generated : ',N_gener
 165       Print*,'number of trials                        : ',N_asked
 166
 167       End
 168 C===================================================================
 169       subroutine Evt3pionsP00(
 170      + alpha_input,iset,
 171      + p_pi_plus,
 172      + p_pi_1_gamma_1,p_pi_1_gamma_2,
 173      + p_pi_2_gamma_1,p_pi_2_gamma_2,
 174      + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm)
 175 c-----------------------------------------------------------------
 176 c       ----------------------------------------------------------
 177 c       --- This is the routine to be called by the Main generator
 178 c           to get the decay of B+ -->-- 3 pions
 179 c       --- AND
 180 c           to get the decay of B- -->-- 3 pions
 181 C       For the sake of clarity, signs refers to B+ decay
 182 c****************************************************************************                                 set to 0.
 183 c       --- Outputs are :
 184 c
 185 c       ---               p_pi_plus      : the four momentum of the pi+
 186 c                                          Then, for the first  pi0
 187 c       ---               p_pi_1_gamma_1 : the four momentum of the first  photon
 188 c       ---               p_pi_1_gamma_2 : the four momentum of the second photon
 189 c                                          Then, for the second pi0
 190 c       ---               p_pi_2_gamma_1 : the four momentum of the first  photon
 191 c       ---               p_pi_2_gamma_2 : the four momentum of the second photon
 192 c
 193 c             Note that : the energy is stored in the fourth component
 194 c                         the values are the ones of the B rest frame
 195 c                         a random rotation has been applied
 196 c
 197 c       ---               Real_B0p   : The real      part of the amplitude of
 198 c                                    the B+    ->- 3 pions decay
 199 c       ---               Imag_B0p   : The imaginary part of the amplitude of
 200 c                                    the B+    ->- 3 pions decay
 201 c
 202 c       ---               Real_B0m   : The real      part of the amplitude of
 203 c                                    the B-    ->- 3 pions decay
 204 c       ---               Imag_B0m   : The imaginary part of the amplitude of
 205 c
 206 c****************************************************************************
 207 c-----------------------------------------------------------------
 208       Implicit none
 209 #include "EvtGenModels/EvtBTo3pi.inc"
 210       Real*8  alpha_input
 211       Integer iset
 212       Real*8  p_pi_plus(4)
 213       Real*8  p_pi_1_gamma_1(4),p_pi_1_gamma_2(4)
 214       Real*8  p_pi_2_gamma_1(4),p_pi_2_gamma_2(4)
 215       Real*8  Real_Bp,Imag_Bp,Real_Bm,Imag_Bm
 216
 217 c Working quantities
 218       Integer i,number
 219       Real*8  p1(5),p2(5),p3(5)
 220       Real*8  Gamma1(5),Gamma2(5),Gamma3(5),Gamma4(5)
 221       Real*8  factor_max,ABp,ABm
 222       Integer ierr
 223         data factor_max/1.D+00/
 224         ierr =0
 225 c-------------------------------------------------------------------
 226       If(iset.eq.0) Then
 227 c-------------------------------------------------------------------
 228 c     this is the normal mode of operation
 229 c     First, generate the kinematics
 230
 231       p1(5)=M_pip**2
 232       p2(5)=M_pi0**2
 233       p3(5)=M_pi0**2
 234
 235
 236  10   continue
 237       call Evtfirst_step_P00(p1,p2,p3)
 238
 239 c     Then, compute the amplitudes
 240
 241       Call EvtCompute_P00(p1,p2,p3,
 242      +    Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
 243         if(ierr.ne.0 ) Go To 10
 244 c-------------------------------------------------------------------
 245       ElseIf(iset.lt.0) Then
 246 c-------------------------------------------------------------------
 247 c     This is an user mode of operation where the kinematics is
 248 c     provided by the user who only wants the corresponding amplitudes
 249 c     to be computed
 250
 251       Do i=1,4
 252       p1(i)= p_pi_plus (i)
 253       p2(i)= p_pi_1_gamma_1 (i) + p_pi_1_gamma_2 (i)
 254       p3(i)= p_pi_2_gamma_1 (i) + p_pi_2_gamma_2 (i)
 255       End Do
 256       p1(5)= M_pip**2
 257       p2(5)= M_pi0**2
 258       p3(5)= M_pi0**2
 259
 260       Call EvtCompute_P00(p1,p2,p3,
 261      +    Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
 262
 263        if(ierr.ne.0) Then
 264         Print*,'the provided kinematics are not physical'
 265         Print*,'ierr=',ierr
 266         Print*,'the program will stop'
 267         Stop
 268        endif
 269 c-------------------------------------------------------------------
 270       ElseIf(iset.gt.0) Then
 271 c-------------------------------------------------------------------
 272 c     This is the pre-run mode of operation where initializations are
 273 c     performed.
 274
 275       factor_max= 0
 276 c      changed by ryd April 24 1998.
 277 c      0.35 is the value of beta? should it not be
 278 c      passed in as an argument????
 279 c      call Evtset_constants(alpha_input,0.362)
 280       call Evtset_constants(alpha_input,0.362D0)
 281       p1(5)=M_pip**2
 282       p2(5)=M_pi0**2
 283       p3(5)=M_pi0**2
 284
 285 c     pre-run
 286       Do number=1,iset
 287
 288   20  continue
 289       call Evtfirst_step_P00(p1,p2,p3)
 290
 291       Call EvtCompute_P00(p1,p2,p3,
 292      +    Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
 293         if(ierr.ne.0) Go To 20
 294       ABp   = Real_Bp   **2 + Imag_Bp   **2
 295       ABm   = Real_Bm   **2 + Imag_Bm   **2
 296
 297       If(ABp.gt.factor_max) factor_max=ABp
 298       If(ABm.gt.factor_max) factor_max=ABm
 299
 300       End Do
 301 c     end of the pre-run
 302
 303       factor_max=1.D+00/Dsqrt(factor_max)
 304
 305 c-------------------------------------------------------------------
 306       End If
 307 c-------------------------------------------------------------------
 308
 309       Real_Bp   =Real_Bp    * factor_max
 310       Imag_Bp   =Imag_Bp    * factor_max
 311       Real_Bm   =Real_Bm    * factor_max
 312       Imag_Bm   =Imag_Bm    * factor_max
 313
 314       if(iset.lt.0) return
 315 c     P1,p2,p3 ---> random rotation in B rest frame
 316
 317       Call EvtRotation(p1,1)
 318       Call EvtRotation(p2,0)
 319       Call EvtRotation(p3,0)
 320
 321 C     Desintegrate the pi_0 s
 322
 323       Call EvtGammaGamma(p2,Gamma1,Gamma2)
 324       Call EvtGammaGamma(p3,Gamma3,Gamma4)
 325
 326 C     Feed the output four vectors
 327
 328       Do i=1,4
 329
 330       p_pi_plus      (i)=p1    (i)
 331       p_pi_1_gamma_1 (i)=Gamma1(i)
 332       p_pi_1_gamma_2 (i)=Gamma2(i)
 333       p_pi_2_gamma_1 (i)=Gamma3(i)
 334       p_pi_2_gamma_2 (i)=Gamma4(i)
 335
 336       End Do
 337
 338       Return
 339
 340       End
 341
 342 c===================================================================
 343       subroutine Evtfirst_step_P00(P1,P2,P3)
 344 c-----------------------------------------------------------------
 345 c       ----------------------------------------------------------
 346 c       --- This routine generates the 5-vectors P1,P2,P3
 347 c       --- Associated respectively with the Pi+ and two Pi0 s
 348 c       ---             P1(1) = Px
 349 c       ---             P1(2) = Py
 350 c       ---             P1(3) = Pz
 351 c       ---             P1(4) = E
 352 c       ---             P1(5) = M**2
 353 c       ----------------------------------------------------------
 354 c       ---     Input Four Vectors
 355 C       ---     Particle [1] is the pi+
 356 C       ---     Particle [2] is the pi0
 357 C       ---     Particle [3] is the pi0
 358 c       ----------------------------------------------------------
 359
 360 c       ----------------------------------------------------------
 361 c       --- commons
 362 c       ----------------------------------------------------------
 363 #include "EvtGenModels/EvtBTo3pi.inc"
 364 c       ----------------------------------------------------------
 365 c       --- Used Functions
 366 c       ----------------------------------------------------------
 367
 368         real*8 evtranf
 369
 370 c       ----------------------------------------------------------
 371 c       --- Variables in Argument
 372 c       ----------------------------------------------------------
 373
 374          real*8 P1(5),P2(5),P3(5)
 375
 376 c       ----------------------------------------------------------
 377 c       --- Local Variables
 378 c       ----------------------------------------------------------
 379
 380         real*8 m12,min_m12, max_m12
 381         real*8 m13,min_m13, max_m13
 382         real*8 m23,min_m23, max_m23
 383         Real*8 cost13,cost12,cost23
 384
 385         real*8 p1mom,p2mom,p3mom
 386
 387         real*8 x, y, z, mass
 388         integer i
 389
 390         Logical Phase_space
 391         data Phase_space/.false./
 392
 393 c       initialize to avoid warning on linux
 394         mass=0.0
 395
 396 c       ----------------------------------------------------------
 397 c       --- Computation
 398 c       ----------------------------------------------------------
 399         max_m12 = M_B**2
 400         min_m12 = P1(5) + P2(5)
 401
 402         max_m13 = M_B**2
 403         min_m13 = P1(5) + P3(5)
 404
 405         max_m23 = M_B**2
 406         min_m23 = P2(5) + P3(5)
 407
 408 100     Continue
 409
 410 c       ----------------------------------------------------------
 411 c       --- Generation of the Mass of the Rho(+)
 412 c       ----------------------------------------------------------
 413
 414         If(.not.Phase_space) Then
 415         y = evtranf()*PI - PI/2.
 416         x = Dtan(y)
 417         mass = x*Gam_rho/2. +Mass_rho
 418         End If
 419
 420 c       ----------------------------------------------------------
 421 c       --- z is the Flag needed to choose between the generation
 422 c       --- of a Rho+ = pi+ pi_0[1] or pi+ pi_0[2]
 423 c       ----------------------------------------------------------
 424
 425         z = evtranf()
 426
 427         if(z.lt..5) Then
 428
 429                 If(Phase_space) Then
 430                 m12 = evtranf()*(max_m12-min_m12)+min_m12
 431                 Else
 432                 m12 = mass**2
 433                 End If
 434
 435         m13 = evtranf()*(max_m13-min_m13)+min_m13
 436         m23 = MB2 - m12 - m13
 437
 438         else
 439
 440                 If(Phase_space) Then
 441                 m13 = evtranf()*(max_m13-min_m13)+min_m13
 442                 Else
 443                 m13 =  mass**2
 444                 End If
 445
 446         m12 = evtranf()*(max_m12-min_m12)+min_m12
 447         m23 = MB2 - m12 - m13
 448
 449         endif
 450
 451 c       ----------------------------------------------------------
 452 c       --- Check that the physics is OK :
 453 c       --- Are the invariant Masses in allowed ranges ?
 454 c       ----------------------------------------------------------
 455
 456         If(m23.lt.min_m23.or.m23.gt.max_m23) Go to 100
 457         If(m13.lt.min_m13.or.m13.gt.max_m13) Go to 100
 458         If(m12.lt.min_m12.or.m12.gt.max_m12) Go to 100
 459
 460 c       ----------------------------------------------------------
 461 c       --- Are the Cosines of the angles between particles
 462 c       --- Between -1 and +1 ?
 463 c       ----------------------------------------------------------
 464
 465         P1(4)=(M_B**2+P1(5)-m23)/(2.*M_B)
 466         P2(4)=(M_B**2+P2(5)-m13)/(2.*M_B)
 467         P3(4)=(M_B**2+P3(5)-m12)/(2.*M_B)
 468
 469         p1mom=p1(4)**2-P1(5)
 470         p2mom=p2(4)**2-P2(5)
 471         p3mom=p3(4)**2-P3(5)
 472         If(p1mom.lt.0) Go to 100
 473         If(p2mom.lt.0) Go to 100
 474         If(p3mom.lt.0) Go to 100
 475         p1mom=Dsqrt(p1mom)
 476         p2mom=Dsqrt(p2mom)
 477         p3mom=Dsqrt(p3mom)
 478
 479         cost13=(2.*p1(4)*p3(4)+P1(5)+p3(5)-m13)/(2.*p1mom*p3mom)
 480         cost12=(2.*p1(4)*p2(4)+P1(5)+p2(5)-m12)/(2.*p1mom*p2mom)
 481         cost23=(2.*p2(4)*p3(4)+P2(5)+p3(5)-m23)/(2.*p2mom*p3mom)
 482         If(Dabs(cost13).gt.1.) Go to 100
 483         If(Dabs(cost12).gt.1.) Go to 100
 484         If(Dabs(cost23).gt.1.) Go to 100
 485
 486 c       ----------------------------------------------------------
 487 c       --- Filling the 5-vectors P1,P2,P3
 488 c       ----------------------------------------------------------
 489
 490         P3(1) = 0
 491         P3(2) = 0
 492         p3(3) = p3mom
 493
 494         P1(3) = p1mom*cost13
 495         P1(1) = p1mom*Dsqrt(1.D+00-cost13**2)
 496         p1(2) = 0.
 497
 498         Do i=1,3
 499         P2(i)=-p1(i)-p3(i)
 500         End do
 501
 502         END
 503
 504
 505 c======================================================================
 506       Subroutine EvtCompute_P00(p1,p2,p3,
 507      +           Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
 508 c-----------------------------------------------------------------------
 509         IMPLICIT None
 510 #include "EvtGenModels/EvtBTo3pi.inc"
 511
 512
 513         Real*8 m12, m13, W12, W13, Wtot
 514         Real*8 evt_gmas
 515         Complex*16 MatBp,MatBm
 516         Real*8  Real_Bp,Imag_Bp,Real_Bm,Imag_Bm
 517         Real*8     p1(5),p2(5),p3(5)
 518         real*8  ASHQ
 519         Data    ASHQ/0.707107 D+00/
 520
 521         Integer ierr,iset
 522         COMPLEX*16 Mat_rhop
 523         Complex*16 BreitWigner
 524
 525         ierr = 0
 526 c       ----------------------------------------------------------------
 527 C       ---     Account for the pole compensation
 528 c       ----------------------------------------------------------------
 529
 530         if(evt_gmas(p1,p2).lt.0.or.evt_gmas(p1,p3)
 531      +          .lt.0.or.evt_gmas(p2,p3).lt.0) Then
 532                 ierr=1
 533                 Print*,'ierr = ',ierr
 534                 return
 535         endif
 536         m12 = sqrt(evt_gmas(p1,p2))
 537         m13 = sqrt(evt_gmas(p1,p3))
 538
 539         W12 = (1./m12)*1./((Mass_rho - m12)**2+(Gam_rho/2.)**2)
 540         W13 = (1./m13)*1./((Mass_rho - m13)**2+(Gam_rho/2.)**2)
 541         if(iset.ge.0) Then
 542         Wtot = 1.D+00/Dsqrt(W12 + W13)
 543             else
 544         Wtot =1.
 545         Endif
 546 c       ----------------------------------------------------------------
 547 C       ---     Compute Breit-Wigners
 548 c       ----------------------------------------------------------------
 549
 550         Mat_rhop = BreitWigner(p1,p2,p3,ierr)
 551      +  + BreitWigner(p1,p3,p2,ierr)
 552
 553 c       ----------------------------------------------------------------
 554 C       ---     Build up the amplitudes
 555 c       ----------------------------------------------------------------
 556 c       The factor ASHQ = 1./sqrt(2) is here just to stick to the notations
 557 c       used by Art Snyder and Helen Quinn. It is irrelevant for the genera-
 558 c       tion of events done here.
 559
 560         MatBp    = Mat_s1 * Mat_rhop * Wtot * ASHQ
 561
 562         MatBm    = Nat_s1 * Mat_rhop * Wtot * ASHQ
 563
 564 c       Pick up the Real and Imaginary parts
 565 c       changed Real to DBLE and Imag DIMAG to make it complie
 566 c       with Absoft and g77 (ryd)
 567
 568
 569         Real_Bp    = DBLE(MatBp)
 570         Imag_Bp    = DIMAG(MatBp)
 571
 572         Real_Bm    = DBLE(MatBm)
 573         Imag_Bm    = DIMAG(MatBm)
 574
 575         Return
 576         End
 577