PHOS/reconstruction/shower_functions_v3.f

   1 c ==================================================================
   2
   3         function coornew(el,er,angl,mmm)
   4         mmnew=mmconvert(mmm)
   5         esum=el+er
   6         ex=el/esum
   7
   8         if(esum.le.0.) then
   9                 write (*,*) ' coornew diagnostic bad esum=',esum
  10                 return
  11         endif
  12
  13         itip=newtipcell(mmnew)
  14         if(itip.eq.1) then
  15                 xx=xcryexgl(ex)
  16         elseif(itip.eq.2) then
  17 c               xx=xcryexcr(ex)
  18                 xx=xcryexcrangl(ex,angl)
  19         endif
  20         coornew=xx
  21         return
  22         end
  23
  24         function xcryexgl(ar)
  25 c FOR LEAD GLASS
  26         data p1,p2 /1.480, 1.600/
  27 c ar - energy fraction right from boundary
  28
  29
  30         if(ar.gt.0.5) then
  31                 a=1.-ar
  32                 s=-1.
  33         else
  34                 a=ar
  35                 s=1.
  36         endif
  37         if(a.le.0.) a=0.0001
  38
  39         y=alog(2.*a)
  40         x=p2*y**2-p1*y
  41 c lead glass for my geant
  42 c       x=x/0.55
  43 c lead glass for my geant 2*2 cells
  44         x=x/0.65
  45 c corrections for crystalls 2*2, systematic is better than 0.25 mm for 3Gev geant showers
  46 c94     x=x*0.54-0.25
  47 c end corrections for crystalls
  48         xcryexgl=s*x
  49         return
  50         end
  51
  52         function xcryexcr(ar)
  53 c FOR CRYSTALS
  54 c       data p1,p2 /1.480, 1.600/
  55 c  very new 26.02.98
  56         data p1,p2 /1.700, 1.150/
  57 c ar - energy fraction right from boundary
  58         if(ar.gt.0.5) then
  59                 a=1.-ar
  60                 s=-1.
  61         else
  62                 a=ar
  63                 s=1.
  64         endif
  65         if(a.le.0.) a=0.0001
  66         y=alog(2.*a)
  67         x=p2*y**2-p1*y
  68         xcryexcr=s*x
  69         return
  70         end
  71
  72         function xcryexcrangl(ar,angl)
  73 c FOR CRYSTALS ALICE parametrisation 23.02.98
  74 c       data p1,p2 /1.480, 1.600/
  75         data p1,p2 /1.700, 1.150/
  76 c ar - energy fraction right from boundary
  77         aangl=abs(angl)
  78 c        write (*,*) ' ar,angl ',ar,angl
  79         if(angl.ge.0.) then
  80                 pp1=p1
  81                 pp2=p2-0.7*aangl
  82         else
  83                 pp1=p1+7.6*aangl
  84                 pp2=p2-0.9*aangl
  85         endif
  86 c        if(angl.ge.0.) then
  87 c                pp1=p1+7.6*aangl
  88 c                pp2=p2-0.9*aangl
  89 c        else
  90 c                pp1=p1
  91 c                pp2=p2-0.7*aangl
  92 c        endif
  93         if(ar.gt.0.5) then
  94                 a=1.-ar
  95                 s=-1.
  96         else
  97                 a=ar
  98                 s=1.
  99         endif
 100         if(a.le.0.) a=0.0001
 101         y=alog(2.*a)
 102         x=pp2*y**2-pp1*y
 103         xcryexcrangl=s*x
 104         return
 105         end
 106 c ==============================================================
 107         real function cumulnew(xc,itip)
 108         real xc
 109         integer itip
 110
 111         if(itip.eq.1) then
 112                 a=cucryexgl(xc)
 113         elseif(itip.eq.2) then
 114                 a=cucryexcr(xc)
 115         endif
 116         cumulnew=a
 117         return
 118         end
 119
 120         real function dcumulnew(xc,itip)
 121         real xc
 122         integer itip
 123
 124         if(itip.eq.1) then
 125                 a=dcucryexgl(xc)
 126         elseif(itip.eq.2) then
 127                 a=dcucryexcr(xc)
 128         endif
 129         dcumulnew=a
 130         return
 131         end
 132
 133         real function cucryexcr(xc)
 134 c fit for ALICE 94 RUN DATA
 135 c       data p1,p2 /2.526, 1.104/
 136 c94     data p1,p2 /2.300, 2.044/
 137 c crystals for my geant
 138 c       data p1,p2 /1.480, 1.600/
 139 c  very new 26.02.98
 140         data p1,p2 /1.700, 1.150/
 141         x=abs(xc)
 142         a=0.5*exp((p1-sqrt(p1**2+4.*p2*x))/2./p2)
 143         if(xc.ge.0.) then
 144                 cucryexcr=1.-a
 145         else
 146                 cucryexcr=a
 147         endif
 148         return
 149         end
 150
 151         real function dcucryexcr(xc)
 152 c fit for ALICE 94 RUN DATA
 153 c       data p1,p2 /2.526, 1.104/
 154 c       data p1,p2 /2.300, 2.044/
 155 c crystals for my geant
 156 c       data p1,p2 /1.480, 1.600/
 157 c  very new 26.02.98
 158         data p1,p2 /1.700, 1.150/
 159         x=abs(xc)
 160
 161         rr=sqrt(p1**2+4.*p2*x)
 162         da=0.5*exp((p1-rr)/2./p2)/rr
 163         dcucryexcr=da
 164
 165         return
 166         end
 167
 168         real function cucryexgl(xc)
 169 c fit for ALICE 94 RUN DATA
 170 c       data p1,p2 /2.526, 1.104/
 171 c94     data p1,p2 /2.300, 2.044/
 172 c crystals for my geant
 173         data p1,p2 /1.480, 1.600/
 174         x=abs(xc)
 175 c lead glass for my geant
 176         x=x*0.55
 177
 178         a=0.5*exp((p1-sqrt(p1**2+4.*p2*x))/2./p2)
 179         if(xc.ge.0.) then
 180                 cucryexgl=1.-a
 181         else
 182                 cucryexgl=a
 183         endif
 184         return
 185         end
 186
 187         real function dcucryexgl(xc)
 188 c fit for ALICE 94 RUN DATA
 189 c       data p1,p2 /2.526, 1.104/
 190 c       data p1,p2 /2.300, 2.044/
 191 c crystals for my geant
 192         data p1,p2 /1.480, 1.600/
 193         x=abs(xc)
 194 c lead glass for my geant
 195         x=x*0.55
 196
 197         rr=sqrt(p1**2+4.*p2*x)
 198         da=0.5*exp((p1-rr)/2./p2)/rr
 199 c lead glass for my geant
 200         da=da*0.55
 201         dcucryexgl=da
 202
 203         return
 204         end
 205 c ==============================================================
 206
 207         real function ampcelnew(e,vx,vc,vb,mmm)
 208         common /comkey/ keykey
 209         real vx(5),vc(3),vb(3)
 210         integer mmm
 211
 212         mmnew=mmconvert(mmm)
 213         itip=newtipcell(mmnew)
 214         xg=vx(1)
 215         yg=vx(2)
 216         xc=vc(1)
 217         yc=vc(2)
 218         if(keykey.eq.0) then
 219                 a=aclnew(xg-xc,yg-yc,vb,itip)
 220         elseif(keykey.eq.1) then
 221 c new with angle
 222                 a=acl3(e,vx,vc,vb)
 223         else
 224                 write (*,*) ' Net takoi keykey !! ',keykey
 225                 a=0.
 226         endif
 227         ampcelnew=a
 228         return
 229         end
 230
 231         real function dampcelnew(nd,e,vx,vc,vb,mmm)
 232         common /comkey/keykey
 233         real vx(5),vc(3),vb(3)
 234         integer mmm
 235
 236         mmnew=mmconvert(mmm)
 237         itip=newtipcell(mmnew)
 238 c old no angle
 239         if(keykey.eq.0) then
 240                 xg=vx(1)
 241                 yg=vx(2)
 242                 xc=vc(1)
 243                 yc=vc(2)
 244            if(nd.eq.1) then
 245                 da=dxaclnew(xg-xc,yg-yc,vb,itip)
 246            elseif(nd.eq.2) then
 247                 da=dxaclnew(yg-yc,xg-xc,vb,itip)
 248            else
 249                 write (*,*) ' net takoi nd ',nd
 250            endif
 251         elseif(keykey.eq.1) then
 252 c new with angle
 253                 da=dacl3(nd,e,vx,vc,vb)
 254         else
 255                 write (*,*) ' net takoi keykey ',keykey
 256                 da=0.
 257         endif
 258         dampcelnew=da
 259
 260         return
 261         end
 262
 263         FUNCTION aclnew(XG,YG,DA,mmn)
 264 C
 265 C     SHOWER FRACTION IN THE CELL
 266 C       XG,YG GAMMA COOR. X=0,Y=0 - CENTER OF THE CELL
 267 C       DA - SIZE OF THE CELL (MM)
 268 C
 269         real da(3)
 270         integer mmn
 271
 272         EXTERNAL a2fnew
 273 C       SX=SIGN(1.,XG)
 274 C       SY=SIGN(1.,YG)    FOR ANGLE
 275         X=-ABS(XG)
 276         Y=-ABS(YG)
 277         Dx=DA(1)/2.
 278         dy=da(2)/2.
 279         XP=X+Dx
 280         YP=Y+Dy
 281         XM=X-Dx
 282         YM=Y-Dy
 283         aclnew=a2fnew(XP,YP,mmn)+a2fnew(XM,YM,mmn)-
 284      -         a2fnew(XP,YM,mmn)-a2fnew(XM,YP,mmn)
 285 c       write (*,*) ' xg,yg,dx,dy,acl ',xg,yg,dx,dy,acl
 286         RETURN
 287         END
 288 C
 289         FUNCTION dxaclnew(XG,YG,DA,mmn)
 290 C
 291 C       DXACL=D(ACL(X,Y,DA))/D(X)
 292 C
 293         integer mmn
 294         EXTERNAL dxa2fnew
 295         SX=SIGN(1.,XG)
 296 C       SY=SIGN(1.,YG)    FOR ANGLE
 297         X=-ABS(XG)
 298         Y=-ABS(YG)
 299         D=DA/2.
 300         XP=X+D
 301         YP=Y+D
 302         XM=X-D
 303         YM=Y-D
 304         A=dxa2fnew(XP,YP,mmn)+dxa2fnew(XM,YM,mmn)-
 305      -    dxa2fnew(XP,YM,mmn)-dxa2fnew(XM,YP,mmn)
 306         dxaclnew=-SX*A
 307         RETURN
 308         END
 309
 310         FUNCTION a2fnew(X,Y,mmn)
 311         integer mmn
 312         EXTERNAL afnew
 313         AX=afnew(X,mmn)
 314         AY=afnew(Y,mmn)
 315         a2fnew=aa2fz(ax,ay)
 316         return
 317         END
 318
 319         FUNCTION dxa2fnew(X,Y,mmn)
 320         integer mmn
 321         EXTERNAL afnew,dafnew
 322         AX=afnew(X,mmn)
 323         AY=afnew(Y,mmn)
 324         dxa2fnew=daa2fz(ax,ay)*dafnew(x,mmn)
 325         return
 326         END
 327
 328         FUNCTION afnew(x,mmn)
 329         afnew=cumulnew(x,mmn)
 330         return
 331         end
 332
 333         FUNCTION dafnew(x,mmn)
 334         dafnew=dcumulnew(x,mmn)
 335         return
 336         end
 337
 338         function sigwlgam0(ewl,egam)
 339 c        write (*,*) ' ewl,egam ',ewl,egam
 340         if(egam.le.0.) then
 341                 sigwlgam0=0.
 342         else
 343                 sqwl=sqrt(ewl)
 344                 fr=ewl/egam
 345                 if(fr.gt.0.8) fr=0.8
 346 c               sigwlgam0=0.36*(1.-fr/0.88)*sqwl
 347         xp=abs(fr-0.05)
 348         sfun=0.5*sin(3.14/0.7*xp)+0.5*sin(3.14/0.76*sqrt(xp*0.7))
 349 c        sfun=2.5*egam*sfun+0.13-0.13*xp
 350 c        fm=(1.+exp(-egam/0.2))
 351 c        if(fm.lt.1.) fm=1.
 352 c        if(fm.gt.2.) fm=2.
 353         fm=2.5
 354         sfun=egam*fm*sfun+0.13-0.13*xp
 355 c        sigwlgam0=sqrt(egam)*sfun*sqwl
 356         sigwlgam0=0.1/sqrt(egam)*sfun*sqwl
 357 c        write (*,*) ' fr sfun ',fr,sfun
 358 c        write (*,*) ' sigwlgam0 ',sigwlgam0,ewl,egam
 359
 360         endif
 361         return
 362         end
 363
 364         function dispwlm(ewl,m)
 365
 366         common /comsigma/sigmaph,sigmapd,sigphsq,sigpdsq
 367
 368
 369         dispph=sigphsq*ewl
 370         disppd=sigpdsq
 371         dispwlm=dispph+disppd
 372
 373         return
 374         end
 375
 376         real function zmidlshower(ee,m)
 377         include 'comgam.for'
 378         zmidlshower = 46.
 379 c        zmidlshower = 50.
 380         return
 381         end
 382
 383