EPOS/epos167/epos-xan-154.f

   1 c---------------------------------------------------------------------
   2       subroutine xiniall
   3 c---------------------------------------------------------------------
   4       include 'epos.inc'
   5       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
   6       parameter (mypara=10,mxpara=10)
   7       logical ilog,icnx,itrevt,idmod
   8       double precision bin,bbin,zcbin,zbbin
   9       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
  10      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
  11      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
  12      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
  13      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
  14      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
  15      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
  16      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
  17      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
  18      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
  19      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
  20       double precision ebin,zebin
  21       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
  22      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
  23
  24       parameter (mxfra=5)
  25       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
  26      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
  27      $     ,emax(mxfra)
  28       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
  29      &     ,multc1
  30       parameter(mxxhis=50)
  31       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
  32      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis)
  33
  34       nhis=0
  35       nfra=0
  36       imulty1=0
  37       do n=0,mxxhis
  38         icorrtrig(n)=0
  39         ihardevent(n)=0
  40         ijetfind1(n)=0
  41         ijetfind2(n)=0
  42       enddo
  43       do n=1,mxhis
  44         xpara(1,n)=0
  45       enddo
  46       ncontrall=0
  47       noerrall=0
  48
  49       end
  50
  51 c---------------------------------------------------------------------
  52       subroutine xini
  53 c---------------------------------------------------------------------
  54 c  called after beginhisto
  55 c---------------------------------------------------------------------
  56       include 'epos.inc'
  57       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
  58       parameter (mypara=10,mxpara=10)
  59       logical ilog,icnx,itrevt,idmod
  60       double precision bin,bbin,zcbin,zbbin
  61       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
  62      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
  63      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
  64      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
  65      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
  66      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
  67      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
  68      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
  69      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
  70      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
  71      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
  72       double precision ebin,zebin
  73       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
  74      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
  75       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
  76      &     ,multc1
  77
  78       parameter (mxfra=5)
  79       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
  80      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
  81      $     ,emax(mxfra)
  82       character line*160,cvar*6
  83       logical go
  84       common/nl/noplin
  85       character*160 cline
  86       common/cjjj/jjj,cline
  87
  88       call utpri('xini  ',ish,ishini,5)
  89
  90       i=1
  91                                 !      iapl=0
  92                                 !      nhis=0
  93       j=jjj     !-1
  94       line=cline
  95                                 !      nfra=1
  96                                 !      ifra(1)=iframe
  97       iapl=0
  98       if(nfra.eq.0)then
  99         nfra=1
 100         ifra(1)=iframe
 101       endif
 102       nhis=nhis+1
 103       if(nhis.gt.mxhis)stop'xini: mxhis too small.       '
 104       noweak(nhis)=0
 105       ionoerr=0
 106 c      newfra=0
 107       indfra=1
 108 c      nepfra=0
 109       inpfra=1
 110  1    call utword(line,i,j,0)
 111       if(line(i:j).eq.'application')then !-----------
 112         call utword(line,i,j,1)
 113         if(line(i:j).eq.'analysis')then
 114                                 !iapl=0
 115                                 !nhis=nhis+1
 116                                 !newfra=0
 117                                 !indfra=1
 118                                 !nepfra=0
 119                                 !inpfra=1
 120         else
 121           iapl=1
 122         endif
 123       elseif(line(i:j).eq.'input')then !-----------
 124         call utword(line,i,j,0)
 125         if(nopen.ge.0)then
 126          nopen=nopen+1
 127          if(nopen.gt.9)stop'too many nested input commands'
 128          open(unit=20+nopen,file=line(i:j),status='old')
 129          if(iprmpt.eq.1)iprmpt=-1
 130         endif
 131       elseif(line(i:j).eq.'runprogram')then !-----------
 132         if(iapl.eq.0)then
 133         else
 134           goto 9999
 135         endif
 136       elseif(line(i:j).eq.'frame'.or.line(i:j).eq.'frame+')then !------
 137         ifp=1
 138         if(line(i:j).eq.'frame+')ifp=2
 139         call utword(line,i,j,1)
 140         if(line(i:j).eq.'total')then
 141           nfp=iframe
 142         elseif(line(i:j).eq.'nucleon-nucleon')then
 143           nfp=11
 144         elseif(line(i:j).eq.'target')then
 145           nfp=12
 146         elseif(line(i:j).eq.'gamma-nucleon')then
 147           nfp=21
 148         elseif(line(i:j).eq.'lab')then
 149           nfp=22
 150         elseif(line(i:j).eq.'breit')then
 151           nfp=23
 152         elseif(line(i:j).eq.'thrust')then
 153           nfp=33
 154         elseif(line(i:j).eq.'sphericity')then
 155           nfp=32
 156         endif
 157         go=.true.
 158         do l=1,nfra
 159           if(ifra(l).eq.nfp)then
 160             inl=l
 161             go=.false.
 162           endif
 163         enddo
 164         if (go) then
 165           nfra=nfra+1
 166           inl=nfra
 167           ifra(nfra)=nfp
 168         endif
 169         if(ifp.eq.1)then
 170           indfra=inl
 171 c          newfra=nfp
 172           ivfra(1,nhis)=indfra
 173           ivfra(2,nhis)=indfra
 174         else
 175           inpfra=inl
 176 c          nepfra=nfp
 177         endif
 178       elseif(line(i:j).eq.'binning')then !-----------
 179         call utword(line,i,j,1)
 180         if(line(i:j).eq.'lin')then
 181           iologb=0
 182           iocnxb=0
 183         elseif(line(i:j).eq.'log')then
 184           iologb=1
 185           iocnxb=0
 186         elseif(line(i:j).eq.'clin')then
 187           iologb=0
 188           iocnxb=1
 189         elseif(line(i:j).eq.'clog')then
 190           iologb=1
 191           iocnxb=1
 192         else
 193           print *, 'what the heck is ',line(i:j),' binning?'
 194           print *, 'I will use the linear (lin) one'
 195         endif
 196       elseif(line(i:j).eq.'setm')then !-----------
 197         if(iapl.eq.0) then
 198           print *,"You should use histogram instead of setm, please"
 199           stop
 200         endif
 201       elseif(line(i:j).eq.'set')then !-----------
 202         call utword(line,i,j,1)
 203         if(line(i:j).eq.'iologb')then
 204           call utword(line,i,j,1)
 205           read(line(i:j),*) iologb
 206         elseif(line(i:j).eq.'iocnxb')then
 207           call utword(line,i,j,1)
 208           read(line(i:j),*) iocnxb
 209         elseif(line(i:j).eq.'etacut')then
 210           call utword(line,i,j,1)
 211           read(line(i:j),*) etacut
 212         elseif(line(i:j).eq.'nemsi')then
 213           call utword(line,i,j,1)
 214           read(line(i:j),*)nemsi
 215         endif
 216       elseif(line(i:j).eq.'xpara')then !-----------
 217         call utword(line,i,j,1)
 218         read(line(i:j),*)ipara
 219         if(ipara.gt.mxpara)stop'mxpara too small.         '
 220         call utword(line,i,j,1)
 221         read(line(i:j),*)val
 222         xpara(ipara,nhis)=val
 223       elseif(line(i:j).eq.'echo')then !-----------
 224         call utword(line,i,j,1)
 225         if(line(i:j).eq.'on')iecho=1
 226         if(line(i:j).eq.'off')iecho=0
 227         if(line(i:j).ne.'on'.and.line(i:j).ne.'off')
 228      *  stop'invalid option'
 229       elseif(line(i:j).eq.'noweak')then !-----------
 230         noweak(nhis)=1
 231       elseif(line(i:j).eq.'histogram')then !-----------
 232         nac(nhis)=1
 233         call utword(line,i,j,1) !xvaria
 234         cvar='      '
 235         cvar=line(i:j)
 236         call xtrans(cvar,inom,ifrnew,nhis)
 237         if(inom.eq.-1)then
 238           if(line(i:i).ge.'0'.and.line(i:i).le.'9')then
 239             inom=298
 240             read(line(i:j),*) sval(1,nhis)
 241           endif
 242         endif
 243         ivar(1,nhis)=inom
 244         if(ifrnew.ne.0)then     !check frame for e+e- event
 245           go=.true.             !shape variables
 246           do l=1,nfra
 247             if(ifra(l).eq.ifrnew)then
 248               indfra=l
 249               go=.false.        !have it already
 250             endif
 251           enddo
 252           if (go) then
 253             nfra=nfra+1
 254             ifra(nfra)=ifrnew
 255             indfra=nfra
 256           endif
 257         endif
 258         call utword(line,i,j,1) !yvaria
 259         cvar='      '
 260         cvar=line(i:j)
 261         call xtrans(cvar,inom,ifrnew,nhis)
 262         ivar(2,nhis)=inom
 263         if(inom.eq.-1)then
 264           if(line(i:i).ge.'0'.and.line(i:i).le.'9')then
 265             inom=299
 266             read(line(i:j),*) sval(2,nhis)
 267           endif
 268         endif
 269         if(inom.eq.-1)ivar(1,nhis)=inom
 270
 271         ivfra(1,nhis)=indfra
 272         ivfra(2,nhis)=indfra
 273
 274         call utword(line,i,j,1) !normation
 275         read(line(i:j),*) inorm(nhis)
 276
 277         call utword(line,i,j,1) !xmin
 278         if(line(i:j).eq.'egy')then
 279          if(engy.gt.0)then
 280           egy=engy
 281          elseif(ecms.gt.0.)then
 282           egy=ecms
 283          elseif(elab.gt.0)then
 284           call idmass(idproj,apj)
 285           call idmass(idtarg,atg)
 286           egy=sqrt( 2*elab*atg+atg**2+apj**2 )
 287          elseif(ekin.gt.0.)then
 288           call idmass(idproj,apj)
 289           call idmass(idtarg,atg)
 290           egy=sqrt( 2*(ekin+apj)*atg+atg**2+apj**2 )
 291          elseif(pnll.gt.0.)then
 292           call idmass(idproj,apj)
 293           call idmass(idtarg,atg)
 294           egy=sqrt( 2*sqrt(pnll**2+apj**2)*atg+atg**2+apj**2 )
 295          else
 296           stop'pb in xini (1).   '
 297          endif
 298          xmin(nhis)=egy-0.001
 299         else
 300          read(line(i:j),*) xmin(nhis)
 301         endif
 302
 303         call utword(line,i,j,1) !xmax
 304         if(line(i:j).eq.'egy')then
 305          if(engy.gt.0)then
 306           egy=engy
 307          elseif(ecms.gt.0.)then
 308           egy=ecms
 309          elseif(elab.gt.0)then
 310           call idmass(idproj,apj)
 311           call idmass(idtarg,atg)
 312           egy=sqrt( 2*elab*atg+atg**2+apj**2 )
 313          elseif(ekin.gt.0.)then
 314           call idmass(idproj,apj)
 315           call idmass(idtarg,atg)
 316           egy=sqrt( 2*(ekin+apj)*atg+atg**2+apj**2 )
 317          elseif(pnll.gt.0.)then
 318           call idmass(idproj,apj)
 319           call idmass(idtarg,atg)
 320           egy=sqrt( 2*sqrt(pnll**2+apj**2)*atg+atg**2+apj**2 )
 321          else
 322           stop'pb in xini (2).   '
 323          endif
 324          xmax(nhis)=egy+0.001
 325         else
 326          read(line(i:j),*) xmax(nhis)
 327         endif
 328
 329         call utword(line,i,j,1) !nbin
 330         read(line(i:j),*) nbin(nhis)
 331         do l=1,nbin(nhis)
 332           bin(l,nac(nhis),nhis)=0.
 333           zcbin(l,nac(nhis),nhis)=0
 334         enddo
 335         lookcontr(nhis)=0
 336         lookcontrx(nhis)=0
 337         inoerr(nhis)=0
 338       elseif(line(i:j).eq.'idcode')then !-----------
 339         call utword(line,i,j,1) !idcode
 340         if(line(i:i+2).eq.'995')stop'xini: idcode 995 not supported'
 341         if(line(i:i+2).eq.'994')stop'xini: idcode 994 not supported'
 342         nidcod(nhis)=nidcod(nhis)+1
 343         read(line(i:j),*) idcod(nidcod(nhis),nhis)
 344         idmod(nidcod(nhis),nhis)=.false.
 345       elseif(line(i:j).eq.'idcode+')then !-----------
 346         stop'xini: idcode+ not supported'
 347         call utword(line,i,j,1) !idcode
 348         if(line(i:i+2).eq.'995')stop'xini: idcode 995 not supported'
 349         if(line(i:i+2).eq.'994')stop'xini: idcode 994 not supported'
 350         nidcod(nhis)=nidcod(nhis)+1
 351         read(line(i:j),*) idcod(nidcod(nhis),nhis)
 352         idmod(nidcod(nhis),nhis)=.true.
 353       elseif(line(i:j).eq.'trigger')then !-----------
 354         call utword(line,i,j,1)
 355         ntc=1
 356         imo=1
 357         ncontr=0
 358         icontrtyp(nhis)=0
 359         if(line(i:j).eq.'or'.or.line(i:j).eq.'contr')then
 360           imo=2
 361           if(line(i:j).eq.'contr')imo=3
 362           call utword(line,i,j,1)
 363           read(line(i:j),*)ztc
 364           ntc=nint(ztc)
 365           call utword(line,i,j,1)
 366           if(imo.eq.3)then
 367             ncontr=ntc
 368             ncontrall=ncontrall+ncontr
 369             if(ncontrall.gt.mxcontr)stop'xini: mxcontr too small.     '
 370             if(ncontr.gt.mxcnt)stop'xini: mxcnt too small.     '
 371             lookcontr(nhis)=ncontrall-ncontr+1
 372             lookcontrx(nhis)=ncontrall
 373             do nb=1,nbin(nhis)
 374               do nn=1,ncontr
 375                 bbin(nb,nac(nhis),lookcontr(nhis)-1+nn)=0.d0
 376                 zbbin(nb,nac(nhis),lookcontr(nhis)-1+nn)=0.d0
 377               enddo
 378             enddo
 379             do nn=1,ncontr
 380                     nccevt(lookcontr(nhis)-1+nn)=0
 381             enddo
 382           endif
 383         endif
 384         do n=1,ntc
 385           if(n.ne.1)call utword(line,i,j,1) !trigger-name
 386           cvar='      '
 387           ifp=1
 388           if(line(j:j).eq.'+')then
 389             cvar=line(i:j-1)
 390             ifp=2
 391           else
 392             cvar=line(i:j)
 393             ifp=1
 394           endif
 395           call xtrans(cvar,inom,ifrnew,nhis)
 396           if(inom.gt.0)then
 397             ntri(nhis)=ntri(nhis)+1
 398             if(ntc.eq.1)then
 399               ntrc(ntri(nhis),nhis)=1
 400             elseif(n.eq.1)then
 401               ntrc(ntri(nhis),nhis)=2
 402             elseif(n.eq.ntc)then
 403               ntrc(ntri(nhis),nhis)=3
 404             else
 405               ntrc(ntri(nhis),nhis)=0
 406             endif
 407             if(imo.eq.3)then
 408               ntrc(ntri(nhis),nhis)=-1
 409               if(n.eq.1)then
 410                 icontrtyp(nhis)=1+inom/100
 411               else
 412                 if(1+inom/100.ne.icontrtyp(nhis))
 413      *               stop'xini: type mismatch'
 414               endif
 415             endif
 416             itri(ntri(nhis),nhis)=inom
 417             if(ifp.eq.1)then
 418               itfra(ntri(nhis),nhis)=indfra
 419             else
 420               itfra(ntri(nhis),nhis)=inpfra
 421             endif
 422             xmitrp(ntri(nhis),nhis)=100.
 423             xmatrp(ntri(nhis),nhis)=100.
 424             call utword(line,i,j,1) !-----------xmin----------
 425             if(line(i:j).eq.'inf')then
 426               xmitri(ntri(nhis),nhis)=1e30
 427             elseif(line(i:j).eq.'-inf')then
 428               xmitri(ntri(nhis),nhis)=-1e30
 429             elseif(line(i:j).eq.'A')then
 430               xmitri(ntri(nhis),nhis)=maproj
 431             elseif(line(i:j).eq.'A+1')then
 432               xmitri(ntri(nhis),nhis)=maproj+1
 433             elseif(line(i:j).eq.'A+B')then
 434               xmitri(ntri(nhis),nhis)=maproj+matarg
 435             elseif(line(i:j).eq.'A+B+1')then
 436               xmitri(ntri(nhis),nhis)=maproj+matarg+1
 437             elseif(line(i:j).eq.'lead')then    !leading particle (neads Standard Variable)
 438               xmitri(ntri(nhis),nhis)=-123456
 439               imulty1=1
 440             else
 441               kk=0
 442               do k=i+1,j-1
 443                 if(line(k:k).eq.'%')kk=k
 444               enddo
 445               if(kk.eq.0)then
 446                 read(line(i:j),*)xmitri(ntri(nhis),nhis)
 447               else
 448                 read(line(i:kk-1),*)xmitrp(ntri(nhis),nhis)
 449                 read(line(kk+1:j),*)xmitri(ntri(nhis),nhis)
 450               endif
 451             endif
 452             call utword(line,i,j,1) !-----------xmax------------
 453             if(line(i:j).eq.'inf')then
 454               xmatri(ntri(nhis),nhis)=1e30
 455             elseif(line(i:j).eq.'-inf')then
 456               xmatri(ntri(nhis),nhis)=-1e30
 457             elseif(line(i:j).eq.'A')then
 458               xmatri(ntri(nhis),nhis)=maproj
 459             elseif(line(i:j).eq.'A+1')then
 460               xmatri(ntri(nhis),nhis)=maproj+1
 461             elseif(line(i:j).eq.'A+B')then
 462               xmatri(ntri(nhis),nhis)=maproj+matarg
 463             elseif(line(i:j).eq.'A+B+1')then
 464               xmatri(ntri(nhis),nhis)=maproj+matarg+1
 465             elseif(line(i:j).eq.'lead')then    !leading particle (neads Standard Variable)
 466               xmatri(ntri(nhis),nhis)=-123456
 467               imulty1=1
 468             else
 469               kk=0
 470               do k=i+1,j-1
 471                 if(line(k:k).eq.'%')kk=k
 472               enddo
 473               if(kk.eq.0)then
 474                 read(line(i:j),*)xmatri(ntri(nhis),nhis)
 475                 xmatrp(ntri(nhis),nhis)=100.
 476               else
 477                 read(line(i:kk-1),*)xmatrp(ntri(nhis),nhis)
 478                 read(line(kk+1:j),*)xmatri(ntri(nhis),nhis)
 479               endif
 480             endif
 481             !---exchange min-max------------------
 482             if(xmitri(ntri(nhis),nhis).gt.xmatri(ntri(nhis),nhis))then
 483               xmatri_save=xmatri(ntri(nhis),nhis)
 484               xmatrp_save=xmatrp(ntri(nhis),nhis)
 485               xmatri(ntri(nhis),nhis)=xmitri(ntri(nhis),nhis)
 486               xmatrp(ntri(nhis),nhis)=xmitrp(ntri(nhis),nhis)
 487               xmitri(ntri(nhis),nhis)=xmatri_save
 488               xmitrp(ntri(nhis),nhis)=xmatrp_save
 489             endif
 490             !-------------------------------------
 491           else
 492             ivar(1,nhis)=-1
 493             call utword(line,i,j,1) !xmin
 494             call utword(line,i,j,1) !xmax
 495           endif
 496         enddo
 497       elseif(line(i:j).eq.'noerrorbut')then !-----------
 498         ionoerr=ionoerr+1
 499         if(ionoerr.gt.2)stop'xini: not more than 2 noerrorbut !   '
 500         noerrall=noerrall+1
 501         if(noerrall.gt.mxhis/2)stop'xini: to many noerrorbut     '
 502
 503         call utword(line,i,j,1) !variable-name
 504         cvar=line(i:j)
 505         call xtrans(cvar,inom,ifrnew,nhis)
 506         if(inom.gt.0)then
 507           if(inom.gt.100)then
 508             write(*,*)'xini: noerrorbut can not be used with :',cvar
 509             stop'xini: error with noerrorbut!'
 510           endif
 511           noerrhis(nhis)=noerrall-ionoerr+1
 512           noerr(noerrhis(nhis),ionoerr)=inom
 513           do nb=1,nbin(nhis)
 514              ebin(nb,nac(nhis),ionoerr-1+noerrhis(nhis))=0.d0
 515                   zebin(nb,nac(nhis),ionoerr-1+noerrhis(nhis))=0.d0
 516           enddo
 517         else
 518           ionoerr=ionoerr-1
 519           noerrall=noerrall-1
 520         endif
 521         inoerr(nhis)=ionoerr
 522       elseif(line(i:j).eq.'write')then !-----------
 523         call utword(line,i,j,1)
 524       elseif(line(i:j).eq.'writearray')then !-----------
 525         call utword(line,i,j,1)
 526         iologb=0
 527         iocnxb=0
 528       elseif(line(i:j).eq.'writehisto')then !-----------
 529         call utword(line,i,j,1)
 530         iologb=0
 531         iocnxb=0
 532       elseif(line(i:j).eq.'endhisto')then   !-----------
 533         ilog(nhis)=.false.
 534         icnx(nhis)=.false.
 535         if(iologb.eq.1)ilog(nhis)=.true.
 536         if(iocnxb.eq.1)icnx(nhis)=.true.
 537         if(ilog(nhis))then
 538           xinc(nhis)=1./log(xmax(nhis)/xmin(nhis))*nbin(nhis)
 539         else
 540           xinc(nhis)=float(nbin(nhis))/(xmax(nhis)-xmin(nhis))
 541         endif
 542         iologb=0
 543         iocnxb=0
 544         jjj=j
 545         cline=line
 546         goto 9999
 547       endif
 548       goto 1
 549
 550  9999 continue
 551       if(ish.ge.5)then
 552         do n=1,nhis
 553           write (ifch,*) ivar(1,n),ivar(2,n),'(',ivfra(1,n),ivfra(2,n)
 554      $         ,')',inorm(n)
 555      $         ,xmin(n),xmax(n),ilog(n),icnx(n)
 556      $         ,nbin(n),(idcod(j,n),j=1,nidcod(n))
 557      $         ,' tri:',ntri(n),(itri(j,n),j=1,ntri(n)),'('
 558      $         ,(itfra(j,n),j=1,ntri(n)),')'
 559      $         ,(xmitri(j,n),j=1,ntri(n)) ,(xmatri(j,n),j=1,ntri(n))
 560         enddo
 561         write (ifch,*) (ifra(j),j=1,nfra)
 562       endif
 563       call utprix('xini  ',ish,ishini,5)
 564       return
 565       end
 566
 567
 568 c---------------------------------------------------------------------
 569       subroutine xana
 570 c---------------------------------------------------------------------
 571       include 'epos.inc'
 572       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
 573       parameter (mypara=10,mxpara=10)
 574       logical ilog,icnx,itrevt,idmod
 575       double precision bin,bbin,zcbin,zbbin
 576       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
 577      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
 578      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
 579      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
 580      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
 581      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
 582      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
 583      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
 584      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
 585      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
 586      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
 587       double precision ebin,zebin
 588       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
 589      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
 590
 591       parameter (mxfra=5)
 592       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
 593      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
 594      $     ,emax(mxfra)
 595       double precision bofra
 596       common/dfra/bofra(5,mxfra)
 597       parameter (ntim=1000)
 598       common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
 599      &     ,idaprt(2,ntim)
 600
 601       double precision pgampr,rgampr
 602       common/cgampr/pgampr(5),rgampr(4)
 603
 604       dimension ten(4,3)
 605       logical go,goo(mxcnt)
 606
 607       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
 608      &     ,multc1
 609       parameter(mxxhis=50)
 610       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
 611      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis)
 612
 613       call utpri('xana  ',ish,ishini,4)
 614
 615       do n=1,nhis
 616       do i=1,mypara
 617         ypara(i,n)=0
 618       enddo
 619       enddo
 620
 621
 622       if(ish.ge.5)write(ifch,*)'frames ...'
 623
 624       if(iappl.eq.6)then
 625         if(mod(iolept/10,10).eq.1) call gakjet(1)
 626         if(mod(iolept/100,10).eq.1) call gakjet(2)
 627       endif
 628
 629       do l=1,nfra
 630         emax(l)=egyevt/2
 631         if(ifra(l).eq.12)emax(l)=sqrt(pnll**2+prom**2)
 632         if(ifra(l).eq.iframe)then
 633           imofra(1,l)=0
 634           imofra(2,l)=0
 635           imofra(3,l)=0
 636         elseif(ifra(l).eq.11.or.ifra(l).eq.12)then
 637           imofra(1,l)=0
 638           imofra(2,l)=0
 639           bofra(1,l)=0d0
 640           bofra(2,l)=0d0
 641           bofra(3,l)=dsinh(dble(yhaha))
 642           bofra(4,l)=dcosh(dble(yhaha))
 643           bofra(5,l)=1d0
 644           if(ifra(l).eq.11.and.iframe.eq.12)then
 645             imofra(3,l)=1       ! target -> NN
 646           elseif(ifra(l).eq.12.and.iframe.eq.11)then
 647             imofra(3,l)=-1      ! NN -> target
 648           else
 649             imofra(3,l)=0       ! not known
 650           endif
 651         elseif(ifra(l).eq.21)then
 652           if(iframe.ne.21)then
 653             print *, 'invalid frame request'
 654             print *, 'choose frame gamma-nucleon for event run'
 655             stop'bye bye'
 656           endif
 657         elseif(ifra(l).eq.22)then
 658           if(iframe.eq.21)then
 659             imofra(1,l)=-1      !'  trafo gN -> lab'
 660             r1fra(1,l)=rgampr(1)
 661             r1fra(2,l)=rgampr(2)
 662             r1fra(3,l)=rgampr(3)
 663             imofra(2,l)=0
 664             imofra(3,l)=-1
 665             bofra(1,l)=pgampr(1)
 666             bofra(2,l)=pgampr(2)
 667             bofra(3,l)=pgampr(3)
 668             bofra(4,l)=pgampr(4)
 669             bofra(5,l)=pgampr(5)
 670           elseif(iframe.eq.22)then
 671                                 ! nothing to do already gN-frame
 672           else
 673             print *, 'invalid frame request'
 674             print *, 'choose frame gamma-nucleon or lab for event run'
 675             stop'bye bye'
 676           endif
 677         elseif(ifra(l).eq.23)then
 678           if(iframe.eq.21)then
 679             imofra(1,l)=0       ! gN -> breit-frame
 680             r1fra(1,l)=rgampr(1)
 681             r1fra(2,l)=rgampr(2)
 682             r1fra(3,l)=rgampr(3)
 683             imofra(2,l)=0
 684             imofra(3,l)=1
 685             bofra(1,l)=0d0
 686             bofra(2,l)=0d0
 687             bofra(3,l)=rgampr(4)
 688             bofra(4,l)=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
 689             bofra(5,l)=sqrt( bofra(4,l)**2-rgampr(4)**2)
 690           elseif(iframe.eq.23)then
 691                                 ! nothing to do already breit-frame
 692           else
 693             print *, 'invalid frame request'
 694             print *, 'choose frame gamma-nucleon or lab for event run'
 695             stop'bye bye'
 696           endif
 697         elseif(ifra(l).eq.33.or.ifra(l).eq.36)then
 698           if(ifra(l).eq.33)then
 699             call gakthru(ten,2)
 700           else
 701             call gakthru(ten,3)
 702           endif
 703           if(ten(4,1).lt.0.)then
 704             imofra(1,l)=0
 705             imofra(2,l)=0
 706             imofra(3,l)=0
 707           else
 708             arox=ten(1,1)
 709             aroy=ten(2,1)
 710             aroz=ten(3,1)
 711             brox=ten(1,2)
 712             broy=ten(2,2)
 713             broz=ten(3,2)
 714             call utrota(1,arox,aroy,aroz,brox,broy,broz)
 715             imofra(1,l)=1
 716             r1fra(1,l)=arox
 717             r1fra(2,l)=aroy
 718             r1fra(3,l)=aroz
 719             imofra(2,l)=1
 720             r2fra(1,l)=brox
 721             r2fra(2,l)=broy
 722             r2fra(3,l)=broz
 723             imofra(3,l)=0       !no boost
 724           endif
 725           bofra(1,l)=dble(ten(4,1)) !usually this is for boosting
 726           bofra(2,l)=dble(ten(4,2)) !I abuse it to store the eigenvalues
 727           bofra(3,l)=dble(ten(4,3)) !
 728         elseif(ifra(l).eq.32.or.ifra(l).eq.34.or.ifra(l).eq.35)then
 729           if(ifra(l).eq.32)then
 730             call gaksphe(ten,2.,2)
 731           elseif(ifra(l).eq.34)then
 732             call gaksphe(ten,1.,2)
 733           else
 734             call gaksphe(ten,2.,3)
 735           endif
 736           if(ten(4,1).lt.0.)then
 737             imofra(1,l)=0
 738             imofra(2,l)=0
 739             imofra(3,l)=0
 740           else
 741             arox=ten(1,1)
 742             aroy=ten(2,1)
 743             aroz=ten(3,1)
 744             brox=ten(1,2)
 745             broy=ten(2,2)
 746             broz=ten(3,2)
 747             call utrota(1,arox,aroy,aroz,brox,broy,broz)
 748             imofra(1,l)=1
 749             r1fra(1,l)=arox
 750             r1fra(2,l)=aroy
 751             r1fra(3,l)=aroz
 752             imofra(2,l)=1
 753             r2fra(1,l)=brox
 754             r2fra(2,l)=broy
 755             r2fra(3,l)=broz
 756             imofra(3,l)=0
 757           endif
 758           bofra(1,l)=dble(ten(4,1))
 759           bofra(2,l)=dble(ten(4,2))
 760           bofra(3,l)=dble(ten(4,3))
 761         endif
 762       enddo
 763
 764       do n=1,nhis
 765         itrevt(n)=.false.
 766         if(ivar(1,n).ge.100.and.ivar(1,n).le.199) sval(1,n)=0.
 767         if(ivar(2,n).ge.100.and.ivar(2,n).le.199) sval(2,n)=0.
 768         if(ivar(1,n).gt.300.and.ivar(1,n).lt.400)then
 769           call xval(n,ivar(1,n),ivfra(1,n),0,x) !initializing of  variables
 770         endif
 771         if(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
 772           call xval(n,ivar(2,n),ivfra(2,n),0,y) !
 773         endif
 774         do j=1,ntri(n)
 775           valtri(j,n)=0.
 776         enddo
 777         do j=1,nbin(n)  !copy bins
 778           bin(j,3-nac(n),n)=bin(j,nac(n),n)
 779           zcbin(j,3-nac(n),n)=zcbin(j,nac(n),n)
 780         enddo
 781         if(lookcontr(n).gt.0)then
 782           do j=1,nbin(n)
 783             do loo=lookcontr(n),lookcontrx(n)
 784                     bbin(j,3-nac(n),loo)=bbin(j,nac(n),loo)
 785                     zbbin(j,3-nac(n),loo)=zbbin(j,nac(n),loo)
 786             enddo
 787           enddo
 788         endif
 789         if(inoerr(n).gt.0)then
 790           do j=1,nbin(n)
 791             do nn=1,inoerr(n)
 792               ebin(j,3-nac(n),nn-1+noerrhis(n))=ebin(j,nac(n),
 793      &                                          nn-1+noerrhis(n))
 794               zebin(j,3-nac(n),nn-1+noerrhis(n))=zebin(j,nac(n),
 795      &                                          nn-1+noerrhis(n))
 796             enddo
 797           enddo
 798         endif
 799       enddo
 800
 801       if(imulty1.eq.1)then
 802         if(ish.ge.5)write(ifch,*)'Calculate standard variables ...'
 803         call StandardVariables
 804       endif
 805       if(ish.ge.5)write(ifch,*)'Call corrtrig ...'
 806       do n=1,icorrtrig(0)
 807         call corrtrig(icorrtrig(n))
 808       enddo
 809       if(ish.ge.5)write(ifch,*)'Call hardevent ...'
 810       do n=1,ihardevent(0)
 811         call hardevent(ihardevent(n))
 812       enddo
 813       if(ish.ge.5)write(ifch,*)'Call jetfind ...'
 814       do n=1,ijetfind1(0)
 815         call jetfind(1,ijetfind1(n))
 816       enddo
 817       do n=1,ijetfind2(0)
 818         call jetfind(2,ijetfind2(n))
 819       enddo
 820
 821 c...........................loop nptl...................................
 822       if(ish.ge.5)write(ifch,*)'Loop nptl ...'
 823       do j=1,nptl
 824         if(iorptl(j).lt.0.or.istptl(j).gt.istmax)goto8
 825         if(ish.ge.5)write(ifch,*)'ptl :',j
 826         call idchrg(idptl(j),ch)
 827         do i=1,nfra
 828           iffra(i)=0            !flag if frame calculated or not
 829         enddo
 830         do n=1,nhis
 831           if(ivar(1,n).eq.-1.or.ivar(2,n).eq.-1)goto 9
 832
 833 c...........check ids
 834           go=nidcod(n).eq.0
 835           do i=1,nidcod(n)
 836             if(istptl(j).eq.0.and.idcod(i,n).eq.9990)then
 837               if(abs(idptl(j)).ge.100
 838      $         .and.abs(idptl(j)).lt.10000) go=.true.
 839             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9970)then
 840               if(abs(ch).gt.0.1.and.abs(idptl(j)).ge.100
 841      $         .and.abs(idptl(j)).lt.10000) go=.true.
 842             elseif(istptl(j).eq.0.and.idcod(i,n).eq.-9960)then
 843               if(ch.lt.-0.1.and.abs(idptl(j)).ge.100
 844      $             .and.abs(idptl(j)).lt.10000)go=.true.
 845             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9960)then
 846               if(ch.gt.0.1.and.abs(idptl(j)).ge.100
 847      $         .and.abs(idptl(j)).lt.10000)go=.true.
 848             elseif((istptl(j).le.1.or.istptl(j).ge.10)
 849      $            .and.idcod(i,n).eq.idptl(j))then
 850               go=.true.
 851             endif
 852           enddo
 853           if(ish.ge.10)write(ifch,*)j,' id,ist',idptl(j),istptl(j),go
 854
 855 c...........check weak decay
 856           if(go)then
 857             if(noweak(n).eq.1)then  !do not consider weak decay products
 858              if(iorptl(j).ne.0)then
 859               idora=abs( idptl(iorptl(j)) )
 860               if(  idora.eq.20   .or.idora.eq.2130
 861      &               .or.idora.eq.2230 .or.idora.eq.1130
 862      &         .or.idora.eq.2330 .or.idora.eq.1330
 863      &         .or.idora.eq.3331 )go=.false.
 864              ! print *, j,n, '   ', idptl(j),idora,go
 865              endif
 866             endif
 867           endif
 868
 869 c...........check triggers
 870           if(go)then
 871             if(ish.ge.7)write(ifch,*)'  check triggers in histogram ',n
 872             ncontr=0
 873             do i=1,ntri(n)
 874               if(ish.ge.7)write(ifch,*)'  trigger variable: ',itri(i,n)
 875               if(itri(i,n).lt.100)then
 876                 call xval(n,itri(i,n),itfra(i,n),j,x)
 877                 if(ntrc(i,n).ne.-1)then
 878                   call triggercondition(i,n,x,go)
 879                 else
 880                   ncontr=ncontr+1
 881                   goo(ncontr)=.true.
 882                   call triggercondition(i,n,x,goo(ncontr))
 883                   if((ivar(1,n).gt.100.and.ivar(1,n).lt.200)
 884      .                 .or.(ivar(2,n).gt.100.and.ivar(2,n).lt.200))then
 885                     print*,'!-----------------------------------------'
 886                     print*,'!  100-199 event variables can not be used'
 887                     print*,'! in connection with "trigger contr ..."  '
 888                     print*,'!-----------------------------------------'
 889                     stop'in xana (1).                      '
 890                   endif
 891                 endif
 892               elseif(itri(i,n).lt.200)then
 893                 if(ntrc(i,n).eq.-1)then
 894                     print*,'!-----------------------------------------'
 895                     print*,'!  100-199 event variables can not be used'
 896                     print*,'! in connection with "trigger contr ..."  '
 897                     print*,'!-----------------------------------------'
 898                     stop'in xana (2).                      '
 899                 endif
 900                 call xval(n,itri(i,n),itfra(i,n),j,x)
 901                 valtri(i,n)=valtri(i,n)+x
 902               endif
 903             enddo
 904           endif
 905
 906 c............fill histogram
 907           if(go)then
 908             if(ish.ge.7)write(ifch,*)'  fill histogram '
 909      &            ,n,ivar(1,n),ivar(2,n),ivfra(2,n)
 910             if(ivar(1,n).lt.100.or.ivar(2,n).lt.100)then
 911               if(ivar(2,n).lt.100)then
 912                 call xval(n,ivar(2,n),ivfra(2,n),j,y)
 913                 sval(2,n)=y
 914               endif
 915               if(ivar(1,n).lt.100)then
 916                 call xval(n,ivar(1,n),ivfra(1,n),j,x)
 917                 if(x.ge.xmin(n).and.x.le.xmax(n))then
 918                   norm3=mod(inorm(n)/100,10)
 919                   if(norm3.eq.1)then
 920                     y=y*x
 921                   elseif(norm3.eq.2.and.ivar(1,n).eq.63.and.x.ne.0.)then
 922                     y=y/(x+pptl(5,j))/2/pi
 923                   elseif(norm3.eq.2.and.ivar(1,n).ne.63.and.x.ne.0.)then
 924                     y=y/x/2/pi
 925                   elseif(norm3.eq.4.and.x.ne.0.)then
 926                     y=y/x**1.5
 927                   elseif(norm3.eq.5.and.x.ne.0.)then
 928                     y=y/x
 929                   elseif(norm3.eq.7.and.x.ne.0.)then
 930                     y=y/x/sqrt(x-pptl(5,j))
 931                   endif
 932                   if(icnx(n))then
 933                     call fillhistoconex(n,x,y,ivfra(2,n),j)   !for conex
 934                   else
 935                     if(ilog(n))then
 936                       nb=1+int(log(x/xmin(n))*xinc(n))
 937                     else
 938                       nb=1+int((x-xmin(n))*xinc(n))
 939                     endif
 940                     bin(nb,nac(n),n)=bin(nb,nac(n),n)+y
 941                     if(ncontr.gt.0)then  !ptl trigger contr
 942                       do nn=1,ncontr
 943                         if(goo(nn))then
 944                              bbin(nb,nac(n),lookcontr(n)-1+nn)=
 945      &                  bbin(nb,nac(n),lookcontr(n)-1+nn)+y
 946                              zbbin(nb,nac(n),lookcontr(n)-1+nn)=
 947      &                  zbbin(nb,nac(n),lookcontr(n)-1+nn)+1
 948                         endif
 949                       enddo
 950                     endif
 951                     if(inoerr(n).gt.0)then
 952                       do nn=1,inoerr(n)
 953                        call xval(n,noerr(noerrhis(n),nn),ivfra(2,n),j,y)
 954                         ebin(nb,nac(n),nn-1+noerrhis(n))=
 955      &                       ebin(nb,nac(n),nn-1+noerrhis(n))+y
 956                         zebin(nb,nac(n),nn-1+noerrhis(n))=
 957      &                       zebin(nb,nac(n),nn-1+noerrhis(n))+1
 958                       enddo
 959                   endif
 960                     zcbin(nb,nac(n),n)=zcbin(nb,nac(n),n)+1
 961                   endif
 962                   itrevt(n)=.true.
 963                 endif
 964               endif
 965             endif
 966             if(ivar(1,n).gt.100.and.ivar(1,n).lt.200)then
 967               call xval(n,ivar(1,n),ivfra(1,n),j,x)
 968               sval(1,n)=sval(1,n)+x
 969             endif
 970             if(ivar(2,n).gt.100.and.ivar(2,n).lt.200)then
 971               call xval(n,ivar(2,n),ivfra(2,n),j,y)
 972               sval(2,n)=sval(2,n)+y
 973             endif
 974             if(ivar(1,n).gt.300.and.ivar(1,n).lt.400)then
 975               call xval(n,ivar(1,n),ivfra(1,n),j,x)
 976             endif
 977             if(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
 978               call xval(n,ivar(2,n),ivfra(2,n),j,y)
 979             endif
 980             if(ish.ge.8)write (ifch,*) '   ---> histo n,x,y:',n,x,y
 981           endif
 982    9      continue
 983         enddo
 984   8     continue
 985       enddo
 986 c...........................end loop nptl...........................
 987
 988
 989       do n=1,nhis
 990       if(ivar(1,n).eq.-1.or.ivar(2,n).eq.-1)goto 99
 991
 992 c........check event triggers
 993
 994        go=.true.
 995         ncontr=0
 996         do i=1,ntri(n)
 997           if(itri(i,n).gt.100)then
 998             if(itri(i,n).lt.200)then
 999               x=valtri(i,n)
1000             else
1001               call xval(n,itri(i,n),itfra(i,n),0,x)
1002             endif
1003             if(ntrc(i,n).ne.-1)then
1004               call triggercondition(i,n,x,go)
1005             else
1006               ncontr=ncontr+1
1007               goo(ncontr)=.true.
1008               call triggercondition(i,n,x,goo(ncontr))
1009             endif
1010           endif
1011         enddo
1012
1013 c........event variables > 200
1014
1015         if(go)then
1016           if(ivar(1,n).gt.100)then
1017             if(ivar(1,n).gt.200.and.ivar(1,n).lt.300)then
1018               call xval(n,ivar(1,n),ivfra(1,n),0,x)
1019             elseif(ivar(1,n).gt.300.and.ivar(1,n).lt.400)then
1020               call xval(n,ivar(1,n),ivfra(1,n),nptl+1,x)
1021             elseif(ivar(1,n).gt.100.and.ivar(1,n).lt.200)then
1022               x=sval(1,n)
1023             else
1024               call xval(n,ivar(1,n),ivfra(1,n),0,x)
1025             endif
1026             if(ivar(2,n).gt.200.and.ivar(2,n).lt.300)then
1027               call xval(n,ivar(2,n),ivfra(2,n),0,y)
1028             elseif(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
1029               call xval(n,ivar(2,n),ivfra(2,n),nptl+1,y)
1030             elseif(ivar(2,n).gt.0.and.ivar(2,n).lt.200)then
1031               y=sval(2,n)
1032             else             !inom>500
1033               call xval(n,ivar(2,n),ivfra(2,n),0,y)
1034             endif
1035 c The following doesn't work for ivar(2,n)<100, since particle number is not defined !
1036 c            if(ivar(2,n).gt.200.and.ivar(2,n).lt.300)then
1037 c              call xval(n,ivar(2,n),ivfra(2,n),0,y)
1038 c            elseif(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
1039 c              call xval(n,ivar(2,n),ivfra(2,n),nptl+1,y)
1040 c            elseif(ivar(2,n).gt.100.and.ivar(2,n).lt.200)then
1041 c              y=sval(2,n)
1042 c            else
1043 c              call xval(n,ivar(2,n),ivfra(2,n),0,y)
1044 c            endif
1045             if(mod(inorm(n)/100,10).eq.1)y=y*x
1046             if(mod(inorm(n)/100,10).eq.2.and.x.ne.0.)y=y/x/2/pi
1047             if(mod(inorm(n)/100,10).eq.4.and.x.ne.0.)y=y/x**1.5
1048             if(mod(inorm(n)/100,10).eq.5.and.x.ne.0.)y=y/x
1049             sval(1,n)=x
1050             sval(2,n)=y
1051             if(ish.ge.6) then
1052               write (ifch,*) 'histo n,x,y:',n,x,y
1053             endif
1054             if(x.ge.xmin(n).and.x.le.xmax(n))then
1055               if(ilog(n))then
1056                 nb=1+int(log(x/xmin(n))*xinc(n))
1057               else
1058                 nb=1+int((x-xmin(n))*xinc(n))
1059               endif
1060               bin(nb,nac(n),n)=bin(nb,nac(n),n)+y
1061               if(ncontr.gt.0)then
1062                 do nn=1,ncontr
1063                   if(goo(nn))
1064      &                    bbin(nb,nac(n),lookcontr(n)-1+nn)=
1065      &             bbin(nb,nac(n),lookcontr(n)-1+nn)+y
1066                 enddo
1067               endif
1068               zcbin(nb,nac(n),n)=zcbin(nb,nac(n),n)+1
1069               itrevt(n)=.true.
1070             endif
1071           endif
1072         endif
1073
1074 c........particle variables
1075
1076         if(go)then
1077           if(ivar(1,n).le.100)then
1078             if(ncontr.gt.0)then  !event trigger contr
1079               do nb=1,nbin(n)
1080                 do nn=1,ncontr
1081                   if(goo(nn))
1082      &                     bbin(nb,nac(n),lookcontr(n)-1+nn)=
1083      &              bbin(nb,nac(n),lookcontr(n)-1+nn)
1084      &              +bin(nb,nac(n),n)-bin(nb,3-nac(n),n)
1085                 enddo
1086               enddo
1087             endif
1088           endif
1089         endif
1090
1091 c............event ok (increase ncevt) or not (take copy)
1092
1093         if(go)then
1094           ncevt(n)=ncevt(n)+1
1095           if(ncontr.gt.0)then
1096             do nn=1,ncontr
1097               loo=lookcontr(n)-1+nn
1098               if(goo(nn))
1099      &        nccevt(loo)=nccevt(loo)+1
1100             enddo
1101           endif
1102         else
1103           nac(n)=3-nac(n)
1104           itrevt(n)=.false.
1105         endif
1106
1107  99   continue
1108       enddo
1109
1110       call utprix('xana  ',ish,ishini,4)
1111       end
1112
1113 c--------------------------------------------------------------------
1114       subroutine triggercondition(i,n,x,go)
1115 c--------------------------------------------------------------------
1116 c ntrc is used to distinguish the different usage of trigger:
1117 c
1118 c    trigger var xmin xmax
1119 c             ntrc=1
1120 c    trigger or n var1 xmin1 xmax1 var2 xmin2 xmax2 ... varn xminn xmaxn
1121 c          1  ntrc=2
1122 c          2  ntrc=0
1123 c              ...
1124 c         n-1 ntrc=0
1125 c          n  ntrc=3
1126 c    trigger contr n var1 xmin1 xmax1 var2 xmin2 xmax2 ... varn xminn xmaxn
1127 c             ntrc=-1
1128 c--------------------------------------------------------------------
1129       include 'epos.inc'
1130       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
1131       parameter (mypara=10,mxpara=10)
1132       logical ilog,icnx,itrevt,idmod
1133       double precision bin,bbin,zcbin,zbbin
1134       common/crvar/idlead
1135       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
1136      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
1137      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
1138      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
1139      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
1140      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
1141      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
1142      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
1143      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
1144      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
1145      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
1146       double precision ebin,zebin
1147       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
1148      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
1149       logical go,gox,ok,goz
1150                 xmn=xmitri(i,n)
1151                 xmx=xmatri(i,n)
1152                 if(xmn.eq.-123456.and.xmx.eq.-123456)then  !for leading part
1153                   xmn=float(idlead)
1154                   xmx=float(idlead)
1155                 endif
1156                 pmn=xmitrp(i,n)
1157                 pmx=xmatrp(i,n)
1158                 if(abs(ntrc(i,n)).eq.1)then
1159                   goz=.true.
1160                   if(pmn.gt.99.999.and.pmx.gt.99.999)then
1161                     if(x.lt.xmn.or.x.gt.xmx)goz=.false.
1162                   else
1163                     if(x.lt.xmn-0.5.or.x.gt.xmx+0.5)goz=.false.
1164                     ok=rangen().le.xmitrp(i,n)/100.
1165                     if(.not.ok.and.x.lt.xmn+0.5)goz=.false.
1166                     ok=rangen().le.xmatrp(i,n)/100.
1167                     if(.not.ok.and.x.gt.xmx-0.5)goz=.false.
1168                   endif
1169                   if(.not.goz)go=.false.
1170                 else
1171                   if(ntrc(i,n).eq.2)gox=.false.
1172                   goz=.true.
1173                   if(pmn.gt.99.999.and.pmx.gt.99.999)then
1174                     if(x.lt.xmn.or.x.gt.xmx)goz=.false.
1175                   else
1176                     if(x.lt.xmn-0.5.or.x.gt.xmx+0.5)goz=.false.
1177                     ok=rangen().le.xmitrp(i,n)/100.
1178                     if(.not.ok.and.x.lt.xmn+0.5)goz=.false.
1179                     ok=rangen().le.xmatrp(i,n)/100.
1180                     if(.not.ok.and.x.gt.xmx-0.5)goz=.false.
1181                   endif
1182                   if(goz)gox=.true.
1183                   if(ntrc(i,n).eq.3.and..not.gox)go=.false.
1184                 endif
1185                 end
1186
1187 c-----------------------------------------------------------------------
1188       subroutine fillhistoconex(n,x,y,lf,j)   !for conex
1189 c-----------------------------------------------------------------------
1190       include 'epos.inc'
1191       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
1192       parameter (mypara=10,mxpara=10)
1193       logical ilog,icnx,itrevt,idmod
1194       double precision bin,bbin,zcbin,zbbin
1195       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
1196      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
1197      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
1198      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
1199      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
1200      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
1201      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
1202      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
1203      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
1204      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
1205      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
1206       double precision ebin,zebin
1207       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
1208      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
1209
1210            if(.not.(mod(inorm(n),10).ne.4
1211      &      .and.mod(inorm(n),10).ne.6
1212      &      .and.mod(inorm(n)/100,10).ne.3))return
1213                     if(ilog(n))then
1214                       c=(xmax(n)/xmin(n))**(1./real(nbin(n)))
1215                       nde=nint(1./log10(c))
1216                       nb=max(1,1+int(log10(x/xmin(n))*nde))
1217                       xmb=xmin(n)*c**(nb-0.5)
1218                       if(x.gt.xmb.and.nb.lt.nbin(n))then
1219                         if(x.gt.xmax(n))
1220      &                    write(ifmt,*)'xana max ?',x,xmax(n),nb
1221                         nbx=1
1222                         xmx=c*xmb
1223                       elseif(x.lt.xmb.and.nb.gt.1)then
1224                         if(x.lt.xmin(n))write(ifmt,*)'xana min ?',x,nb
1225                         nbx=-1
1226                         xmx=xmb/c
1227                       else
1228                         nbx=0
1229                         xmx=0.
1230                       endif
1231                     else
1232                       c=(xmax(n)-xmin(n))/real(nbin(n))
1233                       nb=max(1,1+int((x-xmin(n))/c))
1234                       xmb=xmin(n)+c*(nb-0.5)
1235                       if(x.gt.xmb)then
1236                         nbx=1
1237                         xmx=c+xmb
1238                       elseif(x.lt.xmb)then
1239                         nbx=-1
1240                         xmx=xmb-c
1241                       else
1242                         nbx=0
1243                         xmx=0.
1244                       endif
1245                     endif
1246                     xc=(x-xmx)/(xmb-xmx)
1247                     xc=max(0.,min(1.,xc))
1248                     bin(nb,nac(n),n)=bin(nb,nac(n),n)+xc*y
1249                     if(nbx.ne.0)bin(nb+nbx,nac(n),n)
1250      &                  =bin(nb+nbx,nac(n),n)+(1.-xc)*y
1251                     zcbin(nb,nac(n),n)=zcbin(nb,nac(n),n)+1
1252                     if(inoerr(n).gt.0)then
1253                       do nn=1,inoerr(n)
1254                         call xval(n,noerr(noerrhis(n),nn),lf,j,y2)
1255                         ebin(nb,nac(n),nn-1+noerrhis(n))=
1256      &                       ebin(nb,nac(n),nn-1+noerrhis(n))+y2
1257                         zebin(nb,nac(n),nn-1+noerrhis(n))=
1258      &                       zebin(nb,nac(n),nn-1+noerrhis(n))+1
1259                       enddo
1260                     endif
1261       end
1262
1263 c---------------------------------------------------------------------
1264       subroutine xhis(n)
1265 c---------------------------------------------------------------------
1266       include 'epos.inc'
1267       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
1268       parameter (mypara=10,mxpara=10)
1269       logical ilog,icnx,itrevt,idmod
1270       double precision bin,bbin,zcbin,zbbin
1271       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
1272      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
1273      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
1274      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
1275      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
1276      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
1277      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
1278      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
1279      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
1280      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
1281      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
1282       double precision ebin,zebin
1283       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
1284      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
1285       dimension xx(mxbin)
1286
1287       double precision histoweight
1288       common/chiswei/histoweight
1289       common/cyield/yield
1290       common/csigma/sigma
1291       double precision dcel
1292       common/ems3/dcel,ad
1293       common/geom/rmproj,rmtarg,bmax,bkmx
1294       save cnormx
1295
1296       if(ivar(1,n).eq.-1)then
1297         nrbins=0
1298         goto9999
1299       endif
1300
1301 c.......here normalization.......................................
1302 c           see also   "..........fill histogram"
1303 c.................................................................
1304 c     the norm ( inorm(n) ) is a number hijk which normalizes to:
1305 c
1306 c  k  0:  * 1
1307 c     1:  / number of events
1308 c     2:  / number of triggered events
1309 c     4:  / bin-counts
1310 c     6:  / number of summed bin-counts (yield=1.)
1311 c     7:  uses same normalization as one histo before
1312 c
1313 c  j  0:  * 1
1314 c     1:  / bin-width
1315 c     2:  * sigma_total / bin-width
1316 c     3:  * sigma_diff / bin-width
1317 c
1318 c  i  0:  * 1
1319 c     1:  y => y*x
1320 c     2:  y => y/x/2/pi (modified for mt0)
1321 c     3:  kno-scaling
1322 c     4:  y => y/x**1.5
1323 c     5:  y => y/x
1324 c     6:  y => y*xi (for conex, xi=x of the bin)
1325 c     7:  y => y/x/(x-m)
1326 c
1327 c  h  0: normal
1328 c     1: accumulated
1329 c
1330 c.................................................................
1331       norm1=mod(inorm(n),10)
1332       norm2=mod(inorm(n)/10,10)
1333       norm3=mod(inorm(n)/100,10)
1334       norm4=mod(inorm(n)/1000,10)
1335       nctbin=0
1336       do l=1,nbin(n)
1337         nctbin=nctbin+zcbin(l,nac(n),n)
1338         if(norm1.eq.4.and.zcbin(l,nac(n),n).ne.0d0)then
1339           bin(l,nac(n),n)=bin(l,nac(n),n)/zcbin(l,nac(n),n)
1340           if(lookcontr(n).gt.0)then
1341             do loo=lookcontr(n),lookcontrx(n)
1342               if(zbbin(l,nac(n),loo).ne.0.)
1343      &           bbin(l,nac(n),loo)=bbin(l,nac(n),loo)
1344      &            /zbbin(l,nac(n),loo)
1345             enddo
1346           endif
1347         endif
1348         if(ilog(n))then
1349           xx(l)=xmin(n)*(xmax(n)/xmin(n))**((float(l)-.5)/nbin(n))
1350         else
1351           xx(l)=(float(l)-0.5)*(xmax(n)-xmin(n))/nbin(n)+xmin(n)
1352         endif
1353       enddo
1354       cnorm=1.
1355       if(norm1.eq.1)cnorm=1./float(nevent)
1356       if(norm1.eq.2)then
1357         if(ncevt(n).ne.0)then
1358           cnorm=1./float(ncevt(n))
1359         else
1360           cnorm=0.
1361         endif
1362       endif
1363       if(norm1.eq.6.and.nctbin.ne.0)cnorm=1./float(nctbin)
1364       if(norm1.eq.7)cnorm=cnormx
1365       cnormx=cnorm
1366       if(ntevt.ne.0)
1367      &   sigma=10.*pi*bmax**2.*nevent/ntevt !total (untriggered) sigma
1368       if(norm2.eq.3)then      !differential (triggered) sigma
1369         if(ntevt.ne.0)
1370      &     sigma=10.*pi*bmax**2.*ncevt(n)/ntevt
1371       endif
1372       if(norm3.eq.3)then      !kno
1373         first=0.
1374         secnd=0.
1375         do l=1,nbin(n)
1376           if(nctbin.ne.0)first=first+xx(l)*zcbin(l,nac(n),n)/nctbin
1377           if(nctbin.ne.0)secnd=secnd
1378      $           +xx(l)**2*zcbin(l,nac(n),n)/nctbin
1379         enddo
1380       else
1381         first=1.
1382       endif
1383       if(ilog(n))then
1384         if(norm2.eq.2.or.norm2.eq.3) cnorm=cnorm*sigma
1385       else
1386         if(norm2.ge.1.and.norm2.le.3) cnorm=cnorm*xinc(n)
1387         if(norm2.eq.2.or.norm2.eq.3) cnorm=cnorm*sigma
1388       endif
1389       do l=1,nbin(n)
1390         if(ilog(n).and.norm2.ge.1.and.norm2.le.3)
1391      $      bin(l,nac(n),n) =  bin(l,nac(n),n)
1392      $      /(xmin(n)*exp(float(l)/xinc(n))*(1.-exp(-1./xinc(n))))
1393         bin(l,nac(n),n) =  bin(l,nac(n),n) * cnorm
1394       enddo
1395       f=first
1396       nrbins=nbin(n)
1397       nctbin=0
1398       yield=0.
1399       shft=0
1400        if(nint(xpara(1,n)).eq.999963)shft=xpara(2,n)
1401       do ii=1,nbin(n)
1402         g=1
1403         if(norm3.eq.1.and.xx(ii).ne.0.)g=1./xx(ii)
1404         if(norm3.eq.2)g=2*pi*(xx(ii)+shft)
1405         if(norm3.eq.4)g=xx(ii)**1.5
1406         if(norm3.eq.5)g=xx(ii)
1407         if(norm3.eq.7)g=0
1408         yield=yield+bin(ii,nac(n),n)/xinc(n)*hisfac*f*g
1409       enddo
1410       do l=1,nbin(n)
1411         x=(xx(l)+xshift)      !*xhfact
1412         ar(l,1)=x/f
1413         sigbin=0
1414         if(zcbin(l,nac(n),n).ne.0d0)
1415      *   sigbin=bin(l,nac(n),n)*hisfac*f/sqrt(zcbin(l,nac(n),n))
1416         if(norm4.eq.0.or.l.eq.1)then
1417           ar(l,3)=bin(l,nac(n),n)*hisfac*f
1418           if(lookcontr(n).gt.0)then
1419            do loo=lookcontr(n),lookcontrx(n)
1420              r=1
1421              if(norm1.eq.2.and.nccevt(loo).ne.0.)
1422      *          r=float(ncevt(n))/nccevt(loo)
1423              lo=loo-lookcontr(n)+1
1424              ary(l,lo)=bbin(l,nac(n),loo)*hisfac*f*cnorm*r
1425              if(zbbin(l,nac(n),loo).gt.0.)then
1426                ardy(l,lo)=ary(l,lo)/sqrt(zbbin(l,nac(n),loo))
1427              else
1428                ardy(l,lo)=0
1429              endif
1430              if(norm1.eq.4)ardy(l,lo)=zbbin(l,nac(n),loo)
1431             enddo
1432           endif
1433           if(norm3.eq.6)then   !conex
1434            ar(l,3)=ar(l,3)*xx(l)
1435           endif
1436           ar(l,4)=sigbin
1437         else
1438           ar(l,3)=ar(l-1,3)+bin(l,nac(n),n)*hisfac*f
1439           ar(l,4)=sqrt(ar(l-1,4)**2+sigbin**2)
1440         endif
1441         if(inoerr(n).ge.1)then
1442           if(zebin(l,nac(n),noerrhis(n)).gt.0.d0)then
1443          ar(l,4)=ebin(l,nac(n),noerrhis(n))/zebin(l,nac(n),noerrhis(n))
1444           else
1445             ar(l,4)=0.
1446           endif
1447         endif
1448         if(inoerr(n).eq.2)then
1449           if(zebin(l,nac(n),noerrhis(n)+1).gt.0.d0)then
1450       ar(l,5)=ebin(l,nac(n),noerrhis(n)+1)/zebin(l,nac(n),noerrhis(n)+1)
1451           else
1452             ar(l,5)=0.
1453           endif
1454         endif
1455         if(norm1.eq.4)ar(l,4)=zcbin(l,nac(n),n)
1456       enddo
1457       ionoerr=inoerr(n)
1458       histoweight=dble(ncevt(n))
1459       if(norm1.eq.4)histoweight=0d0
1460
1461  9999 hisfac=1.
1462       xshift=0
1463       end
1464
1465 c-----------------------------------------------------------------------
1466       integer function nsdiff(insdif,now)
1467 c-----------------------------------------------------------------------
1468 c  returns  1 if trigger condition for NSD fulfilled and 0 otherwise
1469 c  for  UA1 (insdif=1) or CDF (insdif=2) or STAR (insdif=3,4)
1470 C  or BRAHMS (insdif=5)
1471 c  now ... noweak(histogram number)
1472 c-----------------------------------------------------------------------
1473       include 'epos.inc'
1474       nsdiff=0
1475            if(insdif.ge.1)then
1476       iii1=0
1477       iii2=0
1478       ipos=0
1479       ineg=0
1480          do npts=1,nptl
1481         if(istptl(npts).ne.0)goto 60
1482         if(   idptl(npts).ne.120 .and.idptl(npts).ne.-120
1483      *   .and.idptl(npts).ne.130 .and.idptl(npts).ne.-130
1484      *   .and.idptl(npts).ne.1120.and.idptl(npts).ne.-1120
1485      *   .and.idptl(npts).ne.1130.and.idptl(npts).ne.-1130
1486      *   .and.idptl(npts).ne.2230.and.idptl(npts).ne.-2230
1487      *   .and.idptl(npts).ne.2330.and.idptl(npts).ne.-2330
1488      *   .and.idptl(npts).ne.3331.and.idptl(npts).ne.-3331)goto 60
1489 c        if(now.eq.1)then  !do not consider weak decay products
1490 c         if(iorptl(npts).ne.0)then
1491 c          idora=abs( idptl(iorptl(npts)) )
1492 c          if(  idora.eq.20   .or.idora.eq.2130
1493 c     &     .or.idora.eq.2230 .or.idora.eq.1130
1494 c     &     .or.idora.eq.2330 .or.idora.eq.1330
1495 c     &     .or.idora.eq.3331 )goto 60
1496 c          endif
1497 c         endif
1498         ppp=sqrt(pptl(1,npts)**2+pptl(2,npts)**2+pptl(3,npts)**2)
1499         if(ppp.gt.abs(pptl(3,npts)))then
1500           yyy=.5*log((ppp+pptl(3,npts))/(ppp-pptl(3,npts)))
1501         else
1502           yyy=sign(100.,pptl(3,npts))
1503         endif
1504         if(insdif.eq.1)then
1505           if(yyy.gt.2.   .and. yyy.lt.5.6)iii1=1
1506           if(yyy.gt.-5.6 .and. yyy.lt.-2.)iii2=1
1507         elseif(insdif.eq.2)then
1508           if(yyy.gt.3.2  .and. yyy.lt.5.9)iii1=1
1509           if(yyy.gt.-5.9 .and. yyy.lt.-3.2)iii2=1
1510           if(yyy.gt.0.   .and. yyy.lt.3.0)ipos=ipos+1
1511           if(yyy.gt.-3.0 .and. yyy.lt.0. )ineg=ineg+1
1512         elseif(insdif.eq.3)then
1513           if(yyy.gt.-5.0 .and. yyy.lt.-3.3 )iii1=1
1514           if(yyy.gt. 3.3 .and. yyy.lt. 5.0 )iii2=1
1515         elseif(insdif.eq.4)then
1516           if(yyy.gt.-5.0 .and. yyy.lt.-3.1 )iii1=1
1517           if(yyy.gt. 3.1 .and. yyy.lt. 5.0 )iii2=1
1518         elseif(insdif.eq.5)then
1519           if(yyy.gt.-5.25 .and. yyy.lt.-3.26 )iii1=1
1520           if(yyy.gt. 3.26 .and. yyy.lt. 5.25 )iii2=1
1521         endif
1522 60      continue
1523          enddo
1524         if(insdif.eq.1)then
1525           if(iii1.eq.1 .and. iii2.eq.1) nsdiff=1
1526         elseif(insdif.eq.2)then
1527           if((iii1.eq.1 .and. iii2.eq.1) .and.
1528      *    ((ipos.ne.0 .and. ineg.ne.0) .and. ipos+ineg.ge.4)) nsdiff=1
1529         elseif(insdif.eq.3.or.insdif.eq.4.or.insdif.eq.5)then
1530           if(iii1.eq.1 .and. iii2.eq.1) nsdiff=1
1531         endif
1532            else
1533       stop'in nsdiff. argument of nsdiff not authorized.        '
1534            endif
1535       end
1536
1537 c----------------------------------------------------------------------
1538       subroutine xtrans(cvar,inom,ifr,n)
1539 c----------------------------------------------------------------------
1540       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
1541      &     ,multc1
1542       parameter(mxxhis=50)
1543       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
1544      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis)
1545
1546       character*6 cvar
1547       ifr=0
1548       if(cvar.eq.'numptl')then
1549         inom=1
1550       elseif(cvar.eq.'npaptl')then
1551         inom=2
1552       elseif(cvar.eq.'npmptl')then
1553         inom=3
1554       elseif(cvar.eq.'ispptl')then
1555         inom=4
1556       elseif(cvar.eq.'rapx')then
1557         inom=5
1558       elseif(cvar.eq.'iptlfr')then
1559         inom=6
1560       elseif(cvar.eq.'rinp')then
1561         inom=7
1562       elseif(cvar.eq.'eco')then
1563         inom=8
1564       elseif(cvar.eq.'absrap')then
1565         inom=12
1566       elseif(cvar.eq.'rap')then
1567         inom=13
1568       elseif(cvar.eq.'xp')then
1569         inom=14
1570       elseif(cvar.eq.'xe')then
1571         inom=15
1572       elseif(cvar.eq.'pt')then
1573         inom=16
1574       elseif(cvar.eq.'p1a')then
1575         inom=17
1576       elseif(cvar.eq.'p2a')then
1577         inom=18
1578       elseif(cvar.eq.'xi')then
1579         inom=19
1580       elseif(cvar.eq.'xf')then
1581         inom=20
1582       elseif(cvar.eq.'t')then
1583         inom=21
1584       elseif(cvar.eq.'rapmi')then
1585         inom=22
1586       elseif(cvar.eq.'eta')then
1587         inom=23
1588       elseif(cvar.eq.'theta')then
1589         inom=24
1590       elseif(cvar.eq.'pt2')then
1591         inom=25
1592       elseif(cvar.eq.'et')then
1593         inom=26
1594       elseif(cvar.eq.'idptl')then
1595         inom=27
1596       elseif(cvar.eq.'istptl')then
1597         inom=28
1598       elseif(cvar.eq.'mass')then
1599         inom=29
1600       elseif(cvar.eq.'idaptl')then
1601         inom=30
1602       elseif(cvar.eq.'egy')then
1603         inom=31
1604       elseif(cvar.eq.'rapwro')then
1605         inom=32
1606       elseif(cvar.eq.'mt')then
1607         inom=33
1608       elseif(cvar.eq.'pplus')then
1609         inom=34
1610       elseif(cvar.eq.'pminus')then
1611         inom=35
1612       elseif(cvar.eq.'p5')then
1613         inom=36
1614       elseif(cvar.eq.'pa')then
1615         inom=37
1616       elseif(cvar.eq.'sob')then
1617         inom=38
1618       elseif(cvar.eq.'idpom')then
1619         inom=39
1620       elseif(cvar.eq.'p3a')then
1621         inom=40
1622       elseif(cvar.eq.'cmass')then
1623         inom=41
1624       elseif(cvar.eq.'arappi')then
1625         inom=42
1626       elseif(cvar.eq.'itsptl')then
1627         inom=50
1628       elseif(cvar.eq.'ityptl')then
1629         inom=51
1630       elseif(cvar.eq.'idoptl')then
1631         inom=52
1632       elseif(cvar.eq.'iptl')then
1633         inom=53
1634       elseif(cvar.eq.'index')then
1635         inom=54
1636       elseif(cvar.eq.'p1')then
1637         inom=55
1638       elseif(cvar.eq.'p2')then
1639         inom=56
1640       elseif(cvar.eq.'p3')then
1641         inom=57
1642       elseif(cvar.eq.'p4')then
1643         inom=58
1644       elseif(cvar.eq.'xg')then
1645         inom=59
1646       elseif(cvar.eq.'ek')then
1647         inom=60
1648       elseif(cvar.eq.'beta')then
1649         inom=61
1650       elseif(cvar.eq.'mt0')then
1651         inom=63
1652       elseif(cvar.eq.'qsqptl')then
1653         inom=64
1654       elseif(cvar.eq.'xelab')then
1655         inom=65
1656       elseif(cvar.eq.'hgtc05')then
1657         inom=66
1658         imulty1=1             !to switch on the calculation of "Standard variable"
1659       elseif(cvar.eq.'hadtyp')then
1660         inom=67
1661         imulty1=1
1662       elseif(cvar.eq.'hgtc1')then
1663         inom=68
1664       elseif(cvar.eq.'x4')then
1665         inom=69
1666       elseif(cvar.eq.'npn')then
1667         inom=70
1668       elseif(cvar.eq.'routp')then
1669         inom=71
1670       elseif(cvar.eq.'hgtc3')then
1671         inom=72
1672         imulty1=1
1673       elseif(cvar.eq.'mu14')then
1674         inom=73
1675         imulty1=1
1676       elseif(cvar.eq.'delphi')then
1677         inom=74
1678         iok=0
1679         !------------------------------------------------------------
1680         !icorrtrig sores the histogram numbers of those histograms which
1681         !use the delphi variable (and therfore require a call corrtrig
1682         !------------------------------------------------------------
1683         do i=1,icorrtrig(0)
1684          if(icorrtrig(i).eq.n)iok=1
1685         enddo
1686         if(iok.eq.0)then
1687           icorrtrig(0)=icorrtrig(0)+1
1688           if(icorrtrig(0).gt.mxxhis)stop'mxxhis too small'
1689         icorrtrig(icorrtrig(0))=n
1690         endif
1691       elseif(cvar.eq.'v2')then
1692         inom=75
1693       elseif(cvar.eq.'pt4')then
1694         inom=76
1695       elseif(cvar.eq.'rin')then
1696         inom=77
1697       elseif(cvar.eq.'mulevt')then
1698         inom=101
1699       elseif(cvar.eq.'etevt')then
1700         inom=102
1701       elseif(cvar.eq.'enevt')then
1702         inom=103
1703       elseif(cvar.eq.'ev6evt')then
1704         inom=104
1705       elseif(cvar.eq.'xenevt')then
1706         inom=105
1707       elseif(cvar.eq.'netevt')then
1708         inom=106
1709       elseif(cvar.eq.'ptevt')then
1710         inom=107
1711       elseif(cvar.eq.'pmxevt')then
1712         inom=108
1713       elseif(cvar.eq.'numevt')then
1714         inom=201
1715       elseif(cvar.eq.'egyevt')then
1716         inom=202
1717       elseif(cvar.eq.'bimevt')then
1718         inom=203
1719       elseif(cvar.eq.'xbjevt')then
1720         inom=204
1721       elseif(cvar.eq.'qsqevt')then
1722         inom=205
1723       elseif(cvar.eq.'yevt')then
1724         inom=206
1725       elseif(cvar.eq.'eloevt')then
1726         inom=207
1727       elseif(cvar.eq.'nd1evt')then
1728         inom=208
1729       elseif(cvar.eq.'nd2evt')then
1730         inom=209
1731       elseif(cvar.eq.'theevt')then
1732         inom=210
1733       elseif(cvar.eq.'nspevt')then
1734         inom=211
1735       elseif(cvar.eq.'nhpevt')then
1736         inom=212
1737       elseif(cvar.eq.'sigtot')then
1738         inom=213
1739       elseif(cvar.eq.'sigela')then
1740         inom=214
1741       elseif(cvar.eq.'sloela')then
1742         inom=215
1743       elseif(cvar.eq.'nrgevt')then
1744         inom=216
1745       elseif(cvar.eq.'qevt')then
1746         inom=217
1747       elseif(cvar.eq.'qtlevt')then
1748         inom=218
1749       elseif(cvar.eq.'threvt')then
1750         inom=220
1751         ifr=33                  !set thrust-frame
1752       elseif(cvar.eq.'omtevt')then
1753         inom=221
1754         ifr=33                  !set thrust-frame
1755       elseif(cvar.eq.'tmaevt')then
1756         inom=222
1757         ifr=33                  !set thrust-frame
1758       elseif(cvar.eq.'tmievt')then
1759         inom=223
1760         ifr=33                  !set thrust-frame
1761       elseif(cvar.eq.'oblevt')then
1762         inom=224
1763         ifr=33                  !set thrust-frame
1764       elseif(cvar.eq.'sphevt')then
1765         inom=230
1766         ifr=32                  !set sph-frame
1767       elseif(cvar.eq.'aplevt')then
1768         inom=231
1769         ifr=32                  !set sph-frame
1770       elseif(cvar.eq.'cpaevt')then
1771         inom=232
1772         ifr=34                  !set sph2-frame
1773       elseif(cvar.eq.'dpaevt')then
1774         inom=233
1775         ifr=34                  !set sph2-frame
1776       elseif(cvar.eq.'npoevt')then
1777         inom=234
1778       elseif(cvar.eq.'npnevt')then
1779         inom=235
1780
1781   !....unused 236, 237
1782
1783       elseif(cvar.eq.'npxevt')then
1784         inom=238
1785       elseif(cvar.eq.'mu1evt')then
1786         inom=240
1787         imulty1=1
1788       elseif(cvar.eq.'muievt')then
1789         inom=241
1790         imulty1=1
1791       elseif(cvar.eq.'hgtevt')then
1792         inom=242
1793         imulty1=1
1794       elseif(cvar.eq.'difevt')then
1795         inom=243
1796       elseif(cvar.eq.'dixevt')then
1797         inom=244
1798       elseif(cvar.eq.'qinevt')then
1799         inom=250
1800       elseif(cvar.eq.'qfievt')then
1801         inom=251
1802       elseif(cvar.eq.'einevt')then
1803         inom=252
1804       elseif(cvar.eq.'efievt')then
1805         inom=253
1806       elseif(cvar.eq.'pinevt')then
1807         inom=254
1808       elseif(cvar.eq.'pfievt')then
1809         inom=255
1810       elseif(cvar.eq.'pxfevt')then    ! leading proton xf in cms
1811         inom=256
1812       elseif(cvar.eq.'pi+xf')then     ! pi+xf: pi+ yield at cms xf>0.01
1813         inom=257
1814       elseif(cvar.eq.'pi-xf')then     ! pi-xf: pi- yield at cms xf>0.01
1815         inom=258
1816       elseif(cvar.eq.'sigcut')then
1817         inom=260
1818       elseif(cvar.eq.'keu')then
1819         inom=261
1820       elseif(cvar.eq.'ked')then
1821         inom=262
1822       elseif(cvar.eq.'kes')then
1823         inom=263
1824       elseif(cvar.eq.'kolevt')then
1825         inom=265
1826       elseif(cvar.eq.'sigsd')then
1827         inom=266
1828       elseif(cvar.eq.'nglevt')then
1829         inom=267
1830       elseif(cvar.eq.'kppevt')then   ! collision numbers per participant
1831         inom=268
1832       elseif(cvar.eq.'npievt')then   ! pion + multiplicity per event
1833         inom=269
1834       elseif(cvar.eq.'np2evt')then   ! pion + multiplicity per participant
1835         inom=270
1836       elseif(cvar.eq.'sigdif'.or.cvar.eq.'sigdifr')then
1837         inom=271
1838       elseif(cvar.eq.'koievt')then
1839         inom=272
1840       elseif(cvar.eq.'ineevt')then
1841         inom=273
1842       elseif(cvar.eq.'elaevt')then
1843         inom=274
1844       elseif(cvar.eq.'itgevt')then
1845         inom=275
1846         iok=0
1847         do i=1,icorrtrig(0)
1848           if(icorrtrig(i).eq.n)iok=1
1849         enddo
1850         if(iok.eq.0)then
1851           icorrtrig(0)=icorrtrig(0)+1
1852           if(icorrtrig(0).gt.mxxhis)stop'mxxhis too small'
1853           icorrtrig(icorrtrig(0))=n
1854         endif
1855       elseif(cvar.eq.'hrdevt')then
1856         inom=276
1857         iok=0
1858         do i=1,ihardevent(0)
1859           if(ihardevent(i).eq.n)iok=1
1860         enddo
1861         if(iok.eq.0)then
1862           ihardevent(0)=ihardevent(0)+1
1863           if(ihardevent(0).gt.mxxhis)stop'mxxhis too small'
1864           ihardevent(ihardevent(0))=n
1865         endif
1866       elseif(cvar(2:6).eq.'j1evt'.or.cvar(2:6).eq.'j2evt')then
1867         iok=0
1868         do i=1,ijetfind1(0)
1869           if(ijetfind1(i).eq.n)iok=1
1870         enddo
1871         if(iok.eq.0)then
1872           ijetfind1(0)=ijetfind1(0)+1
1873           if(ijetfind1(0).gt.mxxhis)stop'mxxhis too small'
1874           ijetfind1(ijetfind1(0))=n
1875         endif
1876         if(cvar.eq.'ej1evt')inom=277
1877         if(cvar.eq.'pj1evt')inom=278
1878         if(cvar(2:6).eq.'j2evt')then
1879           iok=0
1880           do i=1,ijetfind2(0)
1881             if(ijetfind2(i).eq.n)iok=1
1882           enddo
1883           if(iok.eq.0)then
1884             ijetfind2(0)=ijetfind2(0)+1
1885             if(ijetfind2(0).gt.mxxhis)stop'mxxhis too small'
1886             ijetfind2(ijetfind2(0))=n
1887           endif
1888           if(cvar.eq.'ej2evt')inom=279
1889           if(cvar.eq.'pj2evt')inom=280
1890         endif
1891       elseif(cvar.eq.'zppevt')then
1892         inom=281
1893       elseif(cvar.eq.'zptevt')then
1894         inom=282
1895       elseif(cvar.eq.'***not used***')then
1896         inom=283
1897       elseif(cvar.eq.'nd3evt')then
1898         inom=284
1899       elseif(cvar.eq.'nd4evt')then
1900         inom=285
1901       elseif(cvar.eq.'nd5evt')then
1902         inom=287
1903       elseif(cvar.eq.'mubevt')then
1904         inom=286
1905         imulty1=1
1906       elseif(cvar.eq.'aimevt')then
1907         inom=301
1908       elseif(cvar.eq.'wjbevt')then
1909         inom=302
1910       elseif(cvar.eq.'njbevt')then
1911         inom=303
1912       elseif(cvar.eq.'djbevt')then
1913         inom=304
1914       elseif(cvar.eq.'tjbevt')then
1915         inom=305
1916       elseif(cvar.eq.'hjmevt')then
1917         inom=306
1918       elseif(cvar.eq.'ljmevt')then
1919         inom=307
1920       elseif(cvar.eq.'djmevt')then
1921         inom=308
1922       elseif(cvar.eq.'ybal')then
1923         inom=310
1924       elseif(cvar.eq.'yabal')then
1925         inom=310
1926       elseif(cvar.eq.'sigine')then
1927         inom=312
1928       elseif(cvar.eq.'sigiaa')then
1929         inom=313
1930       elseif(cvar.eq.'alpdsf')then
1931         inom=314
1932       elseif(cvar.eq.'alpdsh')then
1933         inom=315
1934       elseif(cvar.eq.'betdsf')then
1935         inom=316
1936       elseif(cvar.eq.'betdsh')then
1937         inom=317
1938       elseif(cvar.eq.'rexdip')then
1939         inom=318
1940       elseif(cvar.eq.'rexdit')then
1941         inom=319
1942       elseif(cvar.eq.'m14evt')then
1943         inom=320
1944       elseif(cvar.eq.'ht3evt')then
1945         inom=321
1946       elseif(cvar.eq.'sigiex')then
1947         inom=322
1948       elseif(cvar.eq.'sigdex')then
1949         inom=323
1950       elseif(cvar.eq.'sigsex')then
1951         inom=324
1952       elseif(cvar.eq.'ox1evt')then
1953         inom=501
1954       elseif(cvar.eq.'ox2evt')then
1955         inom=502
1956       elseif(cvar.eq.'ox3evt')then
1957         inom=503
1958       elseif(cvar.eq.'ox4evt')then
1959         inom=504
1960       else
1961         print *,' '
1962         print *,'              xtrans: unknown variable ',cvar
1963         print *,' '
1964 c       inom=-1
1965         stop
1966       endif
1967       end
1968
1969 c----------------------------------------------------------------------
1970       subroutine xval(n,inom,lf,j,x)
1971 c----------------------------------------------------------------------
1972 c   n ...... histogram index
1973 c   inom ... variable index
1974 c              1-100 particle variables
1975 c              101-200 accumulative event variables
1976 c              > 200 other event variables
1977 c   lf ..... frame index
1978 c   particle index (used for particle variables)
1979 c----------------------------------------------------------------------
1980       include 'epos.inc'
1981       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
1982      &     ,multc1
1983       parameter(mxxhis=50)
1984       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
1985      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis)
1986       common/zeus2/qtl
1987
1988       parameter (ntim=1000)
1989       common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
1990      &     ,idaprt(2,ntim)
1991
1992       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
1993      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
1994
1995       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
1996       parameter (mypara=10,mxpara=10)
1997       logical ilog,icnx,itrevt,idmod
1998       double precision bin,bbin,zcbin,zbbin
1999       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
2000      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
2001      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
2002      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
2003      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
2004      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
2005      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
2006      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
2007      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
2008      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
2009      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
2010       double precision ebin,zebin
2011       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
2012      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
2013       parameter (mxfra=5)
2014       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
2015      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
2016      $     ,emax(mxfra)
2017
2018       double precision bofra
2019       common/dfra/bofra(5,mxfra)
2020       dimension p(5,mxfra),aimuni(10,mxhis),xor(5,mxfra)
2021       save p,aimuni,xor
2022       if(iffra(lf).eq.0.and.j.ne.0)then
2023         do l=1,5
2024           p(l,lf)=pptl(l,j)
2025         enddo
2026         do l=1,4
2027           xor(l,lf)=xorptl(l,j)
2028         enddo
2029         if(imofra(1,lf).ne.0)then
2030           call utrota(imofra(1,lf),r1fra(1,lf),r1fra(2,lf),r1fra(3,lf)
2031      $         ,p(1,lf),p(2,lf),p(3,lf))
2032           call utrota(imofra(1,lf),r1fra(1,lf),r1fra(2,lf),r1fra(3,lf)
2033      $         ,xor(1,lf),xor(2,lf),xor(3,lf))
2034         endif
2035         if(imofra(2,lf).ne.0)then !the x-z exchanged is ok !!
2036           call utrota(imofra(2,lf),r2fra(3,lf),r2fra(2,lf),r2fra(1,lf)
2037      $         ,p(3,lf),p(2,lf),p(1,lf))
2038           call utrota(imofra(2,lf),r2fra(3,lf),r2fra(2,lf),r2fra(1,lf)
2039      $         ,xor(3,lf),xor(2,lf),xor(1,lf))
2040         endif
2041         if(imofra(3,lf).ne.0)then
2042           call utlob4(imofra(3,lf),bofra(1,lf),bofra(2,lf)
2043      $         ,bofra(3,lf) ,bofra(4,lf),bofra(5,lf)
2044      $         ,p(1,lf),p(2,lf),p(3,lf),p(4,lf))
2045           call utlob4(imofra(3,lf),bofra(1,lf),bofra(2,lf)
2046      $         ,bofra(3,lf) ,bofra(4,lf),bofra(5,lf)
2047      $         ,xor(1,lf),xor(2,lf),xor(3,lf),xor(4,lf))
2048         endif
2049         iffra(lf)=1
2050       endif
2051
2052 c--------------------------------- 1 - 100 ----------------------------
2053       if(inom.eq.1)then
2054         x=1.
2055       elseif(inom.eq.2)then
2056         x=isign(1,idptl(j))
2057       elseif(inom.eq.3)then
2058         chrg=0
2059         if(iabs(idptl(j)).le.9999
2060      $       .and.mod(iabs(idptl(j)),10).le.1)
2061      $       call idchrg(idptl(j),chrg)
2062         if(chrg.eq.0.)then
2063           x=0
2064         else
2065           x=int(sign(1.,chrg))
2066         endif
2067       elseif(inom.eq.4)then
2068         iad=abs(idptl(j))
2069         jspin=mod(iad,10)
2070         x=0.
2071         if (iad.ge.100.and.iad.lt.1000) x=1./(1.+2.*jspin)
2072         if (iad.ge.1000.and.iad.lt.9999) x=1./(2.+2*jspin)
2073       elseif(inom.eq.5)then    !'rapx'  !st-rap for string segments only !!!!!!!!!!
2074         x=dezptl(j)
2075       elseif(inom.eq.6)then                                      !'iptlfr'
2076         x=0
2077         if(j.ge.minfra.and.j.le.maxfra)x=1
2078       elseif(inom.eq.7)then                                        !'rinp'
2079         aa=cos(phievt)
2080         bb=sin(phievt)
2081         x=xptl(j)*aa+yptl(j)*bb
2082       elseif(inom.eq.8)then                        !'eco' !engy in comoving frame
2083         x=0
2084         amt=p(5,lf)**2+p(1,lf)**2+p(2,lf)**2
2085         if(amt.gt.0..and.p(4,lf)+abs(p(3,lf)).gt.0.d0)then
2086           amt=sqrt(amt)
2087           rap=sign(1.,p(3,lf))*alog((p(4,lf)+abs(p(3,lf)))/amt)
2088           rapx=dezptl(j)
2089           x=amt*cosh(rap-rapx)
2090         endif
2091       elseif(inom.eq.12)then                                      !'absrap'
2092         amt=p(5,lf)**2+p(1,lf)**2+p(2,lf)**2
2093         if(amt.gt.0..and.p(4,lf)+abs(p(3,lf)).gt.0.d0)then
2094           amt=sqrt(amt)
2095           x=alog((p(4,lf)+abs(p(3,lf)))/amt)
2096         else
2097           x=0.                  !
2098         endif
2099       elseif(inom.eq.13)then    !'rap'
2100         amt=p(5,lf)**2+p(1,lf)**2+p(2,lf)**2
2101         if(amt.gt.0..and.p(4,lf)+abs(p(3,lf)).gt.0.d0)then
2102           amt=sqrt(amt)
2103           x=sign(1.,p(3,lf))*alog((p(4,lf)+abs(p(3,lf)))/amt)
2104         else
2105           x=0.                  !
2106         endif
2107       elseif(inom.eq.14)then                                         !'xp'
2108         x=sqrt(p(3,lf)**2+p(2,lf)**2+p(1,lf)**2)/emax(lf)
2109       elseif(inom.eq.15)then    !'xe'
2110         x=p(4,lf)/emax(lf)
2111       elseif(inom.eq.16)then                                         !'pt'
2112         x=sqrt(p(2,lf)**2+p(1,lf)**2)
2113       elseif(inom.eq.17)then
2114         x=abs(p(1,lf))
2115       elseif(inom.eq.18)then
2116         x=abs(p(2,lf))
2117       elseif(inom.eq.19)then
2118         x=-log(sqrt(p(3,lf)**2+p(2,lf)**2+p(1,lf)**2)/emax(lf))
2119       elseif(inom.eq.20)then                                     !'xf'
2120         m=mod(ifra(lf)/10,10)
2121         if(m.eq.1)then
2122 c          pmax=sqrt((engy/2)**2-prom*2)
2123           pmax=pnullx               !???????????????????
2124           if(ifra(lf).eq.12)pmax=pnll
2125           x=p(3,lf)/pmax
2126         else
2127           x=p(3,lf)/emax(lf)
2128         endif
2129       elseif(inom.eq.21)then
2130         pmax=pmxevt
2131 c        pmax=sqrt((engy/2)**2-prom*2)
2132           pmax=pnullx               !???????????????????
2133         if(ifra(lf).eq.12)pmax=pnll
2134         x=-(amproj**2-2.*sqrt(amproj**2+pmax**2)*p(4,lf)
2135      *      +2.*abs(pmax*p(3,lf))+p(5,lf)**2)
2136       elseif(inom.eq.22)then
2137         amt=sqrt(p(5,lf)**2+p(1,lf)**2+p(2,lf)**2)
2138         if(amt.ne.0.)then
2139           x=-sign(1.,p(3,lf))*alog((p(4,lf)+abs(p(3,lf)))/amt)
2140         else
2141           x=0.                  !
2142         endif
2143       elseif(inom.eq.23)then                                     !'eta'
2144         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2145         x=0
2146         if(p(3,lf).ne.0..and.pt.ne.0.)x=sign(1.,p(3,lf))*
2147      *       alog((sqrt(p(3,lf)**2+pt**2)+abs(p(3,lf)))/pt)
2148       elseif(inom.eq.24)then                                     !'theta'
2149         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2150         x=90
2151         if(p(3,lf).ne.0.)x=atan(pt/p(3,lf))/pi*180.
2152         if(x.lt.0.)x=180.+x
2153       elseif(inom.eq.25)then                                     !'pt2'
2154         x=p(2,lf)**2+p(1,lf)**2
2155       elseif(inom.eq.26)then                                     !'et'
2156         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2157         x=0
2158         eef=p(4,lf)
2159 c        if(idptl(j).ge.1000)eef=eef-prom
2160 c        if(idptl(j).le.-1000)eef=eef+prom
2161         p2=p(3,lf)**2+p(2,lf)**2+p(1,lf)**2
2162         if(p2.ne.0.)x=eef*pt/sqrt(p2)
2163       elseif(inom.eq.27)then                                     !'idptl'
2164         x=idptl(j)
2165       elseif(inom.eq.28)then    !istptl
2166         x=istptl(j)
2167       elseif(inom.eq.29)then    !mass
2168         x=p(5,lf)
2169         if(istptl(j).le.1)call idmass(idptl(j),x)
2170       elseif(inom.eq.30)then    !idaptl
2171         x=abs(idptl(j))
2172       elseif(inom.eq.31)then    !egy
2173         x=egyevt
2174       elseif(inom.eq.32)then    !arapwro
2175         x=0
2176         pt2=p(2,lf)**2+p(1,lf)**2
2177         if(p(3,lf).ne.0.)x=sign(1.,p(3,lf))*
2178      *       alog((sqrt(p(3,lf)**2+pt2+.13957**2)+abs(p(3,lf)))
2179      *       /sqrt(pt2+.13957**2))
2180       elseif(inom.eq.33)then                                  !'mt'
2181         x=sqrt(p(2,lf)**2+p(1,lf)**2+p(5,lf)**2)
2182       elseif(inom.eq.34)then                                  !'pplus'
2183         x=sign(1.,p(3,lf)) * (p(4,lf)+abs(p(3,lf)))
2184       elseif(inom.eq.35)then                                  !'pminus'
2185         x=sign(1.,p(3,lf)) * (p(4,lf)-abs(p(3,lf)))
2186       elseif(inom.eq.36)then                                  !'p5' (mass)
2187         x=p(5,lf)
2188       elseif(inom.eq.37)then                                  !pa
2189         x=sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
2190       elseif(inom.eq.38)then                                  !'pa'
2191         if(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2.ne.0)
2192      *       x=egyevt**2/sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)*p(4,lf)
2193       elseif(inom.eq.39)then                                  !idpom
2194         x=idptl(j)/1000000
2195       elseif(inom.eq.40)then                                  !p3a
2196         x=abs(p(3,lf))
2197       elseif(inom.eq.41)then
2198         cm2=p(4,lf)**2-p(3,lf)**2-p(2,lf)**2-p(1,lf)**2         !cmass
2199         x=sign(sqrt(abs(cm2)),cm2)
2200       elseif(inom.eq.42)then    !arappi
2201         x=0
2202         pt2=p(2,lf)**2+p(1,lf)**2
2203         if(p(3,lf).ne.0.)
2204      *       x=alog((sqrt(p(3,lf)**2+pt2+.13957**2)+abs(p(3,lf)))
2205      *       /sqrt(pt2+.13957**2))
2206       elseif(inom.eq.50)then
2207         x=itsptl(j)
2208       elseif(inom.eq.51)then
2209         x=ityptl(j)
2210       elseif(inom.eq.52)then
2211         x=0.
2212         if(iorptl(j).ne.0) x=idptl(iorptl(j))
2213       elseif(inom.eq.53)then
2214         x=j
2215       elseif(inom.eq.54)then                       !sloela
2216         call idflav(idptl(j),ifl1,ifl2,ifl3,jspin,indx)
2217         x=indx
2218       elseif(inom.eq.55)then                       !p_x
2219         x=p(1,lf)
2220       elseif(inom.eq.56)then                       !p_y
2221         x=p(2,lf)
2222       elseif(inom.eq.57)then                       !p_z
2223         x=p(3,lf)
2224       elseif(inom.eq.58)then                       !'e'
2225         x=p(4,lf)
2226       elseif(inom.eq.59)then                       !E/p_max
2227 c        pmax=sqrt((engy/2)**2-prom*2)
2228           pmax=pnullx               !???????????????????
2229         if(ifra(lf).eq.12)pmax=pnll
2230         x=p(4,lf)/pmax
2231       elseif(inom.eq.60)then                       !e_kin
2232         x=p(4,lf)-p(5,lf)
2233       elseif(inom.eq.61)then                       !'beta'
2234         x=p(3,lf)/p(4,lf)
2235       elseif(inom.eq.63)then                       !'mt0'
2236         x=sqrt(p(2,lf)**2+p(1,lf)**2+p(5,lf)**2)-p(5,lf)
2237       elseif(inom.eq.64)then                       !qsqptl
2238         x=qsqptl(j)
2239       elseif(inom.eq.65)then                       !xelab=Elab/Eolab
2240         x=p(4,lf)/(engy**2/2/prom-prom)
2241         if(x.gt.0.9999999) x=.9999999
2242       elseif(inom.eq.66)then    !hgtc05 ... charged ptl mult |[c]|<0.5
2243         x=multc05
2244       elseif(inom.eq.67)then    !hadtyp ... primary (1) or secondary (2) hadron
2245         if(j.le.nbdky)then
2246           x=1
2247         else
2248           x=2
2249         endif
2250       elseif(inom.eq.68)then    !hgtc1
2251         x=multc1
2252       elseif(inom.eq.69)then                       !'x4'
2253         x=xor(4,lf)
2254       elseif(inom.eq.70)then                       !'npn'
2255         x=npjevt+ntgevt
2256       elseif(inom.eq.71)then                       !'routp'
2257         cc=-sin(phievt)
2258         dd=cos(phievt)
2259         x=xptl(j)*cc+yptl(j)*dd
2260       elseif(inom.eq.72)then    !hgtc3 ... charged ptl mult |eta|<3.15  /6.3
2261         x=multc3/6.3
2262       elseif(inom.eq.73)then    !mu14 ... charged ptl mult |eta|<1  pt>.4
2263         x=multc14
2264       elseif(inom.eq.74)then    !delphi ... azimuthhal correlation
2265         x=10000.
2266         pt=sqrt(p(1,lf)**2+p(2,lf)**2)
2267
2268         if(nint(ypara(1,n)).ne.0.and.j.ne.nint(ypara(1,n)).and.
2269      $           pt.gt.0)then
2270            phi=sign(1.,p(2,lf))*acos(p(1,lf)/pt)
2271                x=phi-ypara(2,n)
2272            if(x.lt.-2.5*pi)then
2273             x=x+4*pi
2274            elseif(x.lt.-0.5*pi)then
2275             x=x+2*pi
2276           elseif(x.gt.3.5*pi)then
2277             x=x-4*pi
2278           elseif(x.gt.1.5*pi)then
2279             x=x-2*pi
2280           endif
2281         endif
2282       elseif(inom.eq.75)then    !'v2'
2283         aa=cos(phievt)
2284         bb=sin(phievt)
2285         cc=-sin(phievt)
2286         dd=cos(phievt)
2287         px=p(1,lf)*aa+p(2,lf)*bb
2288         py=p(1,lf)*cc+p(2,lf)*dd
2289         pt2=p(2,lf)**2+p(1,lf)**2
2290         x=0
2291         if(pt2.gt.0.)x=(px**2-py**2)/pt2
2292       elseif(inom.eq.76)then                                     !'pt4'
2293         x=(p(2,lf)**2+p(1,lf)**2)**2
2294       elseif(inom.eq.77)then                                     !'rin'
2295         x=rinptl(j)
2296
2297
2298 c--------------------------------- 101 - 200 ----------------------------
2299
2300       elseif(inom.eq.101)then           !mulevt
2301         x=1.
2302       elseif(inom.eq.102)then                      !'etevt'
2303         x=0
2304         if(istptl(j).eq.0)then
2305          eef=p(4,lf)
2306          if(idptl(j).ge.1000)eef=eef-prom
2307          if(idptl(j).le.-1000)eef=eef+prom
2308          pp=sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
2309          if(pp.ne.0.)x=eef*sqrt(p(1,lf)**2+p(2,lf)**2)/pp
2310          if(x.ne.x)then
2311            write(ifch,*)x,eef,p(1,lf),p(2,lf),p(3,lf),pp,prom,idptl(j),j
2312            call alist('xan&',1,nptl)
2313            stop 'probleme dans xan'
2314          endif
2315          endif
2316       elseif(inom.eq.103)then
2317         x=p(4,lf)/1000.
2318       elseif(inom.eq.104)then                       !'ev6evt'
2319         x=0
2320         if(istptl(j).eq.0)then
2321          pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2322          eta=0
2323          if(p(3,lf).ne.0..and.pt.ne.0.)eta=sign(1.,p(3,lf))*
2324      *   alog((sqrt(p(3,lf)**2+pt**2)+abs(p(3,lf)))/pt)
2325          if(pt.eq.0.)eta=sign(1e5,p(3,lf))
2326          if(eta.gt.6.0)then
2327           eef=p(4,lf)
2328           if(idptl(j).ge.1000)eef=eef-prom
2329           if(idptl(j).le.-1000)eef=eef+prom
2330           pp=sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
2331           if(pp.ne.0.)x=0.001*eef
2332          endif
2333         endif
2334       elseif(inom.eq.105)then
2335         etot=maproj*emax(lf)+matarg*0.94  !nur richtig fur target frame!!!!!
2336         x=p(4,lf)/etot
2337       elseif(inom.eq.106)then
2338         x=isign(1,idptl(j))
2339       elseif(inom.eq.107)then                       !'ptevt'
2340         x=sqrt(p(2,lf)**2+p(1,lf)**2)
2341       elseif(inom.eq.108)then                       !'pmxevt'
2342         x=pmxevt
2343
2344 c--------------------------------- > 200 ----------------------------
2345
2346       elseif(inom.eq.201)then
2347         x=1.
2348       elseif(inom.eq.202)then
2349         x=egyevt
2350       elseif(inom.eq.203)then
2351         x=bimevt
2352       elseif(inom.eq.204)then                       !'xbjevt'
2353         x=xbjevt
2354       elseif(inom.eq.205)then                       !'qsqevt'
2355         x=qsqevt
2356       elseif(inom.eq.206)then                       !'yevt'
2357         x=qsqevt/xbjevt/engy**2
2358       elseif(inom.eq.207)then                       !'eloevt'
2359         x=qsqevt/4./elepti+elepti*(1.-qsqevt/xbjevt/engy**2)
2360       elseif(inom.eq.208)then                       !nd1evt
2361         x=nsdiff(1,noweak(n))
2362       elseif(inom.eq.209)then                       !'nd2evt'
2363         x=nsdiff(2,noweak(n))
2364       elseif(inom.eq.284)then                       !'nd3evt'
2365         x=nsdiff(3,noweak(n))
2366       elseif(inom.eq.285)then                       !'nd4evt'
2367         x=nsdiff(4,noweak(n))
2368       elseif(inom.eq.287)then                       !'nd5evt'
2369         x=nsdiff(5,noweak(n))
2370       elseif(inom.eq.210)then                       !'theevt'
2371         eloevt=qsqevt/4./elepti+elepti*(1.-qsqevt/xbjevt/engy**2)
2372         x=acos(1-qsqevt/2./elepti/eloevt)/pi*180.
2373       elseif(inom.eq.211)then                       !'nspevt'
2374         x=0
2375         do i=1,nptl
2376          if((istptl(i).eq.30.or.istptl(i).eq.31)
2377      &      .and.int(idptl(i)/1000000).eq.1)x=x+1
2378         enddo
2379       elseif(inom.eq.212)then                       !'nhpevt'
2380         x=0
2381         do i=1,nptl
2382          if((istptl(i).eq.30.or.istptl(i).eq.31)
2383      &      .and.int(idptl(i)/1000000).eq.3)x=x+1
2384         enddo
2385       elseif(inom.eq.213)then                       !'sigtot'
2386         x=sigtot
2387       elseif(inom.eq.214)then                       !'sigela'
2388         x=sigela
2389       elseif(inom.eq.215)then                       !'sloela'
2390         x=sloela
2391       elseif(inom.eq.216)then                       !'nrgevt'
2392         x=0
2393         do i=1,nptl
2394           if(istptl(i).eq.31.and.int(idptl(i)/10000).eq.2)x=x+1
2395         enddo
2396       elseif(inom.eq.217)then                       !qevt
2397         x=sqrt(qsqevt)
2398       elseif(inom.eq.218)then   !qevt
2399         if(iappl.eq.8)then
2400           x=qtl
2401         else
2402           x=pprt(1,5)
2403         endif
2404       elseif(inom.eq.220)then!------------------------------------------
2405         x=sngl(bofra(1,lf))     !thrust
2406       elseif(inom.eq.221)then
2407         x=1.-sngl(bofra(1,lf))  !1-thrust
2408       elseif(inom.eq.222)then
2409         x=sngl(bofra(2,lf))     !major
2410       elseif(inom.eq.223)then
2411         x=sngl(bofra(3,lf))     !minor
2412       elseif(inom.eq.224)then
2413         x=sngl(bofra(2,lf)-bofra(3,lf)) !oblateness
2414       elseif(inom.eq.230)then!------------------------------------------
2415         x=1.5*(1.-sngl(bofra(1,lf))) !spherecity
2416       elseif(inom.eq.231)then
2417         x=1.5*sngl(bofra(3,lf)) !aplanarity
2418       elseif(inom.eq.232)then
2419         x=3.*sngl(bofra(1,lf)*bofra(2,lf)+bofra(1,lf)*bofra(3,lf)
2420      &       +bofra(2,lf)*bofra(3,lf)) !c-parameter
2421       elseif(inom.eq.233)then
2422         x=27.*sngl(bofra(1,lf)*bofra(2,lf)*bofra(3,lf))   !d-parameter
2423       elseif(inom.eq.234)then                       !npoevt
2424         x=0
2425         do i=1,nptl
2426          if(istptl(i).eq.30.or.istptl(i).eq.31)x=x+1
2427         enddo
2428       elseif(inom.eq.235)then                       !npnevt
2429         x=npjevt+ntgevt
2430
2431    !....unused 236, 237
2432
2433       elseif(inom.eq.238)then  !npxevt ... nr of pomerons, including absorbed
2434         x=0
2435         do i=1,nptl
2436          if(istptl(i).eq.30.or.istptl(i).eq.31)x=x+1
2437          if(mod(abs(idptl(i)),100).eq.94)x=x+0.5
2438         enddo
2439       elseif(inom.eq.240)then    !mu1evt ... charged ptl multipl for central rap
2440         x=multy1
2441       elseif(inom.eq.241)then    !muievt ... charged ptl multipl
2442         x=multyi
2443       elseif(inom.eq.242)then    !hgtevt ... charged ptl multipl for central eta
2444         x=multc05
2445       elseif(inom.eq.243)then                       !difevt
2446         npom=0
2447         do i=1,nptl
2448          if(istptl(i).eq.30.or.istptl(i).eq.31)npom=npom+1
2449         enddo
2450         x=0
2451         if(npom.eq.0)x=1
2452       elseif(inom.eq.244)then                       !dixevt
2453         zpom=0
2454         do i=1,nptl
2455          if(istptl(i).eq.30.or.istptl(i).eq.31)zpom=zpom+1
2456          if(mod(abs(idptl(i)),100).eq.94)zpom=zpom+0.5
2457         enddo
2458         x=0
2459         if(abs(zpom).lt.0.001)x=1
2460       elseif(inom.eq.250)then
2461         if(iappl.eq.8)then      !mass in
2462           x=-pptl(5,6)
2463         else
2464           x=pprt(5,3)
2465         endif
2466       elseif(inom.eq.251)then
2467         if(iappl.eq.8)then      !mass out
2468           x=pptl(5,7)
2469         else
2470           x=pprt(5,2)
2471         endif
2472       elseif(inom.eq.252)then
2473         if(iappl.eq.8)then
2474           x=-pptl(4,6)
2475         else
2476           x=pprt(4,2)
2477         endif
2478       elseif(inom.eq.253)then
2479         if(iappl.eq.8)then
2480           x=pptl(4,7)
2481         else
2482           x=pprt(4,3)
2483         endif
2484       elseif(inom.eq.254)then
2485         if(iappl.eq.8)then
2486           x=abs(pptl(3,6))
2487         else
2488           x=abs(pprt(3,2))
2489         endif
2490       elseif(inom.eq.255)then
2491         if(iappl.eq.8)then
2492           x=abs(pptl(3,7))
2493           do l=1,5
2494             p(l,lf)=pptl(l,7)
2495           enddo
2496           if(imofra(1,lf).ne.0)then
2497             call utrota(imofra(1,lf),r1fra(1,lf),r1fra(2,lf),r1fra(3,lf)
2498      $           ,p(1,lf),p(2,lf),p(3,lf))
2499           endif
2500           if(imofra(2,lf).ne.0)then !the x-z exchanged is ok !!
2501             call utrota(imofra(2,lf),r2fra(3,lf),r2fra(2,lf),r2fra(1,lf)
2502      $           ,p(3,lf),p(2,lf),p(1,lf))
2503           endif
2504           if(imofra(3,lf).ne.0)then
2505             call utlob4(imofra(3,lf),bofra(1,lf),bofra(2,lf)
2506      $           ,bofra(3,lf) ,bofra(4,lf),bofra(5,lf)
2507      $           ,p(1,lf),p(2,lf),p(3,lf),p(4,lf))
2508           endif
2509           x=abs(p(3,lf))
2510         else
2511           x=abs(pprt(3,3))
2512         endif
2513       elseif(inom.eq.256)then  !pxfevt: leading proton xf in cms
2514         x=-2
2515 c       pmax=sqrt((engy/2.)**2-prom**2)
2516           pmax=pnullx               !???????????????????
2517         do i=1,nptl
2518           if(idptl(i).eq.1120.and.istptl(i).eq.0)then
2519             if(iframe.eq.11)then
2520               pz=pptl(3,i)
2521             else
2522               amt=sqrt(prom**2+pptl(1,i)**2+pptl(2,i)**2)
2523               rap=alog((pptl(4,i)+pptl(3,i))/amt)
2524      &           -alog((sqrt(pnll**2+engy**2)+pnll)/engy)
2525               pz=amt*sinh(rap)
2526             endif
2527             x=max(x,pz/pmax)
2528           endif
2529         enddo
2530       elseif(inom.eq.257)then  !  pi+xf: pi+ yield at cms xf>0.01
2531         x=0.
2532 c        pmax=sqrt((engy/2)**2-prom*2)
2533           pmax=pnullx               !???????????????????
2534         do i=1,nptl
2535           if(idptl(i).eq.120.and.istptl(i).eq.0)then
2536             if(iframe.eq.11)then
2537               pz=pptl(3,i)
2538             else
2539               amt=sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2)
2540               rap=alog((pptl(4,i)+pptl(3,i))/amt)
2541      &           -alog((sqrt(pnll**2+engy**2)+pnll)/engy)
2542               pz=amt*sinh(rap)
2543             endif
2544             if(pz/pmax.gt.0.01)x=x+1.
2545           endif
2546         enddo
2547       elseif(inom.eq.258)then  !  pi-xf: pi- yield at cms xf>0.01
2548         x=0.
2549 c        pmax=sqrt((engy/2)**2-prom*2)
2550           pmax=pnullx               !???????????????????
2551         do i=1,nptl
2552           if(idptl(i).eq.-120.and.istptl(i).eq.0)then
2553             if(iframe.eq.11)then
2554               pz=pptl(3,i)
2555             else
2556               amt=sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2)
2557               rap=alog((pptl(4,i)+pptl(3,i))/amt)
2558      &           -alog((sqrt(pnll**2+engy**2)+pnll)/engy)
2559               pz=amt*sinh(rap)
2560             endif
2561             if(pz/pmax.gt.0.01)x=x+1.
2562           endif
2563         enddo
2564       elseif(inom.eq.260)then!------------------------------
2565         x=sigcut
2566       elseif(inom.eq.261)then
2567         x=keu
2568       elseif(inom.eq.262)then
2569         x=ked
2570       elseif(inom.eq.263)then
2571         x=kes
2572       elseif(inom.eq.265)then
2573         x=kolevt
2574       elseif(inom.eq.266)then
2575         x=sigsd
2576       elseif(inom.eq.267)then
2577         x=nglevt
2578       elseif(inom.eq.268)then  ! kppevt : collision number per participant
2579         x=kolevt/float(npjevt+ntgevt)
2580       elseif(inom.eq.269)then  ! npievt : pion + multi per event
2581         x=0
2582         do i=1,nptl
2583          if(idptl(i).eq.120)x=x+1
2584         enddo
2585       elseif(inom.eq.270)then  ! np2evt : pion + multi per event
2586         x=0
2587         do i=1,nptl
2588          if(idptl(i).eq.120)x=x+1
2589         enddo
2590         x=x/float(npjevt+ntgevt)
2591       elseif(inom.eq.271)then
2592         x=sigdif
2593       elseif(inom.eq.272)then  !number of inelastic collisions per event
2594         x=koievt
2595       elseif(inom.eq.273)then  ! inelasticity (energy loss of leading particle)
2596         x=0.
2597         do i=maproj+matarg+1,nptl
2598           if(istptl(i).eq.0)then
2599             if((((abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000)
2600      *           .and.idproj.gt.1000).or.(iabs(idptl(i)).gt.100
2601      *           .and.idproj.lt.1000)).and.pptl(4,i)
2602      *           .gt.x.and.pptl(3,i).gt.0.)x=pptl(4,i)
2603           endif
2604         enddo
2605         Eini=pptl(4,1)
2606         if(Eini.gt.0.)x=(Eini-x)/Eini
2607       elseif(inom.eq.274)then  ! elasticity (energy of leading particle)
2608         x=0.
2609         do i=maproj+matarg+1,nptl
2610           if(istptl(i).eq.0)then
2611             if((((abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000)
2612      *           .and.idproj.gt.1000).or.(iabs(idptl(i)).gt.100
2613      *           .and.idproj.lt.1000)).and.pptl(4,i)
2614      *           .gt.x.and.pptl(3,i).gt.0.)x=pptl(4,i)
2615           endif
2616         enddo
2617         Eini=pptl(4,1)
2618         if(Eini.gt.0.)x=x/Eini
2619       elseif(inom.eq.275)then         !'itgevt'
2620         x=0
2621         if(nint(ypara(1,n)).ne.0)x=1
2622       elseif(inom.eq.276)then         !'hrdevt' ......  1 = hard event
2623         x=0
2624         if(nint(ypara(1,n)).ne.0)x=1
2625       elseif(inom.eq.277)then         !'ej1evt' .... et of jet 1
2626         x=0
2627         if(nint(ypara(1,n)).ne.0)
2628      &  x=ypara(2,n)
2629       elseif(inom.eq.278)then         !'pj1evt' .... phi of jet 1
2630         x=1000
2631         if(nint(ypara(1,n)).ne.0)
2632      &  x=ypara(4,n)
2633       elseif(inom.eq.279)then         !'ej2evt' .... et of jet 2
2634         x=0
2635         if(nint(ypara(6,n)).ne.0)
2636      &  x=ypara(7,n)
2637       elseif(inom.eq.280)then         !'pj2evt' .... delta_phi of jet 2 1
2638         x=1000
2639         if(nint(ypara(6,n)).ne.0)then
2640           x=ypara(9,n)-ypara(4,n)
2641            if(x.lt.-2.5*pi)then
2642             x=x+4*pi
2643            elseif(x.lt.-0.5*pi)then
2644             x=x+2*pi
2645           elseif(x.gt.3.5*pi)then
2646             x=x-4*pi
2647           elseif(x.gt.1.5*pi)then
2648             x=x-2*pi
2649           endif
2650         endif
2651       elseif(inom.eq.281)then         !'zppevt'
2652         x=zppevt
2653       elseif(inom.eq.282)then         !'zptevt'
2654         x=zptevt
2655       elseif(inom.eq.283)then
2656         stop '**********not used*********'
2657       elseif(inom.eq.286)then         !'mubevt'
2658         x=multeb
2659       elseif(inom.eq.298)then
2660         x=sval(1,n)
2661       elseif(inom.eq.299)then
2662         x=sval(2,n)
2663       elseif(inom.eq.301)then   !---------------------------------------------
2664         if(j.eq.0)then          !initialize
2665           do l=1,4
2666             aimuni(l,n)=0.
2667           enddo
2668         elseif(j.gt.nptl)then   !final calculation
2669           am2=aimuni(4,n)**2-aimuni(3,n)**2
2670      $         -aimuni(2,n)**2-aimuni(1,n)**2
2671           x=sign(sqrt(abs(am2)),am2)
2672 c          print *, x
2673         else                    !routine work
2674           do l=1,4
2675             aimuni(l,n)=aimuni(l,n)+p(l,lf)
2676           enddo
2677 c          print *, j,(p(l,lf),l=1,5)
2678         endif
2679       elseif(inom.ge.302.and.inom.le.305)then   !-----------------------
2680         if(j.eq.0)then          !initialize
2681           do l=1,4
2682             aimuni(l,n)=0.
2683           enddo
2684         elseif(j.gt.nptl)then   !final calculation
2685           if(inom.eq.302) x=max(aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
2686      $         ,aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n)))
2687           if(inom.eq.303) x=min(aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
2688      $         ,aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n)))
2689           if(inom.eq.304) x=abs(aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
2690      $         -aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n)))
2691           if(inom.eq.305) x=aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
2692      $         +aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n))
2693         else                    !routine work
2694           l=0
2695           if(p(3,lf).lt.0.)l=2
2696           aimuni(1+l,n)=aimuni(1+l,n)+sqrt(p(1,lf)**2+p(2,lf)**2)
2697           aimuni(2+l,n)=aimuni(2+l,n)
2698      $         +sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
2699
2700         endif
2701       elseif(inom.eq.306.or.inom.eq.307.or.inom.eq.308)then !---------
2702         if(j.eq.0)then          !initialize
2703           do ll=1,8
2704             aimuni(ll,n)=0.
2705           enddo
2706         elseif(j.gt.nptl)then   !final calculation
2707           am2a=aimuni(4,n)**2-aimuni(3,n)**2
2708      $         -aimuni(2,n)**2-aimuni(1,n)**2
2709           am2b=aimuni(8,n)**2-aimuni(7,n)**2
2710      $         -aimuni(6,n)**2-aimuni(5,n)**2
2711           if(inom.eq.306)x=(max(0.,am2a,am2b))/engy**2
2712           if(inom.eq.307)x=(max(0.,min(am2a,am2b)))/engy**2
2713           if(inom.eq.308)x=(abs(am2a-am2b))/engy**2
2714         else                    !routine work
2715           ll=0
2716           if(p(3,lf).lt.0.)ll=4
2717           do l=1,4
2718             aimuni(l+ll,n)=aimuni(l+ll,n)+p(l,lf)
2719           enddo
2720         endif
2721       elseif (inom.eq.310.or.inom.eq.311) then !---------
2722         if(j.eq.0)then          !initialize
2723           aimuni(1,n)=0
2724           aimuni(2,n)=0
2725           do i=1,nptl
2726 c            charge=0.
2727              if(istptl(i).eq.0) then
2728                if (idptl(i).eq.idcod(1,n)) aimuni(1,n)=aimuni(1,n)+1.
2729                if (idptl(i).eq.idcod(2,n)) aimuni(2,n)=aimuni(2,n)+1.
2730              endif
2731            enddo
2732         elseif(j.gt.nptl)then   !final calculation
2733           if(aimuni(1,n).eq.0.or.aimuni(2,n).eq.0) then
2734             ncevt(n)=ncevt(n)-1
2735           endif
2736           x=xmin(n)-100.
2737           do i=1,nbin(n)
2738             zcbin(i,nac(n),n)=abs(zcbin(i,nac(n),n))
2739          enddo
2740        else                    !routine work
2741           if( istptl(j).eq.0
2742      $         .and. aimuni(1,n).ne.0. .and. aimuni(2,n).ne.0. ) then
2743             id1=idptl(j)
2744             if(id1.eq.idcod(1,n) .or. id1.eq.idcod(2,n)) then
2745               y1=sign(1.,pptl(3,j))*alog((pptl(4,j)+abs(pptl(3,j)))
2746      *             /sqrt(pptl(5,j)**2+pptl(1,j)**2+pptl(2,j)**2))
2747               do i=1,nptl
2748                 if(i.eq.j .or. istptl(i).ne.0) goto 124
2749                 id2=idptl(i)
2750                 if(id2.eq.idcod(1,n) .or. id2.eq.idcod(2,n)) then
2751                   y2=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
2752      *                 /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2))
2753                   dy=(y2-y1)
2754                   if(inom.eq.311) dy=abs(dy)
2755                   ib=1+int((dy-xmin(n))*xinc(n))
2756                   if(dy.ge.xmin(n).and.dy.le.xmax(n)) then
2757                     if( id1.eq.idcod(1,n) ) then
2758                       if( id2.eq.idcod(2,n) ) then
2759                         bin(ib,nac(n),n)=bin(ib,nac(n),n)+.5/aimuni(2,n)
2760                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)+1
2761                       else
2762                         bin(ib,nac(n),n)=bin(ib,nac(n),n)-.5/aimuni(1,n)
2763                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)-1
2764                       endif
2765                     else        !id1 is idcod(2,n)
2766                       if(id2.eq.idcod(1,n)) then
2767                         bin(ib,nac(n),n)=bin(ib,nac(n),n)+.5/aimuni(1,n)
2768                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)+1
2769                       else
2770                         bin(ib,nac(n),n)=bin(ib,nac(n),n)-.5/aimuni(2,n)
2771                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)-1
2772                       endif
2773                     endif
2774                   endif
2775                 endif
2776  124          enddo
2777             endif
2778           endif
2779         endif
2780       elseif (inom.eq.312) then !---------
2781         x=sigine
2782       elseif (inom.eq.313) then !---------
2783         x=sigineaa
2784       elseif (inom.eq.314) then !---------
2785         x=alpD(idxD0,iclpro,icltar)
2786       elseif (inom.eq.315) then !---------
2787         x=alpD(1,iclpro,icltar)
2788       elseif (inom.eq.316) then !---------
2789         x=betD(idxD0,iclpro,icltar)
2790         if(x.lt.0.)x=-10.*x
2791       elseif (inom.eq.317) then !---------
2792         x=betD(1,iclpro,icltar)
2793       elseif (inom.eq.318) then !---------
2794         x=rexdif(iclpro)
2795       elseif (inom.eq.319) then !---------
2796         x=rexdif(icltar)
2797       elseif(inom.eq.320)then    !m14evt ... multipl |eta|<1, pt>0.4
2798         x=multc14
2799       elseif(inom.eq.321)then    !ht3evt ... height |eta|<3.15
2800         x=multc3/6.3
2801       elseif (inom.eq.322) then !---------
2802         x=sigineex
2803       elseif (inom.eq.323) then !---------
2804         x=sigdifex
2805       elseif (inom.eq.324) then !---------
2806         x=sigsdex
2807       elseif (inom.eq.501) then !---------
2808         x=sval(1,1)
2809       elseif (inom.eq.502) then !---------
2810         x=sval(1,2)
2811       elseif (inom.eq.503) then !---------
2812         x=sval(1,3)
2813       elseif (inom.eq.504) then !---------
2814         x=sval(1,4)
2815       endif                     !---------------------------------------
2816       end
2817
2818 c----------------------------------------------------------------------
2819       subroutine StandardVariables
2820 c----------------------------------------------------------------------
2821       include 'epos.inc'
2822       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
2823      &     ,multc1
2824       common/crvar/idlead
2825       parameter(mxxhis=50)
2826       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
2827      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis)
2828
2829       Emax=0
2830       multy1=0
2831       multc05=0
2832       multc14=0
2833       multc1=0
2834       multyi=0
2835       multc3=0
2836       multeb=0
2837
2838
2839       do i=maproj+matarg+1,nptl
2840 c---multy1---charged ptl multipl for central rap
2841         if(istptl(i).eq.0)then
2842           amt=pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2
2843           pt=pptl(1,i)**2+pptl(2,i)**2
2844           pp=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2845           if(amt.gt.0..and.pptl(4,i).gt.0.)then
2846             amt=sqrt(amt)
2847             rap=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))/amt)
2848           else
2849             rap=1000.
2850           endif
2851           if(pt.gt.0.)then
2852             pt=sqrt(pt)
2853             eta=sign(1.,pptl(3,i))*alog((pp+abs(pptl(3,i)))/pt)
2854           else
2855             eta=1000.
2856           endif
2857           if(abs(idptl(i)).ge.100
2858      $         .and.abs(idptl(i)).lt.10000)then
2859             call idchrg(idptl(i),ch)
2860             if(abs(ch).gt.0.1)then
2861 c---multyi---charged ptl multipl
2862               multyi=multyi+1
2863 c---multy1---charged ptl multipl for central rap
2864               if(abs(rap).le.1.)multy1=multy1+1
2865               if(abs(eta).le.0.5)multc05=multc05+1
2866               if(abs(eta).le.1.and.pt.gt.0.4)multc14=multc14+1
2867               if(abs(eta).le.1)multc1=multc1+1
2868               if(abs(rap).le.3.15)multc3=multc3+1
2869 c---multeb---charged ptl multipl for back rap
2870               if(eta.gt.-3.8.and.eta.lt.-2.8)multeb=multeb+1
2871             endif
2872           endif
2873           if((((abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000)
2874      *         .and.idproj.gt.1000).or.(iabs(idptl(i)).gt.100
2875      *         .and.idproj.lt.1000)).and.pptl(4,i)
2876      *         .gt.Emax.and.pptl(3,i).gt.0.)then
2877             Emax=pptl(4,i)
2878             idlead=i
2879           endif
2880         endif
2881       enddo
2882       end
2883
2884 c----------------------------------------------------------------------
2885       subroutine jetfind(m,n)
2886 c----------------------------------------------------------------------
2887 c   m = 1 ou 2 (two different definitions)
2888 c   n = histogram
2889 c input(jet definition):
2890 c   xpara(1,n) ... 1=used (m=1)
2891 c   xpara(2,n) ... etamin (m=1)
2892 c   xpara(3,n) ... etamax (m=1)
2893 c   xpara(4,n) ... rmax   (m=1) (rmax defining the cone)
2894 c   xpara(5,n) ... ichd   (m=1) (1=charged, 0=all)
2895 c   xpara(6,n) ... 1=used (m=2)
2896 c   xpara(7,n) ... etamin (m=2)
2897 c   xpara(8,n) ... etamax (m=2)
2898 c   xpara(9,n) ... rmax   (m=2)
2899 c   xpara(10,n) .. ichd   (m=2)
2900 c output (jet properties):
2901 c   ypara(1,n) ... 1 (found) or 0 if not  (m=1)
2902 c   ypara(2,n) ... et                     (m=1)
2903 c   ypara(3,n) ... eta of center          (m=1)
2904 c   ypara(4,n) ... phi of center          (m=1)
2905 c   ypara(5,n)
2906 c   ypara(6,n) ... 1 (found) or 0 if not  (m=2)
2907 c   ypara(7,n) ... et                     (m=2)
2908 c   ypara(8,n) ... eta of center          (m=2)
2909 c   ypara(9,n) ... phi of center          (m=2)
2910 c   ypara(10,n)
2911 c----------------------------------------------------------------------
2912
2913       include 'epos.inc'
2914       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
2915       parameter (mypara=10,mxpara=10)
2916       parameter (mxval=5)
2917       real ptx(mxval),lst(mxval),etax(mxval),phix(mxval)
2918       logical ilog,icnx,itrevt,idmod
2919       double precision bin,bbin,zcbin,zbbin
2920       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
2921      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
2922      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
2923      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
2924      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
2925      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
2926      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
2927      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
2928      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
2929      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
2930      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
2931
2932       if(m.ne.1.and.m.ne.2)stop'jetfind: value of m not valid.      '
2933
2934       etamin=           xpara(2+5*(m-1),n)
2935       etamax=           xpara(3+5*(m-1),n)
2936       rmax  =           xpara(4+5*(m-1),n)
2937       ichd  = nint(xpara(5+5*(m-1),n))
2938
2939       ifound=0
2940       do l=1,mxval
2941         ptx(l)=0
2942         lst(l)=0
2943         etax(l)=0
2944         phix(l)=0
2945       enddo
2946
2947 ctp060829      pp1=0
2948 ctp060829      pp2=0
2949 ctp060829      pp3=0
2950
2951       do i=maproj+matarg+1,nptl
2952         iok=0
2953         if(istptl(i).eq.0.and.abs(idptl(i)).lt.10000)iok=1
2954         if(iok.eq.1)call idchrg(idptl(i),ch)
2955         if(ichd.eq.1.and.nint(ch).eq.0)iok=0
2956         if(iok.eq.1)then
2957           p1=pptl(1,i)
2958           p2=pptl(2,i)
2959           p3=pptl(3,i)
2960           pt=sqrt(p1**2+p2**2)
2961                 if(pt.gt.0)then
2962             eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
2963             phi=sign(1.,p2)*acos(p1/pt)
2964           else
2965             eta=10000
2966             phi=0
2967           endif
2968           do k=1,mxval
2969             iok=1
2970             if(m.eq.2)then
2971               dphi=phi-ypara(4,n)
2972               if(dphi.lt.-pi)dphi=dphi+2*pi
2973               if(dphi.gt. pi)dphi=dphi-2*pi
2974               if(abs(dphi).lt.pi/2)iok=0
2975             endif
2976             if(iok.eq.1.and.pt.gt.ptx(k)
2977      &        .and.eta.le.etamax.and.eta.ge.etamin)then
2978               do l=mxval,k+1,-1
2979                ptx(l)=ptx(l-1)
2980                lst(l)=lst(l-1)
2981                etax(l)=etax(l-1)
2982                phix(l)=phix(l-1)
2983               enddo
2984                ptx(k)=pt
2985                lst(k)=i
2986                etax(k)=eta
2987                phix(k)=phi
2988               goto2
2989             endif
2990           enddo
2991   2       continue
2992         endif
2993       enddo
2994
2995       kk=0
2996       etx=0
2997
2998       do k=1,mxval
2999        if(lst(k).ne.0)then
3000
3001         ifound=1
3002         et=0
3003         etaxx=etax(k)
3004         phixx=phix(k)
3005         do j=maproj+matarg+1,nptl
3006           iok=0
3007           if(istptl(j).eq.0.and.abs(idptl(j)).lt.10000)iok=1
3008           if(iok.eq.1)call idchrg(idptl(j),ch)
3009           if(ichd.eq.1.and.nint(ch).eq.0)iok=0
3010           if(iok.eq.1)then
3011             p1=pptl(1,j)
3012             p2=pptl(2,j)
3013             p3=pptl(3,j)
3014             pt=sqrt(p1**2+p2**2)
3015             am=pptl(5,j)
3016                   if(pt.gt.0)then
3017               eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
3018               phi=sign(1.,p2)*acos(p1/pt)
3019             else
3020               eta=-10000
3021               phi=0
3022             endif
3023             if(eta.le.etamax.and.eta.ge.etamin)then
3024               deta=eta-etaxx
3025               dphi=phi-phixx
3026               if(dphi.lt.-pi)dphi=dphi+2*pi
3027               if(dphi.gt. pi)dphi=dphi-2*pi
3028               if(deta**2+dphi**2.lt.rmax**2)et=et+sqrt(pt**2+am**2)
3029             endif
3030           endif
3031         enddo
3032         if(et.gt.etx)then
3033           etx=et
3034           kk=k
3035         endif
3036
3037        endif
3038       enddo
3039
3040       ypara(1+5*(m-1),n)=ifound
3041       ypara(2+5*(m-1),n)=etx
3042       if(kk.gt.0)then
3043        ypara(3+5*(m-1),n)=etax(kk)
3044        ypara(4+5*(m-1),n)=phix(kk)
3045       endif
3046       return
3047       end
3048
3049
3050 c----------------------------------------------------------------------
3051       subroutine hardevent(n)
3052 c----------------------------------------------------------------------
3053 c   n = histogram
3054 c input(jet event conditions):
3055 c   xpara(2,n) ... pt1
3056 c   xpara(3,n) ... pt2
3057 c   xpara(4,n) ... absetamax
3058 c   xpara(5,n) ... rmax        (r=sqrt(deltaeta**2+deltaphi**2))
3059 c   xpara(6,n) ... Et_min
3060 c output (jet event found or not):
3061 c   ypara(1,n) ... 1 (found) or 0 if not
3062 c----------------------------------------------------------------------
3063
3064       include 'epos.inc'
3065       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
3066       parameter (mypara=10,mxpara=10)
3067       logical ilog,icnx,itrevt,idmod
3068       double precision bin,bbin,zcbin,zbbin
3069       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
3070      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
3071      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
3072      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
3073      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
3074      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
3075      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
3076      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
3077      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
3078      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
3079      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
3080
3081       ypara(1,n)=0
3082       do i=maproj+matarg+1,nptl
3083        if(abs(idptl(i)).ge.100.and.abs(idptl(i)).lt.10000.
3084      $  and.istptl(i).eq.0)then
3085         call idchrg(idptl(i),ch)
3086         if(abs(ch).gt.0.1)then
3087           p1=pptl(1,i)
3088           p2=pptl(2,i)
3089           p3=pptl(3,i)
3090           pt=sqrt(p1**2+p2**2)
3091           if(pt.gt.0)then
3092             eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
3093             phi=sign(1.,p2)*acos(p1/pt)
3094           else
3095             eta=10000
3096             phi=0
3097           endif
3098           if(pt.ge.xpara(2,n).and.abs(eta).lt.xpara(4,n))then
3099             et1=pptl(4,i)*pt/sqrt(p3**2+pt**2)
3100             do j=maproj+matarg+1,nptl
3101               if(j.ne.i
3102      $        .and.abs(idptl(j)).ge.100.and.abs(idptl(j)).lt.10000.
3103      $        .and.istptl(j).eq.0)then
3104                 call idchrg(idptl(j),ch)
3105                 if(abs(ch).gt.0.1.and.abs(idptl(j)).ge.100
3106      $          .and.abs(idptl(j)).lt.10000.and.istptl(j).eq.0)then
3107                   p1=pptl(1,j)
3108                   p2=pptl(2,j)
3109                   p3=pptl(3,j)
3110                   pt=sqrt(p1**2+p2**2)
3111                         if(pt.gt.0)then
3112                    etax=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
3113                    phix=sign(1.,p2)*acos(p1/pt)
3114                   else
3115                     etax=-10000
3116                     phix=0
3117                   endif
3118                   if(pt.ge.xpara(3,n).and.abs(etax).lt.xpara(4,n))then
3119                     deta=eta-etax
3120                     dphi=phi-phix
3121                     if(dphi.lt.-pi)dphi=dphi+2*pi
3122                     if(dphi.gt. pi)dphi=dphi-2*pi
3123                     if(deta**2+dphi**2.lt.xpara(5,n)**2)then
3124                     et2=pptl(4,j)*pt/sqrt(p3**2+pt**2)
3125                      if(et1+et2.gt.xpara(6,n))then
3126                       ypara(1,n)=1
3127                       goto1
3128                      endif
3129                     endif
3130                   endif
3131                 endif
3132               endif
3133             enddo
3134           endif
3135         endif
3136        endif
3137       enddo
3138
3139    1  continue
3140       return
3141       end
3142
3143 c----------------------------------------------------------------------
3144       subroutine corrtrig(n)
3145 c----------------------------------------------------------------------
3146 c   n = histogram
3147 c input(trigger conditions):
3148 c   xpara(2,n) ... ptmin
3149 c   xpara(3,n) ... ptmax
3150 c   xpara(4,n) ... etamin
3151 c   xpara(5,n) ... etamax
3152 c   xpara(6,n) ...
3153 c   xpara(7,n) ...
3154 c output (triggered particle (the one with highest pt if there are several)):
3155 c   ypara(1,n) ... iptl or 0 if no particle found
3156 c   ypara(2,n) ... phi of particle
3157 c----------------------------------------------------------------------
3158
3159       include 'epos.inc'
3160       parameter (mxhis=500,mxcontr=300,mxidcd=60,mxtri=20,mxbin=400)
3161       parameter (mypara=10,mxpara=10)
3162       logical ilog,icnx,itrevt,idmod
3163       double precision bin,bbin,zcbin,zbbin
3164       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
3165      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
3166      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
3167      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
3168      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
3169      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
3170      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
3171      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
3172      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
3173      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
3174      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
3175
3176       pt0=xpara(2,n)
3177
3178       do i=maproj+matarg+1,nptl
3179        if(abs(idptl(i)).ge.100.and.abs(idptl(i)).lt.10000.
3180      $  and.istptl(i).eq.0)then
3181         call idchrg(idptl(i),ch)
3182         if(abs(ch).gt.0.1)then
3183           p1=pptl(1,i)
3184           p2=pptl(2,i)
3185           p3=pptl(3,i)
3186           pt=sqrt(p1**2+p2**2)
3187           pt=max(pt,1e-20)
3188           eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
3189           phi=sign(1.,p2)*acos(p1/pt)
3190           if(pt.ge.pt0.and.pt.le.xpara(3,n).and.eta.gt.xpara(4,n)
3191      $      .and.eta.lt.xpara(5,n))then
3192             pt0=pt
3193             ypara(1,n)=i
3194             ypara(2,n)=phi
3195           endif
3196         endif
3197        endif
3198       enddo
3199
3200       return
3201       end
3202
3203 c----------------------------------------------------------------------