1 C--------------------------------------------------------------------------
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.
8 C Copyright Information: See EvtGen/COPYRIGHT
9 C Copyright (C) 1998 Caltech, UCSB
11 C Module: EvtBTo3piMPP.F
15 C Modification history:
17 C DJL/RYD August 11, 1998 Module created
19 C------------------------------------------------------------------------
20 C===================================================================
21 C This package is providing a B -->-- 3pions decay generator
22 C Its is composed of the following subroutines:
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.
30 C===================================================================
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
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
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]
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-
70 C Set the constants values, from the pion mass to the
71 C penguin contributions. It is called by EVT3pions
73 C And other routines which do not deserve comment here.
74 C===================================================================
77 c Common/Pawc/Hmemor(1000000)
78 c Call Hlimit(1000000)
80 c call EvtHowToUse_MPP
81 c call EvtHowToUse_P00
85 subroutine EvtHowToUse_MPP
89 Integer iset,number,j,N_gener,N_asked
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
96 Real*8 m_rho12,m_rho13
104 c run : Simulation of the Anders Ryd Generator
106 Do number=1,N_asked ! weighted events as generated here
109 iset=10000 ! 10^4 events are used to normalize amplitudes
111 iset=0 ! iset must be reset to zero after the first call
114 c Here is the call to EVT3pions_MPP !!!!!!!!!!!!!!!!
115 c communication of data is done by argument only <<<<<<<<
118 + p_pi_minus,p_pi_plus_1,p_pi_plus_2,
119 + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm)
121 C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
130 c a Bm tag => the decay is one from a Bp
135 c a Bp tag => the decay is one from a Bm
139 If(Weight.Gt.evtranf()) Then
140 c----------------------------------------------------------------------------
141 c unweighted event production
142 c----------------------------------------------------------------------------
144 C here is just a Dalitz plot and a few prints
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
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
153 c here is a check that weight_max is one
155 If(Weight.gt.Weight_max) Then
157 Print*,' overweighted event found at weight = ',Weight_max
160 c----------------------------------------------------------------------------
163 c end of the loop over events
166 Print*,'number of unity-weight events generated : ',N_gener
167 Print*,'number of trials : ',N_asked
171 C===================================================================
172 subroutine Evt3pionsMPP(
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
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.
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+
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
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
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
203 c****************************************************************************
204 c-----------------------------------------------------------------
206 #include "EvtGenModels/EvtBTo3pi.inc"
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
214 Real*8 p1(5),p2(5),p3(5)
215 Real*8 factor_max,ABp,ABm
217 data factor_max/1.D+00/
219 c-------------------------------------------------------------------
221 c-------------------------------------------------------------------
222 c this is the normal mode of operation
223 c First, generate the kinematics
230 call Evtfirst_step_MPP(p1,p2,p3)
232 c Then, compute the amplitudes
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
245 p1(i)= p_pi_minus (i)
246 p2(i)= p_pi_plus_1(i)
247 p3(i)= p_pi_plus_2(i)
253 Call EvtCompute_MPP(p1,p2,p3,
254 + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
256 Print*,'the provided kinematics are not physical'
258 Print*,'the program will stop'
261 c-------------------------------------------------------------------
262 ElseIf(iset.gt.0) Then
263 c-------------------------------------------------------------------
264 c This is the pre-run mode of operation where initializations are
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)
281 call Evtfirst_step_MPP(p1,p2,p3)
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
289 If(ABp.gt.factor_max) factor_max=ABp
290 If(ABm.gt.factor_max) factor_max=ABm
295 factor_max=1.D+00/Dsqrt(factor_max)
297 c-------------------------------------------------------------------
299 c-------------------------------------------------------------------
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
307 c P1,p2,p3 ---> random rotation in B rest frame
309 Call EvtRotation(p1,1)
310 Call EvtRotation(p2,0)
311 Call EvtRotation(p3,0)
313 C Feed the output four vectors
317 p_pi_minus (i)=p1 (i)
318 p_pi_plus_1 (i)=p2 (i)
319 p_pi_plus_2 (i)=p3 (i)
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
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 ----------------------------------------------------------
345 c ----------------------------------------------------------
347 c ----------------------------------------------------------
348 #include "EvtGenModels/EvtBTo3pi.inc"
349 c ----------------------------------------------------------
351 c ----------------------------------------------------------
356 c ----------------------------------------------------------
357 c --- Variables in Argument
358 c ----------------------------------------------------------
360 real*8 P1(5),P2(5),P3(5)
362 c ----------------------------------------------------------
363 c --- Local Variables
364 c ----------------------------------------------------------
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
371 real*8 p1mom,p2mom,p3mom
377 data Phase_space/.false./
379 c initialize to avoid warning on linux
382 c ----------------------------------------------------------
384 c ----------------------------------------------------------
386 min_m12 = P1(5) + P2(5)
389 min_m13 = P1(5) + P3(5)
392 min_m23 = P2(5) + P3(5)
396 c ----------------------------------------------------------
397 c --- Generation of the Mass of the Rho(0)
398 c ----------------------------------------------------------
400 If(.not.Phase_space) Then
401 y = evtranf()*PI - PI/2.
403 mass = x*Gam_rho/2. +Mass_rho
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 ----------------------------------------------------------
416 m12 = evtranf()*(max_m12-min_m12)+min_m12
421 m13 = evtranf()*(max_m13-min_m13)+min_m13
422 m23 = MB2 - m12 - m13
427 m13 = evtranf()*(max_m13-min_m13)+min_m13
432 m12 = evtranf()*(max_m12-min_m12)+min_m12
433 m23 = MB2 - m12 - m13
437 c ----------------------------------------------------------
438 c --- Check that the physics is OK :
439 c --- Are the invariant Masses in allowed ranges ?
440 c ----------------------------------------------------------
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
446 c ----------------------------------------------------------
447 c --- Are the Cosines of the angles between particles
448 c --- Between -1 and +1 ?
449 c ----------------------------------------------------------
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)
458 If(p1mom.lt.0) Go to 100
459 If(p2mom.lt.0) Go to 100
460 If(p3mom.lt.0) Go to 100
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
472 c ----------------------------------------------------------
473 c --- Filling the 5-vectors P1,P2,P3
474 c ----------------------------------------------------------
481 P1(1) = p1mom*Dsqrt(1.D+00-cost13**2)
491 c======================================================================
492 Subroutine EvtCompute_MPP(p1,p2,p3,
493 + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
494 c-----------------------------------------------------------------------
496 #include "EvtGenModels/EvtBTo3pi.inc"
499 Real*8 m12, m13, W12, W13, Wtot
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)
506 Data ASHQ/0.707107 D+00/
509 Complex*16 BreitWigner
512 c ----------------------------------------------------------------
513 C --- Account for the pole compensation
514 c ----------------------------------------------------------------
515 if(evt_gmas(p1,p2).lt.0.or.evt_gmas(p1,p3)
521 m12 = sqrt(evt_gmas(p1,p2))
522 m13 = sqrt(evt_gmas(p1,p3))
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)
528 Wtot = 1.D+00/Dsqrt(W12 + W13)
532 c ----------------------------------------------------------------
533 C --- Compute Breit-Wigners
534 c ----------------------------------------------------------------
536 Mat_rhop = BreitWigner(p1,p2,p3,ierr)
537 + + BreitWigner(p1,p3,p2,ierr)
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.
546 MatBp = Mat_s2 * Mat_rhop * Wtot * ASHQ
548 MatBm = Nat_s2 * Mat_rhop * Wtot * ASHQ
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)
554 Real_Bp = DBLE(MatBp)
555 Imag_Bp = DIMAG(MatBp)
557 Real_Bm = DBLE(MatBm)
558 Imag_Bm = DIMAG(MatBm)