PHOS/reconstruction/gamma1_v3.f

   1
   2         subroutine hiwlnew(ng,mg,hidmin,hidmax,prb,mmax)
   3         integer mg(*)
   4         include 'comalwl.for'
   5         include 'comgam.for'
   6         data wlmin /0.02/
   7
   8         hi=0.
   9         himax=0.
  10         hidmax=0.
  11         hidmin=1000000.
  12         if(ng.eq.0) return
  13         nb=0.
  14         do n=1,ng
  15         m=mg(n)
  16         if(wl(m).gt.wlmin) then
  17         nb=nb+1
  18 c        write (*,*)  'swlt(m),swlt0(m)',swlt(m),swlt0(m)
  19         dspe=swlt(m)+swlt0(m)
  20         de=(wl(m)-wlt(m))**2
  21         dde=wl(m)-wlt(m)
  22 c fix by Federico Carminati CERN-EP, just to avoid a crash
  23 c it should be checked by the authors
  24         dspe=max(dspe,1e-11)
  25         hig=de/dspe
  26         hidg=dde/sqrt(dspe)
  27 c       if(wlt(m).eq.0) write (*,*) ' n,m,wl(m),hi ',n,m,wl(m),sqrt(hig)
  28         hi=hi+hig
  29
  30 c       if(hig.gt.himax) then
  31 c               himax=hig
  32 c               mmax=m
  33 c       endif
  34
  35         if(hidg.gt.hidmax) then
  36                 hidmax=hidg
  37                 mmax=m
  38         endif
  39         if(hidg.lt.hidmin) then
  40                 hidmin=hidg
  41         endif
  42
  43         endif
  44         enddo
  45         nbp=nb-3*kg
  46 c        nbp=nb
  47
  48         if(nbp.gt.0) then
  49                 aln=alog(float(nbp))
  50 c               himax0=-0.088+0.548*aln-0.01396*aln**2
  51
  52                 hidmax0=0.318+0.5942*aln-0.0258*aln**2
  53                 hidmax=hidmax-hidmax0
  54                 prb=prob(hi,nbp)
  55         else
  56                 himax0=0.
  57                 hidmax0=0.
  58                 prb=1.
  59         endif
  60 c        write (*,*) ' kg,nb, nbp,prb ',kg,nb,nbp,prb
  61
  62 c       hi=sqrt(h)
  63 c       himax=sqrt(himax)
  64         return
  65         end
  66         subroutine fitgaml(ng,mg,eg,hi,himax,mmax,level)
  67         dimension mg(*),eg(*)
  68         if(level.eq.-1) then
  69                 call fitgam0(ng,mg,eg,hi,himax,mmax)
  70         elseif(level.eq.0) then
  71                 call fitgam1(ng,mg,eg,hi,himax,mmax)
  72         elseif(level.eq.1) then
  73 c               call fitgam2(ng,mg,eg,hi,himax,mmax)
  74                 write (*,*) ' next time '
  75         else
  76                 write (*,*) ' oshibsja parenek? '
  77         endif
  78         return
  79         end
  80
  81         subroutine fitgam0(ng,mg,eg,hi,himax,mmax)
  82         common /error/ier
  83         include 'comalwl.for'
  84         include 'comgam.for'
  85         include 'comwlgen.for'
  86         dimension mg(*),eg(*)
  87
  88 c simplgam
  89         call filet
  90         call filwt
  91         call diven
  92         call ecgam
  93         call clrwt
  94
  95         call filet
  96         call filwt
  97         call diven
  98         call ecgam
  99         call clrwt
 100
 101         call filet
 102         call filwt
 103         call hiwlnew(ng,mg,himin,himax,prb,mmax)
 104         call clrwt
 105
 106         return
 107         end
 108
 109         subroutine fitgam1(ng,mg,eg,prb,himax,mmax)
 110         common /error/ier
 111         include 'comalwl.for'
 112         include 'comgam.for'
 113         include 'comwlgen.for'
 114         dimension mg(*),eg(*)
 115
 116 c       call fitgam(ng,mg,eg)
 117         call again_gmfit
 118         call filet
 119         call filwt
 120         call hiwlnew(ng,mg,himin,himax,prb,mmax)
 121         call clrwt
 122         return
 123         end
 124
 125         subroutine gmplus(ng,mg,eg,gmporog,egamporog,rnenear)
 126         common /error/ier
 127
 128         include 'comgl.for'
 129         include 'comgam.for'
 130         include 'comalwl.for'
 131         include 'comwlgen.for'
 132         common /comwldif/ swldiff(40000),eadd(2000)
 133
 134         dimension mg(*),eg(*)
 135         data dmax /27./
 136 c       data gmporog/0.1/
 137         integer ijdd(2)
 138
 139 c
 140         call filwt
 141         if(ng.gt.2000) then
 142                 write (*,*) ' incr.in gmplus '
 143                 return
 144         endif
 145
 146         do n=1,ng
 147         m=mg(n)
 148         sig=sqrt(swlt(m)+swlt0(m))
 149                 efar=0.
 150                 do 1 kk=1,kg
 151                 dx=abs(x(1,mw(kk))-x(1,m))-(d(1,mw(kk))+d(1,m))/2.
 152                 dy=abs(x(2,mw(kk))-x(2,m))-(d(2,mw(kk))+d(2,m))/2.
 153                 if(dx.le.1..and.dy.le.1.) goto 1
 154                 if(dx.gt.dmax.or.dy.gt.dmax) goto 1
 155                 xd=xw(kk)-x(1,m)
 156                 yd=yw(kk)-x(2,m)
 157                 adefar=e(kk)*aclnew(xd,yd,d(1,m),m)
 158                 efar=efar+adefar
 159 1               continue
 160         sefar=0.15*sqrt(efar)
 161         swldiff(m)=wl(m)-wlt(m)-efar-sqrt(swlt(m)+swlt0(m))-sefar
 162         if(swldiff(m).le.0.) swldiff(m)=0.
 163         enddo
 164         call clrwt
 165 c
 166 c  LET'S TRY TO FINDE NEW GAMMA
 167
 168         do n=1,ng
 169         em=eg(n)
 170         m=mg(n)
 171         adsum=swldiff(m)
 172         do k=1,kg
 173         if(m.eq.mw(k)) goto 2
 174         rrk=sqrt((x(1,m)-xw(k))**2+(x(2,m)-yw(k))**2)
 175         if(rrk.lt.rnenear) goto 2
 176         enddo
 177
 178         if(em.lt.gmporog) goto 2
 179         esumgg1=em
 180         esumgg2=em
 181         do is=2,ns(m)
 182         ms=ks(is,m)
 183         esumgg1=esumgg1+swldiff(ms)
 184         esumgg2=esumgg2+wl(ms)
 185         enddo
 186         if(esumgg1.lt.egamporog) goto 2
 187
 188                 do  3 n1=1,ng
 189                 if(n.ne.n1) then
 190                 m1=mg(n1)
 191                 em1=swldiff(m1)
 192                 dx=abs(x(1,m1)-x(1,m))-(d(1,m)+d(1,m1))/2.
 193                 dy=abs(x(2,m1)-x(2,m))-(d(2,m)+d(2,m1))/2.
 194                 if(dx.le.1..and.dy.le.1.) then
 195                         if(em.gt.em1) then
 196                                 adsum=adsum+em1
 197                         endif
 198                         goto 3
 199                 endif
 200                 if(dx.gt.dmax.or.dy.gt.dmax) goto 3
 201                 xd=dx+d(1,m)/2.
 202                 yd=dy+d(2,m)/2.
 203                 if(xd.le.0.) xd=0.
 204                 if(yd.le.0.) yd=0.
 205                 emax=1.5*em*aclnew(xd,yd,d(1,m),m)
 206                 if(emax.gt.em1) then
 207                         adsum=adsum+em1
 208                 endif
 209                 endif
 210 3               continue
 211 2       eadd(n)=adsum
 212         enddo
 213
 214         nmmm=1
 215         summax=eadd(nmmm)
 216         mmax=mg(nmmm)
 217         emax=eg(nmmm)
 218
 219         do n=2,ng
 220         if(eadd(n).gt.summax) then
 221                 nmmm=n
 222                 summax=eadd(nmmm)
 223                 mmax=mg(nmmm)
 224                 emax=eg(nmmm)
 225         endif
 226         if(eadd(n).eq.summax) then
 227         if(eg(n).gt.emax) then
 228                 nmmm=n
 229                 summax=eadd(nmmm)
 230                 mmax=mg(nmmm)
 231                 emax=eg(nmmm)
 232         endif
 233         endif
 234         enddo
 235
 236 c       write (*,*) ' wl '
 237 c       call prwall(wl)
 238 c       write (*,*) ' wlt '
 239 c       call prwall(wlt)
 240 c       write (*,*) ' swldiff '
 241 c       call prwall(swldiff)
 242
 243         if(emax.gt.gmporog.and.summax.gt.0.) then
 244 c               kg=kg+1
 245 c               e(kg)=emax+summax-swldiff(mmax)
 246 c               mw(kg)=mmax
 247                 ee=emax+summax-swldiff(mmax)
 248                 xx=x(1,mmax)
 249                 yy=x(2,mmax)
 250                 see=ee/4.
 251                 sxx=d(1,mmax)
 252                 syy=d(2,mmax)
 253                 ijdd(1)=1
 254                 ijdd(2)=1
 255         call adgmnew(mmax,ee,xx,yy,see,sxx,syy,ijdd)
 256
 257 c               write (*,*) ' NEW GAMMA e,m',kg,e(kg),mw(kg)
 258 c               write(*,*) ' swldiff kg,mg',kg,mw(1)
 259 c               call prwlm(swldiff,mw(1),5,100.)
 260         else
 261 c               write (*,*) ' NO NEW GAMMA ',kg
 262 c               write(*,*) ' swldiff kg,mg',kg,mg(1)
 263 c               call prwlm(swldiff,mw(1),5,100.)
 264         endif
 265
 266         call vzero(swldiff,nt)
 267
 268         return
 269         end
 270
 271         subroutine filwltgl(ngg,mgg)
 272         include 'comalwl.for'
 273         include 'comgam.for'
 274         include 'comwlgen.for'
 275         common /comdevpar/sigmaph,sigmapd,sigphsq,sigpdsq
 276
 277         integer mgg(*)
 278
 279         if(ngg.eq.0) return
 280
 281         do k=1,kg
 282         do n=1,ngg
 283         mg=mgg(n)
 284
 285         if(mg.gt.0) then
 286         etadd=e(k)*ampcelnew(e(k),raxay(1,k),x(1,mg),d(1,mg),mg)
 287         wlt(mg)=wlt(mg) +etadd
 288 c       swlt(ms)=swlt(ms)+(emimx(2,n,k)-emimx(1,n,k))**2/4.
 289         swlt(mg)=swlt(mg)+
 290      +  (sigwlgam0(etadd,e(k)))**2+sigmphsq*etadd
 291         endif
 292
 293         enddo
 294         enddo
 295
 296         return
 297         end
 298
 299         subroutine clrwltgl(ngg,mgg)
 300         include 'comalwl.for'
 301         include 'comgam.for'
 302         include 'comwlgen.for'
 303
 304         integer mgg(*)
 305
 306         if(ngg.eq.0) return
 307
 308         do n=1,ngg
 309         mg=mgg(n)
 310         if(mg.gt.0) then
 311         wlt(mg)=0.
 312         swlt(mg)=0.
 313         endif
 314         enddo
 315
 316         return
 317         end
 318
 319         subroutine again_gmfit
 320         common /error/ier
 321         include 'comwlgen.for'
 322         include 'comalwl.for'
 323         include 'comgam.for'
 324         common /par/npar,e0,x0,y0,se0,sx0,sy0,Hi,rrr0(5)
 325         real ws(100),pinc(4)
 326
 327
 328 c       call deriv_test
 329
 330         call filet
 331         call filwt
 332         call filwsnew
 333         do k=1,kg
 334         m=mw(k)
 335         e0=e(k)
 336
 337 c       x0=x(1,m)+id(k)*d(1,k)/4.
 338 c       y0=x(2,m)+jd(k)*d(2,k)/4.
 339 c       write (*,*) ' k,xm,ym,id,jd ',x(1,m),x(2,m),id(k),jd(k)
 340
 341         x0=xw(k)
 342         y0=yw(k)
 343
 344 c only for anglx and angly
 345         call ucopy(raxay(1,k),rrr0(1),5)
 346         call ucopy(wsw(1,k),ws,ns(m))
 347
 348 c       write (*,*) ' es k',k,(es(ii,k),ii=1,ns(m))
 349 c       write (*,*) ' ws k',k,(sqrt(1./wsw(ii,k)),ii=1,ns(m))
 350
 351 c       write (*,*) ' e0,x0,y0 before ',e0,x0,y0
 352
 353         se0=0.
 354         sx0=0.
 355         sy0=0.
 356
 357         call gmfit(m,es(1,k),ws)
 358
 359 c       write (*,*) ' e0,x0,y0 after ',e0,x0,y0
 360 c       write (*,*) ' se0,sx0,sy0    ',se0,sx0,sy0
 361
 362
 363         e(k)=e0
 364         xw(k)=x0
 365         yw(k)=y0
 366         raxay(1,k)=x0
 367         raxay(2,k)=y0
 368
 369         sigexy(1,k)=se0
 370         sigexy(2,k)=sx0
 371         sigexy(3,k)=sy0
 372
 373
 374         enddo
 375         call clrwt
 376
 377         return
 378         end
 379
 380         subroutine gmfit(m,eg,weg)
 381 cold    external acl,dxacl
 382 cnew
 383         external ampcelnew,dampcelnew
 384         include 'comalwl.for'
 385         common /error/ier
 386         common /par/npar,e0,x0,y0,se0,sx0,sy0,Hi,rrr0(5)
 387         common /comkey/keykey
 388         real eg(25),weg(25)
 389 cnew
 390 c       real rrr0(5)
 391
 392         real f(4),ff(4,4),de(5,5)
 393         equivalence (aij,f(1)),(dxij,f(2)),(dyij,f(3)),(deij,f(4))
 394         equivalence (ade,ff(1,4)),(dxde,ff(2,4)),(dyde,ff(3,4))
 395         equivalence (aq,ff(1,1)),(dxq,ff(2,2)),(dyq,ff(3,3))
 396         equivalence (adx,ff(1,2)),(ady,ff(1,3)),(dxdy,ff(2,3))
 397         iefc=0
 398         ixfc=0
 399         iyfc=0
 400         nccmx=100
 401 c       write (*,*) ' se0,sx0,sy0 ',se0,sx0,sy0
 402
 403         if(se0.gt.0.) then
 404                 iefc=1
 405                 ec=e0
 406                 wec=1./se0**2
 407         endif
 408         if(sx0.gt.0.) then
 409                 ixfc=1
 410                 xc=x0
 411                 wxc=1./sx0**2
 412         endif
 413         if(sy0.gt.0.) then
 414                 iyfc=1
 415                 yc=y0
 416                 wyc=1./sy0**2
 417         endif
 418         hiold=10000000.
 419         do 10 ncc=1,nccmx
 420
 421         call vzero(f,4)
 422         call vzero(ff,16)
 423         hi=0.
 424         ssww=0.
 425         DO N=1,NS(M)
 426         MS=KS(N,M)
 427         XS=X(1,MS)
 428         YS=x(2,MS)
 429 cold        aij=ACL(X0-XS,Y0-YS,D(1,ms))
 430 cnew
 431         rrr0(1)=x0
 432         rrr0(2)=y0
 433 c       rrr0(3)=z0
 434 c       write (*,*) ' ax0,ay0 ',rrr0(4),rrr0(5)
 435 c       rrr0(4)=ax0
 436 c       rrr0(5)=ay0
 437         aij=ampcelnew(e0,rrr0,x(1,ms),d(1,ms),ms)
 438 cnewend
 439         deij=eg(n)-e0*aij
 440 cold    dxij=dxacl(x0-xs,y0-ys,d(1,ms))
 441 cold    dyij=dxacl(y0-ys,x0-xs,d(1,ms))
 442
 443 cnew
 444         dxij=dampcelnew(1,e0,rrr0,x(1,ms),d(1,ms),ms)
 445         dyij=dampcelnew(2,e0,rrr0,x(1,ms),d(1,ms),ms)
 446 cnewend
 447 c       if(keykey.eq.1) then
 448 c
 449 c       write (*,*) ' deij,dxij,dyij ',deij,dxij,dyij
 450 c
 451 c       endif
 452         do 2 k=1,4
 453         do 2 l=1,4
 454         ff(k,l)=ff(k,l)+f(k)*f(l)*weg(n)
 455 2       continue
 456         ssww=ssww+weg(n)
 457         ENDDO
 458         if(iefc.ne.0) then
 459                 ade=ade+wec*(ec-e0)
 460                 aq=aq+wec
 461         endif
 462         if(ixfc.ne.0) then
 463                 dxde=dxde+wxc*(xc-x0)/e0
 464                 dxq=dxq+wxc/e0**2
 465         endif
 466         if(iyfc.ne.0) then
 467                 dyde=dyde+wyc*(yc-y0)/e0
 468                 dyq=dyq+wyc/e0**2
 469         endif
 470         d11=dxq*dyq-dxdy**2
 471         d12=ady*dxdy-adx*dyq
 472         d13=adx*dxdy-ady*dxq
 473         d22=aq*dyq-ady**2
 474         d23=adx*ady-aq*dxdy
 475         d33=aq*dxq-adx**2
 476         det=aq*d11+adx*d12+ady*d13
 477         if(e0.le.0.) then
 478         write (*,*) ' e0=0 !!!!!!!!!!!!!!!!!!'
 479                 e0=0.
 480                 return
 481         endif
 482         if(det.eq.0.) then
 483         write (*,*) ' det=0 !!!!!!!!!!!!!!!!!!'
 484                 e0=0.
 485                 ier=1
 486                 return
 487         endif
 488         de0=(ade*d11+dxde*d12+dyde*d13)/det
 489         dx0=(ade*d12+dxde*d22+dyde*d23)/(e0*det)
 490         dy0=(ade*d13+dxde*d23+dyde*d33)/(e0*det)
 491
 492
 493         if(abs(de0).gt.3.) de0=3.*de0/abs(de0)
 494         if(abs(dx0).gt.d(1,ms)/3.) dx0=d(1,ms)/3.*dx0/abs(dx0)
 495         if(abs(dy0).gt.d(2,ms)/3.) dy0=d(2,ms)/3.*dy0/abs(dy0)
 496 c       write (*,*) ' de0,dx0,dy0 ',de0,dx0,dy0
 497         e0=e0+de0
 498         x0=x0+dx0
 499         y0=y0+dy0
 500 c e=0!!!!!
 501         if(e0.le.0.03) then
 502 c       write (*,*) ' return e0 !!!',e0
 503 c       write (*,*) ' main ',m,ns(m),x(1,m),x(2,m)
 504 c       write (*,*) ' eg ',(eg(ne),ne=1,ns(m))
 505 c       write (*,*) ' weg ',(sqrt(1./weg(ne)),ne=1,ns(m))
 506         e0fit=e0
 507         e0=0.
 508         DO Ne=1,NS(M)
 509         e0=e0+eg(ne)
 510         ENDDO
 511         if(e0.gt.1000.) write (*,*) ' return e0 from ',e0fit,' to',e0
 512 c       write (*,*) 'nc,e0,x0,y0',ncc,e0,x0,y0
 513         endif
 514         Hi=ff(4,4)
 515         if(iefc.ne.0) Hi=Hi+wec*(ec-e0)**2
 516         if(ixfc.ne.0) Hi=Hi+wxc*(xc-x0)**2
 517         if(iyfc.ne.0) Hi=Hi+wyc*(yc-y0)**2
 518         hi=sqrt(hi)
 519         delthi=(hi-hiold)/hi
 520 c       write (*,*) '       Hi delthi hiold',hi,delthi,hiold
 521         hiold=hi
 522         if(abs(delthi).lt.0.01) goto 11
 523 10      continue
 524 11      continue
 525 c       write (*,*) ' Hi delthi ',hi,delthi
 526         if(d11/det.gt.0.) se0=sqrt(d11/det)
 527         if(d22/det.gt.0.) sx0=sqrt(d22/det/e0)
 528         if(d33/det.gt.0.) sy0=sqrt(d33/det/e0)
 529         if(se0.eq.0) write (*,*) ' se0=0.'
 530         if(sx0.eq.0) write (*,*) ' sx0=0.'
 531         if(sy0.eq.0) write (*,*) ' sy0=0.'
 532 c       write (*,*) ' se0,sx0,sy0 ',se0,sx0,sy0
 533         return
 534         end
 535
 536         subroutine filwsnew
 537         include 'comalwl.for'
 538         include 'comgam.for'
 539         do k=1,kg
 540         m=mw(k)
 541 c       write (*,*) ' new gamma print '
 542         DO N=1,NS(M)
 543         MS=KS(N,M)
 544 c???    wsw(n,k)=swlt(ms)-(emimx(2,n,k)-emimx(1,n,k))**2/4.+swlt0(ms)
 545         wsw(n,k)=swlt(ms)+swlt0(ms)+swltmx(ms)
 546         wsw(n,k)=wsw(n,k)-(emimx(2,n,k)-emimx(1,n,k))**2/4.
 547         wsw(n,k)=1./wsw(n,k)
 548         ENDDO
 549         enddo
 550         return
 551         end
 552
 553
 554
 555
 556