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