]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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: EvtBToKpipi.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 neutral B -->-- K pi pi 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] EVTKpipi | |
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 | |
41 | C routine. | |
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 | |
47 | C [2] compute | |
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. | |
52 | C [3] BreitWigner | |
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. | |
60 | C [4] Set_constants | |
61 | C Set the constants values, from the pion mass to the | |
62 | C penguin contributions. It is called by EVTKpipi | |
63 | C | |
64 | C And other routines which do not deserve comment here. | |
65 | C=================================================================== | |
66 | c Implicit none | |
67 | C Real Hmemor | |
68 | C Common/Pawc/Hmemor(1000000) | |
69 | C Call Hlimit(1000000) | |
70 | C Call EvtHowToUse_Kpipi | |
71 | C Stop | |
72 | C End | |
73 | ||
74 | ||
75 | subroutine EvtHowToUse_Kpipi | |
76 | ||
77 | Implicit none | |
78 | Real*8 alphaCP, betaCP | |
79 | Integer iset,number,j,N_gener,N_asked | |
80 | ||
81 | Real*8 p_K_plus(4),p_pi_minus(4) | |
82 | Real*8 p_gamma_1(4),p_gamma_2(4) | |
83 | ||
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 | |
87 | Real*8 Wrong | |
88 | Real*4 Evtranf,Tag | |
89 | ||
90 | alphaCP = 1.35 | |
91 | betaCP = 0.362 | |
92 | N_gener = 0 | |
93 | N_asked = 100000 | |
94 | ||
95 | weight_max = 1.0 | |
96 | ||
97 | c run : Simulation of the Anders Ryd Generator | |
98 | ||
99 | Do number=1,N_asked ! weighted events as generated here | |
100 | ||
101 | If(number.eq.1) then | |
102 | iset=10000 ! 10^4 events are used to normalize amplitudes | |
103 | Else | |
104 | iset=0 ! iset must be reset to zero after the first call | |
105 | End If | |
106 | ||
107 | c Here is the call to EVTKpipi !!!!!!!!!!!!!!!! | |
108 | c communication of data is done by argument only <<<<<<<< | |
109 | call EVTKpipi( | |
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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | |
115 | c select the Tag | |
116 | Tag =evtranf() | |
117 | c generate acording to the tag | |
118 | If(Tag.gt.0.5) Then | |
119 | ||
120 | c a B0bar tag => the decay is one from a B0 | |
121 | Weight= Real_B0 **2 + Imag_B0 **2 | |
122 | ||
123 | Else | |
124 | ||
125 | c a B0 tag => the decay is one from a B0bar | |
126 | Weight= Real_B0bar **2 + Imag_B0bar **2 | |
127 | ||
128 | End If | |
129 | ||
130 | If(Weight.Gt.evtranf()) Then | |
131 | c---------------------------------------------------------------------------- | |
132 | c unweighted event production | |
133 | c---------------------------------------------------------------------------- | |
134 | N_gener=N_gener+1 | |
135 | C here is just a Dalitz plot and a few prints | |
136 | ||
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 | |
140 | do j=1,3 | |
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 | |
144 | end do | |
145 | ||
146 | c here is the Dalitz plot when assuming the pion mass for the Kaon | |
147 | c the energy is redefined | |
148 | Wrong =0. | |
149 | do j=1,3 | |
150 | Wrong = Wrong + p_K_plus(j)**2 | |
151 | end do | |
152 | Wrong = Dsqrt(Wrong+0.139**2) | |
153 | ||
154 | m_Kstarp=(Wrong +p_gamma_1(4)+p_gamma_2(4))**2 | |
155 | m_Kstar0=(Wrong +p_pi_minus (4)) **2 | |
156 | do j=1,3 | |
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 | |
159 | end do | |
160 | ||
161 | c here is a check that weight_max is one | |
162 | ||
163 | If(Weight.gt.Weight_max) Then | |
164 | Weight_max=Weight | |
165 | Print*,' overweighted event found at weight = ',Weight_max | |
166 | End If | |
167 | ||
168 | c---------------------------------------------------------------------------- | |
169 | End If | |
170 | ||
171 | c end of the loop over events | |
172 | End Do | |
173 | ||
174 | Print*,'number of unity-weight events generated : ',N_gener | |
175 | Print*,'number of trials : ',N_asked | |
176 | ||
177 | End | |
178 | C=================================================================== | |
179 | subroutine EvtKpipi( | |
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 | |
194 | c in itself. | |
195 | c It provides at the same time the CP conjugate decay | |
196 | c B0bar -->-- K- pi+ pi0 | |
197 | c**************************************************************************** | |
198 | c --- Outputs are : | |
199 | 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 | |
204 | c | |
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 | |
208 | c | |
209 | c --- Real_B0 : The real part of the amplitude of the decay | |
210 | c B0 -->-- K+ pi- pi0 | |
211 | c --- Imag_B0 : ... imaginary .................................. | |
212 | c similarly | |
213 | c --- Real_B0bar : The real part of the amplitude of the decay | |
214 | c B0bar -->-- K- pi+ pi0 | |
215 | c --- Imag_B0bar : ... imaginary .................................. | |
216 | ||
217 | c**************************************************************************** | |
218 | c----------------------------------------------------------------- | |
219 | Implicit none | |
220 | #include "EvtGenModels/EvtBTo3pi.inc" | |
221 | Real*8 alpha_input, beta_input | |
222 | Integer iset | |
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 | |
226 | ||
227 | c Working quantities | |
228 | Integer i,number | |
229 | Real*8 p1(5),p2(5),p3(5) | |
230 | Real*8 Gamma1(5),Gamma2(5) | |
231 | Real*8 factor_max,ABp,ABm | |
232 | Integer ierr | |
233 | data factor_max/1.D+00/ | |
234 | ierr =0 | |
235 | c------------------------------------------------------------------- | |
236 | If(iset.eq.0) Then | |
237 | c------------------------------------------------------------------- | |
238 | c this is the normal mode of operation | |
239 | c First, generate the kinematics | |
240 | ||
241 | p1(5)= M_Kp **2 | |
242 | p2(5)= M_pim**2 | |
243 | p3(5)= M_pi0**2 | |
244 | ||
245 | 10 continue | |
246 | call Evtfirst_step_Kpipi(p1,p2,p3) | |
247 | ||
248 | c Then, compute the amplitudes | |
249 | ||
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 | |
258 | c to be computed | |
259 | ||
260 | Do i=1,4 | |
261 | p1(i)= p_K_plus (i) | |
262 | p2(i)= p_pi_minus(i) | |
263 | p3(i)= p_gamma_1 (i) + p_gamma_2 (i) | |
264 | End Do | |
265 | p1(5)= M_Kp **2 | |
266 | p2(5)= M_pim**2 | |
267 | p3(5)= M_pi0**2 | |
268 | ||
269 | Call EvtCompute_Kpipi(p1,p2,p3, | |
270 | + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr) | |
271 | ||
272 | if(ierr.ne.0) Then | |
273 | Print*,'the provided kinematics are not physical' | |
274 | Print*,'ierr=',ierr | |
275 | Print*,'the program will stop' | |
276 | Stop | |
277 | endif | |
278 | c------------------------------------------------------------------- | |
279 | ElseIf(iset.gt.0) Then | |
280 | c------------------------------------------------------------------- | |
281 | c This is the pre-run mode of operation where initializations are | |
282 | c performed. | |
283 | ||
284 | factor_max= 0 | |
285 | call Evtset_constants(alpha_input, beta_input) | |
286 | p1(5)= M_Kp **2 | |
287 | p2(5)= M_pim**2 | |
288 | p3(5)= M_pi0**2 | |
289 | ||
290 | c pre-run | |
291 | Do number=1,iset | |
292 | ||
293 | 20 continue | |
294 | call Evtfirst_step_Kpipi(p1,p2,p3) | |
295 | ||
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 | |
301 | ||
302 | If(ABp.gt.factor_max) factor_max=ABp | |
303 | If(ABm.gt.factor_max) factor_max=ABm | |
304 | ||
305 | End Do | |
306 | c end of the pre-run | |
307 | ||
308 | factor_max=1.D+00/Dsqrt(factor_max) | |
309 | ||
310 | c------------------------------------------------------------------- | |
311 | End If | |
312 | c------------------------------------------------------------------- | |
313 | ||
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 | |
318 | ||
319 | if(iset.lt.0) return | |
320 | c P1,p2,p3 ---> random rotation in B rest frame | |
321 | ||
322 | Call EvtRotation(p1,1) | |
323 | Call EvtRotation(p2,0) | |
324 | Call EvtRotation(p3,0) | |
325 | ||
326 | C Desintegrate the pi_0 s | |
327 | ||
328 | Call EvtGammaGamma(p3,Gamma1,Gamma2) | |
329 | ||
330 | C Feed the output four vectors | |
331 | ||
332 | Do i=1,4 | |
333 | ||
334 | p_K_plus (i)=p1 (i) | |
335 | p_pi_minus(i)=p2 (i) | |
336 | p_gamma_1 (i)=Gamma1(i) | |
337 | p_gamma_2 (i)=Gamma2(i) | |
338 | End Do | |
339 | ||
340 | Return | |
341 | ||
342 | End | |
343 | ||
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 | |
350 | c --- P1(1) = Px | |
351 | c --- P1(2) = Py | |
352 | c --- P1(3) = Pz | |
353 | c --- P1(4) = E | |
354 | c --- P1(5) = M**2 | |
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 ---------------------------------------------------------- | |
361 | ||
362 | c ---------------------------------------------------------- | |
363 | c --- commons | |
364 | c ---------------------------------------------------------- | |
365 | #include "EvtGenModels/EvtBTo3pi.inc" | |
366 | c ---------------------------------------------------------- | |
367 | c --- Used Functions | |
368 | c ---------------------------------------------------------- | |
369 | ||
370 | real evtranf | |
371 | ||
372 | c ---------------------------------------------------------- | |
373 | c --- Variables in Argument | |
374 | c ---------------------------------------------------------- | |
375 | ||
376 | real*8 P1(5),P2(5),P3(5) | |
377 | ||
378 | c ---------------------------------------------------------- | |
379 | c --- Local Variables | |
380 | c ---------------------------------------------------------- | |
381 | ||
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 | |
386 | ||
387 | real*8 p1mom,p2mom,p3mom | |
388 | ||
389 | real*8 x, y, z, mass | |
390 | integer i | |
391 | ||
392 | Logical Phase_space | |
393 | data Phase_space/.false./ | |
394 | ||
395 | c ---------------------------------------------------------- | |
396 | c --- Computation | |
397 | c ---------------------------------------------------------- | |
398 | max_m12 = M_B**2 | |
399 | min_m12 = P1(5) + P2(5) + 2.*Dsqrt(p1(5)*p2(5)) | |
400 | ||
401 | max_m13 = M_B**2 | |
402 | min_m13 = P1(5) + P3(5) + 2.*Dsqrt(p1(5)*p3(5)) | |
403 | ||
404 | max_m23 = M_B**2 | |
405 | min_m23 = P2(5) + P3(5) + 2.*Dsqrt(p2(5)*p3(5)) | |
406 | ||
407 | 100 Continue | |
408 | ||
409 | c ---------------------------------------------------------- | |
410 | c --- Generation of the Mass of the Rho or Kstar | |
411 | c ---------------------------------------------------------- | |
412 | ||
413 | ||
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 ---------------------------------------------------------- | |
418 | ||
419 | z = 3.*evtranf() | |
420 | ||
421 | MC2 = M_B**2 | |
422 | ||
423 | If(z.lt.1.) Then | |
424 | c K*+ production | |
425 | If(Phase_space) Then | |
426 | m13 = evtranf()*(max_m13-min_m13)+min_m13 | |
427 | Else | |
428 | y = evtranf()*PI - PI/2. | |
429 | x = Dtan(y) | |
430 | mass = x*Gam_Kstarp/2. +Mass_Kstarp | |
431 | m13 = mass**2 | |
432 | End If | |
433 | ||
434 | m12 = evtranf()*(max_m12-min_m12)+min_m12 | |
435 | m23 = MC2 - m12 - m13 | |
436 | ||
437 | ElseIf(z.lt.2.) Then | |
438 | c K*0 production | |
439 | If(Phase_space) Then | |
440 | m12 = evtranf()*(max_m12-min_m12)+min_m12 | |
441 | Else | |
442 | y = evtranf()*PI - PI/2. | |
443 | x = Dtan(y) | |
444 | mass = x*Gam_Kstar0/2. +Mass_Kstar0 | |
445 | m12 = mass**2 | |
446 | End If | |
447 | ||
448 | m13 = evtranf()*(max_m13-min_m13)+min_m13 | |
449 | m23 = MC2 - m12 - m13 | |
450 | ||
451 | Else | |
452 | c rho- production | |
453 | If(Phase_space) Then | |
454 | m23 = evtranf()*(max_m23-min_m23)+min_m23 | |
455 | Else | |
456 | y = evtranf()*PI - PI/2. | |
457 | x = Dtan(y) | |
458 | mass = x*Gam_rho/2. +Mass_rho | |
459 | m23 = mass**2 | |
460 | End If | |
461 | ||
462 | m13 = evtranf()*(max_m13-min_m13)+min_m13 | |
463 | m12 = MC2 - m23 - m13 | |
464 | ||
465 | Endif | |
466 | ||
467 | c ---------------------------------------------------------- | |
468 | c --- Check that the physics is OK : | |
469 | c --- Are the invariant Masses in allowed ranges ? | |
470 | c ---------------------------------------------------------- | |
471 | ||
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 | |
475 | ||
476 | c ---------------------------------------------------------- | |
477 | c --- Are the Cosines of the angles between particles | |
478 | c --- Between -1 and +1 ? | |
479 | c ---------------------------------------------------------- | |
480 | ||
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) | |
484 | ||
485 | p1mom=p1(4)**2-P1(5) | |
486 | p2mom=p2(4)**2-P2(5) | |
487 | p3mom=p3(4)**2-P3(5) | |
488 | If(p1mom.lt.0) Go to 100 | |
489 | If(p2mom.lt.0) Go to 100 | |
490 | If(p3mom.lt.0) Go to 100 | |
491 | p1mom=Dsqrt(p1mom) | |
492 | p2mom=Dsqrt(p2mom) | |
493 | p3mom=Dsqrt(p3mom) | |
494 | ||
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 | |
501 | ||
502 | c ---------------------------------------------------------- | |
503 | c --- Filling the 5-vectors P1,P2,P3 | |
504 | c ---------------------------------------------------------- | |
505 | ||
506 | P3(1) = 0 | |
507 | P3(2) = 0 | |
508 | p3(3) = p3mom | |
509 | ||
510 | P1(3) = p1mom*cost13 | |
511 | P1(1) = p1mom*Dsqrt(1.D+00-cost13**2) | |
512 | p1(2) = 0. | |
513 | ||
514 | Do i=1,3 | |
515 | P2(i)=-p1(i)-p3(i) | |
516 | End do | |
517 | ||
518 | END | |
519 | ||
520 | ||
521 | c====================================================================== | |
522 | Subroutine EvtCompute_Kpipi(p1,p2,p3, | |
523 | + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr) | |
524 | c----------------------------------------------------------------------- | |
525 | IMPLICIT None | |
526 | #include "EvtGenModels/EvtBTo3pi.inc" | |
527 | ||
528 | ||
529 | Real*8 m12, m13, m23, W12, W13, W23, Wtot | |
530 | Real*4 evt_gmas | |
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) | |
534 | ||
535 | ||
536 | Integer ierr,iset | |
537 | Complex*16 BrightWagner,BreitWigner | |
538 | ||
539 | ierr = 0 | |
540 | c ---------------------------------------------------------------- | |
541 | C --- Account for the pole compensation | |
542 | c ---------------------------------------------------------------- | |
543 | ||
544 | m12 = evt_gmas(p1,p2) | |
545 | m13 = evt_gmas(p1,p3) | |
546 | m23 = evt_gmas(p2,p3) | |
547 | ||
548 | if(m12.lt.0. .or. m13.lt.0. .or. m23.lt.0.) Then | |
549 | ierr=1 | |
550 | Print*,'ierr = ',ierr | |
551 | return | |
552 | endif | |
553 | ||
554 | m12 = sqrt(m12) | |
555 | m13 = sqrt(m13) | |
556 | m23 = sqrt(m23) | |
557 | ||
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) | |
561 | ||
562 | if(iset.ge.0) Then | |
563 | Wtot = 1.D+00/Dsqrt(W12 + W13 + W23) | |
564 | else | |
565 | Wtot =1. | |
566 | Endif | |
567 | c ---------------------------------------------------------------- | |
568 | C --- Compute Breit-Wigners | |
569 | c ---------------------------------------------------------------- | |
570 | ||
571 | BW13=BrightWagner(p1,p3,p2,Mass_Kstarp,Gam_Kstarp,ierr) | |
572 | If(ierr.ne.0) Return | |
573 | BW12=BrightWagner(p1,p2,p3,Mass_Kstar0,Gam_Kstar0,ierr) | |
574 | If(ierr.ne.0) Return | |
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) | |
578 | If(ierr.ne.0) Return | |
579 | ||
580 | c ---------------------------------------------------------------- | |
581 | c - and | |
582 | C --- Build up the amplitudes | |
583 | c ---------------------------------------------------------------- | |
584 | ||
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 | |
588 | > + MatKstar0 * BW12 | |
589 | > + MatKrho * BW23 | |
590 | c Second, one computes the B0bar decay B0bar ->- K- pi+ pi0 | |
591 | MatB0bar = NatKstarp * BW13 | |
592 | > + NatKstar0 * BW12 | |
593 | > + NatKrho * BW23 | |
594 | ||
595 | ||
596 | c Pick up the Real and Imaginary parts | |
597 | ||
598 | Real_B0 = dreal(MatB0 )*Wtot | |
599 | Imag_B0 = Imag(MatB0 )*Wtot | |
600 | ||
601 | Real_B0bar = dreal(MatB0bar)*Wtot | |
602 | Imag_B0Bar = Imag(MatB0bar)*Wtot | |
603 | ||
604 | Return | |
605 | End | |
606 | c====================================================================== | |
607 | Function BrightWagner(p1,p2,p3,Mass,Width,ierr) | |
608 | c---------------------------------------------------------------------- | |
609 | IMPLICIT None | |
610 | #include "EvtGenModels/EvtBTo3pi.inc" | |
611 | Complex *16 BrightWagner,EvtRBW | |
612 | Integer ierr | |
613 | Logical RelatBW | |
614 | Data RelatBW/.true./ | |
615 | Real*8 Mass,Width,Am2Min | |
616 | ||
617 | c --------------------------------------------------------------- | |
618 | c --- Input Four Vectors | |
619 | c --------------------------------------------------------------- | |
620 | Real*8 p1(5),p2(5),p3(5) | |
621 | ||
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 | |
627 | Integer i | |
628 | Real *8 p1z,p1zcms12,e1cms12,p1cms12 | |
629 | ||
630 | c --------------------------------------------------------------- | |
631 | C --- Boost factor | |
632 | c --------------------------------------------------------------- | |
633 | ||
634 | BrightWagner=Dcmplx(0,0) | |
635 | ||
636 | ierr = 0 | |
637 | E12_2=(p1(4)+p2(4))**2 | |
638 | m12_2=E12_2 | |
639 | Do i=1,3 | |
640 | m12_2=m12_2-(p1(i)+p2(i))**2 | |
641 | End Do | |
642 | Argu = 1.D+00 - m12_2 / E12_2 | |
643 | ||
644 | If(argu.gt.0.) Then | |
645 | beta = Dsqrt(Argu) | |
646 | Else | |
647 | Print *,'Abnormal beta ! Argu = ',Argu | |
648 | ||
649 | Argu = 0. | |
650 | Beta = 0. | |
651 | End If | |
652 | ||
653 | If(m12_2.gt.0.)Then | |
654 | m12 = Dsqrt(m12_2) | |
655 | Else | |
656 | Print *,'Abnormal m12 ! m12_2 = ',m12_2 | |
657 | Print*,'p1 = ',p1 | |
658 | Print*,'p2 = ',p2 | |
659 | Print*,'p3 = ',p3 | |
660 | Stop | |
661 | End if | |
662 | ||
663 | gamma=Dsqrt(E12_2/m12_2) | |
664 | c --------------------------------------------------------------- | |
665 | C --- get the cosine in the B CMS | |
666 | c --------------------------------------------------------------- | |
667 | ||
668 | m13_2=(p1(4)+p3(4))**2 | |
669 | Do i=1,3 | |
670 | m13_2=m13_2-(p1(i)+p3(i))**2 | |
671 | End Do | |
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)) | |
676 | > / | |
677 | > (2.D+00*Dsqrt( (p1(4)**2-p1(5)) * (p3(4)**2-p3(5)) )) | |
678 | ||
679 | c --------------------------------------------------------------- | |
680 | C --- get the costet in the 1-2 CMS | |
681 | c --------------------------------------------------------------- | |
682 | ||
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 | |
688 | ||
689 | c --------------------------------------------------------------- | |
690 | C --- Build the Breit Wigner | |
691 | c --------------------------------------------------------------- | |
692 | ||
693 | If(RelatBW) then | |
694 | Am2Min = p1(5) + p2(5) + 2.*Dsqrt( p1(5)*p2(5) ) | |
695 | BrightWagner = coscms * EvtRBW(m12_2,Mass**2,Width,Am2Min) | |
696 | Else | |
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) | |
702 | End if | |
703 | ||
704 | Return | |
705 | 50 continue | |
706 | ierr = 2 | |
707 | Return | |
708 | End | |
709 | ||
710 | ||
711 | FUNCTION EvtRBW(s,Am2,Gam,Am2Min) | |
712 | ||
713 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
714 | COMPLEX*16 EvtRBW | |
715 | ||
716 | EvtRBW = 0. | |
717 | IF (s.le.Am2Min) RETURN | |
718 | ||
719 | G = Gam* (Am2/s) * ((s-Am2Min)/(Am2-Am2Min))**1.5 | |
720 | D = (Am2-s)**2 + s*G**2 | |
721 | X = Am2*(Am2-s) | |
722 | Y = Am2*SQRT(s)*G | |
723 | ||
724 | EvtRBW = DCMPLX(X/D,Y/D) | |
725 | ||
726 | RETURN | |
727 | END |