Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtBTo3piMPP.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: EvtBTo3piMPP.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 B -->-- 3pions 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] EVT3pions
32 C                 The routine to be called. Note that on the first call
33 C                 some initialization will be made, in particular the
34 C                 computation of a normalization factor designed to help
35 C                 the subsequent time-dependent generation of events.
36 C                 The normalisation done inside EVT3pions is such that
37 C                 at the level of the time implementation, the maximum
38 C                 time-dependant weight is unity : nothing is to be 
39 C                 computed to generate unity-weight events. The exact
40 C                 meaning of the above is made explicit in the HowToUse
41 C                 routine.
42 C [1] first_step_MPP
43 C                 Generation of the kinematics of the 3 pions 
44 C                 It uses the function ranf which is a random number 
45 C                 generator providing an uniform distribution 
46 C                 of Real*4 random number between 0 and 1       
47 C [2] compute
48 C                 return the amplitudes of the B0     -->-- 3pions
49 C                                   and of the B0_bar -->-- 3pions
50 C                 corrected for the generation mechanism
51 C                 The notations used are the ones of the paper of
52 C                 A. Snyder and H. Quinn [Phys.Rev.D48 (1993) 2139] 
53 C [3] BreitWigner
54 C                 compute the Breit-Wigner of the contributing rho s
55 C                 taking into account the cosine term linked to the
56 C                 zero-helicity of the rho s. There is three forms of
57 C                 Breit-Wigners available. The first one is the simple
58 C                 non-relativistic form, while the second and third
59 C                 ones are more involved relativistic expressions which
60 C                 in addition incorporate the contributions of three
61 C                 rho resonances [rho(770:1450:1700)]. The parameters
62 C                 used are the ones resulting from the ALEPH analysis
63 C                 which might be found in CERN-PPE:97013 (submitted to
64 C                 Zeitschrift für Physik C). The two parametrizations
65 C                 of the relativistic Breit-Wigners are the ones of     
66 C                 Kuhn-SantaMaria [default] and of Gounaris-Sakurai.            
67 C                 The default setting is the non-relativistic Breit-
68 C                 Wigner form. 
69 C [4] Set_constants
70 C                 Set the constants values, from the pion mass to the
71 C                 penguin contributions. It is called by EVT3pions
72 C
73 C And other routines which do not deserve comment here.
74 C===================================================================      
75 c      Implicit none     
76 c      Real Hmemor
77 c      Common/Pawc/Hmemor(1000000)
78 c      Call Hlimit(1000000)
79 c      call EvtHowToUse
80 c      call EvtHowToUse_MPP
81 c      call EvtHowToUse_P00
82 c      Stop
83 c      End
84
85       subroutine EvtHowToUse_MPP 
86       
87       Implicit none     
88       Real*8  alpha
89       Integer iset,number,j,N_gener,N_asked
90
91
92       Real*8  p_pi_minus(4),p_pi_plus_1(4),p_pi_plus_2(4)
93       Real*8  Real_Bp,Imag_Bp,Real_Bm,Imag_Bm
94       Real*8  Weight,Weight_max
95       Real*8  ABp,ABm
96       Real*8  m_rho12,m_rho13
97       Real*8  Evtranf,Tag
98       alpha      = 0.4
99       N_gener    = 0
100       N_asked    = 100000
101       
102       weight_max = 1.0
103            
104 c run : Simulation of the Anders Ryd Generator
105
106       Do number=1,N_asked ! weighted events as generated here 
107       
108       If(number.eq.1) then
109       iset=10000          ! 10^4 events are used to normalize amplitudes
110       Else
111       iset=0              ! iset must be reset to zero after the first call 
112       End If
113
114 c Here is the call to EVT3pions_MPP  !!!!!!!!!!!!!!!! 
115 c communication of data is done by argument only <<<<<<<<     
116        call EVT3pionsMPP(
117      + alpha,iset,
118      + p_pi_minus,p_pi_plus_1,p_pi_plus_2,
119      + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm) 
120      
121 C that is it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
122 c select the Tag 
123       Tag   =evtranf()
124 C get the relevant quantities            
125       ABp   = Real_Bp   **2 + Imag_Bp   **2
126       ABm   = Real_Bp   **2 + Imag_Bm   **2
127 c generate acording to the tag 
128       If(Tag.gt.0.5) Then
129       
130 c a Bm tag =>  the decay is one from a Bp      
131       Weight= ABp   
132       
133       Else
134            
135 c a Bp tag =>  the decay is one from a Bm 
136       Weight= ABm
137       
138       End If
139       If(Weight.Gt.evtranf()) Then
140 c----------------------------------------------------------------------------         
141 c unweighted event production
142 c---------------------------------------------------------------------------- 
143       N_gener=N_gener+1  
144 C here is just a Dalitz plot and a few prints
145
146       m_rho12=(p_pi_minus(4)+P_pi_plus_1(4))**2 
147       m_rho13=(p_pi_minus(4)+P_pi_plus_2(4))**2 
148       do j=1,3
149       m_rho12=m_rho12-(p_pi_minus(j)+P_pi_plus_1(j))**2 
150       m_rho13=m_rho13-(p_pi_minus(j)+P_pi_plus_2(j))**2  
151       end do                 
152                      
153 c here is a check that weight_max is one
154
155       If(Weight.gt.Weight_max) Then
156       Weight_max=Weight
157       Print*,' overweighted event found at weight = ',Weight_max           
158       End If
159       
160 c----------------------------------------------------------------------------      
161       End If    
162
163 c end of the loop over events     
164       End Do
165       
166       Print*,'number of unity-weight events generated : ',N_gener
167       Print*,'number of trials                        : ',N_asked
168       
169       
170       End
171 C===================================================================
172       subroutine Evt3pionsMPP(
173      + alpha_input,iset,
174      + p_pi_minus,p_pi_plus_1,p_pi_plus_2,
175      + Real_Bp,Imag_Bp,Real_Bm,Imag_Bm)
176 c-----------------------------------------------------------------
177 c       ----------------------------------------------------------
178 c       --- This is the routine to be called by the Main generator
179 c           to get the decay of B+ -->-- 3 pions
180 c       --- AND
181 c           to get the decay of B- -->-- 3 pions
182 C       For the sake of clarity, signs refers to B+ decay
183 c****************************************************************************                                 set to 0.
184 c       --- Outputs are :
185 c
186 c       ---               p_pi_minus     : the four momentum of the pi- 
187 c       ---               p_pi_plus_1    : the four momentum of the first  pi+
188 c       ---               p_pi_plus_2    : the four momentum of the second pi+   
189 c
190 c             Note that : the energy is stored in the fourth component
191 c                         the values are the ones of the B rest frame
192 c                         a random rotation has been applied 
193
194 c       ---               Real_Bp   : The real      part of the amplitude of 
195 c                                     the B+    ->- 3 pions decay                  
196 c       ---               Imag_Bp   : The imaginary part of the amplitude of 
197 c                                     the B+    ->- 3 pions decay
198
199 c       ---               Real_Bm   : The real      part of the amplitude of 
200 c                                     the B-    ->- 3 pions decay                  
201 c       ---               Imag_Bm   : The imaginary part of the amplitude of 
202 c      
203 c****************************************************************************
204 c-----------------------------------------------------------------      
205       Implicit none
206 #include "EvtGenModels/EvtBTo3pi.inc"
207       Real*8  alpha_input
208       Integer iset
209       Real*8  p_pi_minus(4),p_pi_plus_1(4),p_pi_plus_2(4)
210       Real*8  Real_Bp,Imag_Bp,Real_Bm,Imag_Bm
211       
212 c Working quantities
213       Integer i,number
214       Real*8  p1(5),p2(5),p3(5)
215       Real*8  factor_max,ABp,ABm
216       Integer ierr
217       data factor_max/1.D+00/
218       ierr =0
219 c-------------------------------------------------------------------       
220       If(iset.eq.0) Then
221 c------------------------------------------------------------------- 
222 c     this is the normal mode of operation 
223 c     First, generate the kinematics  
224
225       p1(5)= M_pim**2
226       p2(5)= M_pip**2
227       p3(5)= M_pip**2
228           
229  10   continue  
230       call Evtfirst_step_MPP(p1,p2,p3)
231       
232 c     Then, compute the amplitudes 
233       
234       Call EvtCompute_MPP(p1,p2,p3,
235      +    Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
236       if(ierr.ne.0) Go To 10
237 c-------------------------------------------------------------------       
238       ElseIf(iset.lt.0) Then
239 c-------------------------------------------------------------------
240 c     This is an user mode of operation where the kinematics is
241 c     provided by the user who only wants the corresponding amplitudes
242 c     to be computed
243
244       Do i=1,4
245       p1(i)= p_pi_minus (i)
246       p2(i)= p_pi_plus_1(i) 
247       p3(i)= p_pi_plus_2(i) 
248       End Do
249       p1(5)= M_pim**2
250       p2(5)= M_pip**2
251       p3(5)= M_pip**2
252       
253       Call EvtCompute_MPP(p1,p2,p3,
254      +    Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)
255       if(ierr.ne.0) Then
256         Print*,'the provided kinematics are not physical'
257         Print*,'ierr=',ierr
258         Print*,'the program will stop'
259         Stop
260       endif     
261 c-------------------------------------------------------------------       
262       ElseIf(iset.gt.0) Then
263 c-------------------------------------------------------------------
264 c     This is the pre-run mode of operation where initializations are
265 c     performed.
266  
267       factor_max= 0
268 c       changed by ryd April 27 1998. this causes an 
269 c       unaligned access probelm on OSF.
270 c       (should Beta not be passed as an argument?)    
271 c      call Evtset_constants(alpha_input,0.362)
272       call Evtset_constants(alpha_input,0.362D0)
273       p1(5)=M_pim**2
274       p2(5)=M_pip**2
275       p3(5)=M_pip**2
276       
277 c     pre-run
278       Do number=1,iset
279       
280  20   continue
281       call Evtfirst_step_MPP(p1,p2,p3)
282  
283       Call EvtCompute_MPP(p1,p2,p3,
284      +    Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)     
285        if(ierr.ne.0) Go To 20  
286       ABp   = Real_Bp   **2 + Imag_Bp   **2
287       ABm   = Real_Bm   **2 + Imag_Bm   **2
288       
289       If(ABp.gt.factor_max) factor_max=ABp
290       If(ABm.gt.factor_max) factor_max=ABm
291       
292       End Do
293 c     end of the pre-run 
294
295       factor_max=1.D+00/Dsqrt(factor_max)
296
297 c-------------------------------------------------------------------      
298       End If
299 c------------------------------------------------------------------- 
300     
301       Real_Bp   =Real_Bp    * factor_max
302       Imag_Bp   =Imag_Bp    * factor_max
303       Real_Bm   =Real_Bm    * factor_max
304       Imag_Bm   =Imag_Bm    * factor_max
305        
306       if(iset.lt.0) return      
307 c     P1,p2,p3 ---> random rotation in B rest frame
308
309       Call EvtRotation(p1,1)
310       Call EvtRotation(p2,0)
311       Call EvtRotation(p3,0)
312            
313 C     Feed the output four vectors
314
315       Do i=1,4
316       
317       p_pi_minus     (i)=p1    (i)      
318       p_pi_plus_1    (i)=p2    (i) 
319       p_pi_plus_2    (i)=p3    (i) 
320       
321       End Do
322       
323       Return
324       
325       End
326        
327 c===================================================================
328       subroutine Evtfirst_step_MPP(P1,P2,P3)
329 c-----------------------------------------------------------------
330 c       ----------------------------------------------------------
331 c       --- This routine generates the 5-vectors P1,P2,P3 
332 c       --- Associated respectively with the Pi+ and two Pi0 s 
333 c       ---             P1(1) = Px      
334 c       ---             P1(2) = Py
335 c       ---             P1(3) = Pz
336 c       ---             P1(4) = E
337 c       ---             P1(5) = M**2
338 c       ----------------------------------------------------------
339 c       ---     Input Four Vectors                                                            
340 C       ---     Particle [1] is the pi-                                                       
341 C       ---     Particle [2] is the pi+                                                      
342 C       ---     Particle [3] is the pi+                                                     
343 c       ----------------------------------------------------------
344
345 c       ----------------------------------------------------------              
346 c       --- commons             
347 c       ----------------------------------------------------------      
348 #include "EvtGenModels/EvtBTo3pi.inc"
349 c       ----------------------------------------------------------
350 c       --- Used Functions 
351 c       ----------------------------------------------------------
352
353         real*8 evtranf
354         
355
356 c       ----------------------------------------------------------
357 c       --- Variables in Argument
358 c       ----------------------------------------------------------
359
360          real*8 P1(5),P2(5),P3(5)
361          
362 c       ----------------------------------------------------------
363 c       --- Local Variables
364 c       ----------------------------------------------------------
365         
366         real*8 m12,min_m12, max_m12
367         real*8 m13,min_m13, max_m13
368         real*8 m23,min_m23, max_m23
369         Real*8 cost13,cost12,cost23
370         
371         real*8 p1mom,p2mom,p3mom
372       
373         real*8 x, y, z, mass
374         integer i
375
376         Logical Phase_space     
377         data Phase_space/.false./
378           
379 c       initialize to avoid warning on linux
380         mass=0.0
381
382 c       ----------------------------------------------------------
383 c       --- Computation
384 c       ----------------------------------------------------------
385         max_m12 = M_B**2 
386         min_m12 = P1(5) + P2(5) 
387       
388         max_m13 = M_B**2 
389         min_m13 = P1(5) + P3(5)
390       
391         max_m23 = M_B**2 
392         min_m23 = P2(5) + P3(5)
393         
394 100     Continue
395
396 c       ----------------------------------------------------------      
397 c       --- Generation of the Mass of the Rho(0) 
398 c       ----------------------------------------------------------
399
400         If(.not.Phase_space) Then      
401         y = evtranf()*PI - PI/2.
402         x = Dtan(y)
403         mass = x*Gam_rho/2. +Mass_rho
404         End If
405         
406 c       ----------------------------------------------------------      
407 c       --- z is the Flag needed to choose between the generation 
408 c       --- of a Rho0 = pi- pi+[1] or pi- pi+[2]
409 c       ----------------------------------------------------------
410
411         z = evtranf()
412         
413         if(z.lt..5) Then
414
415                 If(Phase_space) Then
416                 m12 = evtranf()*(max_m12-min_m12)+min_m12
417                 Else
418                 m12 = mass**2
419                 End If  
420                         
421         m13 = evtranf()*(max_m13-min_m13)+min_m13
422         m23 = MB2 - m12 - m13
423         
424         else 
425
426                 If(Phase_space) Then
427                 m13 = evtranf()*(max_m13-min_m13)+min_m13
428                 Else
429                 m13 =  mass**2
430                 End If
431                 
432         m12 = evtranf()*(max_m12-min_m12)+min_m12
433         m23 = MB2 - m12 - m13
434
435         endif
436
437 c       ----------------------------------------------------------      
438 c       --- Check that the physics is OK :
439 c       --- Are the invariant Masses in allowed ranges ?
440 c       ----------------------------------------------------------
441
442         If(m23.lt.min_m23.or.m23.gt.max_m23) Go to 100
443         If(m13.lt.min_m13.or.m13.gt.max_m13) Go to 100
444         If(m12.lt.min_m12.or.m12.gt.max_m12) Go to 100
445         
446 c       ----------------------------------------------------------
447 c       --- Are the Cosines of the angles between particles 
448 c       --- Between -1 and +1 ?
449 c       ----------------------------------------------------------
450       
451         P1(4)=(M_B**2+P1(5)-m23)/(2.*M_B)
452         P2(4)=(M_B**2+P2(5)-m13)/(2.*M_B)
453         P3(4)=(M_B**2+P3(5)-m12)/(2.*M_B)
454         
455         p1mom=p1(4)**2-P1(5)      
456         p2mom=p2(4)**2-P2(5)
457         p3mom=p3(4)**2-P3(5)
458         If(p1mom.lt.0) Go to 100
459         If(p2mom.lt.0) Go to 100
460         If(p3mom.lt.0) Go to 100
461         p1mom=Dsqrt(p1mom)
462         p2mom=Dsqrt(p2mom)
463         p3mom=Dsqrt(p3mom)
464         
465         cost13=(2.*p1(4)*p3(4)+P1(5)+p3(5)-m13)/(2.*p1mom*p3mom)
466         cost12=(2.*p1(4)*p2(4)+P1(5)+p2(5)-m12)/(2.*p1mom*p2mom)      
467         cost23=(2.*p2(4)*p3(4)+P2(5)+p3(5)-m23)/(2.*p2mom*p3mom)
468         If(Dabs(cost13).gt.1.) Go to 100
469         If(Dabs(cost12).gt.1.) Go to 100
470         If(Dabs(cost23).gt.1.) Go to 100
471
472 c       ----------------------------------------------------------
473 c       --- Filling the 5-vectors P1,P2,P3
474 c       ----------------------------------------------------------      
475
476         P3(1) = 0
477         P3(2) = 0
478         p3(3) = p3mom
479
480         P1(3) = p1mom*cost13
481         P1(1) = p1mom*Dsqrt(1.D+00-cost13**2)
482         p1(2) = 0.
483
484         Do i=1,3
485         P2(i)=-p1(i)-p3(i)
486         End do
487         
488         END
489         
490         
491 c======================================================================         
492       Subroutine EvtCompute_MPP(p1,p2,p3,
493      +           Real_Bp,Imag_Bp,Real_Bm,Imag_Bm,iset,ierr)                 
494 c-----------------------------------------------------------------------
495         IMPLICIT None  
496 #include "EvtGenModels/EvtBTo3pi.inc"
497
498                                                                                                  
499         Real*8 m12, m13, W12, W13, Wtot
500         Real*8 evt_gmas
501         Complex*16 MatBp,MatBm
502         Real*8  Real_Bp,Imag_Bp,Real_Bm,Imag_Bm                                                     
503         Real*8     p1(5),p2(5),p3(5)
504         real*8  ASHQ
505         Integer iset,ierr
506         Data    ASHQ/0.707107 D+00/
507                                                       
508         COMPLEX*16 Mat_rhop
509         Complex*16 BreitWigner
510         
511         ierr = 0
512 c       ----------------------------------------------------------------      
513 C       ---     Account for the pole compensation                                                    
514 c       ----------------------------------------------------------------  
515         if(evt_gmas(p1,p2).lt.0.or.evt_gmas(p1,p3)
516      +          .lt.0) Then
517                 ierr=1
518                 return
519         endif   
520
521         m12 = sqrt(evt_gmas(p1,p2))
522         m13 = sqrt(evt_gmas(p1,p3))
523  
524         W12 = (1./m12)*1./((Mass_rho - m12)**2+(Gam_rho/2.)**2)
525         W13 = (1./m13)*1./((Mass_rho - m13)**2+(Gam_rho/2.)**2)
526                      
527         If(iset.ge.0) Then  
528            Wtot = 1.D+00/Dsqrt(W12 + W13)                                        
529         else
530            Wtot = 1.
531         endif
532 c       ----------------------------------------------------------------      
533 C       ---     Compute Breit-Wigners                                                         
534 c       ----------------------------------------------------------------  
535     
536         Mat_rhop = BreitWigner(p1,p2,p3,ierr)
537      +   + BreitWigner(p1,p3,p2,ierr)                                                                               
538                                           
539 c       ----------------------------------------------------------------      
540 C       ---     Build up the amplitudes                                                       
541 c       ---------------------------------------------------------------- 
542 c       The factor ASHQ = 1./sqrt(2) is here just to stick to the notations
543 c       used by Art Snyder and Helen Quinn. It is irrelevant for the genera-
544 c       tion of events done here. 
545                                                                                                                                        
546         MatBp    = Mat_s2 * Mat_rhop * Wtot * ASHQ 
547                                                                                                    
548         MatBm    = Nat_s2 * Mat_rhop * Wtot * ASHQ
549         
550 c       Pick up the Real and Imaginary parts
551 c       changed Real to DBLE and Imag DIMAG to make it complie 
552 c       with Absoft and g77 (ryd)               
553
554         Real_Bp    = DBLE(MatBp)
555         Imag_Bp    = DIMAG(MatBp) 
556                                                    
557         Real_Bm    = DBLE(MatBm)
558         Imag_Bm    = DIMAG(MatBm)
559         
560         Return                                                                    
561         End