← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to models/taudecay_UFO/Fortran/functions.f

  • Committer: olivier Mattelaer
  • Date: 2014-06-11 19:47:48 UTC
  • mfrom: (253.1.38 2.1.2)
  • mto: (253.1.70 2.1.2)
  • mto: This revision was merged to the branch mainline in revision 254.
  • Revision ID: olivier.mattelaer@uclouvain.be-20140611194748-nhxx3210gudd1qfw
merge with 2.1.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
models/taudecay_UFO/Fortran/functions.f
 
1
      
 
2
      subroutine taudecay_param()
 
3
      implicit none
 
4
 
 
5
      DOUBLE PRECISION PIM,PI0
 
6
      COMMON/PARAM_PI/PIM,PI0
 
7
 
 
8
      DOUBLE PRECISION ROM,ROG,ROM1,ROG1
 
9
      COMMON/PARAM_RHO/ROM,ROG,ROM1,ROG1
 
10
 
 
11
      DOUBLE PRECISION a1M,a1G,fpi
 
12
      COMMON/PARAM_A1/a1M,a1G,fpi
 
13
 
 
14
      PIM=0.13957018d0
 
15
      PI0=0.1349766d0
 
16
 
 
17
      ROM=0.77549d0
 
18
      ROG=0.1491d0
 
19
      ROM1=1.465d0
 
20
      ROG1=0.40d0
 
21
 
 
22
      a1M=1.23d0
 
23
      a1G=0.42d0
 
24
      fpi=0.13041d0
 
25
 
 
26
      return
 
27
      end
 
28
 
 
29
      double complex function FFCT2(S)
 
30
      implicit none
 
31
      double complex S,Frho
 
32
      external Frho
 
33
      S=real(S)
 
34
      FFCT2=Frho(S,1)
 
35
      return
 
36
      end
 
37
 
 
38
      double complex function FFCT3(S)
 
39
      implicit none
 
40
      double complex S,Fa1
 
41
      external Fa1
 
42
      S=real(S)
 
43
      FFCT3=Fa1(S,1)
 
44
      return
 
45
      end
 
46
 
 
47
      double complex function FFCT3F1(S)
 
48
      implicit none
 
49
      double complex S,Brho
 
50
      external Brho
 
51
      S=real(S)
 
52
      FFCT3F1=Brho(S,1)
 
53
      return
 
54
      end
 
55
 
 
56
      double complex function FFCT3F0(S)
 
57
      implicit none
 
58
      double complex S,Brho
 
59
      external Brho
 
60
      S=real(S)
 
61
      FFCT3F0=Brho(S,0)
 
62
      return
 
63
      end
 
64
 
 
65
      DOUBLE COMPLEX FUNCTION Frho(S,mode) ! Frho=f2
 
66
      IMPLICIT NONE
 
67
      DOUBLE COMPLEX Brho
 
68
      DOUBLE PRECISION S
 
69
      integer mode
 
70
      EXTERNAL Brho
 
71
      Frho=dsqrt(2d0)*Brho(S,mode)
 
72
      RETURN
 
73
      END
 
74
 
 
75
      DOUBLE COMPLEX FUNCTION Brho(S,mode)
 
76
      IMPLICIT NONE
 
77
      DOUBLE COMPLEX BWrho
 
78
      DOUBLE PRECISION S
 
79
      integer mode
 
80
      EXTERNAL BWrho
 
81
      DOUBLE PRECISION BETA1
 
82
      PARAMETER (BETA1=-0.145d0) !TAUOLA parameters
 
83
      DOUBLE PRECISION ROM,ROG,ROM1,ROG1
 
84
      COMMON/PARAM_RHO/ROM,ROG,ROM1,ROG1
 
85
      !!!PARAMETER (ROM=0.773d0,ROG=0.145d0)   !TAUOLA parameters
 
86
      !!!PARAMETER (ROM1=1.370d0,ROG1=0.510d0) !TAUOLA parameters
 
87
      call taudecay_param()
 
88
      Brho=(BWrho(S,ROM,ROG,mode)+BETA1*BWrho(S,ROM1,ROG1,mode))/(1d0+BETA1)
 
89
      RETURN
 
90
      END
 
91
 
 
92
      DOUBLE COMPLEX FUNCTION BWrho(S,M,G,mode) !Breit-Wigner for the rho mode
 
93
      IMPLICIT NONE
 
94
      DOUBLE PRECISION S,M,G
 
95
      integer mode
 
96
      DOUBLE PRECISION QS,QM,W,GS,pi1,pi2
 
97
      double precision klambda
 
98
      external klambda
 
99
      DOUBLE PRECISION PIM,PI0
 
100
      COMMON/PARAM_PI/PIM,PI0
 
101
      !!!PARAMETER (PIM=0.139d0,PI0=0.139d0)  !TAUOLA parameters
 
102
      call taudecay_param()
 
103
      pi1=pim
 
104
      if (mode.eq.1) then
 
105
         pi2=pi0                ! rho- decay
 
106
      else
 
107
         pi2=pim                ! rho0 decay
 
108
      endif
 
109
      IF (S.GT.(PI1+PI2)**2) THEN
 
110
         QS=dsqrt(klambda(PI1**2/S,PI2**2/S))
 
111
         QM=dsqrt(klambda(PI1**2/M**2,PI2**2/M**2))
 
112
         W=DSQRT(S)
 
113
         GS=G*(W/M)*(QS/QM)**3
 
114
      ELSE
 
115
         GS=0d0
 
116
      ENDIF
 
117
      BWrho=-M**2/DCMPLX(S-M**2,W*GS)
 
118
      RETURN
 
119
      END
 
120
 
 
121
      DOUBLE COMPLEX FUNCTION Fa1(S,mode) ! Fa1=f3
 
