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: EvtBToKpipi.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 neutral B -->-- K pi pi 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 EVTKpipi is such that
37 C at the level of the time implementation, the maximum
38 C time-dependant weight is unity : nothing is to be
39 C computed to generate unity-weight events. The exact
40 C meaning of the above is made explicit in the HowToUse
42 C [1] first_step_Kpipi
43 C Generation of the kinematics of the 3 prongs
44 C It uses the function evtranf which is a random number
45 C generator providing an uniform distribution
46 C of Real*4 random number between 0 and 1
48 C return the amplitudes of the B0 -->-- K+pi-pi0
49 C corrected for the generation mechanism.
50 c Note that this is a Tagging Mode. The CP conjugate
51 c mode (B0bar -->-- K-pi+pi0) is treated at the same time.
53 C compute the Breit-Wigner of the contributing K* and rho
54 C taking into account the cosine term linked to the
55 C zero-helicity of the spin-1 resonances. There is two
56 c forms of Breit-Wigners available. The first one is the
57 c simple non-relativistic form, while the second is the
58 C relativistic expressions.
59 C The default setting is the relativistic one.
61 C Set the constants values, from the pion mass to the
62 C penguin contributions. It is called by EVTKpipi
64 C And other routines which do not deserve comment here.
65 C===================================================================
68 C Common/Pawc/Hmemor(1000000)
69 C Call Hlimit(1000000)
70 C Call EvtHowToUse_Kpipi
75 subroutine EvtHowToUse_Kpipi
78 Real*8 alphaCP, betaCP
79 Integer iset,number,j,N_gener,N_asked
81 Real*8 p_K_plus(4),p_pi_minus(4)
82 Real*8 p_gamma_1(4),p_gamma_2(4)
84 Real*8 Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar
85 Real*8 Weight,Weight_max
86 Real*8 m_Kstarp,m_Kstar0,m_rhom
97 c run : Simulation of the Anders Ryd Generator
99 Do number=1,N_asked ! weighted events as generated here
102 iset=10000 ! 10^4 events are used to normalize amplitudes
104 iset=0 ! iset must be reset to zero after the first call
107 c Here is the call to EVTKpipi !!!!!!!!!!!!!!!!
108 c communication of data is done by argument only <<<<<<<<
110 + alphaCP,betaCP,iset,
111 + p_K_plus,p_pi_minus,
112 + p_gamma_1,p_gamma_2,
113 + Real_B0,Imag_B0,real_B0bar,Imag_B0bar)
114 C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117 c generate acording to the tag
120 c a B0bar tag => the decay is one from a B0
121 Weight= Real_B0 **2 + Imag_B0 **2
125 c a B0 tag => the decay is one from a B0bar
126 Weight= Real_B0bar **2 + Imag_B0bar **2
130 If(Weight.Gt.evtranf()) Then
131 c----------------------------------------------------------------------------
132 c unweighted event production
133 c----------------------------------------------------------------------------
135 C here is just a Dalitz plot and a few prints
137 m_Kstarp=(p_K_plus (4)+p_gamma_1(4)+p_gamma_2(4))**2
138 m_rhom =(p_pi_minus (4)+p_gamma_1(4)+p_gamma_2(4))**2
139 m_Kstar0=(p_K_plus (4)+p_pi_minus (4)) **2
141 m_Kstarp=m_Kstarp -(p_K_plus (j)+p_gamma_1(j)+p_gamma_2(j))**2
142 m_rhom =m_rhom -(p_pi_minus (j)+p_gamma_1(j)+p_gamma_2(j))**2
143 m_Kstar0=m_Kstar0 -(p_K_plus (j)+p_pi_minus (j)) **2
146 c here is the Dalitz plot when assuming the pion mass for the Kaon
147 c the energy is redefined
150 Wrong = Wrong + p_K_plus(j)**2
152 Wrong = Dsqrt(Wrong+0.139**2)
154 m_Kstarp=(Wrong +p_gamma_1(4)+p_gamma_2(4))**2
155 m_Kstar0=(Wrong +p_pi_minus (4)) **2
157 m_Kstarp=m_Kstarp -(p_K_plus (j)+p_gamma_1(j)+p_gamma_2(j))**2
158 m_Kstar0=m_Kstar0 -(p_K_plus (j)+p_pi_minus (j)) **2
161 c here is a check that weight_max is one
163 If(Weight.gt.Weight_max) Then
165 Print*,' overweighted event found at weight = ',Weight_max
168 c----------------------------------------------------------------------------
171 c end of the loop over events
174 Print*,'number of unity-weight events generated : ',N_gener
175 Print*,'number of trials : ',N_asked
178 C===================================================================
180 + alpha_input,beta_input,iset,
181 + p_K_plus,p_pi_minus,
182 + p_gamma_1,p_gamma_2,
183 + Real_B0,Imag_B0,Real_B0bar,Imag_B0bar)
184 c-----------------------------------------------------------------
185 c ----------------------------------------------------------
186 c --- This is the routine to be called by the Main generator
187 c to get the decay of B0 -->-- K+ pi- pi0
188 c The decay proceeeds through three channels:
189 c a) B0 -->-- K*+ pi- ; K*+ -->-- K+ pi0
190 c b) K*0 pi0 ; K*0bar -->-- K+ pi-
191 c c) K- rho+ ; rho+ -->-- pi+ pi0
192 c .) The K0 rho0 channel is not implemented since it does
193 c not lead to Kpipi final state, but it is interesting
195 c It provides at the same time the CP conjugate decay
196 c B0bar -->-- K- pi+ pi0
197 c****************************************************************************
200 c --- p_K_plus : the four momentum of the K+
201 c --- p_pi_minus : ........................ pi-
202 c --- p_gamma_1 : ........................ first gamma of the pi0
203 c --- p_gamma_2 : ........................ second gamma of the pi0
205 c Note that : the energy is stored in the fourth component
206 c the values are the ones of the B rest frame
207 c a random rotation has been applied
209 c --- Real_B0 : The real part of the amplitude of the decay
210 c B0 -->-- K+ pi- pi0
211 c --- Imag_B0 : ... imaginary ..................................
213 c --- Real_B0bar : The real part of the amplitude of the decay
214 c B0bar -->-- K- pi+ pi0
215 c --- Imag_B0bar : ... imaginary ..................................
217 c****************************************************************************
218 c-----------------------------------------------------------------
220 #include "EvtGenModels/EvtBTo3pi.inc"
221 Real*8 alpha_input, beta_input
223 Real*8 p_K_plus(4),p_pi_minus(4)
224 Real*8 p_gamma_1(4),p_gamma_2(4)
225 Real*8 Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar
229 Real*8 p1(5),p2(5),p3(5)
230 Real*8 Gamma1(5),Gamma2(5)
231 Real*8 factor_max,ABp,ABm
233 data factor_max/1.D+00/
235 c-------------------------------------------------------------------
237 c-------------------------------------------------------------------
238 c this is the normal mode of operation
239 c First, generate the kinematics
246 call Evtfirst_step_Kpipi(p1,p2,p3)
248 c Then, compute the amplitudes
250 Call EvtCompute_Kpipi(p1,p2,p3,
251 + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
252 if(ierr.ne.0 ) Go To 10
253 c----------------nedit EvtBto---------------------------------------------------
254 ElseIf(iset.lt.0) Then
255 c-------------------------------------------------------------------
256 c This is an user mode of operation where the kinematics is
257 c provided by the user who only wants the corresponding amplitudes
263 p3(i)= p_gamma_1 (i) + p_gamma_2 (i)
269 Call EvtCompute_Kpipi(p1,p2,p3,
270 + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
273 Print*,'the provided kinematics are not physical'
275 Print*,'the program will stop'
278 c-------------------------------------------------------------------
279 ElseIf(iset.gt.0) Then
280 c-------------------------------------------------------------------
281 c This is the pre-run mode of operation where initializations are
285 call Evtset_constants(alpha_input, beta_input)
294 call Evtfirst_step_Kpipi(p1,p2,p3)
296 Call EvtCompute_Kpipi(p1,p2,p3,
297 + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
298 if(ierr.ne.0) Go To 20
299 ABp = Real_B0 **2 + Imag_B0 **2
300 ABm = Real_B0bar **2 + Imag_B0Bar **2
302 If(ABp.gt.factor_max) factor_max=ABp
303 If(ABm.gt.factor_max) factor_max=ABm
308 factor_max=1.D+00/Dsqrt(factor_max)
310 c-------------------------------------------------------------------
312 c-------------------------------------------------------------------
314 Real_B0 =Real_B0 * factor_max
315 Imag_B0 =Imag_B0 * factor_max
316 Real_B0bar=Real_B0bar * factor_max
317 Imag_B0Bar=Imag_B0Bar * factor_max
320 c P1,p2,p3 ---> random rotation in B rest frame
322 Call EvtRotation(p1,1)
323 Call EvtRotation(p2,0)
324 Call EvtRotation(p3,0)
326 C Desintegrate the pi_0 s
328 Call EvtGammaGamma(p3,Gamma1,Gamma2)
330 C Feed the output four vectors
336 p_gamma_1 (i)=Gamma1(i)
337 p_gamma_2 (i)=Gamma2(i)
344 c===================================================================
345 subroutine Evtfirst_step_Kpipi(P1,P2,P3)
346 c-----------------------------------------------------------------
347 c ----------------------------------------------------------
348 c --- This routine generates the 5-vectors P1,P2,P3
349 c --- Associated respectively with the Pi+ and two Pi0 s
355 c ----------------------------------------------------------
356 c --- Input Four Vectors
357 C --- Particle [1] is the K+
358 C --- Particle [2] is the pi-
359 C --- Particle [3] is the pi0
360 c ----------------------------------------------------------
362 c ----------------------------------------------------------
364 c ----------------------------------------------------------
365 #include "EvtGenModels/EvtBTo3pi.inc"
366 c ----------------------------------------------------------
368 c ----------------------------------------------------------
372 c ----------------------------------------------------------
373 c --- Variables in Argument
374 c ----------------------------------------------------------
376 real*8 P1(5),P2(5),P3(5)
378 c ----------------------------------------------------------
379 c --- Local Variables
380 c ----------------------------------------------------------
382 real*8 m12,min_m12, max_m12
383 real*8 m13,min_m13, max_m13
384 real*8 m23,min_m23, max_m23
385 Real*8 cost13,cost12,cost23
387 real*8 p1mom,p2mom,p3mom
393 data Phase_space/.false./
395 c ----------------------------------------------------------
397 c ----------------------------------------------------------
399 min_m12 = P1(5) + P2(5) + 2.*Dsqrt(p1(5)*p2(5))
402 min_m13 = P1(5) + P3(5) + 2.*Dsqrt(p1(5)*p3(5))
405 min_m23 = P2(5) + P3(5) + 2.*Dsqrt(p2(5)*p3(5))
409 c ----------------------------------------------------------
410 c --- Generation of the Mass of the Rho or Kstar
411 c ----------------------------------------------------------
414 c ----------------------------------------------------------
415 c --- z is the Flag needed to choose between the generation
416 c --- of a K*+, K*0 or rho- resonance
417 c ----------------------------------------------------------
426 m13 = evtranf()*(max_m13-min_m13)+min_m13
428 y = evtranf()*PI - PI/2.
430 mass = x*Gam_Kstarp/2. +Mass_Kstarp
434 m12 = evtranf()*(max_m12-min_m12)+min_m12
435 m23 = MC2 - m12 - m13
440 m12 = evtranf()*(max_m12-min_m12)+min_m12
442 y = evtranf()*PI - PI/2.
444 mass = x*Gam_Kstar0/2. +Mass_Kstar0
448 m13 = evtranf()*(max_m13-min_m13)+min_m13
449 m23 = MC2 - m12 - m13
454 m23 = evtranf()*(max_m23-min_m23)+min_m23
456 y = evtranf()*PI - PI/2.
458 mass = x*Gam_rho/2. +Mass_rho
462 m13 = evtranf()*(max_m13-min_m13)+min_m13
463 m12 = MC2 - m23 - m13
467 c ----------------------------------------------------------
468 c --- Check that the physics is OK :
469 c --- Are the invariant Masses in allowed ranges ?
470 c ----------------------------------------------------------
472 If(m23.lt.min_m23.or.m23.gt.max_m23) Go to 100
473 If(m13.lt.min_m13.or.m13.gt.max_m13) Go to 100
474 If(m12.lt.min_m12.or.m12.gt.max_m12) Go to 100
476 c ----------------------------------------------------------
477 c --- Are the Cosines of the angles between particles
478 c --- Between -1 and +1 ?
479 c ----------------------------------------------------------
481 P1(4)=(M_B**2+P1(5)-m23)/(2.*M_B)
482 P2(4)=(M_B**2+P2(5)-m13)/(2.*M_B)
483 P3(4)=(M_B**2+P3(5)-m12)/(2.*M_B)
488 If(p1mom.lt.0) Go to 100
489 If(p2mom.lt.0) Go to 100
490 If(p3mom.lt.0) Go to 100
495 cost13=(2.*p1(4)*p3(4)+P1(5)+p3(5)-m13)/(2.*p1mom*p3mom)
496 cost12=(2.*p1(4)*p2(4)+P1(5)+p2(5)-m12)/(2.*p1mom*p2mom)
497 cost23=(2.*p2(4)*p3(4)+P2(5)+p3(5)-m23)/(2.*p2mom*p3mom)
498 If(Dabs(cost13).gt.1.) Go to 100
499 If(Dabs(cost12).gt.1.) Go to 100
500 If(Dabs(cost23).gt.1.) Go to 100
502 c ----------------------------------------------------------
503 c --- Filling the 5-vectors P1,P2,P3
504 c ----------------------------------------------------------
511 P1(1) = p1mom*Dsqrt(1.D+00-cost13**2)
521 c======================================================================
522 Subroutine EvtCompute_Kpipi(p1,p2,p3,
523 + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
524 c-----------------------------------------------------------------------
526 #include "EvtGenModels/EvtBTo3pi.inc"
529 Real*8 m12, m13, m23, W12, W13, W23, Wtot
531 Complex*16 MatB0,MatB0bar,BW12,BW13,BW23
532 Real*8 Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar
533 Real*8 p1(5),p2(5),p3(5)
537 Complex*16 BrightWagner,BreitWigner
540 c ----------------------------------------------------------------
541 C --- Account for the pole compensation
542 c ----------------------------------------------------------------
544 m12 = evt_gmas(p1,p2)
545 m13 = evt_gmas(p1,p3)
546 m23 = evt_gmas(p2,p3)
548 if(m12.lt.0. .or. m13.lt.0. .or. m23.lt.0.) Then
550 Print*,'ierr = ',ierr
558 W12 = 1. / (((Mass_Kstar0 - m12)**2+(Gam_Kstar0/2.)**2)*m12)
559 W13 = 1. / (((Mass_Kstarp - m13)**2+(Gam_Kstarp/2.)**2)*m13)
560 W23 = 1. / (((Mass_rho - m23)**2+(Gam_rho /2.)**2)*m23)
563 Wtot = 1.D+00/Dsqrt(W12 + W13 + W23)
567 c ----------------------------------------------------------------
568 C --- Compute Breit-Wigners
569 c ----------------------------------------------------------------
571 BW13=BrightWagner(p1,p3,p2,Mass_Kstarp,Gam_Kstarp,ierr)
573 BW12=BrightWagner(p1,p2,p3,Mass_Kstar0,Gam_Kstar0,ierr)
575 c If the rho is to be treated on the same footing as K* ==> use the line below
576 c BW23=BrightWagner(p2,p3,p1,Mass_Rho ,Gam_Rho ,ierr)
577 BW23=BreitWigner(p2,p3,p1,ierr)
580 c ----------------------------------------------------------------
582 C --- Build up the amplitudes
583 c ----------------------------------------------------------------
585 c Here come the relative amplitudes of the three decay channels
586 c First, one computes the B0 decay B0 ->- K+ pi- pi0
587 MatB0 = MatKstarp * BW13
590 c Second, one computes the B0bar decay B0bar ->- K- pi+ pi0
591 MatB0bar = NatKstarp * BW13
596 c Pick up the Real and Imaginary parts
598 Real_B0 = dreal(MatB0 )*Wtot
599 Imag_B0 = Imag(MatB0 )*Wtot
601 Real_B0bar = dreal(MatB0bar)*Wtot
602 Imag_B0Bar = Imag(MatB0bar)*Wtot
606 c======================================================================
607 Function BrightWagner(p1,p2,p3,Mass,Width,ierr)
608 c----------------------------------------------------------------------
610 #include "EvtGenModels/EvtBTo3pi.inc"
611 Complex *16 BrightWagner,EvtRBW
615 Real*8 Mass,Width,Am2Min
617 c ---------------------------------------------------------------
618 c --- Input Four Vectors
619 c ---------------------------------------------------------------
620 Real*8 p1(5),p2(5),p3(5)
622 c ---------------------------------------------------------------
623 C --- intermediate variables
624 c ---------------------------------------------------------------
625 Real*8 E12_2,m12_2,beta,gamma,argu,m13_2,costet,coscms,m12
626 Real*8 Factor,Num_real,Num_imag
628 Real *8 p1z,p1zcms12,e1cms12,p1cms12
630 c ---------------------------------------------------------------
632 c ---------------------------------------------------------------
634 BrightWagner=Dcmplx(0,0)
637 E12_2=(p1(4)+p2(4))**2
640 m12_2=m12_2-(p1(i)+p2(i))**2
642 Argu = 1.D+00 - m12_2 / E12_2
647 Print *,'Abnormal beta ! Argu = ',Argu
656 Print *,'Abnormal m12 ! m12_2 = ',m12_2
663 gamma=Dsqrt(E12_2/m12_2)
664 c ---------------------------------------------------------------
665 C --- get the cosine in the B CMS
666 c ---------------------------------------------------------------
668 m13_2=(p1(4)+p3(4))**2
670 m13_2=m13_2-(p1(i)+p3(i))**2
672 if(m13_2.lt.0) Go To 50
673 if((p1(4)**2-p1(5)).lt.0) Go To 50
674 if((p3(4)**2-p3(5)).lt.0) Go To 50
675 costet= (2.D+00*p1(4)*p3(4)-m13_2+p1(5)+p3(5))
677 > (2.D+00*Dsqrt( (p1(4)**2-p1(5)) * (p3(4)**2-p3(5)) ))
679 c ---------------------------------------------------------------
680 C --- get the costet in the 1-2 CMS
681 c ---------------------------------------------------------------
683 p1z=dsqrt(P1(4)**2-p1(5))*costet
684 p1zcms12=gamma*(p1z+beta*P1(4))
685 e1cms12 =gamma*(p1(4)+beta*p1z)
686 p1cms12 =Dsqrt(e1cms12**2-p1(5))
687 coscms=p1zcms12/p1cms12
689 c ---------------------------------------------------------------
690 C --- Build the Breit Wigner
691 c ---------------------------------------------------------------
694 Am2Min = p1(5) + p2(5) + 2.*Dsqrt( p1(5)*p2(5) )
695 BrightWagner = coscms * EvtRBW(m12_2,Mass**2,Width,Am2Min)
697 Factor = 2.D+00* ( (Mass-m12)**2+(0.5D+00*Width)**2 )
698 Factor = coscms * Width / Factor
699 Num_real= (Mass-m12) * Factor
700 Num_imag= 0.5D+00*Width * Factor
701 BrightWagner=Dcmplx(Num_real,Num_imag)
711 FUNCTION EvtRBW(s,Am2,Gam,Am2Min)
713 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
717 IF (s.le.Am2Min) RETURN
719 G = Gam* (Am2/s) * ((s-Am2Min)/(Am2-Am2Min))**1.5
720 D = (Am2-s)**2 + s*G**2
724 EvtRBW = DCMPLX(X/D,Y/D)