]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenModels/EvtBToKpipi.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding task for subtracting background after jet finding, used for all clustering...
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBToKpipi.F
CommitLineData
da0e9ce3 1C--------------------------------------------------------------------------
2C
3C Environment:
4C This software is part of the EvtGen package developed jointly
5C for the BaBar and CLEO collaborations. If you use all or part
6C of it, please give an appropriate acknowledgement.
7C
8C Copyright Information: See EvtGen/COPYRIGHT
9C Copyright (C) 1998 Caltech, UCSB
10C
11C Module: EvtBToKpipi.F
12C
13C Description:
14C
15C Modification history:
16C
17C DJL/RYD August 11, 1998 Module created
18C
19C------------------------------------------------------------------------
20C===================================================================
21C This package is providing a neutral B -->-- K pi pi decay generator
22C Its is composed of the following subroutines:
23C
24C [*] HowToUse
25C This is an How To Use routine where one may find the
26C implementation of the time dependance: That is to
27C say that it shows how the output of the routine is
28C supposed to be used in the mind of the authors.
29C
30C===================================================================
31C [0] EVTKpipi
32C The routine to be called. Note that on the first call
33C some initialization will be made, in particular the
34C computation of a normalization factor designed to help
35C the subsequent time-dependent generation of events.
36C The normalisation done inside EVTKpipi is such that
37C at the level of the time implementation, the maximum
38C time-dependant weight is unity : nothing is to be
39C computed to generate unity-weight events. The exact
40C meaning of the above is made explicit in the HowToUse
41C routine.
42C [1] first_step_Kpipi
43C Generation of the kinematics of the 3 prongs
44C It uses the function evtranf which is a random number
45C generator providing an uniform distribution
46C of Real*4 random number between 0 and 1
47C [2] compute
48C return the amplitudes of the B0 -->-- K+pi-pi0
49C corrected for the generation mechanism.
50c Note that this is a Tagging Mode. The CP conjugate
51c mode (B0bar -->-- K-pi+pi0) is treated at the same time.
52C [3] BreitWigner
53C compute the Breit-Wigner of the contributing K* and rho
54C taking into account the cosine term linked to the
55C zero-helicity of the spin-1 resonances. There is two
56c forms of Breit-Wigners available. The first one is the
57c simple non-relativistic form, while the second is the
58C relativistic expressions.
59C The default setting is the relativistic one.
60C [4] Set_constants
61C Set the constants values, from the pion mass to the
62C penguin contributions. It is called by EVTKpipi
63C
64C And other routines which do not deserve comment here.
65C===================================================================
66c Implicit none
67C Real Hmemor
68C Common/Pawc/Hmemor(1000000)
69C Call Hlimit(1000000)
70C Call EvtHowToUse_Kpipi
71C Stop
72C 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
97c 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
107c Here is the call to EVTKpipi !!!!!!!!!!!!!!!!
108c 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)
114C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
115c select the Tag
116 Tag =evtranf()
117c generate acording to the tag
118 If(Tag.gt.0.5) Then
119
120c a B0bar tag => the decay is one from a B0
121 Weight= Real_B0 **2 + Imag_B0 **2
122
123 Else
124
125c 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
131c----------------------------------------------------------------------------
132c unweighted event production
133c----------------------------------------------------------------------------
134 N_gener=N_gener+1
135C 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
146c here is the Dalitz plot when assuming the pion mass for the Kaon
147c 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
161c 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
168c----------------------------------------------------------------------------
169 End If
170
171c 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
178C===================================================================
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)
184c-----------------------------------------------------------------
185c ----------------------------------------------------------
186c --- This is the routine to be called by the Main generator
187c to get the decay of B0 -->-- K+ pi- pi0
188c The decay proceeeds through three channels:
189c a) B0 -->-- K*+ pi- ; K*+ -->-- K+ pi0
190c b) K*0 pi0 ; K*0bar -->-- K+ pi-
191c c) K- rho+ ; rho+ -->-- pi+ pi0
192c .) The K0 rho0 channel is not implemented since it does
193c not lead to Kpipi final state, but it is interesting
194c in itself.
195c It provides at the same time the CP conjugate decay
196c B0bar -->-- K- pi+ pi0
197c****************************************************************************
198c --- Outputs are :
199c
200c --- p_K_plus : the four momentum of the K+
201c --- p_pi_minus : ........................ pi-
202c --- p_gamma_1 : ........................ first gamma of the pi0
203c --- p_gamma_2 : ........................ second gamma of the pi0
204c
205c Note that : the energy is stored in the fourth component
206c the values are the ones of the B rest frame
207c a random rotation has been applied
208c
209c --- Real_B0 : The real part of the amplitude of the decay
210c B0 -->-- K+ pi- pi0
211c --- Imag_B0 : ... imaginary ..................................
212c similarly
213c --- Real_B0bar : The real part of the amplitude of the decay
214c B0bar -->-- K- pi+ pi0
215c --- Imag_B0bar : ... imaginary ..................................
216
217c****************************************************************************
218c-----------------------------------------------------------------
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
227c 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
235c-------------------------------------------------------------------
236 If(iset.eq.0) Then
237c-------------------------------------------------------------------
238c this is the normal mode of operation
239c 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
248c 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
253c----------------nedit EvtBto---------------------------------------------------
254 ElseIf(iset.lt.0) Then
255c-------------------------------------------------------------------
256c This is an user mode of operation where the kinematics is
257c provided by the user who only wants the corresponding amplitudes
258c 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
278c-------------------------------------------------------------------
279 ElseIf(iset.gt.0) Then
280c-------------------------------------------------------------------
281c This is the pre-run mode of operation where initializations are
282c 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
290c 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
306c end of the pre-run
307
308 factor_max=1.D+00/Dsqrt(factor_max)
309
310c-------------------------------------------------------------------
311 End If
312c-------------------------------------------------------------------
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
320c 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
326C Desintegrate the pi_0 s
327
328 Call EvtGammaGamma(p3,Gamma1,Gamma2)
329
330C 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
344c===================================================================
345 subroutine Evtfirst_step_Kpipi(P1,P2,P3)
346c-----------------------------------------------------------------
347c ----------------------------------------------------------
348c --- This routine generates the 5-vectors P1,P2,P3
349c --- Associated respectively with the Pi+ and two Pi0 s
350c --- P1(1) = Px
351c --- P1(2) = Py
352c --- P1(3) = Pz
353c --- P1(4) = E
354c --- P1(5) = M**2
355c ----------------------------------------------------------
356c --- Input Four Vectors
357C --- Particle [1] is the K+
358C --- Particle [2] is the pi-
359C --- Particle [3] is the pi0
360c ----------------------------------------------------------
361
362c ----------------------------------------------------------
363c --- commons
364c ----------------------------------------------------------
365#include "EvtGenModels/EvtBTo3pi.inc"
366c ----------------------------------------------------------
367c --- Used Functions
368c ----------------------------------------------------------
369
370 real evtranf
371
372c ----------------------------------------------------------
373c --- Variables in Argument
374c ----------------------------------------------------------
375
376 real*8 P1(5),P2(5),P3(5)
377
378c ----------------------------------------------------------
379c --- Local Variables
380c ----------------------------------------------------------
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
395c ----------------------------------------------------------
396c --- Computation
397c ----------------------------------------------------------
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
407100 Continue
408
409c ----------------------------------------------------------
410c --- Generation of the Mass of the Rho or Kstar
411c ----------------------------------------------------------
412
413
414c ----------------------------------------------------------
415c --- z is the Flag needed to choose between the generation
416c --- of a K*+, K*0 or rho- resonance
417c ----------------------------------------------------------
418
419 z = 3.*evtranf()
420
421 MC2 = M_B**2
422
423 If(z.lt.1.) Then
424c 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
438c 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
452c 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
467c ----------------------------------------------------------
468c --- Check that the physics is OK :
469c --- Are the invariant Masses in allowed ranges ?
470c ----------------------------------------------------------
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
476c ----------------------------------------------------------
477c --- Are the Cosines of the angles between particles
478c --- Between -1 and +1 ?
479c ----------------------------------------------------------
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
502c ----------------------------------------------------------
503c --- Filling the 5-vectors P1,P2,P3
504c ----------------------------------------------------------
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
521c======================================================================
522 Subroutine EvtCompute_Kpipi(p1,p2,p3,
523 + Real_B0,Imag_B0,Real_B0bar,Imag_B0Bar,iset,ierr)
524c-----------------------------------------------------------------------
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
540c ----------------------------------------------------------------
541C --- Account for the pole compensation
542c ----------------------------------------------------------------
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
567c ----------------------------------------------------------------
568C --- Compute Breit-Wigners
569c ----------------------------------------------------------------
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
575c If the rho is to be treated on the same footing as K* ==> use the line below
576c BW23=BrightWagner(p2,p3,p1,Mass_Rho ,Gam_Rho ,ierr)
577 BW23=BreitWigner(p2,p3,p1,ierr)
578 If(ierr.ne.0) Return
579
580c ----------------------------------------------------------------
581c - and
582C --- Build up the amplitudes
583c ----------------------------------------------------------------
584
585c Here come the relative amplitudes of the three decay channels
586c First, one computes the B0 decay B0 ->- K+ pi- pi0
587 MatB0 = MatKstarp * BW13
588 > + MatKstar0 * BW12
589 > + MatKrho * BW23
590c Second, one computes the B0bar decay B0bar ->- K- pi+ pi0
591 MatB0bar = NatKstarp * BW13
592 > + NatKstar0 * BW12
593 > + NatKrho * BW23
594
595
596c 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
606c======================================================================
607 Function BrightWagner(p1,p2,p3,Mass,Width,ierr)
608c----------------------------------------------------------------------
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
617c ---------------------------------------------------------------
618c --- Input Four Vectors
619c ---------------------------------------------------------------
620 Real*8 p1(5),p2(5),p3(5)
621
622c ---------------------------------------------------------------
623C --- intermediate variables
624c ---------------------------------------------------------------
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
630c ---------------------------------------------------------------
631C --- Boost factor
632c ---------------------------------------------------------------
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)
664c ---------------------------------------------------------------
665C --- get the cosine in the B CMS
666c ---------------------------------------------------------------
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
679c ---------------------------------------------------------------
680C --- get the costet in the 1-2 CMS
681c ---------------------------------------------------------------
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
689c ---------------------------------------------------------------
690C --- Build the Breit Wigner
691c ---------------------------------------------------------------
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