PHOS/reconstruction/hisgamnew_v3.f

   1
   2         subroutine hisgamnew
   3
   4         include 'comgam.for'
   5
   6         nh0=110000
   7         if(ifl.eq.0) then
   8         ngmax=200
   9         egmax=5.
  10         etmax=150.
  11         ptmax=5.
  12         smax=1500.
  13         fngmax=float(ngmax)
  14 c               call hbook1(nh0+1,' ngam$',ngmax,0.,fngmax,0)
  15 c               call hbook1(nh0+2,' egam$',100,0.,egmax,0)
  16 c               call hbook1(nh0+3,' etot$',150,0.,etmax,0)
  17 c               call hbook1(nh0+4,'  pt  $',100,0.,ptmax,0)
  18 c       call hbook2(nh0+5,'  x*y  $',50,-smax,smax,50,-smax,smax,0)
  19 c               call hbook1(nh0+6,'  x  $',500,-smax,smax,0)
  20 c               call hbook1(nh0+7,'  y  $',500,-smax,smax,0)
  21 c               call hbook1(nh0+8,'  ptalice  $',100,0.,ptmax,0)
  22         ifl=1
  23         endif
  24
  25 c        call hf1(nh0+1,float(kg),1.)
  26         etot=0.
  27         do k=1,kg
  28              etot=etot+e(k)
  29         pt=sqrt(p4(1,k)**2+p4(2,k)**2)
  30         ptalice=sqrt(p4(1,k)**2+p4(3,k)**2)
  31 c             call hf1(nh0+2,e(k),1.)
  32 c             call hf2(nh0+5,xw(k),yw(k),1.)
  33 c             call hf1(nh0+6,xw(k),1.)
  34 c             call hf1(nh0+7,yw(k),1.)
  35 c             call hf1(nh0+4,pt,1.)
  36 c             call hf1(nh0+8,ptalice,1.)
  37          enddo
  38 c        call hf1(nh0+3,etot,1.)
  39
  40         return
  41         end
  42
  43         subroutine p4gamnew
  44         common /comgeom/igeomflag
  45
  46         include 'comgam.for'
  47         include 'comwlgen.for'
  48         external xgamvert,ygamvert,zgamvert
  49
  50         do k=1,kg
  51         iv=igamvert(k)
  52         xv=xgamvert(iv)
  53         yv=ygamvert(iv)
  54         zv=zgamvert(iv)
  55         igd=igdev(k)
  56         if(igeomflag.eq.2) then
  57                 rg=zgdev(k)+xyzwall(3,igd)
  58                 rr=4600.
  59                 fi=xw(k)/rr
  60                 xg=rg*sin(fi)-xv
  61                 zg=rg*cos(fi)-zv
  62                 yg=yw(k)-yv
  63         else
  64         xg=xw(k)-xv
  65         yg=yw(k)-yv
  66 c        write (*,*) ' k zgdev xyzwall zv',k,zgdev(k),xyzwall(3,igd),zv
  67         zg=zgdev(k)+xyzwall(3,igd)-zv
  68         endif
  69         rr=sqrt(xg**2+yg**2+zg**2)
  70         p4(1,k)=e(k)*xg/rr
  71         p4(2,k)=e(k)*yg/rr
  72         p4(3,k)=e(k)*zg/rr
  73         p4(4,k)=e(k)
  74         enddo
  75
  76         return
  77         end
  78
  79         real function xgamvert(ivert)
  80         include 'comgam.for'
  81         if(ivert.le.0) then
  82                 xgamvert=0.
  83         elseif(ivert.gt.ntotvertx) then
  84                 xgamvert=0.
  85                 write (*,*) ' wrong vertex number '
  86         else
  87                 xgamvert=gamvertex(1,ivert)
  88         endif
  89         return
  90         end
  91
  92         real function ygamvert(ivert)
  93         include 'comgam.for'
  94         if(ivert.le.0) then
  95                 ygamvert=0.
  96         elseif(ivert.gt.ntotvertx) then
  97                 ygamvert=0.
  98                 write (*,*) ' wrong vertex number '
  99         else
 100                 ygamvert=gamvertex(2,ivert)
 101         endif
 102         return
 103         end
 104
 105         real function zgamvert(ivert)
 106         include 'comgam.for'
 107         if(ivert.le.0) then
 108                 zgamvert=0.
 109         elseif(ivert.gt.ntotvertx) then
 110                 zgamvert=0.
 111                 write (*,*) ' wrong vertex number '
 112         else
 113                 zgamvert=gamvertex(3,ivert)
 114         endif
 115         return
 116         end
 117
 118         subroutine clrcompart
 119         include 'compart.for'
 120         nbpart=0
 121         return
 122         end
 123
 124         subroutine addp4tocompart(ncrad)
 125
 126         common /comgeom/igeomflag
 127         include 'event_format.inc'
 128         include 'comgam.for'
 129         include 'compart.for'
 130
 131         data ifl/0/
 132         real fir(10)
 133         data fir/10*0./
 134
 135         if(ifl.eq.0) then
 136         fi0=22.*44./4600.
 137         fi0=atan(fi0)
 138         nmmmax=crystals_matrix_amount_PHOS
 139         do i=1,nmmmax
 140 cnew
 141         fir(i)=fi0*float(nmmmax-1)-2.*fi0*float(i-1)
 142         enddo
 143         ifl=1
 144         endif
 145
 146         do k=1,kg
 147         if(nbpart+1.gt.nbpartmax) then
 148                 write (*,*) ' incr nbpartmax ',nbpartmax
 149                 return
 150         endif
 151         nbpart=nbpart+1
 152
 153         if(igeomflag.eq.1) then
 154                 call rotrefr(p4(1,k),ppart(1,nbpart),fir(ncrad))
 155         else
 156                 call rotrefc(p4(1,k),ppart(1,nbpart))
 157         endif
 158
 159         enddo
 160         return
 161         end
 162
 163         subroutine extrp4wrf(ncrad)
 164
 165         common /comgeom/igeomflag
 166         include 'event_format.inc'
 167         include 'comgam.for'
 168         include 'compart.for'
 169         include 'comggen.for'
 170
 171         data ifl/0/
 172         real fir(10)
 173         data fir/10*0./
 174         integer ijdd(2)
 175         real p4temp(4)
 176
 177
 178         if(ifl.eq.0) then
 179         fi0=22.*44./4600.
 180         fi0=atan(fi0)
 181         nmmmax=crystals_matrix_amount_PHOS
 182         do i=1,nmmmax
 183         fir(i)=fi0*float(nmmmax-1)-2.*fi0*float(i-1)
 184         enddo
 185         ifl=1
 186         endif
 187
 188 c        do n=1,nbpart
 189 c        if(igeomflag.eq.1) then
 190 c                call rotrefr_a(p4temp,ppart(1,n),fir(ncrad))
 191 c        else
 192 c                call rotrefc_a(p4temp,ppart(1,n))
 193 c        endif
 194 c
 195 c
 196 c        call xydetwlnew(xx,yy,p4temp)
 197 c       call mdetwlnew(mmww,xx,yy,ijdd,1)
 198 c
 199 c        enddo
 200
 201         do n=1,ngamgen
 202         if(igeomflag.eq.1) then
 203                 call rotrefr_a(p4temp,g(1,n),fir(ncrad))
 204         else
 205                 call rotrefc_a(p4temp,g(1,n))
 206         endif
 207
 208
 209         call xydetwlnew(xx,yy,p4temp)
 210         write (*,*) ' xx,yy ',xx,yy,xgamgen(n),ygamgen(n)
 211 c       call mdetwlnew(mmww,xx,yy,ijdd,1)
 212
 213        enddo
 214
 215         return
 216         end
 217
 218         subroutine calcgamgennew(ncrad)
 219
 220         common /comgeom/igeomflag
 221         include 'event_format.inc'
 222         include 'comgam.for'
 223         include 'compart.for'
 224         include 'comggen.for'
 225
 226         data ifl/0/
 227         real fir(10)
 228         data fir/10*0./
 229         integer ijdd(2)
 230         real p4temp(4)
 231
 232
 233         if(ifl.eq.0) then
 234         fi0=22.*44./4600.
 235         fi0=atan(fi0)
 236         nmmmax=crystals_matrix_amount_PHOS
 237         do i=1,nmmmax
 238         fir(i)=fi0*float(nmmmax-1)-2.*fi0*float(i-1)
 239         enddo
 240         ifl=1
 241         endif
 242
 243         do n=1,ngamgen
 244         call ucopy(g(1,n),p4ggen(1,n),3)
 245         p4ggen(4,n)=sqrt(g(1,n)**2+g(2,n)**2+g(3,n)**2)
 246         egamgen(n)=p4ggen(4,n)
 247         if(igeomflag.eq.1) then
 248                 call rotrefr_a(p4temp,p4ggen(1,n),fir(ncrad))
 249         else
 250                 call rotrefc_a(p4temp,p4ggen(1,n))
 251         endif
 252
 253
 254         call xydetwlnew(xx,yy,p4temp)
 255         xgamgen(n)=xx
 256         ygamgen(n)=yy
 257
 258        enddo
 259
 260         do n=1,ngamgen
 261         rmingen(n)=10000.
 262         kmingen(n)=n
 263         enddo
 264
 265         do n1=1,ngamgen-1
 266         do n2=n1+1,ngamgen
 267         rx=xgamgen(n1)-xgamgen(n2)
 268         ry=ygamgen(n1)-ygamgen(n2)
 269         rr=sqrt(rx**2+ry**2)
 270         if(rr.lt.rmingen(n1)) then
 271                 rmingen(n1)=rr
 272                 kmingen(n1)=n2
 273         endif
 274         if(rr.lt.rmingen(n2)) then
 275                 rmingen(n2)=rr
 276                 kmingen(n2)=n1
 277         endif
 278
 279         enddo
 280         enddo
 281
 282         return
 283         end
 284
 285         subroutine xydetwlnew(x,y,p4)
 286         common /comgeom/igeomflag
 287         real x,y,p4(4)
 288
 289 c       zc=4600.+46.
 290         zc=4600.+40.
 291         if(igeomflag.eq.2) then
 292 c cylindric case
 293                 rcentr=4600.
 294                 fi1=atan(p4(1)/p4(3))
 295                 g12=sqrt(p4(3)**2+p4(1)**2)
 296                 x=fi1*rcentr
 297                 y=p4(2)/g12*zc
 298         elseif(igeomflag.eq.1) then
 299 c rect case
 300                 x=p4(1)/p4(3)*zc
 301                 y=p4(2)/p4(3)*zc
 302         else
 303                 write (*,*) ' check geometry flag '
 304                 stop
 305         endif
 306         return
 307         end
 308
 309
 310
 311         subroutine rotrefr(p1,p2,fi)
 312         real p1(4),p2(4)
 313
 314         cosfi=cos(fi)
 315         sinfi=sin(fi)
 316         xp=p1(3)
 317         yp=p1(1)
 318         p2(1)=xp*cosfi+yp*sinfi
 319         p2(2)=yp*cosfi-xp*sinfi
 320         p2(3)=p1(2)
 321         p2(4)=p1(4)
 322         return
 323         end
 324
 325         subroutine rotrefc(p1,p2)
 326         real p1(4),p2(4)
 327         p2(1)=p1(3)
 328         p2(2)=p1(1)
 329         p2(3)=p1(2)
 330         p2(4)=p1(4)
 331         return
 332         end
 333
 334         subroutine rotrefr_a(p1,p2,fi)
 335         real p1(4),p2(4)
 336
 337         cosfi=cos(fi)
 338         sinfi=sin(fi)
 339
 340         p1(3)=p2(1)*cosfi-p2(2)*sinfi
 341         p1(1)=p2(2)*cosfi+p2(1)*sinfi
 342
 343         p1(2)=p2(3)
 344         p1(4)=p2(4)
 345         return
 346         end
 347
 348         subroutine rotrefc_a(p1,p2)
 349         real p1(4),p2(4)
 350         p1(3)=p2(1)
 351         p1(1)=p2(2)
 352         p1(2)=p2(3)
 353         p1(4)=p2(4)
 354         return
 355         end
 356
 357
 358