Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBToKpipi.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: 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
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*8  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
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*8 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*8 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