122
      IMPLICIT NONE
 
123
      DOUBLE COMPLEX BWa1
 
124
      DOUBLE PRECISION S
 
125
      integer mode
 
126
      EXTERNAL BWa1
 
127
      DOUBLE PRECISION a1M,a1G,fpi
 
128
      COMMON/PARAM_A1/a1M,a1G,fpi
 
129
      call taudecay_param()
 
130
      Fa1=4d0/3d0/fpi *BWa1(S,a1M,a1G,mode)
 
131
      RETURN
 
132
      END
 
133
 
 
134
      DOUBLE COMPLEX FUNCTION BWa1(S,M,G,mode) !Breit-Wigner for the a1 mode
 
135
      IMPLICIT NONE
 
136
      DOUBLE PRECISION S,M,G
 
137
      integer mode
 
138
      DOUBLE PRECISION GS,GFUN,W,pi3
 
139
      DOUBLE PRECISION PIM,PI0
 
140
      COMMON/PARAM_PI/PIM,PI0
 
141
      !!!PARAMETER (PIM=0.139d0,PI0=0.139d0) !TAUOLA parameters
 
142
      call taudecay_param()
 
143
      W=DSQRT(S)
 
144
      if (mode.eq.1) then
 
145
         pi3=2d0*pi0+pim        ! rho- decay
 
146
      else
 
147
         pi3=3d0*pim            ! rho0 decay
 
148
      endif      
 
149
      IF (S.GT.pi3**2) THEN
 
150
        GS=G*(W/M)*GFUN(S,mode)/GFUN(M**2,mode)
 
151
      ELSE
 
152
        GS=0d0
 
153
      ENDIF
 
154
      BWa1=-M**2/DCMPLX(S-M**2,W*GS)
 
155
      RETURN
 
156
      END
 
157
 
 
158
      DOUBLE PRECISION FUNCTION GFUN(QKWA,mode) ! running width for a1
 
159
      IMPLICIT NONE
 
160
      DOUBLE PRECISION QKWA
 
161
      integer mode
 
162
      double precision pi1,pi3
 
163
      DOUBLE PRECISION PIM,PI0
 
164
      COMMON/PARAM_PI/PIM,PI0
 
165
      DOUBLE PRECISION ROM,ROG,ROM1,ROG1
 
166
      COMMON/PARAM_RHO/ROM,ROG,ROM1,ROG1
 
167
      call taudecay_param()
 
168
      if (mode.eq.1) then
 
169
         pi3=2d0*pi0+pim        ! rho- decay
 
170
         pi1=pi0
 
171
      else
 
172
         pi3=3d0*pim            ! rho0 decay
 
173
         pi1=pim
 
174
      endif        
 
175
      IF (QKWA.LT.(ROM+pi1)**2) THEN
 
176
        GFUN=4.1d0/QKWA*(QKWA-pi3**2)**3*(1d0-3.3d0*(QKWA-pi3**2)+5.8d0*(QKWA-pi3**2)**2)
 
177
      ELSE
 
178
        GFUN=1.623d0+10.38d0/QKWA-9.32d0/QKWA**2+0.65d0/QKWA**3
 
179
      ENDIF
 
180
      END
 
181
 
 
182
      double precision function klambda(a,b)
 
183
      implicit none 
 
184
      double precision a,b,c
 
185
      klambda=1d0+a**2+b**2-2d0*(a*b+b+a)
 
186
      return
 
187
      end
 
188
 
 
189
      DOUBLE PRECISION FUNCTION pSumDot(P1,P2,dsign) ! invariant mass
 
190
      IMPLICIT NONE
 
191
      double precision p1(0:3),p2(0:3),dsign
 
192
      integer i
 
193
      double precision ptot(0:3)
 
194
      double precision pdot
 
195
      external pdot
 
196
      do i=0,3
 
197
        ptot(i)=p1(i)+dsign*p2(i)
 
198
      enddo
 
199
      pSumDot = pdot(ptot,ptot)
 
200
      RETURN
 
201
      END
 
202
 
 
203
      double precision function pdot(p1,p2) !4-vector dot product
 
204
      implicit none
 
205
      double precision p1(0:3),p2(0:3)
 
206
      pdot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
 
207
      if(dabs(pdot).lt.1d-6)then ! solve numerical problem 
 
208
         pdot=0d0
 
209
      endif
 
210
      return
 
211
      end
 
212
 
 
213
      subroutine psum(p1,p2, q) !4-vector sum
 
214
      implicit none
 
215
      double precision p1(0:3),p2(0:3), q(0:3)
 
216
      q(0)=p1(0)+p2(0)
 
217
      q(1)=p1(1)+p2(1)
 
218
      q(2)=p1(2)+p2(2)
 
219
      q(3)=p1(3)+p2(3)
 
220
      return
 
221
      end
 
222
 
 
223
      subroutine psub(p1,p2, q) !4-vector subtract
 
224
      implicit none
 
225
      double precision p1(0:3),p2(0:3), q(0:3)
 
226
      q(0)=p1(0)-p2(0)
 
227
      q(1)=p1(1)-p2(1)
 
228
      q(2)=p1(2)-p2(2)
 
229
      q(3)=p1(3)-p2(3)
 
230
      return
 
231
      end
 
232
 
 
233
      subroutine psum3(p1,p2,p3, q) !4-vector sum
 
234
      implicit none
 
235
      double precision p1(0:3),p2(0:3),p3(0:3), q(0:3)
 
236
      q(0)=p1(0)+p2(0)+p3(0)
 
237
      q(1)=p1(1)+p2(1)+p3(1)
 
238
      q(2)=p1(2)+p2(2)+p3(2)
 
239
      q(3)=p1(3)+p2(3)+p3(3)
 
240
      return
 
241
      end