Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:29

0001 c---------------------------------------------------------------------
0002       subroutine xiniall
0003 c---------------------------------------------------------------------
0004       include 'epos.inc'
0005       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
0006       parameter (mypara=10,mxpara=10)
0007       logical ilog,icnx,itrevt,idmod
0008       double precision bin,bbin,zcbin,zbbin
0009       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
0010      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
0011      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
0012      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
0013      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
0014      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
0015      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
0016      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
0017      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
0018      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
0019      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
0020       double precision ebin,zebin
0021       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
0022      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
0023 
0024       parameter (mxfra=5)
0025       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
0026      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
0027      $     ,emax(mxfra)
0028       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
0029      &     ,multc1,multc83,multc24,multc25,rapgap,ipairs1,xsi
0030       parameter(mxxhis=70)
0031       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
0032      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis),imux(0:mxxhis)
0033      &,ifastjet(0:mxxhis),ijetevent(0:mxxhis),icaltrig(0:mxxhis)
0034 
0035       nhis=0
0036       nfra=0
0037       imulty1=0
0038       ipairs1=0
0039       do n=0,mxxhis
0040         icorrtrig(n)=0
0041         ihardevent(n)=0
0042         ijetfind1(n)=0
0043         ijetfind2(n)=0
0044         ifastjet(n)=0
0045         ijetevent(n)=0
0046         imux(n)=0
0047         icaltrig(n)=0
0048       enddo
0049       do n=1,mxhis
0050         do m=1,mxpara
0051           xpara(m,n)=0
0052         enddo
0053       enddo
0054       ncontrall=0
0055       noerrall=0
0056 
0057       end
0058 
0059 c---------------------------------------------------------------------
0060       subroutine xini
0061 c---------------------------------------------------------------------
0062 c  called after beginhisto
0063 c---------------------------------------------------------------------
0064       include 'epos.inc'
0065       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
0066       parameter (mypara=10,mxpara=10)
0067       logical ilog,icnx,itrevt,idmod
0068       double precision bin,bbin,zcbin,zbbin
0069       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
0070      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
0071      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
0072      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
0073      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
0074      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
0075      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
0076      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
0077      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
0078      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
0079      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
0080       double precision ebin,zebin
0081       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
0082      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
0083       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
0084      &     ,multc1,multc83,multc24,multc25,rapgap,ipairs1,xsi
0085       parameter(mxxhis=70)
0086       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
0087      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis),imux(0:mxxhis)
0088      &,ifastjet(0:mxxhis),ijetevent(0:mxxhis),icaltrig(0:mxxhis)
0089 
0090       parameter (mxfra=5)
0091       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
0092      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
0093      $     ,emax(mxfra)
0094       character line*1000,cvar*6
0095       logical go
0096       common/nl/noplin  /cnnnhis/nnnhis
0097       character*1000 cline
0098       common/cjjj/jjj,cline
0099 
0100       call utpri('xini  ',ish,ishini,5)
0101 
0102       i=1
0103                                 !      iapl=0
0104                                 !      nhis=0
0105       j=jjj     !-1
0106       line=cline
0107                                 !      nfra=1
0108                                 !      ifra(1)=iframe
0109       iapl=0
0110       if(nfra.eq.0)then
0111         nfra=1
0112         ifra(1)=iframe
0113       endif
0114       nhis=nhis+1
0115       nnnhis=nhis
0116       if(nhis.gt.mxhis)stop'xini: mxhis too small.       '
0117       noweak(nhis)=0
0118       ionoerr=0
0119 c      newfra=0
0120       indfra=1
0121 c      nepfra=0
0122       inpfra=1
0123  1    call utword(line,i,j,0)
0124       if(line(i:j).eq.'application')then !-----------
0125         call utword(line,i,j,1)
0126         if(line(i:j).eq.'analysis')then
0127                                 !iapl=0
0128                                 !nhis=nhis+1
0129                                 !newfra=0
0130                                 !indfra=1
0131                                 !nepfra=0
0132                                 !inpfra=1
0133         else
0134           iapl=1
0135         endif
0136       elseif(line(i:j).eq.'input')then !-----------
0137         call utword(line,i,j,0)
0138         if(nopen.ge.0)then
0139          nopen=nopen+1
0140          if(nopen.gt.9)stop'too many nested input commands'
0141          open(unit=20+nopen,file=line(i:j),status='old')
0142          if(iprmpt.eq.1)iprmpt=-1
0143         endif
0144       elseif(line(i:j).eq.'runprogram')then !-----------
0145         if(iapl.eq.0)then
0146         else
0147           goto 9999
0148         endif
0149       elseif(line(i:j).eq.'frame'.or.line(i:j).eq.'frame+')then !------
0150         ifp=1
0151         if(line(i:j).eq.'frame+')ifp=2
0152         call utword(line,i,j,1)
0153         if(line(i:j).eq.'total')then
0154           nfp=iframe
0155         elseif(line(i:j).eq.'nucleon-nucleon')then
0156           nfp=11
0157         elseif(line(i:j).eq.'target')then
0158           nfp=12
0159         elseif(line(i:j).eq.'gamma-nucleon')then
0160           nfp=21
0161         elseif(line(i:j).eq.'lab')then
0162           nfp=22
0163         elseif(line(i:j).eq.'breit')then
0164           nfp=23
0165         elseif(line(i:j).eq.'thrust')then
0166           nfp=33
0167         elseif(line(i:j).eq.'sphericity')then
0168           nfp=32
0169         else
0170           nfp=0
0171           call utstop("Wrong frame in xini !&")
0172         endif
0173         go=.true.
0174         inl=0
0175         do l=1,nfra
0176           if(ifra(l).eq.nfp)then
0177             inl=l
0178             go=.false.
0179           endif
0180         enddo
0181         if (go) then
0182           nfra=nfra+1
0183           inl=nfra
0184           ifra(nfra)=nfp
0185         endif
0186         if(ifp.eq.1)then
0187           indfra=inl
0188 c          newfra=nfp
0189           ivfra(1,nhis)=indfra
0190           ivfra(2,nhis)=indfra
0191         else
0192           inpfra=inl
0193 c          nepfra=nfp
0194         endif
0195       elseif(line(i:j).eq.'binning')then !-----------
0196         call utword(line,i,j,1)
0197         if(line(i:j).eq.'lin')then
0198           iologb=0
0199           iocnxb=0
0200         elseif(line(i:j).eq.'log')then
0201           iologb=1
0202           iocnxb=0
0203         elseif(line(i:j).eq.'clin')then
0204           iologb=0
0205           iocnxb=1
0206         elseif(line(i:j).eq.'clog')then
0207           iologb=1
0208           iocnxb=1
0209         else
0210           print *, 'what the heck is ',line(i:j),' binning?'
0211           print *, 'I will use the linear (lin) one'
0212         endif
0213       elseif(line(i:j).eq.'setm')then !-----------
0214         if(iapl.eq.0) then
0215           print *,"You should use histogram instead of setm, please"
0216           stop
0217         endif
0218       elseif(line(i:j).eq.'set')then !-----------
0219         call utword(line,i,j,1)
0220         if(line(i:j).eq.'iologb')then
0221           call utword(line,i,j,1)
0222           read(line(i:j),*) iologb
0223         elseif(line(i:j).eq.'iocnxb')then
0224           call utword(line,i,j,1)
0225           read(line(i:j),*) iocnxb
0226         elseif(line(i:j).eq.'etacut')then
0227           call utword(line,i,j,1)
0228           read(line(i:j),*) etacut
0229         elseif(line(i:j).eq.'nemsi')then
0230           call utword(line,i,j,1)
0231           read(line(i:j),*)nemsi
0232         endif
0233       elseif(line(i:j).eq.'xpara')then !-----------
0234         call utword(line,i,j,1)
0235         read(line(i:j),*)ipara
0236         if(ipara.gt.mxpara)stop'mxpara too small.         '
0237         call utword(line,i,j,1)
0238         read(line(i:j),*)val
0239         xpara(ipara,nhis)=val
0240       elseif(line(i:j).eq.'xparas')then !-----------
0241         call utword(line,i,j,1)
0242         read(line(i:j),*)ipara
0243         if(ipara.gt.mxpara)stop'mxpara too small.'
0244         do ii=1,ipara
0245         call utword(line,i,j,1)
0246         read(line(i:j),*)val
0247         xpara(ii,nhis)=val
0248         enddo
0249       elseif(line(i:j).eq.'echo')then !-----------
0250         call utword(line,i,j,1)
0251         if(line(i:j).eq.'on')iecho=1
0252         if(line(i:j).eq.'off')iecho=0
0253         if(line(i:j).ne.'on'.and.line(i:j).ne.'off')
0254      *  stop'invalid option'
0255       elseif(line(i:j).eq.'noweak')then !-----------
0256         noweak(nhis)=1
0257       elseif(line(i:j).eq.'histogram'
0258      *       .or.line(i:j).eq.'hi')then !-----------
0259         nac(nhis)=1
0260         call utword(line,i,j,1) !xvaria
0261         cvar='      '
0262         cvar=line(i:j)
0263         call xtrans(cvar,inom,ifrnew,nhis)
0264         if(inom.eq.-1)then
0265           if(line(i:i).ge.'0'.and.line(i:i).le.'9')then
0266             inom=298
0267             read(line(i:j),*) sval(1,nhis)
0268           endif
0269         endif
0270         ivar(1,nhis)=inom
0271         if(ifrnew.ne.0)then     !check frame for e+e- event
0272           go=.true.             !shape variables
0273           do l=1,nfra
0274             if(ifra(l).eq.ifrnew)then
0275               indfra=l
0276               go=.false.        !have it already
0277             endif
0278           enddo
0279           if (go) then
0280             nfra=nfra+1
0281             ifra(nfra)=ifrnew
0282             indfra=nfra
0283           endif
0284         endif
0285         call utword(line,i,j,1) !yvaria
0286         cvar='      '
0287         cvar=line(i:j)
0288         call xtrans(cvar,inom,ifrnew,nhis)
0289         ivar(2,nhis)=inom
0290         if(inom.eq.-1)then
0291           if(line(i:i).ge.'0'.and.line(i:i).le.'9')then
0292             inom=299
0293             read(line(i:j),*) sval(2,nhis)
0294           endif
0295         endif
0296         if(inom.eq.-1)ivar(1,nhis)=inom
0297 
0298         ivfra(1,nhis)=indfra
0299         ivfra(2,nhis)=indfra
0300 
0301         call utword(line,i,j,1) !normation
0302         read(line(i:j),*) inorm(nhis)
0303 
0304         call utword(line,i,j,1) !xmin
0305         if(line(i:j).eq.'egy')then
0306          if(engy.gt.0)then
0307           egy=engy
0308          elseif(ecms.gt.0.)then
0309           egy=ecms
0310          elseif(elab.gt.0)then
0311           call idmass(idproj,apj)
0312           call idmass(idtarg,atg)
0313           egy=sqrt( 2*elab*atg+atg**2+apj**2 )
0314          elseif(ekin.gt.0.)then
0315           call idmass(idproj,apj)
0316           call idmass(idtarg,atg)
0317           egy=sqrt( 2*(ekin+apj)*atg+atg**2+apj**2 )
0318          elseif(pnll.gt.0.)then
0319           call idmass(idproj,apj)
0320           call idmass(idtarg,atg)
0321           egy=sqrt( 2*sqrt(pnll**2+apj**2)*atg+atg**2+apj**2 )
0322          else
0323           stop'pb in xini (1).   '
0324          endif
0325          xmin(nhis)=egy-0.5
0326         else
0327          read(line(i:j),*) xmin(nhis)
0328         endif
0329 
0330         call utword(line,i,j,1) !xmax
0331         if(line(i:j).eq.'egy')then
0332          if(engy.gt.0)then
0333           egy=engy
0334          elseif(ecms.gt.0.)then
0335           egy=ecms
0336          elseif(elab.gt.0)then
0337           call idmass(idproj,apj)
0338           call idmass(idtarg,atg)
0339           egy=sqrt( 2*elab*atg+atg**2+apj**2 )
0340          elseif(ekin.gt.0.)then
0341           call idmass(idproj,apj)
0342           call idmass(idtarg,atg)
0343           egy=sqrt( 2*(ekin+apj)*atg+atg**2+apj**2 )
0344          elseif(pnll.gt.0.)then
0345           call idmass(idproj,apj)
0346           call idmass(idtarg,atg)
0347           egy=sqrt( 2*sqrt(pnll**2+apj**2)*atg+atg**2+apj**2 )
0348          else
0349           stop'pb in xini (2).   '
0350          endif
0351          xmax(nhis)=egy+0.5
0352         else
0353          read(line(i:j),*) xmax(nhis)
0354         endif
0355 
0356         call utword(line,i,j,1) !nbin
0357         read(line(i:j),*) nbin(nhis)
0358         do l=1,nbin(nhis)
0359           bin(l,nac(nhis),nhis)=0.
0360           zcbin(l,nac(nhis),nhis)=0
0361         enddo
0362         lookcontr(nhis)=0
0363         lookcontrx(nhis)=0
0364         inoerr(nhis)=0
0365       elseif(line(i:j).eq.'idcode')then !-----------
0366         call utword(line,i,j,1) !idcode
0367         if(line(i:i+2).eq.'995')stop'xini: idcode 995 not supported'
0368         if(line(i:i+2).eq.'994')stop'xini: idcode 994 not supported'
0369         nidcod(nhis)=nidcod(nhis)+1
0370         read(line(i:j),*) idcod(nidcod(nhis),nhis)
0371         idmod(nidcod(nhis),nhis)=.false.
0372       elseif(line(i:j).eq.'idcode+')then !-----------
0373         stop'xini: idcode+ not supported'
0374         call utword(line,i,j,1) !idcode
0375         if(line(i:i+2).eq.'995')stop'xini: idcode 995 not supported'
0376         if(line(i:i+2).eq.'994')stop'xini: idcode 994 not supported'
0377         nidcod(nhis)=nidcod(nhis)+1
0378         read(line(i:j),*) idcod(nidcod(nhis),nhis)
0379         idmod(nidcod(nhis),nhis)=.true.
0380       elseif(line(i:j).eq.'trigger')then !-----------
0381         call utword(line,i,j,1)
0382         ntc=1
0383         imo=1
0384         ncontr=0
0385         icontrtyp(nhis)=0
0386         if(line(i:j).eq.'or'.or.line(i:j).eq.'contr')then
0387           imo=2
0388           if(line(i:j).eq.'contr')imo=3
0389           call utword(line,i,j,1)
0390           read(line(i:j),*)ztc
0391           ntc=nint(ztc)
0392           call utword(line,i,j,1)
0393           if(imo.eq.3)then
0394             ncontr=ntc
0395             ncontrall=ncontrall+ncontr
0396             if(ncontrall.gt.mxcontr)stop'xini: mxcontr too small.     '
0397             if(ncontr.gt.mxcnt)stop'xini: mxcnt too small.     '
0398             lookcontr(nhis)=ncontrall-ncontr+1
0399             lookcontrx(nhis)=ncontrall
0400             do nb=1,nbin(nhis)
0401               do nn=1,ncontr
0402                 bbin(nb,nac(nhis),lookcontr(nhis)-1+nn)=0.d0
0403                 zbbin(nb,nac(nhis),lookcontr(nhis)-1+nn)=0.d0
0404               enddo
0405             enddo
0406             do nn=1,ncontr
0407                     nccevt(lookcontr(nhis)-1+nn)=0
0408             enddo
0409           endif
0410         endif
0411         do n=1,ntc
0412           if(n.ne.1)call utword(line,i,j,1) !trigger-name
0413           cvar='      '
0414           ifp=1
0415           if(line(j:j).eq.'+')then
0416             cvar=line(i:j-1)
0417             ifp=2
0418           else
0419             cvar=line(i:j)
0420             ifp=1
0421           endif
0422           call xtrans(cvar,inom,ifrnew,nhis)
0423           if(inom.gt.0)then
0424             ntri(nhis)=ntri(nhis)+1
0425             if(ntc.eq.1)then
0426               ntrc(ntri(nhis),nhis)=1
0427             elseif(n.eq.1)then
0428               ntrc(ntri(nhis),nhis)=2
0429             elseif(n.eq.ntc)then
0430               ntrc(ntri(nhis),nhis)=3
0431             else
0432               ntrc(ntri(nhis),nhis)=0
0433             endif
0434             if(imo.eq.3)then
0435               ntrc(ntri(nhis),nhis)=-1
0436               if(n.eq.1)then
0437                 icontrtyp(nhis)=1+inom/100
0438               else
0439                 if(1+inom/100.ne.icontrtyp(nhis))
0440      *               stop'xini: type mismatch'
0441               endif
0442             endif
0443             itri(ntri(nhis),nhis)=inom
0444             if(ifp.eq.1)then
0445               itfra(ntri(nhis),nhis)=indfra
0446             else
0447               itfra(ntri(nhis),nhis)=inpfra
0448             endif
0449             xmitrp(ntri(nhis),nhis)=100.
0450             xmatrp(ntri(nhis),nhis)=100.
0451             call utword(line,i,j,1) !-----------xmin----------
0452             if(line(i:j).eq.'inf')then
0453               xmitri(ntri(nhis),nhis)=1e30
0454             elseif(line(i:j).eq.'-inf')then
0455               xmitri(ntri(nhis),nhis)=-1e30
0456             elseif(line(i:j).eq.'A')then
0457               xmitri(ntri(nhis),nhis)=maproj
0458             elseif(line(i:j).eq.'A+1')then
0459               xmitri(ntri(nhis),nhis)=maproj+1
0460             elseif(line(i:j).eq.'A+B')then
0461               xmitri(ntri(nhis),nhis)=maproj+matarg
0462             elseif(line(i:j).eq.'A+B+1')then
0463               xmitri(ntri(nhis),nhis)=maproj+matarg+1
0464             elseif(line(i:j).eq.'lead')then    !leading particle (neads Standard Variable)
0465               xmitri(ntri(nhis),nhis)=-123456
0466               imulty1=1
0467             elseif(line(i:j).eq.'jet')then    !jet from fastjet
0468               xmitri(ntri(nhis),nhis)=nhis*100
0469               iok=0
0470               do i=1,ifastjet(0)
0471                 if(ifastjet(i).eq.nhis)iok=1
0472               enddo
0473               if(iok.eq.0)then
0474                 ifastjet(0)=ifastjet(0)+1
0475                 if(ifastjet(0).gt.mxxhis)stop'mxxhis too small'
0476                 ifastjet(ifastjet(0))=nhis
0477               endif
0478             else
0479               kk=0
0480               do k=i+1,j-1
0481                 if(line(k:k).eq.'%')kk=k
0482               enddo
0483               if(kk.eq.0)then
0484                 read(line(i:j),*)xmitri(ntri(nhis),nhis)
0485               else
0486                 read(line(i:kk-1),*)xmitrp(ntri(nhis),nhis)
0487                 read(line(kk+1:j),*)xmitri(ntri(nhis),nhis)
0488               endif
0489             endif
0490             call utword(line,i,j,1) !-----------xmax------------
0491             if(line(i:j).eq.'inf')then
0492               xmatri(ntri(nhis),nhis)=1e30
0493             elseif(line(i:j).eq.'-inf')then
0494               xmatri(ntri(nhis),nhis)=-1e30
0495             elseif(line(i:j).eq.'A')then
0496               xmatri(ntri(nhis),nhis)=maproj
0497             elseif(line(i:j).eq.'A+1')then
0498               xmatri(ntri(nhis),nhis)=maproj+1
0499             elseif(line(i:j).eq.'A+B')then
0500               xmatri(ntri(nhis),nhis)=maproj+matarg
0501             elseif(line(i:j).eq.'A+B+1')then
0502               xmatri(ntri(nhis),nhis)=maproj+matarg+1
0503             elseif(line(i:j).eq.'lead')then    !leading particle (neads Standard Variable)
0504               xmatri(ntri(nhis),nhis)=-123456
0505               imulty1=1
0506             elseif(line(i:j).eq.'jet')then    !jet form fastjet
0507               xmatri(ntri(nhis),nhis)=nhis*100
0508             else
0509               kk=0
0510               do k=i+1,j-1
0511                 if(line(k:k).eq.'%')kk=k
0512               enddo
0513               if(kk.eq.0)then
0514                 read(line(i:j),*)xmatri(ntri(nhis),nhis)
0515                 xmatrp(ntri(nhis),nhis)=100.
0516               else
0517                 read(line(i:kk-1),*)xmatrp(ntri(nhis),nhis)
0518                 read(line(kk+1:j),*)xmatri(ntri(nhis),nhis)
0519               endif
0520             endif
0521             !---exchange min-max------------------
0522             if(xmitri(ntri(nhis),nhis).gt.xmatri(ntri(nhis),nhis))then
0523               xmatri_save=xmatri(ntri(nhis),nhis)
0524               xmatrp_save=xmatrp(ntri(nhis),nhis)
0525               xmatri(ntri(nhis),nhis)=xmitri(ntri(nhis),nhis)
0526               xmatrp(ntri(nhis),nhis)=xmitrp(ntri(nhis),nhis)
0527               xmitri(ntri(nhis),nhis)=xmatri_save
0528               xmitrp(ntri(nhis),nhis)=xmatrp_save
0529             endif
0530             !-------------------------------------
0531           else
0532             ivar(1,nhis)=-1
0533             call utword(line,i,j,1) !xmin
0534             call utword(line,i,j,1) !xmax
0535           endif
0536         enddo
0537       elseif(line(i:j).eq.'noerrorbut')then !-----------
0538         ionoerr=ionoerr+1
0539         if(ionoerr.gt.2)stop'xini: not more than 2 noerrorbut !   '
0540         noerrall=noerrall+1
0541         if(noerrall.gt.mxhis/2)stop'xini: to many noerrorbut     '
0542 
0543         call utword(line,i,j,1) !variable-name
0544         cvar=line(i:j)
0545         call xtrans(cvar,inom,ifrnew,nhis)
0546         if(inom.gt.0)then
0547           if(inom.gt.100)then
0548             write(*,*)'xini: noerrorbut can not be used with :',cvar
0549             stop'xini: error with noerrorbut!'
0550           endif
0551           noerrhis(nhis)=noerrall-ionoerr+1
0552           noerr(noerrhis(nhis),ionoerr)=inom
0553           do nb=1,nbin(nhis)
0554              ebin(nb,nac(nhis),ionoerr-1+noerrhis(nhis))=0.d0
0555                   zebin(nb,nac(nhis),ionoerr-1+noerrhis(nhis))=0.d0
0556           enddo
0557         else
0558           ionoerr=ionoerr-1
0559           noerrall=noerrall-1
0560         endif
0561         inoerr(nhis)=ionoerr
0562       elseif(line(i:j).eq.'write')then !-----------
0563         call utword(line,i,j,1)
0564       elseif(line(i:j).eq.'writearray')then !-----------
0565         call utword(line,i,j,1)
0566         iologb=0
0567         iocnxb=0
0568       elseif(line(i:j).eq.'writehisto')then !-----------
0569         call utword(line,i,j,1)
0570         iologb=0
0571         iocnxb=0
0572       elseif(line(i:j).eq.'endhisto'
0573      .      .or.line(i:j).eq.'eh')then   !-----------
0574         ilog(nhis)=.false.
0575         icnx(nhis)=.false.
0576         if(iologb.eq.1)ilog(nhis)=.true.
0577         if(iocnxb.eq.1)icnx(nhis)=.true.
0578         if(ilog(nhis))then
0579           xinc(nhis)=1./log(xmax(nhis)/xmin(nhis))*nbin(nhis)
0580         else
0581           xinc(nhis)=float(nbin(nhis))/(xmax(nhis)-xmin(nhis))
0582         endif
0583         iologb=0
0584         iocnxb=0
0585         jjj=j
0586         cline=line
0587         goto 9999
0588       endif
0589       goto 1
0590 
0591  9999 continue
0592       if(ish.ge.5)then
0593         do n=1,nhis
0594           write (ifch,*) n,': ',ivar(1,n),ivar(2,n)
0595      $         ,'(',ivfra(1,n),ivfra(2,n)
0596      $         ,')',inorm(n)
0597      $         ,xmin(n),xmax(n),ilog(n),icnx(n)
0598      $         ,nbin(n),(idcod(j,n),j=1,nidcod(n))
0599      $         ,' tri:',ntri(n),(itri(j,n),j=1,ntri(n)),'('
0600      $         ,(itfra(j,n),j=1,ntri(n)),')'
0601      $         ,(xmitri(j,n),j=1,ntri(n)) ,(xmatri(j,n),j=1,ntri(n))
0602         enddo
0603         write (ifch,*) (ifra(j),j=1,nfra)
0604       endif
0605       call utprix('xini  ',ish,ishini,5)
0606       return
0607       end
0608 
0609 
0610 c---------------------------------------------------------------------
0611       subroutine xana
0612 c---------------------------------------------------------------------
0613       include 'epos.inc'
0614       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
0615       parameter (mypara=10,mxpara=10)
0616       logical ilog,icnx,itrevt,idmod
0617       double precision bin,bbin,zcbin,zbbin
0618       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
0619      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
0620      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
0621      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
0622      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
0623      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
0624      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
0625      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
0626      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
0627      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
0628      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
0629       double precision ebin,zebin
0630       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
0631      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
0632 
0633       parameter (mxfra=5)
0634       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
0635      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
0636      $     ,emax(mxfra)
0637       double precision bofra
0638       common/dfra/bofra(5,mxfra)
0639       parameter (ntim=1000)
0640       common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0641      &,iorprt(ntim),jorprt(ntim),nprtj
0642 
0643       double precision pgampr,rgampr
0644       common/cgampr/pgampr(5),rgampr(4)
0645 
0646       common/photrans/phoele(4),ebeam
0647 
0648       dimension ten(4,3)
0649       logical go,goo(mxcnt),cont
0650 
0651       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
0652      &     ,multc1,multc83,multc24,multc25,rapgap,ipairs1,xsi
0653       parameter(mxxhis=70)
0654       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
0655      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis),imux(0:mxxhis)
0656      &,ifastjet(0:mxxhis),ijetevent(0:mxxhis),icaltrig(0:mxxhis)
0657 
0658       call utpri('xana  ',ish,ishini,4)
0659 
0660       if(ish.ge.2)then
0661           call alist('fill histograms&',0,0)
0662       endif
0663 
0664       do n=1,nhis
0665       do i=1,mypara
0666         ypara(i,n)=0
0667       enddo
0668       enddo
0669 
0670 
0671       if(ish.ge.5)write(ifch,*)'frames ...'
0672 
0673       if(iappl.eq.6)then
0674         if(mod(iolept/10,10).eq.1) call gakjet(1)
0675         if(mod(iolept/100,10).eq.1) call gakjet(2)
0676       endif
0677 
0678       do l=1,nfra
0679         emax(l)=egyevt/2
0680         if(ifra(l).eq.12)emax(l)=sqrt(pnll**2+prom**2)
0681         if(ifra(l).eq.iframe)then
0682           if(iappl.eq.1.and.iframe.eq.22)emax(l)=ebeam
0683           imofra(1,l)=0
0684           imofra(2,l)=0
0685           imofra(3,l)=0
0686         elseif(ifra(l).eq.11.or.ifra(l).eq.12)then
0687           imofra(1,l)=0
0688           imofra(2,l)=0
0689           bofra(1,l)=0d0
0690           bofra(2,l)=0d0
0691           bofra(3,l)=dsinh(dble(yhaha))
0692           bofra(4,l)=dcosh(dble(yhaha))
0693           bofra(5,l)=1d0
0694           if(ifra(l).eq.11.and.iframe.eq.12)then
0695             imofra(3,l)=1       ! target -> NN
0696           elseif(ifra(l).eq.12.and.iframe.eq.11)then
0697             imofra(3,l)=-1      ! NN -> target
0698           else
0699             imofra(3,l)=0       ! not known
0700           endif
0701         elseif(ifra(l).eq.21)then
0702           if(iframe.ne.21)then
0703             print *, 'invalid frame request'
0704             print *, 'choose frame gamma-nucleon for event run'
0705             stop'bye bye'
0706           endif
0707         elseif(ifra(l).eq.22)then
0708           if(iappl.eq.1)emax(l)=ebeam
0709           if(iframe.eq.21)then
0710             imofra(1,l)=-1      !'  trafo gN -> lab'
0711             imofra(1,l)=0
0712             r1fra(1,l)=rgampr(1)
0713             r1fra(2,l)=rgampr(2)
0714             r1fra(3,l)=rgampr(3)
0715             imofra(2,l)=0
0716             if(iappl.eq.1)then
0717               imofra(3,l)=-2
0718               bofra(1,l)=dsinh(dble(yhaha)) !used for first boost in targ frame
0719               bofra(2,l)=dcosh(dble(yhaha)) !here : pgampr(1)=pgampr(2)=0.
0720               bofra(3,l)=pgampr(3)
0721               bofra(4,l)=pgampr(4)
0722               bofra(5,l)=pgampr(5)
0723             else
0724               imofra(3,l)=-1
0725               bofra(1,l)=pgampr(1)
0726               bofra(2,l)=pgampr(2)
0727               bofra(3,l)=pgampr(3)
0728               bofra(4,l)=pgampr(4)
0729               bofra(5,l)=pgampr(5)
0730             endif
0731           elseif(iframe.eq.22)then
0732                                 ! nothing to do already gN-frame
0733           else
0734             print *, 'invalid frame request'
0735             print *, 'choose frame gamma-nucleon or lab for event run'
0736             stop'bye bye'
0737           endif
0738         elseif(ifra(l).eq.23)then
0739           if(iframe.eq.21)then
0740             imofra(1,l)=0       ! gN -> breit-frame
0741             r1fra(1,l)=rgampr(1)
0742             r1fra(2,l)=rgampr(2)
0743             r1fra(3,l)=rgampr(3)
0744             imofra(2,l)=0
0745             imofra(3,l)=1
0746             bofra(1,l)=0d0
0747             bofra(2,l)=0d0
0748             bofra(3,l)=rgampr(4)
0749             bofra(4,l)=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
0750             bofra(5,l)=sqrt( bofra(4,l)**2-rgampr(4)**2)
0751           elseif(iframe.eq.23)then
0752                                 ! nothing to do already breit-frame
0753           else
0754             print *, 'invalid frame request'
0755             print *, 'choose frame gamma-nucleon or lab for event run'
0756             stop'bye bye'
0757           endif
0758         elseif(ifra(l).eq.33.or.ifra(l).eq.36)then
0759           if(ifra(l).eq.33)then
0760             call gakthru(ten,2)
0761           else
0762             call gakthru(ten,3)
0763           endif
0764           if(ten(4,1).lt.0.)then
0765             imofra(1,l)=0
0766             imofra(2,l)=0
0767             imofra(3,l)=0
0768           else
0769             arox=ten(1,1)
0770             aroy=ten(2,1)
0771             aroz=ten(3,1)
0772             brox=ten(1,2)
0773             broy=ten(2,2)
0774             broz=ten(3,2)
0775             call utrota(1,arox,aroy,aroz,brox,broy,broz)
0776             imofra(1,l)=1
0777             r1fra(1,l)=arox
0778             r1fra(2,l)=aroy
0779             r1fra(3,l)=aroz
0780             imofra(2,l)=1
0781             r2fra(1,l)=brox
0782             r2fra(2,l)=broy
0783             r2fra(3,l)=broz
0784             imofra(3,l)=0       !no boost
0785           endif
0786           bofra(1,l)=dble(ten(4,1)) !usually this is for boosting
0787           bofra(2,l)=dble(ten(4,2)) !I abuse it to store the eigenvalues
0788           bofra(3,l)=dble(ten(4,3)) !
0789         elseif(ifra(l).eq.32.or.ifra(l).eq.34.or.ifra(l).eq.35)then
0790           if(ifra(l).eq.32)then
0791             call gaksphe(ten,2.,2)
0792           elseif(ifra(l).eq.34)then
0793             call gaksphe(ten,1.,2)
0794           else
0795             call gaksphe(ten,2.,3)
0796           endif
0797           if(ten(4,1).lt.0.)then
0798             imofra(1,l)=0
0799             imofra(2,l)=0
0800             imofra(3,l)=0
0801           else
0802             arox=ten(1,1)
0803             aroy=ten(2,1)
0804             aroz=ten(3,1)
0805             brox=ten(1,2)
0806             broy=ten(2,2)
0807             broz=ten(3,2)
0808             call utrota(1,arox,aroy,aroz,brox,broy,broz)
0809             imofra(1,l)=1
0810             r1fra(1,l)=arox
0811             r1fra(2,l)=aroy
0812             r1fra(3,l)=aroz
0813             imofra(2,l)=1
0814             r2fra(1,l)=brox
0815             r2fra(2,l)=broy
0816             r2fra(3,l)=broz
0817             imofra(3,l)=0
0818           endif
0819           bofra(1,l)=dble(ten(4,1))
0820           bofra(2,l)=dble(ten(4,2))
0821           bofra(3,l)=dble(ten(4,3))
0822         endif
0823       enddo
0824 
0825       do n=1,nhis
0826         itrevt(n)=.false.
0827         if(ivar(1,n).ge.100.and.ivar(1,n).le.199) sval(1,n)=0.
0828         if(ivar(2,n).ge.100.and.ivar(2,n).le.199) sval(2,n)=0.
0829         if(ivar(1,n).gt.300.and.ivar(1,n).lt.400)then
0830           call xval(n,ivar(1,n),ivfra(1,n),0,x) !initializing of  variables
0831         endif
0832         if(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
0833           call xval(n,ivar(2,n),ivfra(2,n),0,y) !
0834         endif
0835         do j=1,ntri(n)
0836           valtri(j,n)=0.
0837         enddo
0838         do j=1,nbin(n)  !copy bins
0839           bin(j,3-nac(n),n)=bin(j,nac(n),n)
0840           zcbin(j,3-nac(n),n)=zcbin(j,nac(n),n)
0841         enddo
0842         if(lookcontr(n).gt.0)then
0843           do j=1,nbin(n)
0844             do loo=lookcontr(n),lookcontrx(n)
0845                     bbin(j,3-nac(n),loo)=bbin(j,nac(n),loo)
0846                     zbbin(j,3-nac(n),loo)=zbbin(j,nac(n),loo)
0847             enddo
0848           enddo
0849         endif
0850         if(inoerr(n).gt.0)then
0851           do j=1,nbin(n)
0852             do nn=1,inoerr(n)
0853               ebin(j,3-nac(n),nn-1+noerrhis(n))=ebin(j,nac(n),
0854      &                                          nn-1+noerrhis(n))
0855               zebin(j,3-nac(n),nn-1+noerrhis(n))=zebin(j,nac(n),
0856      &                                          nn-1+noerrhis(n))
0857             enddo
0858           enddo
0859         endif
0860       enddo
0861 
0862       if(imulty1.eq.1)then
0863         if(ish.ge.5)write(ifch,*)'Calculate standard variables ...'
0864         call StandardVariables
0865       endif
0866       if(ipairs1.eq.1)then
0867         if(ish.ge.5)write(ifch,*)'Calculate pair variables ...'
0868         call PairVariables
0869       endif
0870       if(ish.ge.5)write(ifch,*)'Call corrtrig ...'
0871       do n=1,icorrtrig(0)
0872         call corrtrig(icorrtrig(n))
0873       enddo
0874       if(ish.ge.5)write(ifch,*)'Call hardevent ...'
0875       do n=1,ihardevent(0)
0876         call hardevent(ihardevent(n))
0877       enddo
0878       if(ish.ge.5)write(ifch,*)'Call mux ...'
0879       do n=1,imux(0)
0880         call mux(imux(n))
0881       enddo
0882       if(ish.ge.5)write(ifch,*)'Call caltrig ...'
0883       do n=1,icaltrig(0)
0884         call caltrig(icaltrig(n))
0885       enddo
0886       if(ish.ge.5)write(ifch,*)'Call jetfind ...'
0887       do n=1,ijetfind1(0)
0888         call jetfind(1,ijetfind1(n))
0889       enddo
0890       do n=1,ijetfind2(0)
0891         call jetfind(2,ijetfind2(n))
0892       enddo
0893       if(ish.ge.5)write(ifch,*)'Call fastjet ...'
0894       do n=1,ifastjet(0)
0895         call fastjet(ifastjet(n))
0896       enddo
0897 c don't change order here : jetevent should be called always after fastjet !
0898       if(ish.ge.5)write(ifch,*)'Call jetevent ...'
0899       do n=1,ijetevent(0)
0900         call jetevent(ijetevent(n))
0901       enddo
0902 
0903 
0904 c...........................loop nptl...................................
0905       ncontr=0
0906       if(ish.ge.5)write(ifch,*)'Loop nptl ...'
0907       do j=1,nptl
0908         if(iorptl(j).lt.0.or.(istptl(j).lt.100.and.istptl(j).gt.istmax))
0909      &                                                            goto 8
0910         if(ish.ge.5)write(ifch,*)'ptl :',j
0911         call idchrg(idptl(j),ch)
0912         do i=1,nfra
0913           iffra(i)=0            !flag if frame calculated or not
0914         enddo
0915         do n=1,nhis
0916           if(ivar(1,n).eq.-1.or.ivar(2,n).eq.-1)goto 9
0917           if(ivar(1,n).ge.200.and.ivar(2,n).ge.200)goto 9 !skip particle loop if event variables
0918 
0919 c...........check ids
0920           go=nidcod(n).eq.0
0921           do i=1,nidcod(n)
0922             if(istptl(j).eq.0.and.idcod(i,n).eq.10000)then    !all final particle
0923               go=.true.     
0924             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9995)then    !all particles but nuclei
0925               if(abs(idptl(j)).lt.10000) go=.true.     
0926             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9990)then  !all hadrons
0927               if((abs(idptl(j)).ge.100.or.abs(idptl(j)).eq.20)
0928      $         .and.abs(idptl(j)).lt.10000) go=.true.      
0929             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9985)then !neutral particles
0930               if(abs(ch).lt.0.1.and.abs(idptl(j)).lt.10000) go=.true. 
0931             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9980)then !charged particles
0932               if(abs(ch).gt.0.1.and.abs(idptl(j)).lt.10000) go=.true. 
0933             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9975)then  !neutral hadrons
0934               if(abs(ch).lt.0.1.and.(abs(idptl(j)).ge.100
0935      $ .or.abs(idptl(j)).eq.20).and.abs(idptl(j)).lt.10000) go=.true.
0936             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9970)then  !charged hadrons
0937               if(abs(ch).gt.0.1.and.abs(idptl(j)).ge.100
0938      $         .and.abs(idptl(j)).lt.10000) go=.true.
0939             elseif(istptl(j).eq.0.and.idcod(i,n).eq.-9960)then !negative hadrons
0940               if(ch.lt.-0.1.and.abs(idptl(j)).ge.100
0941      $             .and.abs(idptl(j)).lt.10000)go=.true.
0942             elseif(istptl(j).eq.0.and.idcod(i,n).eq.9960)then  !positive hadrons
0943               if(ch.gt.0.1.and.abs(idptl(j)).ge.100
0944      $         .and.abs(idptl(j)).lt.10000)go=.true.
0945             elseif((istptl(j).le.1.or.istptl(j).ge.10)
0946      $            .and.ityptl(n).ne.61
0947      $            .and.idcod(i,n).eq.idptl(j))then
0948               go=.true.
0949             elseif(istptl(j).gt.100.and.idcod(i,n).eq.9999)then   !jets from fastjet
0950               go=.true.
0951             endif
0952           enddo
0953           if(ish.ge.10)write(ifch,*)j,' id,ist',idptl(j),istptl(j),go
0954 
0955 c...........check weak decay
0956           if(go)then
0957             if(noweak(n).eq.1)then  !do not consider weak decay products
0958              if(iorptl(j).ne.0)then
0959               idora=abs( idptl(iorptl(j)) )
0960               if(  idora.eq.20   .or.idora.eq.2130
0961      &               .or.idora.eq.2230 .or.idora.eq.1130
0962      &         .or.idora.eq.2330 .or.idora.eq.1330
0963      &         .or.idora.eq.3331 )go=.false.
0964              ! print *, j,n, '   ', idptl(j),idora,go
0965              endif
0966             endif
0967           endif
0968 
0969 c...........check triggers
0970           if(go)then
0971             if(ish.ge.7)write(ifch,*)'  check triggers in histogram ',n
0972             ncontr=0
0973             do i=1,ntri(n)
0974               if(ish.ge.7)write(ifch,*)'  trigger variable: ',itri(i,n)
0975               if(itri(i,n).lt.100)then
0976                 call xval(n,itri(i,n),itfra(i,n),j,x)
0977                 if(ntrc(i,n).ne.-1)then
0978                   call triggercondition(i,n,x,go)
0979                 else
0980                   ncontr=ncontr+1
0981                   goo(ncontr)=.true.
0982                   call triggercondition(i,n,x,goo(ncontr))
0983                   if((ivar(1,n).gt.100.and.ivar(1,n).lt.200)
0984      .                 .or.(ivar(2,n).gt.100.and.ivar(2,n).lt.200))then
0985                     print*,'!-----------------------------------------'
0986                     print*,'!  100-199 event variables can not be used'
0987                     print*,'! in connection with "trigger contr ..."  '
0988                     print*,'!-----------------------------------------'
0989                     stop'in xana (1).                      '
0990                   endif
0991                 endif
0992               elseif(itri(i,n).lt.200)then
0993                 if(ntrc(i,n).eq.-1)then
0994                     print*,'!-----------------------------------------'
0995                     print*,'!  100-199 event variables can not be used'
0996                     print*,'! in connection with "trigger contr ..."  '
0997                     print*,'!-----------------------------------------'
0998                     stop'in xana (2).                      '
0999                 endif
1000                 call xval(n,itri(i,n),itfra(i,n),j,x)
1001                 valtri(i,n)=valtri(i,n)+x
1002               endif
1003             enddo
1004           endif
1005 
1006 c............fill histogram
1007           if(go)then
1008             if(ish.ge.6)write(ifch,*)'  fill histogram '
1009      &            ,n,ivar(1,n),ivar(2,n),ivfra(2,n)
1010             cont=.false.
1011             if(ivar(1,n).lt.100.or.ivar(2,n).lt.100)then
1012               if(ivar(2,n).lt.100)then
1013                 call xval(n,ivar(2,n),ivfra(2,n),j,y)
1014                 sval(2,n)=y
1015                 cont=.true.
1016               endif
1017               if(ivar(1,n).lt.100)then
1018                 call xval(n,ivar(1,n),ivfra(1,n),j,x)
1019                 if(x.ge.xmin(n).and.x.le.xmax(n))then
1020                   norm3=mod(inorm(n)/100,10)
1021                   if(norm3.eq.1)then
1022                     y=y*x
1023                   elseif(norm3.eq.2.and.ivar(1,n).eq.63.and.x.ne.0.)then
1024                     y=y/(x+pptl(5,j))/2/pi
1025                   elseif(norm3.eq.2.and.ivar(1,n).ne.63.and.x.ne.0.)then
1026                     y=y/x/2/pi
1027                   elseif(norm3.eq.4.and.x.ne.0.)then
1028                     y=y/x**1.5
1029                   elseif(norm3.eq.5.and.x.ne.0.)then
1030                     y=y/x
1031                   elseif(norm3.eq.7.and.x.ne.0.)then
1032                     y=y/x/sqrt(x-pptl(5,j))
1033                   endif
1034                   if(icnx(n))then
1035                     call fillhistoconex(n,x,y,ivfra(2,n),j)   !for conex
1036                   else
1037                     if(ilog(n))then
1038                       nb=1+int(log(x/xmin(n))*xinc(n))
1039                     else
1040                       nb=1+int((x-xmin(n))*xinc(n))
1041                     endif
1042                     bin(nb,nac(n),n)=bin(nb,nac(n),n)+y
1043                    if(ncontr.gt.0)then  !ptl trigger contr
1044                       do nn=1,ncontr
1045                         if(goo(nn))then
1046                              bbin(nb,nac(n),lookcontr(n)-1+nn)=
1047      &                  bbin(nb,nac(n),lookcontr(n)-1+nn)+y
1048                              zbbin(nb,nac(n),lookcontr(n)-1+nn)=
1049      &                  zbbin(nb,nac(n),lookcontr(n)-1+nn)+1
1050                         endif
1051                       enddo
1052                     endif
1053                     if(inoerr(n).gt.0)then
1054                       do nn=1,inoerr(n)
1055                        call xval(n,noerr(noerrhis(n),nn),ivfra(2,n),j,y)
1056                         ebin(nb,nac(n),nn-1+noerrhis(n))=
1057      &                       ebin(nb,nac(n),nn-1+noerrhis(n))+y
1058                         zebin(nb,nac(n),nn-1+noerrhis(n))=
1059      &                       zebin(nb,nac(n),nn-1+noerrhis(n))+1
1060                       enddo
1061                     endif
1062                     zcbin(nb,nac(n),n)=zcbin(nb,nac(n),n)+1
1063                   endif
1064                   itrevt(n)=.true.
1065                 endif
1066               endif
1067             endif
1068             if(ivar(1,n).gt.100.and.ivar(1,n).lt.200)then
1069               call xval(n,ivar(1,n),ivfra(1,n),j,x)
1070               sval(1,n)=sval(1,n)+x
1071             endif
1072             if(ivar(2,n).gt.100.and.ivar(2,n).lt.200)then
1073               call xval(n,ivar(2,n),ivfra(2,n),j,y)
1074               sval(2,n)=sval(2,n)+y
1075               cont=.true.
1076             endif
1077             if(ivar(1,n).gt.300.and.ivar(1,n).lt.400)then
1078               call xval(n,ivar(1,n),ivfra(1,n),j,x)
1079             endif
1080             if(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
1081               call xval(n,ivar(2,n),ivfra(2,n),j,y)
1082               cont=.true.
1083             endif
1084             if(ish.ge.6.and.cont)write (ifch,*)
1085      *                           '   ---> histo n,x,y:',n,x,y
1086           endif
1087    9      continue
1088         enddo
1089   8     continue
1090       enddo
1091 c...........................end loop nptl...........................
1092 
1093       do n=1,nhis
1094       if(ivar(1,n).eq.-1.or.ivar(2,n).eq.-1)goto 99
1095 
1096 c........check event triggers
1097 
1098       if(ish.ge.7)write(ifch,*)'  check event triggers in histogram ',n
1099        go=.true.
1100         ncontr=0
1101         do i=1,ntri(n)
1102           if(itri(i,n).gt.100)then
1103             if(itri(i,n).lt.200)then
1104               x=valtri(i,n)
1105             else
1106               call xval(n,itri(i,n),itfra(i,n),0,x)
1107             endif
1108             if(ntrc(i,n).ne.-1)then
1109               call triggercondition(i,n,x,go)
1110             else
1111               ncontr=ncontr+1
1112               goo(ncontr)=.true.
1113               call triggercondition(i,n,x,goo(ncontr))
1114             endif
1115           endif
1116         enddo
1117 
1118 c........event variables > 200
1119 
1120         if(go)then
1121           if(ivar(1,n).gt.100)then
1122             if(ivar(1,n).gt.200.and.ivar(1,n).lt.300)then
1123               call xval(n,ivar(1,n),ivfra(1,n),0,x)
1124             elseif(ivar(1,n).gt.300.and.ivar(1,n).lt.400)then
1125               call xval(n,ivar(1,n),ivfra(1,n),nptl+1,x)
1126             elseif(ivar(1,n).gt.100.and.ivar(1,n).lt.200)then
1127               x=sval(1,n)
1128             else
1129               call xval(n,ivar(1,n),ivfra(1,n),0,x)
1130             endif
1131             if(ivar(2,n).gt.200.and.ivar(2,n).lt.300)then
1132               call xval(n,ivar(2,n),ivfra(2,n),0,y)
1133             elseif(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
1134               call xval(n,ivar(2,n),ivfra(2,n),nptl+1,y)
1135             elseif(ivar(2,n).gt.0.and.ivar(2,n).lt.200)then
1136               y=sval(2,n)
1137             else             !inom>500
1138               call xval(n,ivar(2,n),ivfra(2,n),0,y)
1139             endif
1140 c The following doesn't work for ivar(2,n)<100, since particle number is not defined !
1141 c            if(ivar(2,n).gt.200.and.ivar(2,n).lt.300)then
1142 c              call xval(n,ivar(2,n),ivfra(2,n),0,y)
1143 c            elseif(ivar(2,n).gt.300.and.ivar(2,n).lt.400)then
1144 c              call xval(n,ivar(2,n),ivfra(2,n),nptl+1,y)
1145 c            elseif(ivar(2,n).gt.100.and.ivar(2,n).lt.200)then
1146 c              y=sval(2,n)
1147 c            else
1148 c              call xval(n,ivar(2,n),ivfra(2,n),0,y)
1149 c            endif
1150             if(mod(inorm(n)/100,10).eq.1)y=y*x
1151             if(mod(inorm(n)/100,10).eq.2.and.x.ne.0.)y=y/x/2/pi
1152             if(mod(inorm(n)/100,10).eq.4.and.x.ne.0.)y=y/x**1.5
1153             if(mod(inorm(n)/100,10).eq.5.and.x.ne.0.)y=y/x
1154             sval(1,n)=x
1155             sval(2,n)=y
1156             if(ish.ge.6)write (ifch,*) 'histo n,x,y:',n,x,y
1157             if(x.ge.xmin(n).and.x.le.xmax(n))then
1158               if(ilog(n))then
1159                 nb=1+int(log(x/xmin(n))*xinc(n))
1160               else
1161                 nb=1+int((x-xmin(n))*xinc(n))
1162               endif
1163               bin(nb,nac(n),n)=bin(nb,nac(n),n)+y
1164               if(ncontr.gt.0)then
1165                 do nn=1,ncontr
1166                   if(goo(nn))
1167      &                    bbin(nb,nac(n),lookcontr(n)-1+nn)=
1168      &             bbin(nb,nac(n),lookcontr(n)-1+nn)+y
1169                 enddo
1170               endif
1171               zcbin(nb,nac(n),n)=zcbin(nb,nac(n),n)+1
1172               itrevt(n)=.true.
1173             endif
1174           endif
1175         endif
1176 
1177 c........particle variables
1178 
1179         if(go)then
1180           if(ivar(1,n).le.100)then
1181             if(ncontr.gt.0)then  !event trigger contr
1182               do nb=1,nbin(n)
1183                 do nn=1,ncontr
1184                   if(goo(nn))
1185      &                     bbin(nb,nac(n),lookcontr(n)-1+nn)=
1186      &              bbin(nb,nac(n),lookcontr(n)-1+nn)
1187      &              +bin(nb,nac(n),n)-bin(nb,3-nac(n),n)
1188                 enddo
1189               enddo
1190             endif
1191           endif
1192         endif
1193 
1194 c............event ok (increase ncevt) or not (take copy)
1195 
1196         if(go)then
1197           ncevt(n)=ncevt(n)+1
1198           if(ncontr.gt.0)then
1199             do nn=1,ncontr
1200               loo=lookcontr(n)-1+nn
1201               if(goo(nn))
1202      &        nccevt(loo)=nccevt(loo)+1
1203             enddo
1204           endif
1205         else
1206           if(ish.ge.6)write (ifch,*) 'event rejected for histo',n
1207           nac(n)=3-nac(n)
1208           itrevt(n)=.false.
1209         endif
1210 
1211  99   continue
1212       enddo
1213 
1214       call utprix('xana  ',ish,ishini,4)
1215       end
1216 
1217 c--------------------------------------------------------------------
1218       subroutine triggercondition(i,n,x,go)
1219 c--------------------------------------------------------------------
1220 c ntrc is used to distinguish the different usage of trigger:
1221 c
1222 c    trigger var xmin xmax
1223 c             ntrc=1
1224 c    trigger or n var1 xmin1 xmax1 var2 xmin2 xmax2 ... varn xminn xmaxn
1225 c          1  ntrc=2
1226 c          2  ntrc=0
1227 c              ...
1228 c         n-1 ntrc=0
1229 c          n  ntrc=3
1230 c    trigger contr n var1 xmin1 xmax1 var2 xmin2 xmax2 ... varn xminn xmaxn
1231 c             ntrc=-1
1232 c--------------------------------------------------------------------
1233       include 'epos.inc'
1234       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
1235       parameter (mypara=10,mxpara=10)
1236       logical ilog,icnx,itrevt,idmod
1237       double precision bin,bbin,zcbin,zbbin
1238       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
1239      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
1240      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
1241      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
1242      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
1243      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
1244      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
1245      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
1246      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
1247      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
1248      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
1249       double precision ebin,zebin
1250       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
1251      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
1252       logical go,gox,ok,goz
1253                 xmn=xmitri(i,n)
1254                 xmx=xmatri(i,n)
1255                 if(xmn.eq.-123456.and.xmx.eq.-123456)then  !for leading part
1256                   xmn=float(idlead)
1257                   xmx=float(idlead)
1258                 endif
1259                 pmn=xmitrp(i,n)
1260                 pmx=xmatrp(i,n)
1261                 if(abs(ntrc(i,n)).eq.1)then
1262                   goz=.true.
1263                   if(pmn.gt.99.999.and.pmx.gt.99.999)then
1264                     if(x.lt.xmn.or.x.gt.xmx)goz=.false.
1265                   else
1266                     if(x.lt.xmn-0.5.or.x.gt.xmx+0.5)goz=.false.
1267                     ok=rangen().le.xmitrp(i,n)/100.
1268                     if(.not.ok.and.x.lt.xmn+0.5)goz=.false.
1269                     ok=rangen().le.xmatrp(i,n)/100.
1270                     if(.not.ok.and.x.gt.xmx-0.5)goz=.false.
1271                   endif
1272                   if(.not.goz)go=.false.
1273                 else
1274                   if(ntrc(i,n).eq.2)gox=.false.
1275                   goz=.true.
1276                   if(pmn.gt.99.999.and.pmx.gt.99.999)then
1277                     if(x.lt.xmn.or.x.gt.xmx)goz=.false.
1278                   else
1279                     if(x.lt.xmn-0.5.or.x.gt.xmx+0.5)goz=.false.
1280                     ok=rangen().le.xmitrp(i,n)/100.
1281                     if(.not.ok.and.x.lt.xmn+0.5)goz=.false.
1282                     ok=rangen().le.xmatrp(i,n)/100.
1283                     if(.not.ok.and.x.gt.xmx-0.5)goz=.false.
1284                   endif
1285                   if(goz)gox=.true.
1286                   if(ntrc(i,n).eq.3.and..not.gox)go=.false.
1287                 endif
1288                 if(ish.ge.9)write(ifch,*)'trigger conditions '
1289      &                                   ,i,n,xmn,x,xmx,go
1290                 end
1291 
1292 c-----------------------------------------------------------------------
1293       subroutine fillhistoconex(n,x,y,lf,j)   !for conex
1294 c-----------------------------------------------------------------------
1295       include 'epos.inc'
1296       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
1297       parameter (mypara=10,mxpara=10)
1298       logical ilog,icnx,itrevt,idmod
1299       double precision bin,bbin,zcbin,zbbin
1300       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
1301      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
1302      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
1303      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
1304      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
1305      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
1306      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
1307      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
1308      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
1309      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
1310      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
1311       double precision ebin,zebin
1312       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
1313      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
1314 
1315            if(.not.(mod(inorm(n),10).ne.4
1316      &      .and.mod(inorm(n),10).ne.6
1317      &      .and.mod(inorm(n)/100,10).ne.3))return
1318                     if(ilog(n))then
1319                       c=(xmax(n)/xmin(n))**(1./real(nbin(n)))
1320                       nde=nint(1./log10(c))
1321                       nb=max(1,1+int(log10(x/xmin(n))*nde))
1322                       xmb=xmin(n)*c**(nb-0.5)
1323                       if(x.gt.xmb.and.nb.lt.nbin(n))then
1324                         if(x.gt.xmax(n))
1325      &                    write(ifmt,*)'xana max ?',x,xmax(n),nb
1326                         nbx=1
1327                         xmx=c*xmb
1328                       elseif(x.lt.xmb.and.nb.gt.1)then
1329                         if(x.lt.xmin(n))write(ifmt,*)'xana min ?',x,nb
1330                         nbx=-1
1331                         xmx=xmb/c
1332                       else
1333                         nbx=0
1334                         xmx=0.
1335                       endif
1336                     else
1337                       c=(xmax(n)-xmin(n))/real(nbin(n))
1338                       nb=max(1,1+int((x-xmin(n))/c))
1339                       xmb=xmin(n)+c*(nb-0.5)
1340                       if(x.gt.xmb)then
1341                         nbx=1
1342                         xmx=c+xmb
1343                       elseif(x.lt.xmb)then
1344                         nbx=-1
1345                         xmx=xmb-c
1346                       else
1347                         nbx=0
1348                         xmx=0.
1349                       endif
1350                     endif
1351                     xc=(x-xmx)/(xmb-xmx)
1352                     xc=max(0.,min(1.,xc))
1353                     bin(nb,nac(n),n)=bin(nb,nac(n),n)+xc*y
1354                     if(nbx.ne.0)bin(nb+nbx,nac(n),n)
1355      &                  =bin(nb+nbx,nac(n),n)+(1.-xc)*y
1356                     zcbin(nb,nac(n),n)=zcbin(nb,nac(n),n)+1
1357                     if(inoerr(n).gt.0)then
1358                       do nn=1,inoerr(n)
1359                         call xval(n,noerr(noerrhis(n),nn),lf,j,y2)
1360                         ebin(nb,nac(n),nn-1+noerrhis(n))=
1361      &                       ebin(nb,nac(n),nn-1+noerrhis(n))+y2
1362                         zebin(nb,nac(n),nn-1+noerrhis(n))=
1363      &                       zebin(nb,nac(n),nn-1+noerrhis(n))+1
1364                       enddo
1365                     endif
1366       end
1367 
1368 c---------------------------------------------------------------------
1369       subroutine xhis(n)
1370 c---------------------------------------------------------------------
1371       include 'epos.inc'
1372       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
1373       parameter (mypara=10,mxpara=10)
1374       logical ilog,icnx,itrevt,idmod
1375       double precision bin,bbin,zcbin,zbbin
1376       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
1377      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
1378      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
1379      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
1380      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
1381      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
1382      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
1383      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
1384      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
1385      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
1386      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
1387       double precision ebin,zebin
1388       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
1389      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
1390       dimension xx(mxbin)
1391 
1392       double precision histoweight
1393       common/chiswei/histoweight
1394       common/cyield/yield
1395       common/csigma/sigma
1396       double precision dcel
1397       common/ems3/dcel,ad
1398       common/geom/rmproj,rmtarg,bmax,bkmx
1399       save cnormx
1400 
1401       if(ivar(1,n).eq.-1)then
1402         nrbins=0
1403         goto 9999
1404       endif
1405 
1406 c.......here normalization.......................................
1407 c           see also   "..........fill histogram"
1408 c.................................................................
1409 c     the norm ( inorm(n) ) is a number hijk which normalizes to:
1410 c
1411 c  k  0:  * 1
1412 c     1:  / number of events
1413 c     2:  / number of triggered events
1414 c     4:  / bin-counts
1415 c     5:  / bin sum
1416 c     6:  / number of summed bin-counts (yield=1.)
1417 c     7:  uses same normalization as one histo before
1418 c
1419 c  j  0:  * 1
1420 c     1:  / bin-width
1421 c     2:  * sigma_total / bin-width
1422 c     3:  * sigma_diff / bin-width
1423 c
1424 c  i  0:  * 1
1425 c     1:  y => y*x
1426 c     2:  y => y/x/2/pi (modified for mt0)
1427 c     3:  kno-scaling
1428 c     4:  y => y/x**1.5
1429 c     5:  y => y/x
1430 c     6:  y => y*xi (for conex, xi=x of the bin)
1431 c     7:  y => y/x/(x-m)
1432 c
1433 c  h  0: normal
1434 c     1: accumulated
1435 c
1436 c.................................................................
1437 
1438       norm1=mod(inorm(n),10)
1439       norm2=mod(inorm(n)/10,10)
1440       norm3=mod(inorm(n)/100,10)
1441       norm4=mod(inorm(n)/1000,10)
1442       nctbin=0
1443       sumbin=0
1444       do l=1,nbin(n)
1445         nctbin=nctbin+zcbin(l,nac(n),n)
1446         sumbin=sumbin+bin(l,nac(n),n)
1447         if(norm1.eq.4.and.zcbin(l,nac(n),n).ne.0d0)then
1448           bin(l,nac(n),n)=bin(l,nac(n),n)/zcbin(l,nac(n),n)
1449           if(lookcontr(n).gt.0)then
1450             do loo=lookcontr(n),lookcontrx(n)
1451               if(zbbin(l,nac(n),loo).ne.0.)
1452      &           bbin(l,nac(n),loo)=bbin(l,nac(n),loo)
1453      &            /zbbin(l,nac(n),loo)
1454             enddo
1455           endif
1456         endif
1457         if(ilog(n))then
1458           xx(l)=xmin(n)*(xmax(n)/xmin(n))**((float(l)-.5)/nbin(n))
1459         else
1460           xx(l)=(float(l)-0.5)*(xmax(n)-xmin(n))/nbin(n)+xmin(n)
1461         endif
1462       enddo
1463       cnorm=1.
1464       if(norm1.eq.1)cnorm=1./float(nevent)
1465       if(norm1.eq.2)then
1466         if(ncevt(n).ne.0)then
1467           cnorm=1./float(ncevt(n))
1468         else
1469           cnorm=0.
1470         endif
1471       endif
1472       if(norm1.eq.5.and.sumbin.ne.0.)cnorm=1./sumbin
1473       if(norm1.eq.6.and.nctbin.ne.0)cnorm=1./float(nctbin)
1474       if(norm1.eq.7)cnorm=cnormx
1475       cnormx=cnorm
1476       if(ntevt.ne.0)
1477      &   sigma=10.*pi*bmax**2.*nevent/ntevt !total (untriggered) sigma
1478       if(norm2.eq.3)then      !differential (triggered) sigma
1479         if(ntevt.ne.0)
1480      &     sigma=10.*pi*bmax**2.*ncevt(n)/ntevt
1481       endif
1482       if(norm3.eq.3)then      !kno
1483         first=0.
1484         secnd=0.
1485         do l=1,nbin(n)
1486           if(nctbin.ne.0)first=first+xx(l)*zcbin(l,nac(n),n)/nctbin
1487           if(nctbin.ne.0)secnd=secnd
1488      $           +xx(l)**2*zcbin(l,nac(n),n)/nctbin
1489         enddo
1490       else
1491         first=1.
1492       endif
1493       if(ilog(n))then
1494         if(norm2.eq.2.or.norm2.eq.3) cnorm=cnorm*sigma
1495       else
1496         if(norm2.ge.1.and.norm2.le.3) cnorm=cnorm*xinc(n)
1497         if(norm2.eq.2.or.norm2.eq.3) cnorm=cnorm*sigma
1498       endif
1499       do l=1,nbin(n)
1500         bnorm=0.
1501         if(ilog(n).and.norm2.ge.1.and.norm2.le.3)then
1502           bnorm=1./(xmin(n)*exp(float(l)/xinc(n))*(1.-exp(-1./xinc(n))))
1503           bin(l,nac(n),n) =  bin(l,nac(n),n) * bnorm
1504         endif
1505         bin(l,nac(n),n) =  bin(l,nac(n),n) * cnorm
1506         if(lookcontr(n).gt.0)then
1507           if(ilog(n).and.norm2.ge.1.and.norm2.le.3)then
1508             do loo=lookcontr(n),lookcontrx(n)
1509               bbin(l,nac(n),loo)=bbin(l,nac(n),loo) * bnorm
1510             enddo
1511           endif
1512         endif
1513       enddo
1514       f=first
1515       nrbins=nbin(n)
1516       nctbin=0
1517       yield=0.
1518       shft=0
1519        if(nint(xpara(1,n)).eq.999963)shft=xpara(2,n)
1520       do ii=1,nbin(n)
1521         g=1
1522         if(norm3.eq.1.and.xx(ii).ne.0.)g=1./xx(ii)
1523         if(norm3.eq.2)g=2*pi*(xx(ii)+shft)
1524         if(norm3.eq.4)g=xx(ii)**1.5
1525         if(norm3.eq.5)g=xx(ii)
1526         if(norm3.eq.7)g=0
1527         yield=yield+bin(ii,nac(n),n)/xinc(n)*hisfac*f*g
1528       enddo
1529       do l=1,nbin(n)
1530         x=(xx(l)+xshift)      !*xhfact
1531         ar(l,1)=x/f
1532         sigbin=0
1533         if(zcbin(l,nac(n),n).ne.0d0)
1534      *   sigbin=bin(l,nac(n),n)*hisfac*f/sqrt(zcbin(l,nac(n),n))
1535         if(norm4.eq.0.or.l.eq.1)then
1536           ar(l,3)=bin(l,nac(n),n)*hisfac*f
1537           if(lookcontr(n).gt.0)then
1538            do loo=lookcontr(n),lookcontrx(n)
1539              r=1
1540              if(norm1.eq.2.and.nccevt(loo).ne.0.)
1541      *          r=float(ncevt(n))/nccevt(loo)
1542              lo=loo-lookcontr(n)+1
1543              ary(l,lo)=bbin(l,nac(n),loo)*hisfac*f*cnorm*r
1544              if(zbbin(l,nac(n),loo).gt.0.)then
1545                ardy(l,lo)=ary(l,lo)/sqrt(zbbin(l,nac(n),loo))
1546              else
1547                ardy(l,lo)=0
1548              endif
1549              if(norm1.eq.4)ardy(l,lo)=zbbin(l,nac(n),loo)
1550             enddo
1551           endif
1552           if(norm3.eq.6)then   !conex
1553            ar(l,3)=ar(l,3)*xx(l)
1554           endif
1555           ar(l,4)=sigbin
1556         else
1557           ar(l,3)=ar(l-1,3)+bin(l,nac(n),n)*hisfac*f
1558           ar(l,4)=sqrt(ar(l-1,4)**2+sigbin**2)
1559         endif
1560         if(inoerr(n).ge.1)then
1561           if(zebin(l,nac(n),noerrhis(n)).gt.0.d0)then
1562          ar(l,4)=ebin(l,nac(n),noerrhis(n))/zebin(l,nac(n),noerrhis(n))
1563           else
1564             ar(l,4)=0.
1565           endif
1566         endif
1567         if(inoerr(n).eq.2)then
1568           if(zebin(l,nac(n),noerrhis(n)+1).gt.0.d0)then
1569       ar(l,5)=ebin(l,nac(n),noerrhis(n)+1)/zebin(l,nac(n),noerrhis(n)+1)
1570           else
1571             ar(l,5)=0.
1572           endif
1573         endif
1574         if(norm1.eq.4)ar(l,4)=zcbin(l,nac(n),n)
1575       enddo
1576       ionoerr=inoerr(n)
1577       histoweight=dble(ncevt(n))
1578       if(norm1.eq.1)histoweight=dble(nevent)
1579       if(norm1.eq.4)histoweight=0d0
1580 
1581  9999 hisfac=1.
1582       xshift=0
1583       end
1584 
1585 c-----------------------------------------------------------------------
1586       integer function nsdiff(insdif,now)
1587 c-----------------------------------------------------------------------
1588 c  returns  1 if trigger condition for NSD fulfilled and 0 otherwise
1589 c  for  UA1 (insdif=0) or UA5 (insdif=1) or CDF (insdif=2) or STAR (insdif=3,4)
1590 C  or BRAHMS (insdif=5) or NA61 (insdif=6) or CMS (insdif=7)
1591 c  or ATLAS (insdif=8) or ALICE (insdif=9, 10 and 11) 
1592 c  or CMS hadron level (insdif=12) or CMS hadron level double sided (insdif=13)
1593 c  now ... noweak(histogram number) (obsolete)
1594 c-----------------------------------------------------------------------
1595       include 'epos.inc'
1596       integer ncevt,nsdi(0:20)
1597       logical cont
1598       data ncevt/1/
1599       save nsdi,ncevt
1600 c initialization for each event
1601       if(ncevt.eq.nrevt)then
1602         ncevt=ncevt+1
1603         do i=0,20
1604           nsdi(i)=-1
1605         enddo
1606       endif
1607       nsdiff=0
1608            if(insdif.ge.0)then
1609          if(nsdi(insdif).lt.0)then
1610       iii1=0
1611       iii2=0
1612       iii3=0
1613       ipos=0
1614       ineg=0
1615          do npts=1,nptl
1616         if(istptl(npts).ne.0)goto 60
1617           cont=   idptl(npts).ne.120 .and.idptl(npts).ne.-120
1618      *   .and.idptl(npts).ne.130 .and.idptl(npts).ne.-130
1619      *   .and.idptl(npts).ne.1120.and.idptl(npts).ne.-1120
1620      *   .and.idptl(npts).ne.1130.and.idptl(npts).ne.-1130
1621      *   .and.idptl(npts).ne.2230.and.idptl(npts).ne.-2230
1622      *   .and.idptl(npts).ne.2330.and.idptl(npts).ne.-2330
1623      *   .and.idptl(npts).ne.3331.and.idptl(npts).ne.-3331
1624         if(insdif.ne.7.and.insdif.ne.14)then
1625           if(cont)goto 60
1626         endif
1627 c just to avoid warning with gfortran when "now" is not used
1628         nowdum=now
1629 c        if(now.eq.1)then  !do not consider weak decay products
1630 c         if(iorptl(npts).ne.0)then
1631 c          idora=abs( idptl(iorptl(npts)) )
1632 c          if(  idora.eq.20   .or.idora.eq.2130
1633 c     &     .or.idora.eq.2230 .or.idora.eq.1130
1634 c     &     .or.idora.eq.2330 .or.idora.eq.1330
1635 c     &     .or.idora.eq.3331 )goto 60
1636 c          endif
1637 c         endif
1638         pz=pptl(3,npts)
1639         pt=sqrt(pptl(2,npts)**2+pptl(1,npts)**2)
1640         ppp=sqrt(pz**2+pt**2)
1641         Etot=pptl(4,npts)
1642         if(ppp.gt.abs(pz))then
1643           yyy=.5*log((ppp+pz)/(ppp-pz))
1644         else
1645           yyy=sign(100.,pz)
1646         endif
1647         if(insdif.eq.0)then
1648           if(yyy.gt.1.5  .and. yyy.lt.5.5)iii1=1
1649           if(yyy.gt.-5.5 .and. yyy.lt.-1.5)iii2=1
1650         elseif(insdif.eq.1)then
1651           if(yyy.gt.2.   .and. yyy.lt.5.6)iii1=1
1652           if(yyy.gt.-5.6 .and. yyy.lt.-2.)iii2=1
1653         elseif(insdif.eq.2)then
1654           if(yyy.gt.3.2  .and. yyy.lt.5.9)iii1=1
1655           if(yyy.gt.-5.9 .and. yyy.lt.-3.2)iii2=1
1656           if(yyy.gt.0.   .and. yyy.lt.3.0)ipos=ipos+1
1657           if(yyy.gt.-3.0 .and. yyy.lt.0. )ineg=ineg+1
1658         elseif(insdif.eq.3)then
1659           if(yyy.gt.-5.0 .and. yyy.lt.-3.3 )iii1=1
1660           if(yyy.gt. 3.3 .and. yyy.lt. 5.0 )iii2=1
1661         elseif(insdif.eq.4)then
1662           if(yyy.gt.-5.0 .and. yyy.lt.-3.1 )iii1=1
1663           if(yyy.gt. 3.1 .and. yyy.lt. 5.0 )iii2=1
1664         elseif(insdif.eq.5)then
1665           if(yyy.gt.-5.25 .and. yyy.lt.-3.26 )iii1=1
1666           if(yyy.gt. 3.26 .and. yyy.lt. 5.25 )iii2=1
1667         elseif(insdif.eq.6)then  !NA61 trigger if NO charged particle with theta<5.26 mrad
1668           if(pptl(3,npts).gt.0..and.yyy.lt.100.)then
1669            theta=sqrt(pptl(1,npts)**2+pptl(2,npts)**2)/pptl(3,npts)
1670            if(theta.lt.5.26e-3)iii1=1
1671           endif
1672         elseif(insdif.eq.7)then   !CMS NSD corrected using PYTHIA (2010)
1673           if(yyy.gt.-5.2 .and. yyy.lt.-2.9 .and. Etot .gt.3.)iii1=1
1674           if(yyy.gt. 2.9 .and. yyy.lt. 5.2 .and. Etot .gt.3.)iii2=1
1675           if(yyy.gt. -2.5 .and. yyy.lt. 2.5 .and. pt .gt.0.2 .and. cont)
1676      &                                                       iii3=1
1677         elseif(insdif.eq.8)then   !ATLAS
1678           if(yyy.gt.-2.5 .and. yyy.lt.2.5.and.pt.gt.0.5)iii1=1
1679           iii2=1
1680         elseif(insdif.eq.9)then   !ALICE 900 GeV
1681           if(yyy.gt.-3.7 .and. yyy.lt.-1.7 )iii1=1
1682           if(yyy.gt. 2.8 .and. yyy.lt. 5.1 )iii2=1
1683         elseif(insdif.eq.10)then   !ALICE 2.36 TeV
1684           if(yyy.gt.-2 .and. yyy.lt.2 )iii1=1
1685           iii2=1
1686         elseif(insdif.eq.11)then   !ALICE Inel>0
1687           if(yyy.gt.-1 .and. yyy.lt.1 )iii1=1
1688           iii2=1
1689         elseif(insdif.eq.12)then   !CMS hadron level NSD trigger (2011)
1690           if(yyy.gt.-4.4 .and. yyy.lt.-3.9 )iii1=1
1691           if(yyy.gt. 3.9 .and. yyy.lt. 4.4 )iii2=1
1692         elseif(insdif.eq.13)then   !CMS hadron level doubl sided trigger (2012)
1693           if(yyy.gt.-5. .and. yyy.lt.-3. .and. Etot .gt. 3. )iii1=1
1694           if(yyy.gt. 3. .and. yyy.lt. 5. .and. Etot .gt. 3. )iii2=1
1695         elseif(insdif.eq.14)then   !CMS hadron level single sided trigger (HF 2012)
1696           if(yyy.gt.-4.9 .and. yyy.lt.-2.9 .and. Etot .gt. 5. )iii1=1
1697           if(yyy.gt. 2.9 .and. yyy.lt. 4.9 .and. Etot .gt. 5. )iii2=1
1698         endif
1699 60      continue
1700          enddo
1701         if(insdif.le.1)then
1702           if(iii1.eq.1 .and. iii2.eq.1) nsdiff=1
1703         elseif(insdif.eq.2)then
1704           if((iii1.eq.1 .and. iii2.eq.1) .and.
1705      *    ((ipos.ne.0 .and. ineg.ne.0) .and. ipos+ineg.ge.4)) nsdiff=1
1706         elseif(insdif.eq.14)then
1707           if(iii1.eq.1 .or. iii2.eq.1) nsdiff=1
1708         elseif(insdif.eq.3.or.insdif.eq.4
1709      *         .or.insdif.eq.5.or.insdif.eq.8.or.insdif.ge.9)then
1710           if(iii1.eq.1 .and. iii2.eq.1) nsdiff=1
1711         elseif(insdif.eq.6)then
1712           if(iii1.eq.0 .and. iii2.eq.0)then
1713             nsdiff=1
1714           endif
1715         elseif(insdif.eq.7)then
1716           if(iii1.eq.1 .and. iii2.eq.1 .and.iii3.eq.1)then
1717             nsdiff=1
1718           endif
1719         endif
1720         nsdi(insdif)=nsdiff
1721         else
1722           nsdiff=nsdi(insdif)
1723         endif
1724            else
1725       stop'in nsdiff. argument of nsdiff not authorized.        '
1726            endif
1727       end
1728 
1729 c-----------------------------------------------------------------------
1730       integer function isdiff(isdif)
1731 c-----------------------------------------------------------------------
1732 c  returns  1 if trigger condition for single diff fulfilled and 0 otherwise
1733 c  for  UA4 Mult distri (isdif=1) or UA4 xsection (isdif=2) 
1734 c  or CDF SD (isdif=3) or CDF DPE (isdif=4) or CDF min bias (for DD) 
1735 c  (isdif=5)
1736 c-----------------------------------------------------------------------
1737       include 'epos.inc'
1738       isdiff=0
1739            if(isdif.ge.1)then
1740       iii0=0
1741       iii1=0
1742       iii2=0
1743       iii3=0
1744       iii4=0
1745       Et1=0.
1746       Et2=0.
1747       do npts=1,nptl
1748         if(istptl(npts).ne.0)goto 60
1749         if(   abs(idptl(npts)).ne.120
1750      *   .and.abs(idptl(npts)).ne.130
1751      *   .and.abs(idptl(npts)).ne.1120
1752      *   .and.abs(idptl(npts)).ne.1130
1753      *   .and.abs(idptl(npts)).ne.2230
1754      *   .and.abs(idptl(npts)).ne.2330
1755      *   .and.abs(idptl(npts)).ne.3331)goto 60
1756         ppt=pptl(1,npts)**2+pptl(2,npts)**2
1757         ppp=sqrt(ppt+pptl(3,npts)**2)
1758         ppt=sqrt(ppt)
1759         yyy=0.
1760         if(pptl(3,npts).ne.0..and.ppt.ne.0.)yyy=sign(1.,pptl(3,npts))*
1761      *   log((ppp+abs(pptl(3,npts)))/ppt)
1762 c        if(ppp.gt.abs(pptl(3,npts)))then
1763 c          yyy=.5*log((ppp+pptl(3,npts))/(ppp-pptl(3,npts)))
1764 c        else
1765 c          yyy=sign(100.,pptl(3,npts))
1766 c        endif
1767         if(isdif.le.2)yyy=-sign(1.,float(ilprtg))*yyy   !trigger on antiproton (target : ilprt=-1)
1768 c        if(idptl(npts).eq.-1120)then
1769           if(abs(pptl(3,npts)).gt.0.)then
1770             theta=sign(1.,float(ilprtg))*ppt/pptl(3,npts)
1771             if((isdif.le.2.and.theta.gt.2.5e-3.and.theta.lt.4.5e-3).or.
1772      *         (isdif.gt.2.and.theta.gt.0.2e-3.and.theta.lt.1.2e-3))then
1773               iii0=iii0+1
1774 c              write(ifch,*)'la',
1775 c              print *,'la',ilprtg,yyy,iii1,iii2
1776 c     &      ,npts,idptl(npts),ppt,pptl(3,npts),theta,ityptl(npts)
1777             endif
1778 c          endif
1779         endif
1780         if(isdif.eq.1)then
1781           if(yyy.gt.2.5   .and. yyy.lt.5.6)iii1=1
1782         elseif(isdif.eq.2)then
1783           if(yyy.gt.3.  .and. yyy.lt.5.6)iii1=1
1784           if(yyy.gt.-5.6 .and. yyy.lt.-4.4 )iii2=1
1785         elseif(isdif.eq.3)then
1786           if(yyy.gt.2.4  .and. yyy.lt.5.9)iii1=1
1787           if(yyy.gt.-4.2 .and. yyy.lt.1.1)iii3=1
1788           if(yyy.gt.-5.9 .and. yyy.lt.-2.4)iii2=1
1789           if(yyy.gt.-1.1 .and. yyy.lt.4.2)iii4=1
1790         elseif(isdif.eq.4)then
1791           if(ilprtg.eq.-1)then  !antiproton = target
1792             if(yyy.gt.2.4  .and. yyy.lt.5.9)iii1=1
1793             if(yyy.gt.-5.9 .and. yyy.lt.-3.2)iii2=iii2+1
1794           else                  !antiproton = projectile
1795             if(yyy.gt.-5.9 .and. yyy.lt.-2.4)iii1=1
1796             if(yyy.gt. 3.2 .and. yyy.lt.5.9)iii2=iii2+1
1797           endif
1798         elseif(isdif.eq.5)then
1799           if(yyy.gt.3.2  .and. yyy.lt.5.9)iii1=1
1800           if(yyy.gt.-5.9 .and. yyy.lt.-3.2)iii2=1
1801           if(abs(yyy).lt.2.4)Et1=Et1+ppt
1802           if(Et1.gt.0.2)iii3=1
1803           if(abs(yyy).gt.2.2 .and. abs(yyy).lt.4.2 )Et2=Et2+ppt
1804           if(Et2.gt.1.)iii4=1
1805         endif
1806 60      continue
1807       enddo
1808       if(isdif.eq.1)then
1809         if(iii0.eq.1 .and. iii1.eq.1) isdiff=1
1810       elseif(isdif.eq.2)then
1811         if(iii0.eq.1 .and. iii1.eq.1 .and. iii2.ne.1) isdiff=1
1812       elseif(isdif.eq.3)then
1813         if( (iii1.ne.1 .and. iii3.eq.1) .or.
1814      &      (iii2.ne.1 .and. iii4.eq.1)       ) isdiff=1
1815       elseif(isdif.eq.4)then
1816         if(iii0.eq.1 .and. iii1.ne.1 .and. iii2.le.6) isdiff=1
1817 c        if(isdiff.eq.1)then
1818 c      do npts=1,nptl
1819 c        if(istptl(npts).ne.0)goto 80
1820 c        if(   abs(idptl(npts)).ne.120
1821 c     *   .and.abs(idptl(npts)).ne.130
1822 c     *   .and.abs(idptl(npts)).ne.1120
1823 c     *   .and.abs(idptl(npts)).ne.1130
1824 c     *   .and.abs(idptl(npts)).ne.2230
1825 c     *   .and.abs(idptl(npts)).ne.2330
1826 c     *   .and.abs(idptl(npts)).ne.3331)goto 80
1827 c        ppt=pptl(1,npts)**2+pptl(2,npts)**2
1828 c        ppp=sqrt(ppt+pptl(3,npts)**2)
1829 c        ppt=sqrt(ppt)
1830 c        yyy=0.
1831 c        if(pptl(3,npts).ne.0..and.ppt.ne.0.)yyy=sign(1.,pptl(3,npts))*
1832 c     *   log((ppp+abs(pptl(3,npts)))/ppt)
1833 c        print *,nrevt,yyy,idptl(npts),ityptl(npts)
1834 c80      continue
1835 c        enddo
1836 c          print *,'dpe',iii0
1837 c        endif
1838       elseif(isdif.eq.5)then
1839         if(iii1+iii2+iii3+iii4.eq.4) isdiff=1
1840       else
1841         stop'in sdiff. argument of sdiff not authorized.        '
1842       endif
1843              endif
1844       end
1845 
1846 c----------------------------------------------------------------------
1847       subroutine xtrans(cvar,inom,ifr,n)
1848 c----------------------------------------------------------------------
1849       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
1850      &     ,multc1,multc83,multc24,multc25,rapgap,ipairs1,xsi
1851       parameter(mxxhis=70)
1852       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
1853      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis),imux(0:mxxhis)
1854      &,ifastjet(0:mxxhis),ijetevent(0:mxxhis),icaltrig(0:mxxhis)
1855 
1856       character*6 cvar
1857       ifr=0
1858       if(cvar.eq.'numptl')then
1859         inom=1
1860       elseif(cvar.eq.'npaptl')then
1861         inom=2
1862       elseif(cvar.eq.'npmptl')then
1863         inom=3
1864       elseif(cvar.eq.'ispptl')then
1865         inom=4
1866       elseif(cvar.eq.'rapx')then
1867         inom=5
1868       elseif(cvar.eq.'iptlfr')then
1869         inom=6
1870       elseif(cvar.eq.'rinp')then
1871         inom=7
1872       elseif(cvar.eq.'eco')then
1873         inom=8
1874       elseif(cvar.eq.'tau')then
1875         inom=9
1876       elseif(cvar.eq.'ctr')then
1877         inom=10
1878       elseif(cvar.eq.'v2np')then
1879         inom=11
1880         imulty1=1             !to switch on the calculation of "Standard variable"
1881       elseif(cvar.eq.'absrap')then
1882         inom=12
1883       elseif(cvar.eq.'rap')then
1884         inom=13
1885       elseif(cvar.eq.'xp')then
1886         inom=14
1887       elseif(cvar.eq.'xe')then
1888         inom=15
1889       elseif(cvar.eq.'pt')then
1890         inom=16
1891       elseif(cvar.eq.'p1a')then
1892         inom=17
1893       elseif(cvar.eq.'p2a')then
1894         inom=18
1895       elseif(cvar.eq.'xi')then
1896         inom=19
1897       elseif(cvar.eq.'xf')then
1898         inom=20
1899       elseif(cvar.eq.'t')then
1900         inom=21
1901       elseif(cvar.eq.'rapmi')then
1902         inom=22
1903       elseif(cvar.eq.'eta')then
1904         inom=23
1905       elseif(cvar.eq.'theta')then
1906         inom=24
1907       elseif(cvar.eq.'pt2')then
1908         inom=25
1909       elseif(cvar.eq.'et')then
1910         inom=26
1911       elseif(cvar.eq.'idptl')then
1912         inom=27
1913       elseif(cvar.eq.'istptl')then
1914         inom=28
1915       elseif(cvar.eq.'mass')then
1916         inom=29
1917       elseif(cvar.eq.'idaptl')then
1918         inom=30
1919       elseif(cvar.eq.'egy')then
1920         inom=31
1921       elseif(cvar.eq.'rapwro')then
1922         inom=32
1923       elseif(cvar.eq.'mt')then
1924         inom=33
1925       elseif(cvar.eq.'pplus')then
1926         inom=34
1927       elseif(cvar.eq.'pminus')then
1928         inom=35
1929       elseif(cvar.eq.'p5')then
1930         inom=36
1931       elseif(cvar.eq.'pa')then
1932         inom=37
1933       elseif(cvar.eq.'sob')then
1934         inom=38
1935       elseif(cvar.eq.'idpom')then
1936         inom=39
1937       elseif(cvar.eq.'p3a')then
1938         inom=40
1939       elseif(cvar.eq.'cmass')then
1940         inom=41
1941       elseif(cvar.eq.'arappi')then
1942         inom=42
1943       elseif(cvar.eq.'itsptl')then
1944         inom=50
1945       elseif(cvar.eq.'ityptl')then
1946         inom=51
1947       elseif(cvar.eq.'idoptl')then
1948         inom=52
1949       elseif(cvar.eq.'iptl')then
1950         inom=53
1951       elseif(cvar.eq.'index')then
1952         inom=54
1953       elseif(cvar.eq.'p1')then
1954         inom=55
1955       elseif(cvar.eq.'p2')then
1956         inom=56
1957       elseif(cvar.eq.'p3')then
1958         inom=57
1959       elseif(cvar.eq.'p4')then
1960         inom=58
1961       elseif(cvar.eq.'xg')then
1962         inom=59
1963       elseif(cvar.eq.'ek')then
1964         inom=60
1965       elseif(cvar.eq.'beta')then
1966         inom=61
1967       elseif(cvar.eq.'mt0')then
1968         inom=63
1969       elseif(cvar.eq.'qsqptl')then
1970         inom=64
1971       elseif(cvar.eq.'xelab')then
1972         inom=65
1973       elseif(cvar.eq.'hgtc05')then
1974         inom=66
1975         imulty1=1             !to switch on the calculation of "Standard variable"
1976       elseif(cvar.eq.'hadtyp')then
1977         inom=67
1978         imulty1=1
1979       elseif(cvar.eq.'hgtc1')then
1980         inom=68
1981         imulty1=1
1982       elseif(cvar.eq.'x4')then
1983         inom=69
1984       elseif(cvar.eq.'npn')then
1985         inom=70
1986       elseif(cvar.eq.'routp')then
1987         inom=71
1988       elseif(cvar.eq.'hgtc3')then
1989         inom=72
1990         imulty1=1
1991       elseif(cvar.eq.'mu14')then
1992         inom=73
1993         imulty1=1
1994       elseif(cvar.eq.'delphi')then
1995         inom=74
1996         iok=0
1997         !------------------------------------------------------------
1998         !icorrtrig stores the histogram numbers of those histograms which
1999         !use the delphi variable (and therfore require a call corrtrig
2000         !------------------------------------------------------------
2001         do i=1,icorrtrig(0)
2002          if(icorrtrig(i).eq.n)iok=1
2003         enddo
2004         if(iok.eq.0)then
2005           icorrtrig(0)=icorrtrig(0)+1
2006           if(icorrtrig(0).gt.mxxhis)stop'mxxhis too small'
2007         icorrtrig(icorrtrig(0))=n
2008         endif
2009       elseif(cvar.eq.'v2')then
2010         inom=75
2011       elseif(cvar.eq.'pt4')then
2012         inom=76
2013       elseif(cvar.eq.'rin')then
2014         inom=77
2015       elseif(cvar.eq.'theh1p')then
2016         inom=78
2017       elseif(cvar.eq.'theh1t')then
2018         inom=79
2019       elseif(cvar.eq.'phi')then
2020         inom=80
2021       elseif(cvar.eq.'isoft')then
2022         inom=81
2023       elseif(cvar.eq.'mux')then
2024         inom=82
2025         imux(0)=imux(0)+1
2026         imux(imux(0))=n
2027       elseif(cvar.eq.'v4')then
2028         inom=83
2029       elseif(cvar.eq.'x3')then
2030         inom=84
2031       elseif(cvar.eq.'jorptl')then
2032         inom=85
2033       elseif(cvar.eq.'ptlead')then
2034         inom=86
2035         iok=0
2036         !------------------------------------------------------------
2037         !icorrtrig stores the histogram numbers of those histograms which
2038         !use the ptlead variable (and therfore require a call corrtrig
2039         !------------------------------------------------------------
2040         do i=1,icorrtrig(0)
2041          if(icorrtrig(i).eq.n)iok=1
2042         enddo
2043         if(iok.eq.0)then
2044           icorrtrig(0)=icorrtrig(0)+1
2045           if(icorrtrig(0).gt.mxxhis)stop'mxxhis too small'
2046         icorrtrig(icorrtrig(0))=n
2047         endif
2048       elseif(cvar.eq.'mu25')then
2049         inom=87
2050         imulty1=1
2051       elseif(cvar.eq.'pai')then
2052         inom=88
2053         ipairs1=1
2054       elseif(cvar.eq.'co2')then
2055         inom=89
2056         ipairs1=1
2057       elseif(cvar.eq.'co3')then
2058         inom=90
2059         ipairs1=1
2060       elseif(cvar.eq.'rad')then
2061         inom=91
2062       elseif(cvar.eq.'abseta')then
2063         inom=92
2064       elseif(cvar.eq.'phiexp')then
2065         inom=93
2066       elseif(cvar.eq.'mu24')then
2067         inom=94
2068         imulty1=1
2069       elseif(cvar.eq.'mulevt')then
2070         inom=101
2071       elseif(cvar.eq.'etevt')then
2072         inom=102
2073       elseif(cvar.eq.'enevt')then
2074         inom=103
2075       elseif(cvar.eq.'ev6evt')then
2076         inom=104
2077       elseif(cvar.eq.'xenevt')then
2078         inom=105
2079       elseif(cvar.eq.'netevt')then
2080         inom=106
2081       elseif(cvar.eq.'ptevt')then
2082         inom=107
2083       elseif(cvar.eq.'pmxevt')then
2084         inom=108
2085       elseif(cvar.eq.'numevt')then
2086         inom=201
2087       elseif(cvar.eq.'egyevt')then
2088         inom=202
2089       elseif(cvar.eq.'bimevt')then
2090         inom=203
2091       elseif(cvar.eq.'xbjevt')then
2092         inom=204
2093       elseif(cvar.eq.'qsqevt')then
2094         inom=205
2095       elseif(cvar.eq.'yevt')then
2096         inom=206
2097       elseif(cvar.eq.'eloevt')then
2098         inom=207
2099       elseif(cvar.eq.'nd1evt')then
2100         inom=208
2101       elseif(cvar.eq.'nd2evt')then
2102         inom=209
2103       elseif(cvar.eq.'theevt')then
2104         inom=210
2105       elseif(cvar.eq.'nspevt')then
2106         inom=211
2107       elseif(cvar.eq.'nhpevt')then
2108         inom=212
2109       elseif(cvar.eq.'sigtot')then
2110         inom=213
2111       elseif(cvar.eq.'sigela')then
2112         inom=214
2113       elseif(cvar.eq.'sloela')then
2114         inom=215
2115       elseif(cvar.eq.'nrgevt')then
2116         inom=216
2117       elseif(cvar.eq.'qevt')then
2118         inom=217
2119       elseif(cvar.eq.'qtlevt')then
2120         inom=218
2121       elseif(cvar.eq.'nd0evt')then
2122         inom=219
2123       elseif(cvar.eq.'threvt')then
2124         inom=220
2125         ifr=33                  !set thrust-frame
2126       elseif(cvar.eq.'omtevt')then
2127         inom=221
2128         ifr=33                  !set thrust-frame
2129       elseif(cvar.eq.'tmaevt')then
2130         inom=222
2131         ifr=33                  !set thrust-frame
2132       elseif(cvar.eq.'tmievt')then
2133         inom=223
2134         ifr=33                  !set thrust-frame
2135       elseif(cvar.eq.'oblevt')then
2136         inom=224
2137         ifr=33                  !set thrust-frame
2138       elseif(cvar.eq.'sphevt')then
2139         inom=230
2140         ifr=32                  !set sph-frame
2141       elseif(cvar.eq.'aplevt')then
2142         inom=231
2143         ifr=32                  !set sph-frame
2144       elseif(cvar.eq.'cpaevt')then
2145         inom=232
2146         ifr=34                  !set sph2-frame
2147       elseif(cvar.eq.'dpaevt')then
2148         inom=233
2149         ifr=34                  !set sph2-frame
2150       elseif(cvar.eq.'npoevt')then
2151         inom=234
2152       elseif(cvar.eq.'npnevt')then
2153         inom=235
2154       elseif(cvar.eq.'ikoevt')then
2155         inom=236
2156       elseif(cvar.eq.'iktevt')then
2157         inom=237
2158       elseif(cvar.eq.'npxevt')then
2159         inom=238
2160       elseif(cvar.eq.'nd6evt')then
2161         inom=239
2162       elseif(cvar.eq.'mu1evt')then
2163         inom=240
2164         imulty1=1
2165       elseif(cvar.eq.'muievt')then
2166         inom=241
2167         imulty1=1
2168       elseif(cvar.eq.'hgtevt')then
2169         inom=242
2170         imulty1=1
2171       elseif(cvar.eq.'difevt')then
2172         inom=243
2173       elseif(cvar.eq.'dixevt')then
2174         inom=244
2175       elseif(cvar.eq.'nd7evt')then
2176         inom=245
2177       elseif(cvar.eq.'nd8evt')then
2178         inom=246
2179       elseif(cvar.eq.'nd9evt')then
2180         inom=247
2181       elseif(cvar.eq.'ndaevt')then
2182         inom=248
2183       elseif(cvar.eq.'ndbevt')then
2184         inom=249
2185       elseif(cvar.eq.'qinevt')then
2186         inom=250
2187       elseif(cvar.eq.'qfievt')then
2188         inom=251
2189       elseif(cvar.eq.'einevt')then
2190         inom=252
2191       elseif(cvar.eq.'efievt')then
2192         inom=253
2193       elseif(cvar.eq.'pinevt')then
2194         inom=254
2195       elseif(cvar.eq.'pfievt')then
2196         inom=255
2197       elseif(cvar.eq.'pxfevt')then    ! leading proton xf in cms
2198         inom=256
2199       elseif(cvar.eq.'pi+xf')then     ! pi+xf: pi+ yield at cms xf>0.01
2200         inom=257
2201       elseif(cvar.eq.'pi-xf')then     ! pi-xf: pi- yield at cms xf>0.01
2202         inom=258
2203       elseif(cvar.eq.'sigcut')then
2204         inom=260
2205       elseif(cvar.eq.'keu')then
2206         inom=261
2207       elseif(cvar.eq.'ked')then
2208         inom=262
2209       elseif(cvar.eq.'kes')then
2210         inom=263
2211       elseif(cvar.eq.'kolevt')then
2212         inom=265
2213       elseif(cvar.eq.'sigsd')then
2214         inom=266
2215       elseif(cvar.eq.'nglevt')then
2216         inom=267
2217       elseif(cvar.eq.'kppevt')then   ! collision numbers per participant
2218         inom=268
2219       elseif(cvar.eq.'npievt')then   ! pion + multiplicity per event
2220         inom=269
2221       elseif(cvar.eq.'np2evt')then   ! pion + multiplicity per participant
2222         inom=270
2223       elseif(cvar.eq.'sigdif'.or.cvar.eq.'sigdifr')then
2224         inom=271
2225       elseif(cvar.eq.'koievt')then
2226         inom=272
2227       elseif(cvar.eq.'ineevt')then
2228         inom=273
2229       elseif(cvar.eq.'elaevt')then
2230         inom=274
2231       elseif(cvar.eq.'itgevt')then
2232         inom=275
2233         iok=0
2234         do i=1,icorrtrig(0)
2235           if(icorrtrig(i).eq.n)iok=1
2236         enddo
2237         if(iok.eq.0)then
2238           icorrtrig(0)=icorrtrig(0)+1
2239           if(icorrtrig(0).gt.mxxhis)stop'mxxhis too small'
2240           icorrtrig(icorrtrig(0))=n
2241         endif
2242       elseif(cvar.eq.'hrdevt')then
2243         inom=276
2244         iok=0
2245         do i=1,ihardevent(0)
2246           if(ihardevent(i).eq.n)iok=1
2247         enddo
2248         if(iok.eq.0)then
2249           ihardevent(0)=ihardevent(0)+1
2250           if(ihardevent(0).gt.mxxhis)stop'mxxhis too small'
2251           ihardevent(ihardevent(0))=n
2252         endif
2253       elseif(cvar(2:6).eq.'j1evt'.or.cvar(2:6).eq.'j2evt')then
2254         iok=0
2255         do i=1,ijetfind1(0)
2256           if(ijetfind1(i).eq.n)iok=1
2257         enddo
2258         if(iok.eq.0)then
2259           ijetfind1(0)=ijetfind1(0)+1
2260           if(ijetfind1(0).gt.mxxhis)stop'mxxhis too small'
2261           ijetfind1(ijetfind1(0))=n
2262         endif
2263         if(cvar.eq.'ej1evt')inom=277
2264         if(cvar.eq.'pj1evt')inom=278
2265         if(cvar(2:6).eq.'j2evt')then
2266           iok=0
2267           do i=1,ijetfind2(0)
2268             if(ijetfind2(i).eq.n)iok=1
2269           enddo
2270           if(iok.eq.0)then
2271             ijetfind2(0)=ijetfind2(0)+1
2272             if(ijetfind2(0).gt.mxxhis)stop'mxxhis too small'
2273             ijetfind2(ijetfind2(0))=n
2274           endif
2275           if(cvar.eq.'ej2evt')inom=279
2276           if(cvar.eq.'pj2evt')inom=280
2277         endif
2278       elseif(cvar.eq.'zppevt')then
2279         inom=281
2280       elseif(cvar.eq.'zptevt')then
2281         inom=282
2282       elseif(cvar.eq.'***not used***')then
2283         inom=283
2284       elseif(cvar.eq.'nd3evt')then
2285         inom=284
2286       elseif(cvar.eq.'nd4evt')then
2287         inom=285
2288       elseif(cvar.eq.'mubevt')then
2289         inom=286
2290         imulty1=1
2291       elseif(cvar.eq.'nd5evt')then
2292         inom=287
2293       elseif(cvar.eq.'ekievt')then
2294         inom=288
2295       elseif(cvar.eq.'sd1evt')then
2296         inom=289
2297       elseif(cvar.eq.'sd2evt')then
2298         inom=290
2299       elseif(cvar.eq.'mdevt')then
2300         inom=291
2301         imulty1=1     !to switch on the calculation of "Standard variable"
2302       elseif(cvar.eq.'m2devt')then
2303         inom=292
2304         imulty1=1
2305       elseif(cvar.eq.'tdevt')then
2306         inom=293
2307         imulty1=1
2308       elseif(cvar.eq.'ndpevt')then
2309         inom=294
2310       elseif(cvar.eq.'rapgap')then
2311         inom=295
2312         imulty1=1
2313       elseif(cvar.eq.'ng1evt')then
2314         inom=296
2315       elseif(cvar.eq.'r21evt')then
2316         inom=297
2317       elseif(cvar.eq.'aimevt')then
2318         inom=301
2319       elseif(cvar.eq.'wjbevt')then
2320         inom=302
2321       elseif(cvar.eq.'njbevt')then
2322         inom=303
2323       elseif(cvar.eq.'djbevt')then
2324         inom=304
2325       elseif(cvar.eq.'tjbevt')then
2326         inom=305
2327       elseif(cvar.eq.'hjmevt')then
2328         inom=306
2329       elseif(cvar.eq.'ljmevt')then
2330         inom=307
2331       elseif(cvar.eq.'djmevt')then
2332         inom=308
2333       elseif(cvar.eq.'ybal')then
2334         inom=310
2335       elseif(cvar.eq.'yabal')then
2336         inom=310
2337       elseif(cvar.eq.'sigine')then
2338         inom=312
2339       elseif(cvar.eq.'sigiaa')then
2340         inom=313
2341       elseif(cvar.eq.'alpdsf')then
2342         inom=314
2343       elseif(cvar.eq.'alpdsh')then
2344         inom=315
2345       elseif(cvar.eq.'betdsf')then
2346         inom=316
2347       elseif(cvar.eq.'betdsh')then
2348         inom=317
2349       elseif(cvar.eq.'rexdip')then
2350         inom=318
2351       elseif(cvar.eq.'rexdit')then
2352         inom=319
2353       elseif(cvar.eq.'m14evt')then
2354         inom=320
2355         imulty1=1
2356       elseif(cvar.eq.'ht3evt')then
2357         inom=321
2358       elseif(cvar.eq.'sigiex')then
2359         inom=322
2360       elseif(cvar.eq.'sigdex')then
2361         inom=323
2362       elseif(cvar.eq.'sigsex')then
2363         inom=324
2364       elseif(cvar.eq.'ekievt')then
2365         inom=325
2366       elseif(cvar.eq.'sigcaa')then
2367         inom=326
2368       elseif(cvar.eq.'sigtaa')then
2369         inom=327
2370       elseif(cvar.eq.'xkappa')then
2371         inom=328
2372       elseif(cvar.eq.'gamdsf')then
2373         inom=329
2374       elseif(cvar.eq.'gamdsh')then
2375         inom=330
2376       elseif(cvar.eq.'deldsf')then
2377         inom=331
2378       elseif(cvar.eq.'deldsh')then
2379         inom=332
2380       elseif(cvar.eq.'nd6evt')then
2381         inom=333
2382       elseif(cvar.eq.'muxevt')then
2383         inom=334
2384         imux(0)=imux(0)+1
2385         imux(imux(0))=n
2386       elseif(cvar.eq.'typevt')then
2387         inom=335
2388       elseif(cvar.eq.'m25evt')then
2389         inom=339
2390         imulty1=1
2391       elseif(cvar.eq.'segevt')then
2392         inom=340
2393       elseif(cvar.eq.'ielevt')then
2394         inom=341
2395       elseif(cvar.eq.'mc1evt')then
2396         inom=342
2397         imulty1=1
2398       elseif(cvar.eq.'sdcdf')then
2399         inom=343
2400       elseif(cvar.eq.'dpecdf')then
2401         inom=344
2402       elseif(cvar.eq.'ddcdf')then
2403         inom=345
2404       elseif(cvar.eq.'phievt')then
2405         inom=346
2406       elseif(cvar.eq.'ndcevt')then
2407         inom=347
2408       elseif(cvar.eq.'jetevt')then
2409         inom=348
2410 c activate jetevent trigger
2411         iok=0
2412         do i=1,ijetevent(0)
2413           if(ijetevent(i).eq.n)iok=1
2414         enddo
2415         if(iok.eq.0)then
2416           ijetevent(0)=ijetevent(0)+1
2417           if(ijetevent(0).gt.mxxhis)stop'mxxhis too small'
2418           ijetevent(ijetevent(0))=n
2419         endif
2420 c activate fastjet
2421         iok=0
2422         do i=1,ifastjet(0)
2423           if(ifastjet(i).eq.n)iok=1
2424         enddo
2425         if(iok.eq.0)then
2426           ifastjet(0)=ifastjet(0)+1
2427           if(ifastjet(0).gt.mxxhis)stop'mxxhis too small'
2428           ifastjet(ifastjet(0))=n
2429         endif
2430       elseif(cvar.eq.'epszer')then
2431         inom=349
2432       elseif(cvar.eq.'xsievt')then
2433         inom=350
2434         imulty1=1
2435       elseif(cvar.eq.'xsicms')then
2436         inom=351
2437       elseif(cvar.eq.'calevt')then
2438         inom=352
2439         icaltrig(0)=icaltrig(0)+1
2440         icaltrig(icaltrig(0))=n
2441       elseif(cvar.eq.'fgpevt')then
2442         inom=353
2443         icaltrig(0)=icaltrig(0)+1
2444         icaltrig(icaltrig(0))=n
2445       elseif(cvar.eq.'bgpevt')then
2446         inom=354
2447         icaltrig(0)=icaltrig(0)+1
2448         icaltrig(icaltrig(0))=n
2449       elseif(cvar.eq.'gapevt')then
2450         inom=355
2451         icaltrig(0)=icaltrig(0)+1
2452         icaltrig(icaltrig(0))=n
2453       elseif(cvar.eq.'sigdd')then
2454         inom=356
2455       elseif(cvar.eq.'styevt')then
2456         inom=357
2457       elseif(cvar.eq.'ndsevt')then
2458         inom=358
2459       elseif(cvar.eq.'m24evt')then
2460         inom=359
2461         imulty1=1
2462       elseif(cvar.eq.'ndhevt')then
2463         inom=360
2464       elseif(cvar.eq.'ox1evt')then
2465         inom=501
2466       elseif(cvar.eq.'ox2evt')then
2467         inom=502
2468       elseif(cvar.eq.'ox3evt')then
2469         inom=503
2470       elseif(cvar.eq.'ox4evt')then
2471         inom=504
2472       elseif(cvar.eq.'eglevt')then  ! eccentricity
2473         inom=505
2474       elseif(cvar.eq.'fglevt')then  ! eccentricity_part
2475         inom=506
2476       elseif(cvar.eq.'rglevt')then  ! ratio ng2 / ng1
2477         inom=507
2478       elseif(cvar.eq.'sglevt')then  ! area S
2479         inom=508
2480       elseif(cvar.eq.'ptrevt')then
2481         inom=509
2482       elseif(cvar.eq.'rr2evt')then
2483         inom=510
2484         imulty1=1             !to switch on the calculation of "Standard variable"
2485       elseif(cvar.eq.'perevt')then
2486         inom=511
2487       elseif(cvar.eq.'paievt')then
2488         inom=512
2489         ipairs1=1
2490       elseif(cvar.eq.'co2evt')then
2491         inom=513
2492         ipairs1=1
2493       elseif(cvar.eq.'co3evt')then
2494         inom=514
2495         ipairs1=1
2496       else
2497         print *,' '
2498         print *,'              xtrans: unknown variable ',cvar
2499         print *,' '
2500 c       inom=-1
2501         stop
2502       endif
2503       end
2504 
2505 c----------------------------------------------------------------------
2506       subroutine xval(n,inom,lf,j,x)
2507 c----------------------------------------------------------------------
2508 c   n ...... histogram index
2509 c   inom ... variable index
2510 c              1-100 particle variables
2511 c              101-200 accumulative event variables
2512 c              > 200 other event variables
2513 c   lf ..... frame index
2514 c   particle index (used for particle variables)
2515 c----------------------------------------------------------------------
2516       include 'epos.inc'
2517       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
2518      &     ,multc1,multc83,multc24,multc25,rapgap,ipairs1,xsi
2519       parameter(mxxhis=70)
2520       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
2521      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis),imux(0:mxxhis)
2522      &,ifastjet(0:mxxhis),ijetevent(0:mxxhis),icaltrig(0:mxxhis)
2523       common/zeus2/qtl
2524 
2525       parameter (ntim=1000)
2526       common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
2527      &,iorprt(ntim),jorprt(ntim),nprtj
2528 
2529       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
2530      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
2531       common/cpairs/paievt,co2evt,co3evt
2532 
2533       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
2534       parameter (mypara=10,mxpara=10)
2535       logical ilog,icnx,itrevt,idmod
2536       double precision bin,bbin,zcbin,zbbin
2537       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
2538      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
2539      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
2540      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
2541      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
2542      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
2543      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
2544      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
2545      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
2546      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
2547      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
2548       double precision ebin,zebin
2549       common/errbins/ebin(mxbin,2,mxhis/2),zebin(mxbin,2,mxhis/2),
2550      $inoerr(mxhis),noerr(mxhis/2,2),noerrhis(mxhis/2),noerrall
2551       parameter (mxfra=5)
2552       common/pfra/nfra,ifra(mxfra),ivfra(2,mxhis),itfra(mxtri,mxhis)
2553      $     ,imofra(3,mxfra),iffra(mxfra),r1fra(3,mxfra),r2fra(3,mxfra)
2554      $     ,emax(mxfra)
2555       common/cphi2/phi2pos,phi2neg
2556       common/cen/ncentr  
2557       parameter(nbkbin=100)
2558       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
2559       common/cpairs2/paipi(40),co2pi(40),co3pi(40)
2560      .              ,paipr(40),co2pr(40),co3pr(40)
2561      .            ,ipaipi(40),ico2pi(40),ico3pi(40)
2562      .            ,ipaipr(40),ico2pr(40),ico3pr(40)
2563      .            ,maxpt,delpt
2564 
2565       double precision bofra,bofra1,bofra2,bofra3,bofra4,bofra5,xd
2566       common/dfra/bofra(5,mxfra)
2567       dimension p(5,mxfra),aimuni(10,mxhis),xor(5,mxfra)
2568       common/cranphi/ranphi
2569       save p,aimuni,xor
2570       phinll=phievt+ranphi
2571       if(phinll.lt.-pi)phinll=phinll+2*pi
2572       if(phinll.gt.pi)phinll=phinll-2*pi
2573       if(iffra(lf).eq.0.and.j.ne.0)then
2574         do l=1,5
2575           p(l,lf)=pptl(l,j)
2576         enddo
2577         do l=1,4
2578           xor(l,lf)=xorptl(l,j)
2579         enddo
2580         if(imofra(1,lf).ne.0)then
2581           call utrota(imofra(1,lf),r1fra(1,lf),r1fra(2,lf),r1fra(3,lf)
2582      $         ,p(1,lf),p(2,lf),p(3,lf))
2583           call utrota(imofra(1,lf),r1fra(1,lf),r1fra(2,lf),r1fra(3,lf)
2584      $         ,xor(1,lf),xor(2,lf),xor(3,lf))
2585         endif
2586         if(imofra(2,lf).ne.0)then !the x-z exchanged is ok !!
2587           call utrota(imofra(2,lf),r2fra(3,lf),r2fra(2,lf),r2fra(1,lf)
2588      $         ,p(3,lf),p(2,lf),p(1,lf))
2589           call utrota(imofra(2,lf),r2fra(3,lf),r2fra(2,lf),r2fra(1,lf)
2590      $         ,xor(3,lf),xor(2,lf),xor(1,lf))
2591         endif
2592         if(imofra(3,lf).ne.0)then
2593           imof3=sign(1,imofra(3,lf))
2594           if(abs(imofra(3,lf)).gt.1)then
2595             bofra1=0d0
2596             bofra2=0d0
2597             bofra5=1d0
2598             call utlob5(imof3*yhaha
2599      $                 ,p(1,lf),p(2,lf),p(3,lf),p(4,lf),p(5,lf))
2600             bofra3=bofra(1,lf)
2601             bofra4=bofra(2,lf)
2602             call utlob4(imof3,bofra1,bofra2,bofra3,bofra4,bofra5
2603      $         ,xor(1,lf),xor(2,lf),xor(3,lf),xor(4,lf))
2604             bofra3=bofra(3,lf)
2605             bofra4=bofra(4,lf)
2606             bofra5=bofra(5,lf)
2607           else
2608             bofra1=bofra(1,lf)
2609             bofra2=bofra(2,lf)
2610             bofra3=bofra(3,lf)
2611             bofra4=bofra(4,lf)
2612             bofra5=bofra(5,lf)
2613           endif
2614           call utlob4(imof3,bofra1,bofra2,bofra3,bofra4,bofra5
2615      $         ,p(1,lf),p(2,lf),p(3,lf),p(4,lf))
2616           call utlob4(imof3,bofra1,bofra2,bofra3,bofra4,bofra5
2617      $         ,xor(1,lf),xor(2,lf),xor(3,lf),xor(4,lf))
2618         endif
2619         iffra(lf)=1
2620       endif
2621 
2622 c--------------------------------- 1 - 100 ----------------------------
2623       if(inom.eq.1)then
2624         x=1.
2625       elseif(inom.eq.2)then
2626         x=isign(1,idptl(j))
2627       elseif(inom.eq.3)then
2628         chrg=0
2629         if(iabs(idptl(j)).le.9999
2630      $       .and.mod(iabs(idptl(j)),10).le.1)
2631      $       call idchrg(idptl(j),chrg)
2632         if(chrg.eq.0.)then
2633           x=0
2634         else
2635           x=int(sign(1.,chrg))
2636         endif
2637       elseif(inom.eq.4)then
2638         iad=abs(idptl(j))
2639         jspin=mod(iad,10)
2640         x=0.
2641         if (iad.ge.100.and.iad.lt.1000) x=1./(1.+2.*jspin)
2642         if (iad.ge.1000.and.iad.lt.9999) x=1./(2.+2*jspin)
2643       elseif(inom.eq.5)then    !'rapx'  !st-rap for string segments only !!!!!!!!!!
2644         x=dezptl(j)
2645       elseif(inom.eq.6)then                                      !'iptlfr'
2646         x=0
2647         if(j.ge.minfra.and.j.le.maxfra)x=1
2648       elseif(inom.eq.7)then                                        !'rinp'
2649         aa=cos(phinll)
2650         bb=sin(phinll)
2651         x=xptl(j)*aa+yptl(j)*bb
2652       elseif(inom.eq.8)then                        !'eco' !engy in comoving frame
2653         x=0
2654         amt=p(5,lf)**2+p(1,lf)**2+p(2,lf)**2
2655         if(amt.gt.0..and.p(4,lf)+abs(p(3,lf)).gt.0.d0)then
2656           amt=sqrt(amt)
2657           rap=sign(1.,p(3,lf))*alog((p(4,lf)+abs(p(3,lf)))/amt)
2658           rapx=dezptl(j)
2659           x=amt*cosh(rap-rapx)
2660         endif
2661       elseif(inom.eq.9)then                                       !'tau'
2662         x=-999999
2663         !if(iorptl(j).ne.0)then
2664         ! jo=iorptl(j)
2665          dt=xorptl(4,j)  !-xorptl(4,jo)
2666          dz=xorptl(3,j)  !-xorptl(3,jo)
2667          x2=dt**2-dz**2
2668          if(x2.gt.0.)x=sqrt(x2)
2669         !endif
2670       elseif(inom.eq.10)then                                       !'ctr'
2671 c        x=ctrevt
2672       elseif(inom.eq.11)then                                       !'v2np'
2673         phi=polar( p(1,lf) , p(2,lf) )
2674         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2675         eta=0
2676         if(p(3,lf).ne.0..and.pt.ne.0.)eta=sign(1.,p(3,lf))*
2677      *       alog((sqrt(p(3,lf)**2+pt**2)+abs(p(3,lf)))/pt)
2678         if(eta.gt.0)then
2679         phi2=phi2neg
2680         else
2681         phi2=phi2pos
2682         endif
2683         x=cos(2*(phi-phi2))
2684       elseif(inom.eq.12)then                                      !'absrap'
2685         amt=p(5,lf)**2+p(1,lf)**2+p(2,lf)**2
2686         if(amt.gt.0..and.p(4,lf)+abs(p(3,lf)).gt.0.d0)then
2687           amt=sqrt(amt)
2688           x=alog((p(4,lf)+abs(p(3,lf)))/amt)
2689         else
2690           x=0.                  !
2691         endif
2692       elseif(inom.eq.13)then    !'rap'
2693         amt=p(5,lf)**2+p(1,lf)**2+p(2,lf)**2
2694         if(amt.gt.0..and.p(4,lf)+abs(p(3,lf)).gt.0.d0)then  !not correct if particles off-shell
2695 c        if(amt.gt.0..and.p(4,lf)+p(3,lf).gt.0.d0)then    !not correct and assymetric if particles off-shell
2696           amt=sqrt(amt)
2697           x=sign(1.,p(3,lf))*log((p(4,lf)+abs(p(3,lf)))/amt)  !not correct if particles off-shell
2698 c          x=alog((p(4,lf)+p(3,lf))/amt)    !not correct and assymetric if particles off-shell
2699 c          x=0.5*alog((p(4,lf)+p(3,lf))/(p(4,lf)-p(3,lf)))  !always correct but numerically unstable
2700 c          if(abs(x).lt.0.05.and.idptl(j).eq.120)
2701 c     &    write(ifch,*)'pion epo',p(4,lf),p(3,lf),amt,x
2702         else
2703           x=0.                  !
2704         endif
2705       elseif(inom.eq.14)then                                         !'xp'
2706         x=sqrt(p(3,lf)**2+p(2,lf)**2+p(1,lf)**2)/emax(lf)
2707       elseif(inom.eq.15)then                                         !'xe'
2708         x=min(1.,p(4,lf)/emax(lf))
2709       elseif(inom.eq.16)then                                         !'pt'
2710         x=sqrt(p(2,lf)**2+p(1,lf)**2)
2711       elseif(inom.eq.17)then
2712         x=abs(p(1,lf))
2713       elseif(inom.eq.18)then
2714         x=abs(p(2,lf))
2715       elseif(inom.eq.19)then
2716         x=-log(sqrt(p(3,lf)**2+p(2,lf)**2+p(1,lf)**2)/emax(lf))
2717       elseif(inom.eq.20)then                                     !'xf'
2718         m=mod(ifra(lf)/10,10)
2719         if(m.eq.1.or.noebin.lt.0)then
2720 c          pmax=sqrt((engy/2)**2-prom*2)
2721           pmax=pnullx               !???????????????????
2722           if(mod(ifra(lf),10).eq.2)pmax=pnll
2723           x=p(3,lf)/pmax
2724         else
2725           x=p(3,lf)/emax(lf)
2726         endif
2727 c        if(x.gt.0.95.and.idptl(j).eq.1220)then
2728 c          write(ifch,'(a,d25.15)')'ici !!!!!!!!!',seedc
2729 c          stop
2730 c        endif
2731       elseif(inom.eq.21)then
2732 c        pmax=pmxevt
2733 c        pmax=sqrt((engy/2)**2-prom*2)
2734         pmax=pnullx             !???????????????????
2735         if(mod(ifra(lf),10).eq.2)pmax=pnll
2736         x=-(amproj**2-2.*sqrt(amproj**2+pmax**2)*p(4,lf)
2737      *      +2.*abs(pmax*p(3,lf))+p(5,lf)**2)
2738       elseif(inom.eq.22)then
2739         amt=sqrt(p(5,lf)**2+p(1,lf)**2+p(2,lf)**2)
2740         if(amt.ne.0.)then
2741           x=-sign(1.,p(3,lf))*alog((p(4,lf)+abs(p(3,lf)))/amt)
2742         else
2743           x=0.                  !
2744         endif
2745       elseif(inom.eq.23)then                                     !'eta'
2746         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2747         if(p(3,lf).eq.0.)then
2748           x=0.
2749         elseif(pt.ne.0.)then
2750           x=sign(1.,p(3,lf))*
2751      *       alog((sqrt(p(3,lf)**2+pt**2)+abs(p(3,lf)))/pt)
2752         else
2753           x=sign(1000.,p(3,lf))
2754         endif
2755       elseif(inom.eq.24)then                                     !'theta (deg)'
2756         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2757         x=90
2758         if(p(3,lf).ne.0.)x=atan(pt/p(3,lf))/pi*180.
2759         if(x.lt.0.)x=180.+x
2760       elseif(inom.eq.25)then                                     !'pt2'
2761         x=p(2,lf)**2+p(1,lf)**2
2762       elseif(inom.eq.26)then                                     !'et'
2763         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2764         x=0
2765         eef=p(4,lf)
2766 c        if(idptl(j).ge.1000)eef=eef-prom
2767 c        if(idptl(j).le.-1000)eef=eef+prom
2768         p2=p(3,lf)**2+p(2,lf)**2+p(1,lf)**2
2769         if(p2.ne.0.)x=eef*pt/sqrt(p2)
2770       elseif(inom.eq.27)then                                     !'idptl'
2771         x=idptl(j)
2772       elseif(inom.eq.28)then    !istptl
2773         x=istptl(j)
2774       elseif(inom.eq.29)then    !mass
2775         x=p(5,lf)
2776         if(istptl(j).le.1)call idmass(idptl(j),x)
2777       elseif(inom.eq.30)then    !idaptl
2778         x=abs(idptl(j))
2779       elseif(inom.eq.31)then    !egy
2780         x=egyevt
2781       elseif(inom.eq.32)then    !arapwro
2782         x=0
2783         pt2=p(2,lf)**2+p(1,lf)**2
2784         if(p(3,lf).ne.0.)x=sign(1.,p(3,lf))*
2785      *       alog((sqrt(p(3,lf)**2+pt2+.13957**2)+abs(p(3,lf)))
2786      *       /sqrt(pt2+.13957**2))
2787       elseif(inom.eq.33)then                                  !'mt'
2788         x=sqrt(p(2,lf)**2+p(1,lf)**2+p(5,lf)**2)
2789       elseif(inom.eq.34)then                                  !'pplus'
2790         x=sign(1.,p(3,lf)) * (p(4,lf)+abs(p(3,lf)))
2791       elseif(inom.eq.35)then                                  !'pminus'
2792         x=sign(1.,p(3,lf)) * (p(4,lf)-abs(p(3,lf)))
2793       elseif(inom.eq.36)then                                  !'p5' (mass)
2794         x=p(5,lf)
2795       elseif(inom.eq.37)then                                  !pa
2796         x=sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
2797       elseif(inom.eq.38)then                                  !'pa'
2798         if(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2.ne.0)
2799      *       x=egyevt**2/sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)*p(4,lf)
2800       elseif(inom.eq.39)then                                  !idpom
2801         x=idptl(j)/1000000
2802       elseif(inom.eq.40)then                                  !p3a
2803         x=abs(p(3,lf))
2804       elseif(inom.eq.41)then
2805         cm2=p(4,lf)**2-p(3,lf)**2-p(2,lf)**2-p(1,lf)**2         !cmass
2806         x=sign(sqrt(abs(cm2)),cm2)
2807       elseif(inom.eq.42)then    !arappi
2808         x=0
2809         pt2=p(2,lf)**2+p(1,lf)**2
2810         if(p(3,lf).ne.0.)
2811      *       x=alog((sqrt(p(3,lf)**2+pt2+.13957**2)+abs(p(3,lf)))
2812      *       /sqrt(pt2+.13957**2))
2813       elseif(inom.eq.50)then
2814         x=itsptl(j)
2815       elseif(inom.eq.51)then
2816         x=ityptl(j)
2817       elseif(inom.eq.52)then                      !'idoptl'
2818         x=0.
2819         if(iorptl(j).ne.0) x=idptl(iorptl(j))
2820       elseif(inom.eq.53)then
2821         x=j
2822       elseif(inom.eq.54)then                       !'sloela'
2823         call idflav(idptl(j),ifl1,ifl2,ifl3,jspin,indx)
2824         x=indx
2825       elseif(inom.eq.55)then                       !'p1'
2826         x=p(1,lf)
2827       elseif(inom.eq.56)then                       !'p2'
2828         x=p(2,lf)
2829       elseif(inom.eq.57)then                       !'p3'
2830         x=p(3,lf)
2831       elseif(inom.eq.58)then                       !'p4'
2832         x=p(4,lf)
2833       elseif(inom.eq.59)then                       !E/p_max
2834 c        pmax=sqrt((engy/2)**2-prom*2)
2835           pmax=pnullx               !???????????????????
2836         if(mod(ifra(lf),10).eq.2)pmax=pnll
2837         x=p(4,lf)/pmax
2838       elseif(inom.eq.60)then                       !'ek'
2839         x=p(4,lf)-p(5,lf)
2840       elseif(inom.eq.61)then                       !'beta'
2841         x=p(3,lf)/p(4,lf)
2842       elseif(inom.eq.63)then                       !'mt0'
2843         x=sqrt(p(2,lf)**2+p(1,lf)**2+p(5,lf)**2)-p(5,lf)
2844       elseif(inom.eq.64)then                       !qsqptl
2845         x=qsqptl(j)
2846       elseif(inom.eq.65)then                       !xelab=Elab/Eolab
2847         x=p(4,lf)/(ecms**2/2/prom-prom)
2848         if(x.gt.0.9999999) x=.9999999
2849       elseif(inom.eq.66)then    !'hgtc05' ... charged ptl mult |[c]|<0.5
2850         x=multc05
2851       elseif(inom.eq.67)then    !'hadtyp' ... primary (1) or secondary (2) hadron
2852         if(j.le.nbdky)then
2853           x=1
2854         else
2855           x=2
2856         endif
2857       elseif(inom.eq.68)then    !'hgtc1'
2858         x=multc1
2859       elseif(inom.eq.69)then                       !'x4'
2860         x=xor(4,lf)
2861       elseif(inom.eq.70)then                       !'npn'
2862         x=npjevt+ntgevt
2863       elseif(inom.eq.71)then                       !'routp'
2864         cc=-sin(phinll)
2865         dd=cos(phinll)
2866         x=xptl(j)*cc+yptl(j)*dd
2867       elseif(inom.eq.72)then    !'hgtc3' ... charged ptl mult |eta|<3.15  /6.3
2868         x=multc3/6.3
2869       elseif(inom.eq.73)then    !'mu14' ... charged ptl mult |eta|<1  pt>.4
2870         x=multc14
2871       elseif(inom.eq.74)then    !'delphi' ... azimuthhal correlation
2872         x=10000.
2873         pt=sqrt(p(1,lf)**2+p(2,lf)**2)
2874 
2875         if(nint(ypara(1,n)).ne.0.and.j.ne.nint(ypara(1,n)).and.
2876      $           pt.gt.0)then
2877            phi=sign(1.,p(2,lf))*acos(p(1,lf)/pt)
2878                x=phi-ypara(2,n)
2879            phiz= ypara(3,n)
2880            if   (x.lt.(-2+phiz)*pi)then
2881             x=x+4*pi
2882            elseif(x.lt.(0+phiz)*pi)then
2883             x=x+2*pi
2884            elseif(x.gt.(4+phiz)*pi)then
2885             x=x-4*pi
2886            elseif(x.gt.(2+phiz)*pi)then
2887             x=x-2*pi
2888            endif
2889         endif
2890       elseif(inom.eq.75)then    !'v2'
2891 c        if(iranphi.ne.1)stop'\n\n ERROR 29062010b \n\n'
2892         aa=cos(phinll)
2893         bb=sin(phinll)
2894         cc=-sin(phinll)
2895         dd=cos(phinll)
2896         px=p(1,lf)*aa+p(2,lf)*bb
2897         py=p(1,lf)*cc+p(2,lf)*dd
2898         pt2=p(2,lf)**2+p(1,lf)**2
2899         x=0
2900         if(pt2.gt.0.)x=(px**2-py**2)/pt2
2901       elseif(inom.eq.76)then                                     !'pt4'
2902         x=(p(2,lf)**2+p(1,lf)**2)**2
2903       elseif(inom.eq.77)then                                     !'rin'
2904         x=rinptl(j)
2905       elseif(inom.eq.78)then              !'theta for H1 (rad) for (proj side)'
2906         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2907         p1=p(1,lf)
2908         p2=p(2,lf)
2909         p3=p(3,lf)
2910         p4=p(4,lf)
2911         p5=p(5,lf)
2912         if(abs(p3).gt.1e-5)then
2913 c put the particle in the projectile frame
2914           call utlob5(yhaha,p1,p2,p3,p4,p5)
2915 c put the particle in a frame where the projectile (proton) has 820 GeV (HERA)
2916           call utlob4(-1,0d0,0d0,819.99946d0,820.d0,0.938d0,p1,p2,p3,p4)
2917           x=atan(pt/p3)
2918           if(x.lt.0.)x=pi+x
2919         else
2920           x=0.5*pi
2921         endif
2922       elseif(inom.eq.79)then             !'theta for H1 (rad) (for target side)'
2923         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2924         p1=p(1,lf)
2925         p2=p(2,lf)
2926         p3=p(3,lf)
2927         p4=p(4,lf)
2928         p5=p(5,lf)
2929         if(abs(p3).gt.1e-5)then
2930 c put the particle in the projectile frame
2931           call utlob5(-yhaha,p1,p2,p3,p4,p5)
2932 c put the particle in a frame where the projectile (proton) has 820 GeV (HERA)
2933           call utlob4(-1,0d0,0d0,-819.99946d0,820d0,0.938d0,p1,p2,p3,p4)
2934           x=atan(pt/p3)
2935           if(x.gt.0.)x=pi-x
2936         else
2937           x=0.5*pi
2938         endif
2939       elseif(inom.eq.80)then    !'phi'
2940         x=1000
2941         pt=sqrt(p(1,lf)**2+p(2,lf)**2)
2942         if(pt.gt.0.)then
2943            phi=sign(1.,p(2,lf))*acos(p(1,lf)/pt)
2944            x=phi-phinll
2945         endif
2946         if(x.lt.-pi)x=x+2*pi
2947         if(x.gt.pi)x=x-2*pi
2948       elseif(inom.eq.81)then  !'isoft'
2949         x=0
2950         it=ityptl(j)
2951         if(it.ge.20.and.it.le.29)x=1
2952         if(it.ge.40.and.it.le.60)x=1
2953       elseif(inom.eq.82)then  !'mux' ... charged ptl mult 
2954         x= ypara(1,n)        
2955       elseif(inom.eq.83)then    !'v4'
2956         aa=cos(phinll)
2957         bb=sin(phinll)
2958         cc=-sin(phinll)
2959         dd=cos(phinll)
2960         px=p(1,lf)*aa+p(2,lf)*bb
2961         py=p(1,lf)*cc+p(2,lf)*dd
2962         pt2=p(2,lf)**2+p(1,lf)**2
2963         x=0
2964         if(pt2.gt.0.)x=px**2/pt2   !cos**2
2965         x=8*x**2-8*x+1
2966       elseif(inom.eq.84)then                       !'x3'
2967         x=xor(3,lf)
2968       elseif(inom.eq.85)then   !jorptl
2969         x=jorptl(j)
2970       elseif(inom.eq.86)then    !'ptlead' ... pt of particle with higher pt
2971         x=0.
2972         if(nint(ypara(1,n)).ne.0)x=ypara(4,n)
2973       elseif(inom.eq.87)then    !'mu25' ... charged ptl mult |eta|<2.5  pt>.5
2974         x=multc25
2975       elseif(inom.eq.88)then    !'pai' 
2976         x=0
2977         if(bimevt.gt.xpara(1,n).and.bimevt.lt.xpara(2,n))then
2978           ida=idptl(j)
2979           pt=sqrt(p(2,lf)**2+p(1,lf)**2)
2980           ipt=pt/delpt+1
2981           if(ipt.ge.1.and.ipt.le.maxpt)then
2982             if(ida.eq.120)then
2983               if(ipaipi(ipt).eq.0)then 
2984                 x=paipi(ipt)
2985                 ipaipi(ipt)=1
2986               endif
2987             elseif(ida.eq.1120)then
2988               if(ipaipr(ipt).eq.0)then 
2989                 x=paipr(ipt)
2990                 !print*,'++++++++++',ipt,ipaipr(ipt),x
2991                 ipaipr(ipt)=1
2992               endif
2993             endif
2994           endif 
2995         endif 
2996       elseif(inom.eq.89)then    !'co2' 
2997         x=0
2998         if(bimevt.gt.xpara(1,n).and.bimevt.lt.xpara(2,n))then
2999           ida=idptl(j)
3000           pt=sqrt(p(2,lf)**2+p(1,lf)**2)
3001           ipt=pt/delpt+1
3002           if(ipt.ge.1.and.ipt.le.maxpt)then
3003             if(ida.eq.120)then
3004               if(ico2pi(ipt).eq.0)then 
3005                 x=co2pi(ipt)
3006                 ico2pi(ipt)=1
3007               endif
3008             elseif(ida.eq.1120)then
3009               if(ico2pr(ipt).eq.0)then 
3010                 x=co2pr(ipt)
3011                 ico2pr(ipt)=1
3012               endif
3013             endif
3014           endif  
3015         endif 
3016       elseif(inom.eq.90)then    !'co3' 
3017         x=0
3018         if(bimevt.gt.xpara(1,n).and.bimevt.lt.xpara(2,n))then
3019           ida=idptl(j)
3020           pt=sqrt(p(2,lf)**2+p(1,lf)**2)
3021           ipt=pt/delpt+1
3022           if(ipt.ge.1.and.ipt.le.maxpt)then
3023             if(ida.eq.120)then
3024               if(ico3pi(ipt).eq.0)then 
3025                 x=co3pi(ipt)
3026                 ico3pi(ipt)=1
3027               endif
3028             elseif(ida.eq.1120)then
3029               if(ico3pr(ipt).eq.0)then 
3030                 x=co3pr(ipt)
3031                 ico3pr(ipt)=1
3032               endif
3033             endif
3034           endif  
3035         endif  
3036       elseif(inom.eq.91)then                       !'rad'
3037         x=0.001       !unit is fm !
3038         x1=xor(1,lf)
3039         x2=xor(2,lf)
3040         xd=dble(x1)**2+dble(x2)**2
3041         if(xd.gt.0.d0.and.xd.eq.xd)x=sqrt(xd)
3042       elseif(inom.eq.92)then                                     !'abseta'
3043         pt=sqrt(p(2,lf)**2+p(1,lf)**2)
3044         if(p(3,lf).eq.0.)then
3045           x=0.
3046         elseif(pt.ne.0.)then
3047           x=sign(1.,p(3,lf))*
3048      *       alog((sqrt(p(3,lf)**2+pt**2)+abs(p(3,lf)))/pt)
3049         else
3050           x=sign(1000.,p(3,lf))
3051         endif
3052 c        pp=sqrt(p(2,lf)**2+p(1,lf)**2+p(3,lf)**2)
3053 c        x=0
3054 c        if(pp-p(3,lf).gt.0..and.pp+p(3,lf).gt.0.)x=
3055 c     *       0.5*log((pp+p(3,lf))/(pp-p(3,lf)))
3056          x=abs(x)
3057       elseif(inom.eq.93)then                                     !'phiexp'
3058         x=1000
3059         pt=sqrt(p(1,lf)**2+p(2,lf)**2)
3060         if(pt.gt.0.)x=sign(1.,p(2,lf))*acos(p(1,lf)/pt)
3061         if(x.lt.-pi)x=x+2*pi
3062         if(x.gt.pi)x=x-2*pi
3063       elseif(inom.eq.94)then    !'mu24' ... charged ptl mult |eta|<2.4 (CMS)
3064         x=multc24
3065 c--------------------------------- 101 - 200 ----------------------------
3066 
3067       elseif(inom.eq.101)then           !mulevt
3068         x=1.
3069       elseif(inom.eq.102)then                      !'etevt'
3070         x=0
3071         if(istptl(j).eq.0)then
3072          eef=p(4,lf)
3073          if(maproj.gt.1.or.matarg.gt.1)then
3074            if(idptl(j).ge.1000)eef=eef-prom
3075            if(idptl(j).le.-1000)eef=eef+prom
3076          endif
3077          pp=sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
3078          if(pp.ne.0.)x=eef*sqrt(p(1,lf)**2+p(2,lf)**2)/pp
3079          if(x.ne.x)then
3080            write(ifch,*)x,eef,p(1,lf),p(2,lf),p(3,lf),pp,prom,idptl(j),j
3081            call alist('xan&',1,nptl)
3082            stop 'probleme dans xan'
3083          endif
3084          endif
3085       elseif(inom.eq.103)then
3086         x=p(4,lf)/1000.
3087       elseif(inom.eq.104)then                       !'ev6evt'
3088         x=0
3089         if(istptl(j).eq.0)then
3090          pt=sqrt(p(2,lf)**2+p(1,lf)**2)
3091          eta=0
3092          if(p(3,lf).ne.0..and.pt.ne.0.)eta=sign(1.,p(3,lf))*
3093      *   alog((sqrt(p(3,lf)**2+pt**2)+abs(p(3,lf)))/pt)
3094          if(pt.eq.0.)eta=sign(1e5,p(3,lf))
3095          if(eta.gt.6.0)then
3096           eef=p(4,lf)
3097           if(idptl(j).ge.1000)eef=eef-prom
3098           if(idptl(j).le.-1000)eef=eef+prom
3099           pp=sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
3100           if(pp.ne.0.)x=0.001*eef
3101          endif
3102         endif
3103       elseif(inom.eq.105)then
3104         !etot=maproj*emax(lf)+matarg*0.94  !nur richtig fur target frame!!!!!
3105         etot=maproj*emax(lf)+matarg*emax(lf)
3106         x=p(4,lf)/etot
3107       elseif(inom.eq.106)then
3108         x=isign(1,idptl(j))
3109       elseif(inom.eq.107)then                       !'ptevt'
3110         x=sqrt(p(2,lf)**2+p(1,lf)**2)
3111       elseif(inom.eq.108)then                       !'pmxevt'
3112         x=pmxevt
3113 
3114 c--------------------------------- > 200 ----------------------------
3115 
3116       elseif(inom.eq.201)then
3117         x=1.
3118       elseif(inom.eq.202)then
3119         x=egyevt
3120       elseif(inom.eq.203)then
3121         x=bimevt
3122       elseif(inom.eq.204)then                       !'xbjevt'
3123         x=xbjevt
3124       elseif(inom.eq.205)then                       !'qsqevt'
3125         x=qsqevt
3126       elseif(inom.eq.206)then                       !'yevt'
3127         x=qsqevt/xbjevt/ecms**2
3128       elseif(inom.eq.207)then                       !'eloevt'
3129 c        x=qsqevt/4./elepti+elepti*(1.-qsqevt/xbjevt/ecms**2)
3130         x=elepto
3131       elseif(inom.eq.208)then                       !nd1evt
3132         x=nsdiff(1,noweak(n))
3133       elseif(inom.eq.209)then                       !'nd2evt'
3134         x=nsdiff(2,noweak(n))
3135       elseif(inom.eq.210)then                       !'theevt'
3136 c        eloevt=qsqevt/4./elepti+elepti*(1.-qsqevt/xbjevt/ecms**2)
3137 c        x=acos(1-qsqevt/2./elepti/eloevt)/pi*180.
3138         x=acos(1-qsqevt/2./elepti/elepto)/pi*180.
3139       elseif(inom.eq.211)then                       !'nspevt'
3140         x=0
3141         do i=1,nptl
3142          if((istptl(i).eq.30.or.istptl(i).eq.31)
3143      &      .and.int(idptl(i)/1000000).eq.1)x=x+1
3144         enddo
3145       elseif(inom.eq.212)then                       !'nhpevt'
3146         x=0
3147         do i=1,nptl
3148          if((istptl(i).eq.30.or.istptl(i).eq.31)
3149      &      .and.int(idptl(i)/1000000).eq.3)x=x+1
3150         enddo
3151       elseif(inom.eq.213)then                       !'sigtot'
3152         x=sigtot
3153       elseif(inom.eq.214)then                       !'sigela'
3154         x=sigela
3155       elseif(inom.eq.215)then                       !'sloela'
3156         x=sloela
3157       elseif(inom.eq.216)then                       !'nrgevt'
3158         x=0
3159         do i=1,nptl
3160           if(istptl(i).eq.31.and.int(idptl(i)/10000).eq.2)x=x+1
3161         enddo
3162       elseif(inom.eq.217)then                       !qevt
3163         x=sqrt(qsqevt)
3164       elseif(inom.eq.218)then   !qevt
3165         if(iappl.eq.8)then
3166           x=qtl
3167         else
3168           x=pprt(1,5)
3169         endif
3170       elseif(inom.eq.219)then                       !'nd0evt'  UA1
3171         x=nsdiff(0,noweak(n))
3172       elseif(inom.eq.220)then!------------------------------------------
3173         x=sngl(bofra(1,lf))     !thrust
3174       elseif(inom.eq.221)then
3175         x=1.-sngl(bofra(1,lf))  !1-thrust
3176       elseif(inom.eq.222)then
3177         x=sngl(bofra(2,lf))     !major
3178       elseif(inom.eq.223)then
3179         x=sngl(bofra(3,lf))     !minor
3180       elseif(inom.eq.224)then
3181         x=sngl(bofra(2,lf)-bofra(3,lf)) !oblateness
3182       elseif(inom.eq.230)then!------------------------------------------
3183         x=1.5*(1.-sngl(bofra(1,lf))) !spherecity
3184       elseif(inom.eq.231)then
3185         x=1.5*sngl(bofra(3,lf)) !aplanarity
3186       elseif(inom.eq.232)then
3187         x=3.*sngl(bofra(1,lf)*bofra(2,lf)+bofra(1,lf)*bofra(3,lf)
3188      &       +bofra(2,lf)*bofra(3,lf)) !c-parameter
3189       elseif(inom.eq.233)then
3190         x=27.*sngl(bofra(1,lf)*bofra(2,lf)*bofra(3,lf))   !d-parameter
3191       elseif(inom.eq.234)then                       !npoevt
3192         x=0
3193         do i=1,nptl
3194          if(istptl(i).eq.30.or.istptl(i).eq.31)x=x+1
3195         enddo
3196       elseif(inom.eq.235)then                       !npnevt
3197         x=npjevt+ntgevt   !npnevt
3198       elseif(inom.eq.236)then                       !ikoevt
3199         x=ikoevt
3200       elseif(inom.eq.237)then                       !iktevt
3201 c        x=zkotest
3202       elseif(inom.eq.238)then  !npxevt ... nr of pomerons, including absorbed
3203         x=0
3204         do i=1,nptl
3205          if(istptl(i).eq.30.or.istptl(i).eq.31)x=x+1
3206          if(mod(abs(idptl(i)),100).eq.94)x=x+0.5
3207         enddo
3208       elseif(inom.eq.239)then                       !'nd6evt'
3209         x=nsdiff(6,noweak(n))
3210       elseif(inom.eq.240)then    !mu1evt ... charged ptl multipl for central rap
3211         x=multy1
3212       elseif(inom.eq.241)then    !muievt ... charged ptl multipl
3213         x=multyi
3214       elseif(inom.eq.242)then    !hgtevt ... charged ptl multipl for central eta
3215         x=multc05
3216       elseif(inom.eq.243)then                       !difevt
3217         npom=0
3218         do i=1,nptl
3219          if(istptl(i).eq.30.or.istptl(i).eq.31)npom=npom+1
3220         enddo
3221         x=0
3222         if(npom.eq.0)x=1
3223       elseif(inom.eq.244)then                       !dixevt
3224         zpom=0
3225         do i=1,nptl
3226          if(istptl(i).eq.30.or.istptl(i).eq.31)zpom=zpom+1
3227          if(mod(abs(idptl(i)),100).eq.94)zpom=zpom+0.5
3228         enddo
3229         x=0
3230         if(abs(zpom).lt.0.001)x=1
3231       elseif(inom.eq.245)then                       !'nd7evt' CMS NSD
3232         x=nsdiff(7,noweak(n))
3233       elseif(inom.eq.246)then                       !'nd8evt'  ATLAS
3234         x=nsdiff(8,noweak(n))
3235       elseif(inom.eq.247)then                       !'nd9evt'  ALICE 900 GeV
3236         x=nsdiff(9,noweak(n))
3237       elseif(inom.eq.248)then                       !'ndaevt'  ALICE 2.36 TeV
3238         x=nsdiff(10,noweak(n))
3239       elseif(inom.eq.249)then                       !'ndbevt'  ALICE Inel > 0
3240         x=nsdiff(11,noweak(n))
3241       elseif(inom.eq.250)then
3242         if(iappl.eq.8)then      !mass in
3243           x=-pptl(5,6)
3244         else
3245           x=pprt(5,3)
3246         endif
3247       elseif(inom.eq.251)then
3248         if(iappl.eq.8)then      !mass out
3249           x=pptl(5,7)
3250         else
3251           x=pprt(5,2)
3252         endif
3253       elseif(inom.eq.252)then
3254         if(iappl.eq.8)then
3255           x=-pptl(4,6)
3256         else
3257           x=pprt(4,2)
3258         endif
3259       elseif(inom.eq.253)then
3260         if(iappl.eq.8)then
3261           x=pptl(4,7)
3262         else
3263           x=pprt(4,3)
3264         endif
3265       elseif(inom.eq.254)then
3266         if(iappl.eq.8)then
3267           x=abs(pptl(3,6))
3268         else
3269           x=abs(pprt(3,2))
3270         endif
3271       elseif(inom.eq.255)then
3272         if(iappl.eq.8)then
3273           x=abs(pptl(3,7))
3274           do l=1,5
3275             p(l,lf)=pptl(l,7)
3276           enddo
3277           if(imofra(1,lf).ne.0)then
3278             call utrota(imofra(1,lf),r1fra(1,lf),r1fra(2,lf),r1fra(3,lf)
3279      $           ,p(1,lf),p(2,lf),p(3,lf))
3280           endif
3281           if(imofra(2,lf).ne.0)then !the x-z exchanged is ok !!
3282             call utrota(imofra(2,lf),r2fra(3,lf),r2fra(2,lf),r2fra(1,lf)
3283      $           ,p(3,lf),p(2,lf),p(1,lf))
3284           endif
3285           if(imofra(3,lf).ne.0)then
3286             call utlob4(imofra(3,lf),bofra(1,lf),bofra(2,lf)
3287      $           ,bofra(3,lf) ,bofra(4,lf),bofra(5,lf)
3288      $           ,p(1,lf),p(2,lf),p(3,lf),p(4,lf))
3289           endif
3290           x=abs(p(3,lf))
3291         else
3292           x=abs(pprt(3,3))
3293         endif
3294       elseif(inom.eq.256)then  !pxfevt: leading proton xf in cms
3295         x=-2
3296 c       pmax=sqrt((ecms/2.)**2-prom**2)
3297           pmax=pnullx               !???????????????????
3298         do i=1,nptl
3299           if(idptl(i).eq.1120.and.istptl(i).eq.0)then
3300             if(iframe.eq.11)then
3301               pz=pptl(3,i)
3302             else
3303               amt=sqrt(prom**2+pptl(1,i)**2+pptl(2,i)**2)
3304               rap=alog((pptl(4,i)+pptl(3,i))/amt)
3305      &           -alog((sqrt(pnll**2+ecms**2)+pnll)/ecms)
3306               pz=amt*sinh(rap)
3307             endif
3308             x=max(x,pz/pmax)
3309           endif
3310         enddo
3311       elseif(inom.eq.257)then  !  pi+xf: pi+ yield at cms xf>0.01
3312         x=0.
3313 c        pmax=sqrt((ecms/2)**2-prom*2)
3314           pmax=pnullx               !???????????????????
3315         do i=1,nptl
3316           if(idptl(i).eq.120.and.istptl(i).eq.0)then
3317             if(iframe.eq.11)then
3318               pz=pptl(3,i)
3319             else
3320               amt=sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2)
3321               rap=alog((pptl(4,i)+pptl(3,i))/amt)
3322      &           -alog((sqrt(pnll**2+ecms**2)+pnll)/ecms)
3323               pz=amt*sinh(rap)
3324             endif
3325             if(pz/pmax.gt.0.01)x=x+1.
3326           endif
3327         enddo
3328       elseif(inom.eq.258)then  !  pi-xf: pi- yield at cms xf>0.01
3329         x=0.
3330 c        pmax=sqrt((ecms/2)**2-prom*2)
3331           pmax=pnullx               !???????????????????
3332         do i=1,nptl
3333           if(idptl(i).eq.-120.and.istptl(i).eq.0)then
3334             if(iframe.eq.11)then
3335               pz=pptl(3,i)
3336             else
3337               amt=sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2)
3338               rap=alog((pptl(4,i)+pptl(3,i))/amt)
3339      &           -alog((sqrt(pnll**2+ecms**2)+pnll)/ecms)
3340               pz=amt*sinh(rap)
3341             endif
3342             if(pz/pmax.gt.0.01)x=x+1.
3343           endif
3344         enddo
3345       elseif(inom.eq.260)then!------------------------------
3346         x=sigcut
3347       elseif(inom.eq.261)then
3348         x=keu
3349       elseif(inom.eq.262)then
3350         x=ked
3351       elseif(inom.eq.263)then
3352         x=kes
3353       elseif(inom.eq.265)then
3354         x=kolevt
3355       elseif(inom.eq.266)then
3356         x=sigsd
3357       elseif(inom.eq.267)then
3358         x=nglevt
3359       elseif(inom.eq.268)then  ! kppevt : collision number per participant
3360         x=kolevt/float(npjevt+ntgevt)
3361       elseif(inom.eq.269)then  ! npievt : pion + multi per event
3362         x=0
3363         do i=1,nptl
3364          if(idptl(i).eq.120)x=x+1
3365         enddo
3366       elseif(inom.eq.270)then  ! np2evt : pion + multi per event
3367         x=0
3368         do i=1,nptl
3369          if(idptl(i).eq.120)x=x+1
3370         enddo
3371         x=x/float(npjevt+ntgevt)
3372       elseif(inom.eq.271)then
3373         x=sigdif
3374       elseif(inom.eq.272)then  !number of inelastic collisions per event
3375         x=koievt
3376       elseif(inom.eq.273)then  ! inelasticity (energy loss of leading particle)
3377         x=0.
3378         do i=maproj+matarg+1,nptl
3379           if(istptl(i).eq.0)then
3380             if((((abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000)
3381      *           .and.idproj.gt.1000).or.(iabs(idptl(i)).gt.100
3382      *           .and.idproj.lt.1000)).and.pptl(4,i)
3383      *           .gt.x.and.pptl(3,i).gt.0.)x=pptl(4,i)
3384           endif
3385         enddo
3386         Eini=pptl(4,1)
3387         if(Eini.gt.0.)x=(Eini-x)/Eini
3388       elseif(inom.eq.274)then  ! elasticity (energy of leading particle)
3389         x=0.
3390         do i=maproj+matarg+1,nptl
3391           if(istptl(i).eq.0)then
3392             if((((abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000)
3393      *           .and.idproj.gt.1000).or.(iabs(idptl(i)).gt.100
3394      *           .and.idproj.lt.1000)).and.pptl(4,i)
3395      *           .gt.x.and.pptl(3,i).gt.0.)x=pptl(4,i)
3396           endif
3397         enddo
3398         Eini=pptl(4,1)
3399         if(Eini.gt.0.)x=x/Eini
3400       elseif(inom.eq.275)then         !'itgevt'
3401         x=0
3402         if(nint(ypara(1,n)).ne.0)x=1
3403       elseif(inom.eq.276)then         !'hrdevt' ......  1 = hard event
3404         x=0
3405         if(nint(ypara(1,n)).ne.0)x=1
3406       elseif(inom.eq.277)then         !'ej1evt' .... et of jet 1
3407         x=0
3408         if(nint(ypara(1,n)).ne.0)
3409      &  x=ypara(2,n)
3410       elseif(inom.eq.278)then         !'pj1evt' .... phi of jet 1
3411         x=1000
3412         if(nint(ypara(1,n)).ne.0)
3413      &  x=ypara(4,n)
3414       elseif(inom.eq.279)then         !'ej2evt' .... et of jet 2
3415         x=0
3416         if(nint(ypara(6,n)).ne.0)
3417      &  x=ypara(7,n)
3418       elseif(inom.eq.280)then         !'pj2evt' .... delta_phi of jet 2 1
3419         x=1000
3420         if(nint(ypara(6,n)).ne.0)then
3421           x=ypara(9,n)-ypara(4,n)
3422            if(x.lt.-2.5*pi)then
3423             x=x+4*pi
3424            elseif(x.lt.-0.5*pi)then
3425             x=x+2*pi
3426           elseif(x.gt.3.5*pi)then
3427             x=x-4*pi
3428           elseif(x.gt.1.5*pi)then
3429             x=x-2*pi
3430           endif
3431         endif
3432       elseif(inom.eq.281)then         !'zppevt'
3433         x=zppevt
3434       elseif(inom.eq.282)then         !'zptevt'
3435         x=zptevt
3436       elseif(inom.eq.283)then
3437         stop '**********not used*********'
3438       elseif(inom.eq.284)then                       !'nd3evt'
3439         x=nsdiff(3,noweak(n))
3440       elseif(inom.eq.285)then                       !'nd4evt'
3441         x=nsdiff(4,noweak(n))
3442       elseif(inom.eq.286)then         !'mubevt'
3443         x=multeb
3444       elseif(inom.eq.287)then                       !'nd5evt'
3445         x=nsdiff(5,noweak(n))
3446       elseif(inom.eq.288)then
3447         x=ekievt
3448       elseif(inom.eq.289)then                       !'diffmevt'
3449         x=isdiff(1)
3450       elseif(inom.eq.290)then                       !'diffxevt'
3451         x=isdiff(2)
3452       elseif(inom.eq.291.or.inom.eq.292)then  ! mass of produced system (inelasticity of leading particle )
3453         x=0.
3454         i=idlead
3455         if(i.gt.0)then
3456           pmax=pnullx
3457           if(mod(ifra(lf),10).eq.2)pmax=pnll
3458           x=abs(pptl(3,i)/pmax)
3459           x=(1.-x)*engy*engy
3460           if(inom.eq.291.and.x.gt.0.)x=sqrt(x)
3461 c         write(ifch,*)'ici',i,idptl(i),x,abs(pptl(3,i)/pmax)
3462         endif
3463       elseif(inom.eq.293)then  ! tdevt : -t of leading particle
3464         x=0.
3465         i=idlead
3466         if(i.gt.0)then
3467           pmax=pnullx
3468           if(mod(ifra(lf),10).eq.2)pmax=pnll
3469 c        xxx=(amproj**2-2.*sqrt(amproj**2+pmax**2)*pptl(4,i)
3470 c     *      +2.*abs(pmax*pptl(4,i))+pptl(5,i)**2)
3471           ppt=sqrt(pptl(1,i)**2+pptl(2,i)**2)
3472           if(abs(pptl(3,i)).gt.0.)then
3473             theta=atan(ppt/pptl(3,i))
3474           else
3475             theta=pi/2.
3476           endif
3477           x=abs(pptl(3,i)/pmax)
3478 c -t definition of UA4 (Phys Lett B136,217)
3479           x=pptl(5,i)**2*(1.-x)**2/x+2*x*pmax*pmax*(1.-cos(theta))
3480 c         write(*,*)'ici',i,idptl(i),theta,x,xxx
3481         endif
3482       elseif(inom.eq.294)then          !'ndpevt' pomeron from diffraction
3483         x=0
3484         do i=1,nptl
3485          if((istptl(i).eq.30.or.istptl(i).eq.31)
3486      &      .and.mod(ityptl(i),10).eq.5)x=x+1
3487         enddo
3488       elseif(inom.eq.295)then          !'rapgap' rapidity gap
3489         x=rapgap
3490       elseif(inom.eq.296)then        !'ng1evt'
3491         x=ng1evt
3492       elseif(inom.eq.297)then        !'r21evt'
3493         x=0
3494         if(ng1evt.ne.0)x=ng2evt/float(ng1evt)
3495       elseif(inom.eq.298)then
3496         x=sval(1,n)
3497       elseif(inom.eq.299)then
3498         x=sval(2,n)
3499       elseif(inom.eq.301)then   !---------------------------------------------
3500         if(j.eq.0)then          !initialize
3501           do l=1,4
3502             aimuni(l,n)=0.
3503           enddo
3504         elseif(j.gt.nptl)then   !final calculation
3505           am2=aimuni(4,n)**2-aimuni(3,n)**2
3506      $         -aimuni(2,n)**2-aimuni(1,n)**2
3507           x=sign(sqrt(abs(am2)),am2)
3508 c          print *, x
3509         else                    !routine work
3510           do l=1,4
3511             aimuni(l,n)=aimuni(l,n)+p(l,lf)
3512           enddo
3513 c          print *, j,(p(l,lf),l=1,5)
3514         endif
3515       elseif(inom.ge.302.and.inom.le.305)then   !-----------------------
3516         if(j.eq.0)then          !initialize
3517           do l=1,4
3518             aimuni(l,n)=0.
3519           enddo
3520         elseif(j.gt.nptl)then   !final calculation
3521           if(inom.eq.302) x=max(aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
3522      $         ,aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n)))
3523           if(inom.eq.303) x=min(aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
3524      $         ,aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n)))
3525           if(inom.eq.304) x=abs(aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
3526      $         -aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n)))
3527           if(inom.eq.305) x=aimuni(1,n)/2/(aimuni(2,n)+aimuni(4,n))
3528      $         +aimuni(3,n)/2/(aimuni(2,n)+aimuni(4,n))
3529         else                    !routine work
3530           l=0
3531           if(p(3,lf).lt.0.)l=2
3532           aimuni(1+l,n)=aimuni(1+l,n)+sqrt(p(1,lf)**2+p(2,lf)**2)
3533           aimuni(2+l,n)=aimuni(2+l,n)
3534      $         +sqrt(p(1,lf)**2+p(2,lf)**2+p(3,lf)**2)
3535 
3536         endif
3537       elseif(inom.eq.306.or.inom.eq.307.or.inom.eq.308)then !---------
3538         if(j.eq.0)then          !initialize
3539           do ll=1,8
3540             aimuni(ll,n)=0.
3541           enddo
3542         elseif(j.gt.nptl)then   !final calculation
3543           am2a=aimuni(4,n)**2-aimuni(3,n)**2
3544      $         -aimuni(2,n)**2-aimuni(1,n)**2
3545           am2b=aimuni(8,n)**2-aimuni(7,n)**2
3546      $         -aimuni(6,n)**2-aimuni(5,n)**2
3547           if(inom.eq.306)x=(max(0.,am2a,am2b))/engy**2
3548           if(inom.eq.307)x=(max(0.,min(am2a,am2b)))/engy**2
3549           if(inom.eq.308)x=(abs(am2a-am2b))/engy**2
3550         else                    !routine work
3551           ll=0
3552           if(p(3,lf).lt.0.)ll=4
3553           do l=1,4
3554             aimuni(l+ll,n)=aimuni(l+ll,n)+p(l,lf)
3555           enddo
3556         endif
3557       elseif (inom.eq.310.or.inom.eq.311) then !---------
3558         if(j.eq.0)then          !initialize
3559           aimuni(1,n)=0
3560           aimuni(2,n)=0
3561           do i=1,nptl
3562 c            charge=0.
3563              if(istptl(i).eq.0) then
3564                if (idptl(i).eq.idcod(1,n)) aimuni(1,n)=aimuni(1,n)+1.
3565                if (idptl(i).eq.idcod(2,n)) aimuni(2,n)=aimuni(2,n)+1.
3566              endif
3567            enddo
3568         elseif(j.gt.nptl)then   !final calculation
3569           if(aimuni(1,n).eq.0.or.aimuni(2,n).eq.0) then
3570             ncevt(n)=ncevt(n)-1
3571           endif
3572           x=xmin(n)-100.
3573           do i=1,nbin(n)
3574             zcbin(i,nac(n),n)=abs(zcbin(i,nac(n),n))
3575          enddo
3576        else                    !routine work
3577           if( istptl(j).eq.0
3578      $         .and. aimuni(1,n).ne.0. .and. aimuni(2,n).ne.0. ) then
3579             id1=idptl(j)
3580             if(id1.eq.idcod(1,n) .or. id1.eq.idcod(2,n)) then
3581               y1=sign(1.,pptl(3,j))*alog((pptl(4,j)+abs(pptl(3,j)))
3582      *             /sqrt(pptl(5,j)**2+pptl(1,j)**2+pptl(2,j)**2))
3583               do i=1,nptl
3584                 if(i.eq.j .or. istptl(i).ne.0) goto 124
3585                 id2=idptl(i)
3586                 if(id2.eq.idcod(1,n) .or. id2.eq.idcod(2,n)) then
3587                   y2=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
3588      *                 /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2))
3589                   dy=(y2-y1)
3590                   if(inom.eq.311) dy=abs(dy)
3591                   ib=1+int((dy-xmin(n))*xinc(n))
3592                   if(dy.ge.xmin(n).and.dy.le.xmax(n)) then
3593                     if( id1.eq.idcod(1,n) ) then
3594                       if( id2.eq.idcod(2,n) ) then
3595                         bin(ib,nac(n),n)=bin(ib,nac(n),n)+.5/aimuni(2,n)
3596                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)+1
3597                       else
3598                         bin(ib,nac(n),n)=bin(ib,nac(n),n)-.5/aimuni(1,n)
3599                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)-1
3600                       endif
3601                     else        !id1 is idcod(2,n)
3602                       if(id2.eq.idcod(1,n)) then
3603                         bin(ib,nac(n),n)=bin(ib,nac(n),n)+.5/aimuni(1,n)
3604                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)+1
3605                       else
3606                         bin(ib,nac(n),n)=bin(ib,nac(n),n)-.5/aimuni(2,n)
3607                         zcbin(ib,nac(n),n)=zcbin(ib,nac(n),n)-1
3608                       endif
3609                     endif
3610                   endif
3611                 endif
3612  124            continue 
3613               enddo
3614             endif
3615           endif
3616         endif
3617       elseif (inom.eq.312) then !---------
3618         x=sigine
3619       elseif (inom.eq.313) then !---------
3620         x=sigineaa
3621       elseif (inom.eq.314) then !---------
3622         x=alpD(idxD0,iclpro,icltar)
3623       elseif (inom.eq.315) then !---------
3624         x=alpD(1,iclpro,icltar)
3625       elseif (inom.eq.316) then !---------
3626         x=betD(idxD0,iclpro,icltar)
3627         if(x.lt.0.)x=-10.*x
3628       elseif (inom.eq.317) then !---------
3629         x=betD(1,iclpro,icltar)
3630       elseif (inom.eq.318) then !---------
3631         x=rexdif(iclpro)
3632       elseif (inom.eq.319) then !---------
3633         x=rexdif(icltar)
3634       elseif(inom.eq.320)then    !m14evt ... multipl |eta|<1, pt>0.4
3635         x=multc14
3636       elseif(inom.eq.321)then    !ht3evt ... height |eta|<3.15
3637         x=multc3/6.3
3638       elseif (inom.eq.322) then !---------
3639         x=sigineex
3640       elseif (inom.eq.323) then !---------
3641         x=sigdifex
3642       elseif (inom.eq.324) then !---------
3643         x=sigsdex
3644       elseif (inom.eq.325) then !---------
3645         x=ekin
3646       elseif (inom.eq.326) then !---------
3647         x=sigcutaa
3648       elseif (inom.eq.327) then !---------
3649         x=sigtotaa
3650       elseif (inom.eq.328) then !---------
3651         x=xkappafit(iclegy,iclpro,icltar,1)
3652       elseif (inom.eq.329) then !---------
3653         x=gamD(idxD0,iclpro,icltar)
3654       elseif (inom.eq.330) then !---------
3655         x=gamD(1,iclpro,icltar)
3656       elseif (inom.eq.331) then !---------
3657         x=delD(idxD0,iclpro,icltar)
3658       elseif (inom.eq.332) then !---------
3659         x=delD(1,iclpro,icltar)
3660       elseif(inom.eq.333)then                       !'nd6evt'
3661         x=nsdiff(6,noweak(n))
3662       elseif(inom.eq.334)then !'muxevt' ... multipl 
3663         x= ypara(1,n)         
3664       elseif(inom.eq.335)then                       
3665         x=abs(nint(typevt))                        !ND(1), DD(2), or SD(3)
3666       elseif(inom.eq.339)then    !m25evt ... multipl |eta|<2.5, pt>0.5
3667         x=multc25
3668       elseif(inom.eq.340)then    !'segevt' ... segment multiplicity
3669 c        x=segevt
3670       elseif(inom.eq.341)then    !'ielevt'
3671         x=nsdiff(11,noweak(n))
3672       elseif(inom.eq.342)then                       !'mc1evt' charged particle
3673         x=multc1                                    ! mult for |eta|<1
3674       elseif(inom.eq.343)then   !CDF SD trigger 'sdcdf'
3675         x=isdiff(3)
3676       elseif(inom.eq.344)then   !CDF DPE trigger 'dpecdf'
3677         x=isdiff(4)
3678       elseif(inom.eq.345)then   !CDF DD trigger 'ddcdf'
3679         x=isdiff(5)
3680       elseif(inom.eq.346)then   !'phievt'
3681         x=phievt
3682       elseif(inom.eq.347)then   !'ndcevt'  CMS hadron level NSD (2011)
3683         x=nsdiff(12,noweak(n))
3684       elseif(inom.eq.348)then   !'jetevt' ......  1 = jet event
3685         x=0
3686         if(nint(ypara(1,n)).ne.0)x=1
3687       elseif(inom.eq.349)then   !'epszero (Z for pp)'
3688         x=epszero
3689       elseif(inom.eq.350)then  !xsievt:  xsi = (M_X^2/s) 
3690 c (where M_X = sqrt{ (sum E)^2 - (sum vec-p)^2 } with sum E and sum vec-p 
3691 c being resp. the sum of the energies and the sum of the 3-momenta of the 
3692 c generated particles in the event, excluding the proton with the largest 
3693 c laboratory momentum)
3694         x=xsi
3695 cc xsievt:  xsi = 1-xF_leading
3696 c        x=-2
3697 cc       pmax=sqrt((ecms/2.)**2-prom**2)
3698 c          pmax=pnullx               !???????????????????
3699 c        do i=1,nptl
3700 c          if(abs(idptl(i)).gt.100.and.istptl(i).eq.0)then
3701 c            if(iframe.eq.11)then
3702 c              pz=pptl(3,i)
3703 c            else
3704 c              amt=sqrt(prom**2+pptl(1,i)**2+pptl(2,i)**2)
3705 c              rap=alog((pptl(4,i)+pptl(3,i))/amt)
3706 c     &           -alog((sqrt(pnll**2+ecms**2)+pnll)/ecms)
3707 c              pz=amt*sinh(rap)
3708 c            endif
3709 c            x=max(x,abs(pz/pmax))
3710 c          endif
3711 c        enddo
3712 c        x=max(0.,1.-x)
3713       elseif(inom.eq.351)then  !xsicms: CMS determination of xsi=M2_X/s
3714                                !using xsi=Sum [(E+pz)_i/sqrt(s)]      where i runs on every reconstructed particles (= |eta|<4.9 charged + neutral)
3715         x=-2
3716 c       pmax=sqrt((ecms/2.)**2-prom**2)
3717         pmax=pnullx               !???????????????????
3718         Ef=0.
3719         Eb=0.
3720         Pf=0.
3721         Pb=0.
3722         Esum=0.
3723         Psum=0.
3724         do i=1,nptl
3725           if(istptl(i).eq.0)then
3726             pt=sqrt(pptl(2,i)**2+pptl(1,i)**2)
3727             pz=pptl(3,i)
3728             eta=0.
3729             if(abs(pz).gt.0..and.pt.gt.0.)eta=sign(1.,pz)*
3730      *   log((sqrt(pz**2+pt**2)+abs(pz))/pt)
3731             if(pt.eq.0.)eta=sign(1e5,pz)
3732             if(eta.ge.4.9)then
3733               Ef=Ef+pptl(4,i)
3734               Pf=Pf+pptl(3,i)
3735             elseif(eta.le.-4.9)then
3736               Eb=Eb+pptl(4,i)
3737               Pb=Pb+pptl(3,i)
3738             else
3739               Esum=Esum+pptl(4,i)
3740               Psum=Psum+pptl(3,i)
3741             endif
3742 c        write(ifch,*)'ici',i,eta,Ef,Eb,Esum,Psum,Ef-Pf,Eb+Pb
3743           endif
3744         enddo
3745         if(Ef.ge.Eb)then
3746           x=Esum+Psum+Eb+Pb
3747         else
3748           x=Esum-Psum+Ef-Pf
3749         endif
3750         x=max(0.,x/ecms)
3751 c        write(ifmt,*)'ici',x,min(Esum-Psum,Esum+Psum)/ecms
3752       elseif(inom.eq.352)then   !'calevt' ......  energy in eta range
3753         x= ypara(1,n)      
3754       elseif(inom.eq.353)then   !'fgpevt' ......  max forward rapidity gap in eta range
3755         x= ypara(2,n)      
3756       elseif(inom.eq.354)then   !'bgpevt' ......  max backward rapidity gap in eta range
3757         x= ypara(3,n)      
3758       elseif(inom.eq.355)then   !'gapevt' ......  max backward rapidity gap in eta range
3759         x= ypara(4,n)      
3760       elseif(inom.eq.356)then
3761         x=sigdd
3762       elseif(inom.eq.357)then                       
3763         x=nint(typevt)            !ND(1), DD(2), or SD(3) with fusion < 0
3764       elseif(inom.eq.358)then   !'ndsevt'  CMS hadron level double sided trigger (2012)
3765         x=nsdiff(13,noweak(n))
3766       elseif(inom.eq.359)then    !m24evt ... multipl |eta|<2.4 (CMS)
3767         x=multc24
3768       elseif(inom.eq.360)then   !'ndhevt'  CMS hadron level single sided trigger (HF 2012)
3769         x=nsdiff(14,noweak(n))
3770       elseif (inom.eq.501) then !---------
3771         x=sval(1,1)
3772       elseif (inom.eq.502) then !---------
3773         x=sval(1,2)
3774       elseif (inom.eq.503) then !---------
3775         x=sval(1,3)
3776       elseif (inom.eq.504) then !---------
3777         x=sval(1,4)
3778       elseif(inom.eq.505)then    !'eglevt'  eccentricity
3779         x=eglevt
3780       elseif(inom.eq.506)then    !'fglevt'  eccentricity_part
3781         x=fglevt
3782       elseif(inom.eq.507)then    !'rglevt'  ratio ng2 / ng1
3783         x=0
3784         if(ng1evt.ne.0)
3785      .  x=ng2evt/float(ng1evt)
3786       elseif(inom.eq.508)then    !'sglevt'  area S
3787         x=sglevt
3788       elseif(inom.eq.509)then                       !'ptrevt'
3789         x=0
3790         do i=maproj+matarg+1,minfra
3791         if(istptl(i).eq.25)then
3792           pt=sqrt(pptl(1,i)**2+pptl(2,i)**2)
3793           x=max(x,pt)
3794         endif
3795         enddo
3796       elseif(inom.eq.510)then                       !'rr2evt'
3797         x=cos(2*(phi2neg-phi2pos))
3798         !write(ifmt,*)'+++++ rr2evt +++++ ',x,phi2neg,phi2pos
3799       elseif(inom.eq.511)then  !'perevt'
3800         call getJKNcentr
3801         x=(ncentr-0.5)*5
3802       elseif(inom.eq.512)then  !'paievt'
3803         x=paievt  
3804       elseif(inom.eq.513)then  !'co2evt'
3805         x=co2evt      
3806       elseif(inom.eq.514)then  !'co3evt'
3807         x=co3evt
3808       endif                     !---------------------------------------
3809       end
3810 
3811 c----------------------------------------------------------------------
3812       subroutine mux(n)
3813 c----------------------------------------------------------------------
3814       ! input
3815       !   n = histogram number
3816       !   xpara(1,n) ... etamin
3817       !   xpara(2,n) ... etamax
3818       !   xpara(3,n) ... ptmin
3819       !   xpara(4,n) ... ptmax
3820       !   xpara(5,n) ... factor
3821       !   xpara(6,n) ... divisor
3822       !   xpara(7,n) ... absolute value of eta (1)
3823       !
3824       ! output
3825       !   ypara(1,n) ... multiplicity
3826       !--------------------------------------------------------------
3827       include 'epos.inc'
3828       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
3829       parameter (mypara=10,mxpara=10)
3830       logical ilog,icnx,itrevt,idmod
3831       double precision bin,bbin,zcbin,zbbin
3832       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
3833      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
3834      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
3835      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
3836      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
3837      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
3838      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
3839      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
3840      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
3841      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
3842      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
3843       mul=0
3844       do i=maproj+matarg+1,nptl
3845         if(istptl(i).eq.0)then
3846           pt=pptl(1,i)**2+pptl(2,i)**2
3847           pp=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
3848           if(pt.gt.0.)then
3849             pt=sqrt(pt)
3850             eta=sign(1.,pptl(3,i))*alog((pp+abs(pptl(3,i)))/pt)
3851           else
3852             eta=1000.
3853           endif
3854           if(xpara(7,n).ge.1.)eta=abs(eta)
3855           if(abs(idptl(i)).ge.100
3856      $   .and.abs(idptl(i)).lt.10000)then
3857             call idchrg(idptl(i),ch)
3858             if(abs(ch).gt.0.1)then
3859               if(     eta.ge.xpara(1,n)
3860      *           .and.eta.le.xpara(2,n)
3861      *           .and.pt .gt.xpara(3,n)
3862      *           .and.pt .lt.xpara(4,n)   )mul=mul+1
3863             endif
3864           endif
3865         endif
3866       enddo
3867       ypara(1,n)=mul*xpara(5,n)/xpara(6,n) 
3868       !print*,'+++++++',n,mul
3869       !. ,xpara(1,n),xpara(2,n),xpara(3,n)
3870       !. ,xpara(4,n),xpara(5,n),xpara(6,n)
3871       end
3872 
3873 c----------------------------------------------------------------------
3874       subroutine PairVariables
3875 c----------------------------------------------------------------------
3876       include 'epos.inc'
3877       parameter (maxpv=20000)
3878       real phixx(maxpv),etaxx(maxpv),ptxx(maxpv)
3879       integer idxx(maxpv)
3880       common/cpairs/paievt,co2evt,co3evt
3881       common/cpairs2/paipi(40),co2pi(40),co3pi(40)
3882      .              ,paipr(40),co2pr(40),co3pr(40)
3883      .            ,ipaipi(40),ico2pi(40),ico3pi(40)
3884      .            ,ipaipr(40),ico2pr(40),ico3pr(40)
3885      .            ,maxpt,delpt
3886       maxpt=40
3887       delpt=8./maxpt
3888       do j=1,maxpt
3889         ipaipi(j)=0
3890         ico2pi(j)=0
3891         ico3pi(j)=0
3892         ipaipr(j)=0
3893         ico2pr(j)=0
3894         ico3pr(j)=0
3895         paipi(j)=0.
3896         co2pi(j)=0.
3897         co3pi(j)=0.
3898         paipr(j)=0.
3899         co2pr(j)=0.
3900         co3pr(j)=0.
3901       enddo
3902       paievt=0
3903       co2evt=0
3904       co3evt=0
3905       ii=0
3906       do i=maproj+matarg+1,nptl
3907        if(istptl(i).eq.0)then
3908           px=pptl(1,i)
3909           py=pptl(2,i)
3910           pt=pptl(1,i)**2+pptl(2,i)**2
3911           if(pt.gt.0.)then
3912             pt=sqrt(pt)
3913             theta=atan(pt/pptl(3,i))
3914             if(theta.lt.0.)theta=theta+pi
3915             eta=-log(tan(theta*0.5))
3916           else
3917             eta=1000.
3918           endif
3919           if(abs(eta).lt.0.8)then
3920             call idchrg(idptl(i),ch)
3921             if(abs(ch).gt.0.1)then
3922               ii=ii+1
3923               if(ii.gt.maxpv)then
3924                 write(ifmt,*)
3925      .          '***** ERROR: PairVariables: arrays too small'
3926                  stop'\n\n PairVariables: arrays too small \n\n'
3927               endif
3928               phixx(ii)=polar(px,py)
3929               etaxx(ii)=eta
3930               ptxx(ii)=pt
3931               idxx(ii)=idptl(i)
3932               !print*,'+++++',ii,pt,idxx(ii)
3933             endif    
3934           endif  
3935         endif
3936       enddo      
3937       do m=1,ii
3938       do n=1,ii 
3939         if(    (etaxx(m).lt.-0.5 .and. etaxx(n).gt.0.5)
3940      .    .or. (etaxx(n).lt.-0.5 .and. etaxx(m).gt.0.5) )then
3941           if(abs(etaxx(m)-etaxx(n)).lt.0.999)
3942      .    stop'\n\n ERROR 04112011\n\n' 
3943           !m.ne.n automatic in this case
3944           paievt=paievt+1
3945           co2mn=cos(2*(phixx(m)-phixx(n)))
3946           co3mn=cos(3*(phixx(m)-phixx(n)))
3947           co2evt=co2evt+co2mn
3948           co3evt=co3evt+co3mn
3949           !iam=abs(idxx(m))
3950           !im=ptxx(m)/delpt+1
3951           ian=abs(idxx(n))
3952           in=ptxx(n)/delpt+1
3953           if(in.ge.1.and.in.le.maxpt)then
3954             if(ian.eq.120)then
3955               paipi(in)=paipi(in)+1
3956               co2pi(in)=co2pi(in)+co2mn
3957               co3pi(in)=co3pi(in)+co3mn
3958               !print*,'+++pion++',in, idxx(n),ptxx(n)
3959             elseif(ian.eq.1120)then
3960               paipr(in)=paipr(in)+1
3961               co2pr(in)=co2pr(in)+co2mn
3962               co3pr(in)=co3pr(in)+co3mn
3963               !print*,'+++p++',in, idxx(n),ptxx(n)
3964             endif
3965           endif
3966         endif
3967       enddo
3968       enddo
3969       !print*,'+++++ pairs +++++',ii,paievt,co2evt,co3evt
3970       sum1=0
3971       sum2=0
3972       sum3=0
3973       do j=1,min(20,maxpt)
3974          !write(*,'(a,i5,2(f7.0,2f7.3,3x))') '+++++++',j,
3975          !.   paipi(j),co2pi(j),co3pi(j),paipr(j),co2pr(j),co3pr(j)
3976          sum1=sum1+paipi(j)+paipr(j)
3977          sum2=sum2+co2pi(j)+co2pr(j)
3978          sum3=sum3+co3pi(j)+co3pr(j)
3979       enddo
3980       !print*,'+++++ sum +++++     ',sum1,sum2,sum3
3981       end
3982 
3983 c----------------------------------------------------------------------
3984       subroutine StandardVariables
3985 c----------------------------------------------------------------------
3986       include 'epos.inc'
3987       common/stavar/multc05,multy1,multc14,multyi,multc3,imulty1,multeb
3988      &     ,multc1,multc83,multc24,multc25,rapgap,ipairs1,xsi
3989       parameter(mxxhis=70)
3990       common/varhis/icorrtrig(0:mxxhis),ihardevent(0:mxxhis)
3991      &,ijetfind1(0:mxxhis),ijetfind2(0:mxxhis),imux(0:mxxhis)
3992      &,ifastjet(0:mxxhis),ijetevent(0:mxxhis),icaltrig(0:mxxhis)
3993       logical CDF
3994       common/cphi2/phi2pos,phi2neg
3995       double precision sumE,sumP(4)
3996 
3997       Emax=0.
3998       Pmax=0.
3999       multy1=0
4000       multc05=0
4001       multc14=0
4002       multc24=0
4003       multc25=0
4004       multc1=0
4005       multyi=0
4006       multc3=0
4007       multeb=0
4008       multc83=0
4009       rapgap=0.
4010       etamn=-1000.
4011       etamx=1000.
4012       avcos2p=0
4013       avsin2p=0
4014       avcos2n=0
4015       avsin2n=0
4016       sumE=0d0
4017       sumP(1)=0d0
4018       sumP(2)=0d0
4019       sumP(3)=0d0
4020       sumP(4)=0d0
4021       imax=0
4022 
4023       do i=maproj+matarg+1,nptl
4024         if(istptl(i).eq.0)then
4025           amt=pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2
4026           px=pptl(1,i)
4027           py=pptl(2,i)
4028           pt=pptl(1,i)**2+pptl(2,i)**2
4029           pp=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
4030           et=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(5,i)**2)
4031           if(amt.gt.0..and.pptl(4,i).gt.0.)then
4032             amt=sqrt(amt)
4033             rap=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))/amt)
4034           else
4035             rap=1000.
4036           endif
4037           if(pt.gt.0.)then
4038             pt=sqrt(pt)
4039             theta=atan(pt/pptl(3,i))
4040             if(theta.lt.0.)theta=theta+pi
4041             et=sin(theta)*pptl(4,i)
4042             eta=-log(tan(theta*0.5))
4043           else
4044             eta=1000.
4045           endif
4046           if(eta.gt.0.)then
4047             a=polar(px,py)
4048             avcos2p=avcos2p+cos(2*a)
4049             avsin2p=avsin2p+sin(2*a)
4050           else
4051             a=polar(px,py)
4052             avcos2n=avcos2n+cos(2*a)
4053             avsin2n=avsin2n+sin(2*a)
4054           endif
4055           sumE=sumE+dble(pptl(4,i))
4056           sumP(1)=sumP(1)+dble(pptl(1,i))
4057           sumP(2)=sumP(2)+dble(pptl(2,i))
4058           sumP(3)=sumP(3)+dble(pptl(3,i))
4059           if(idptl(i).eq.idproj.and.pp.gt.Pmax)then
4060             imax=i
4061             Pmax=pp
4062           endif
4063           if(abs(idptl(i)).ge.100
4064      $         .and.abs(idptl(i)).lt.10000)then
4065             call idchrg(idptl(i),ch)
4066             CDF=.false.
4067             if(abs(ch).gt.0.1)then
4068 c---multyi---charged ptl multipl
4069               multyi=multyi+1
4070 c---multy1---charged ptl multipl for central rap
4071               if(abs(rap).le.1.)multy1=multy1+1
4072               if(abs(eta).le.0.5)multc05=multc05+1
4073               if(abs(eta).le.1.and.pt.gt.0.4)multc14=multc14+1
4074               if(abs(eta).le.0.8.and.pt.gt.0.3
4075      *        .and.pt.lt.4.)multc83=multc83+1
4076               if(abs(eta).lt.2.4)multc24=multc24+1
4077               if(abs(eta).le.2.5.and.pt.gt.0.5)multc25=multc25+1
4078               if(abs(eta).le.1)multc1=multc1+1
4079               if(abs(rap).le.3.15)multc3=multc3+1
4080 c---multeb---charged ptl multipl for back rap
4081               if(eta.gt.-3.8.and.eta.lt.-2.8)multeb=multeb+1
4082               if(abs(eta).lt.1.2.and.pt.gt.0.3)then
4083                 CDF=.true.      !CDF CTC acceptance
4084               elseif(abs(eta).gt.3.2.and.abs(eta).lt.5.9)then
4085                 CDF=.true.      !CDF BBC acceptance
4086               endif
4087             endif
4088             if(abs(eta).lt.2.4.and.et.gt.0.2)then
4089               CDF=.true.     !CDF central and plug calorimeters acceptance
4090             elseif(abs(eta).gt.2.2.and.abs(eta).lt.4.2.and.et.gt.1.)then
4091               CDF=.true.     !CDF forward calorimeters acceptance
4092             endif
4093             if(CDF)then
4094               if(eta.le.0)etamn=max(etamn,eta)
4095               if(eta.ge.0)etamx=min(etamx,eta)
4096             endif
4097           endif
4098           if(ilprtg.eq.1)then
4099             if((((abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000)
4100      *         .and.abs(idproj).gt.1000).or.(iabs(idptl(i)).gt.100
4101      *         .and.abs(idproj).lt.1000)).and.pptl(4,i)
4102      *         .gt.Emax.and.pptl(3,i).gt.0.)then
4103               Emax=pptl(4,i)
4104               idlead=i
4105             endif
4106           else
4107             if(abs(idptl(i)).gt.1000.and.abs(idptl(i)).lt.10000
4108      *        .and.pptl(4,i).gt.Emax.and.pptl(3,i).lt.0.)then
4109               Emax=pptl(4,i)
4110               idlead=i
4111             endif
4112           endif
4113         endif
4114       enddo
4115       if(imax.gt.0)then
4116         sumE=sumE-dble(pptl(4,imax))
4117         sumP(1)=sumP(1)-dble(pptl(1,imax))
4118         sumP(2)=sumP(2)-dble(pptl(2,imax))
4119         sumP(3)=sumP(3)-dble(pptl(3,imax))
4120       endif
4121       sumP(4)=sqrt(sumP(1)**2+sumP(2)**2+sumP(3)**2)
4122       xsi=sngl((sumE+sumP(4))*(sumE-sumP(4))/dble(engy)**2)
4123       rapgap=etamx-etamn
4124       if(rapgap.gt.100)rapgap=-1.    !not defined
4125       phi2pos=polar(avcos2p,avsin2p)/2.
4126       phi2neg=polar(avcos2n,avsin2n)/2.
4127       !write(ifmt,*)'+++++ phi2pos phi2neg +++++ '
4128       !. ,phi2pos, phi2neg
4129       end
4130 
4131 c----------------------------------------------------------------------
4132       subroutine jetfind(m,n)
4133 c----------------------------------------------------------------------
4134 c   m = 1 ou 2 (two different definitions)
4135 c   n = histogram
4136 c input(jet definition):
4137 c   xpara(1,n) ... output (m=1) (0=et, 1=pt)
4138 c   xpara(2,n) ... etamin (m=1)
4139 c   xpara(3,n) ... etamax (m=1)
4140 c   xpara(4,n) ... rmax   (m=1) (rmax defining the cone)
4141 c   xpara(5,n) ... ichd   (m=1) (1=charged, 0=all)
4142 c   xpara(6,n) ... output (m=2) (0=et, 1=pt)
4143 c   xpara(7,n) ... etamin (m=2)
4144 c   xpara(8,n) ... etamax (m=2)
4145 c   xpara(9,n) ... rmax   (m=2)
4146 c   xpara(10,n) .. ichd   (m=2)
4147 c output (jet properties):
4148 c   ypara(1,n) ... 1 (found) or 0 if not  (m=1)
4149 c   ypara(2,n) ... et or pt               (m=1)
4150 c   ypara(3,n) ... eta of center          (m=1)
4151 c   ypara(4,n) ... phi of center          (m=1)
4152 c   ypara(5,n)
4153 c   ypara(6,n) ... 1 (found) or 0 if not  (m=2)
4154 c   ypara(7,n) ... et or pt               (m=2)
4155 c   ypara(8,n) ... eta of center          (m=2)
4156 c   ypara(9,n) ... phi of center          (m=2)
4157 c   ypara(10,n)
4158 c----------------------------------------------------------------------
4159 
4160       include 'epos.inc'
4161       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
4162       parameter (mypara=10,mxpara=10)
4163       logical ilog,icnx,itrevt,idmod
4164       parameter (mxval=5)
4165       real ptx(mxval),lst(mxval),etax(mxval),phix(mxval)
4166       double precision bin,bbin,zcbin,zbbin
4167       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
4168      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
4169      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
4170      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
4171      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
4172      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
4173      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
4174      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
4175      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
4176      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
4177      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
4178 
4179       if(m.ne.1.and.m.ne.2)stop'jetfind: value of m not valid.      '
4180 
4181       ipt   = nint(xpara(1+5*(m-1),n))
4182       etamin=           xpara(2+5*(m-1),n)
4183       etamax=           xpara(3+5*(m-1),n)
4184       rmax  =           xpara(4+5*(m-1),n)
4185       ichd  = nint(xpara(5+5*(m-1),n))
4186 
4187       ifound=0
4188       do l=1,mxval
4189         ptx(l)=0
4190         lst(l)=0
4191         etax(l)=0
4192         phix(l)=0
4193       enddo
4194 
4195 ctp060829      pp1=0
4196 ctp060829      pp2=0
4197 ctp060829      pp3=0
4198 
4199       do i=maproj+matarg+1,nptl
4200         iok=0
4201         if(istptl(i).eq.0.and.abs(idptl(i)).lt.10000)iok=1
4202         if(iok.eq.1)call idchrg(idptl(i),ch)
4203         if(ichd.eq.1.and.nint(ch).eq.0)iok=0
4204         if(iok.eq.1)then
4205           p1=pptl(1,i)
4206           p2=pptl(2,i)
4207           p3=pptl(3,i)
4208           pt=sqrt(p1**2+p2**2)
4209                 if(pt.gt.0)then
4210             eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
4211             phi=sign(1.,p2)*acos(p1/pt)
4212           else
4213             eta=10000
4214             phi=0
4215           endif
4216           do k=1,mxval
4217             iok=1
4218             if(m.eq.2)then
4219               dphi=phi-ypara(4,n)
4220               if(dphi.lt.-pi)dphi=dphi+2*pi
4221               if(dphi.gt. pi)dphi=dphi-2*pi
4222               if(abs(dphi).lt.pi/2)iok=0
4223             endif
4224             if(iok.eq.1.and.pt.gt.ptx(k)
4225      &        .and.eta.le.etamax.and.eta.ge.etamin)then
4226               do l=mxval,k+1,-1
4227                ptx(l)=ptx(l-1)
4228                lst(l)=lst(l-1)
4229                etax(l)=etax(l-1)
4230                phix(l)=phix(l-1)
4231               enddo
4232                ptx(k)=pt
4233                lst(k)=i
4234                etax(k)=eta
4235                phix(k)=phi
4236               goto2
4237             endif
4238           enddo
4239   2       continue
4240         endif
4241       enddo
4242 
4243       kk=0
4244       etx=0
4245 
4246       do k=1,mxval
4247        if(lst(k).ne.0)then
4248 
4249         ifound=1
4250         et=0
4251         etaxx=etax(k)
4252         phixx=phix(k)
4253         do j=maproj+matarg+1,nptl
4254           iok=0
4255           if(istptl(j).eq.0.and.abs(idptl(j)).lt.10000)iok=1
4256           if(iok.eq.1)call idchrg(idptl(j),ch)
4257           if(ichd.eq.1.and.nint(ch).eq.0)iok=0
4258           if(iok.eq.1)then
4259             p1=pptl(1,j)
4260             p2=pptl(2,j)
4261             p3=pptl(3,j)
4262             pt=sqrt(p1**2+p2**2)
4263             am=pptl(5,j)
4264                   if(pt.gt.0)then
4265               eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
4266               phi=sign(1.,p2)*acos(p1/pt)
4267             else
4268               eta=-10000
4269               phi=0
4270             endif
4271             if(eta.le.etamax.and.eta.ge.etamin)then
4272               deta=eta-etaxx
4273               dphi=phi-phixx
4274               if(dphi.lt.-pi)dphi=dphi+2*pi
4275               if(dphi.gt. pi)dphi=dphi-2*pi
4276               if(deta**2+dphi**2.lt.rmax**2)then
4277                 if(ipt.eq.0)then           !output is et
4278                   et=et+sqrt(pt**2+am**2)
4279                 else                       !output is pt
4280                   et=et+pt
4281                 endif
4282               endif
4283             endif
4284           endif
4285         enddo
4286         if(et.gt.etx)then
4287           etx=et
4288           kk=k
4289         endif
4290 
4291        endif
4292       enddo
4293 
4294       ypara(1+5*(m-1),n)=ifound
4295       ypara(2+5*(m-1),n)=etx
4296       if(kk.gt.0)then
4297        ypara(3+5*(m-1),n)=etax(kk)
4298        ypara(4+5*(m-1),n)=phix(kk)
4299       endif
4300       return
4301       end
4302 
4303 
4304 c----------------------------------------------------------------------
4305       subroutine fastjet(n)
4306 c----------------------------------------------------------------------
4307 c   n = histogram (to define istptl=100*n of id=9999)
4308 c input(jet definition):
4309 c   xpara(1,n) ... algorithm (1.0=kt, 0.0=Cam/Aachen,  -1.0 = anti-kt)
4310 c   xpara(2,n) ... rapmin (for particles used to define jet)
4311 c   xpara(3,n) ... rapmax (for particles used to define jet)
4312 c   xpara(4,n) ... rmax   (rmax defining the cone)
4313 c   xpara(5,n) ... ichd   (1=charged, 0=all)
4314 c output : new particles (jets four momentum) in particle list triggered
4315 c          by : trigger istptl jet jet
4316 c          can be used as usual particle to plot pt, phi, etc ...
4317 c----------------------------------------------------------------------
4318 
4319       include 'epos.inc'
4320       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
4321       parameter (mypara=10,mxpara=10)
4322       parameter (mxval=10000)
4323       double precision p(4,mxval), rmax, algo,jets(4,mxval)
4324       logical ilog,icnx,itrevt,idmod
4325       double precision bin,bbin,zcbin,zbbin
4326       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
4327      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
4328      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
4329      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
4330      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
4331      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
4332      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
4333      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
4334      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
4335      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
4336      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
4337 
4338       algo  = dble(xpara(1,n))
4339       rapmin=      xpara(2,n)
4340       rapmax=      xpara(3,n)
4341       rmax  = dble(xpara(4,n))
4342       ichd  = nint(xpara(5,n))
4343 
4344       npart=0
4345       do i=maproj+matarg+1,nptl
4346         if(istptl(i).eq.0.and.abs(idptl(i)).lt.10000)then
4347           call idchrg(idptl(i),ch)
4348           if(nint(ch).ne.0.or.ichd.eq.0)then
4349             p1=pptl(1,i)
4350             p2=pptl(2,i)
4351             p3=pptl(3,i)
4352             p4=pptl(4,i)
4353             p5=pptl(5,i)
4354             amt=p5**2+p1**2+p2**2
4355             if(amt.gt.0..and.p4+abs(p3).gt.0.)then 
4356               amt=sqrt(amt)
4357               rap=sign(1.,p3)*log((p4+abs(p3))/amt) !not correct if particles off-shell
4358 c             rap=0.5*log((p4+p3)/(p4-p3))  !always correct but numerically unstable
4359             else
4360               rap=0.         
4361             endif
4362             if(rap.ge.rapmin.and.rap.le.rapmax)then  !particle used for jet
4363               npart=npart+1
4364               if(npart.gt.mxval)then
4365         write(ifmt,*)'Too many particles (mxval) for Fastjet ! skip ...'
4366                 return
4367               endif
4368               p(1,npart)=dble(p1)
4369               p(2,npart)=dble(p2)
4370               p(3,npart)=dble(p3)
4371               p(4,npart)=dble(p4)
4372             endif
4373           endif
4374         endif
4375       enddo
4376 
4377       if(npart.gt.0)then
4378 c.....run the clustering with a pp generalised-kt sequential recombination alg
4379         call fastjetppgenkt(p,npart,rmax,algo,jets,njets) ! ... now you have the jets
4380         if(njets.gt.0)then
4381           nptl0=nptl+1
4382           do i=1,njets
4383             nptl=nptl+1
4384             istptl(nptl)=100*n     !used for trigger
4385             idptl(nptl)=9999       !jet ID
4386             pptl(1,nptl)=sngl(jets(1,i))  !jet px momentum
4387             pptl(2,nptl)=sngl(jets(2,i))  !jet py momentum
4388             pptl(3,nptl)=sngl(jets(3,i))  !jet pz momentum
4389             pptl(4,nptl)=sngl(jets(4,i))  !jet E  momentum
4390             p5x=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2)  !jet Et
4391             p5=0
4392             p5= pptl(4,nptl)**2
4393      .      -pptl(1,nptl)**2-pptl(2,nptl)**2-pptl(3,nptl)**2
4394             if(p5.gt.0)p5=sqrt(p5)
4395             pptl(5,nptl)=p5
4396             !print*,'+++++',p5,p5x
4397 c            write(ifch,*)'fastjet',i,njets,nptl,idptl(nptl),istptl(nptl)
4398 c     &         ,sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2)
4399           enddo
4400           if(ish.ge.5)call alist('list after fastjet&',nptl0,nptl)
4401         endif
4402       endif
4403 
4404 
4405       return
4406       end
4407 
4408 
4409 c----------------------------------------------------------------------
4410       subroutine jetevent(n)
4411 c----------------------------------------------------------------------
4412 c   n = histogram
4413 c input(jet event conditions) using particles (jets) found with fastjet:
4414 c common parameters with fastjet
4415 c   xpara(1,n) ... algorithm (1.0=kt, 0.0=Cam/Aachen,  -1.0 = anti-kt)
4416 c   xpara(2,n) ... rapmin (for particles used to define jet)
4417 c   xpara(3,n) ... rapmax (for particles used to define jet)
4418 c   xpara(4,n) ... rmax   (rmax defining the cone)
4419 c   xpara(5,n) ... ichd   (1=charged, 0=all)
4420 c   xpara(6,n) ... Et_min for event selection
4421 c   xpara(7,n) ... number of jets needed for event selection
4422 c   xpara(8,n) ... delta phi needed for dijet event selection (if xpara(7,n)=2)
4423 c                  if 0, not used, >0 uses delta_phi-pi, <0 uses delta_phi
4424 c output (jet event found or not):
4425 c   ypara(1,n) ... 1 (found) or 0 if not
4426 c----------------------------------------------------------------------
4427 
4428       include 'epos.inc'
4429       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
4430       parameter (mypara=10,mxpara=10)
4431       logical ilog,icnx,itrevt,idmod
4432       double precision bin,bbin,zcbin,zbbin
4433       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
4434      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
4435      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
4436      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
4437      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
4438      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
4439      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
4440      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
4441      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
4442      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
4443      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
4444       dimension inumj(10000)
4445 
4446       ypara(1,n)=0
4447       njet=nint(xpara(7,n))
4448       numj=0
4449 c count the number of jets with Et>Et_min in the event
4450       do i=maproj+matarg+1,nptl
4451        if(istptl(i).eq.100*n.and.pptl(5,i).ge.xpara(6,n))then
4452          numj=numj+1
4453          if(numj.le.10000)then
4454 c save position of jet in particle list
4455            inumj(numj)=i
4456          else
4457            write(ifmt,*)"Too many jets in jetevent, last are skipped!"
4458          endif
4459        endif
4460       enddo
4461 c      write(ifch,*)"jetevent",numj,abs(xpara(8,n)),inumj(1),inumj(2)
4462 c if enough jets, analyse them
4463       if(numj.ge.njet)then
4464 c dijet selection
4465         if(abs(xpara(8,n)).gt.0)then
4466 c test delta phi for 2 higher pt jets (fastjet provide jet list ordered by pt)
4467           j=inumj(1)
4468           p1=pptl(1,j)
4469           p2=pptl(2,j)
4470           pt=pptl(5,j)
4471           if(pt.gt.0)then
4472             phi1=sign(1.,p2)*acos(p1/pt)
4473           else
4474             phi1=0.
4475           endif
4476           j=inumj(2)
4477           p1=pptl(1,j)
4478           p2=pptl(2,j)
4479           pt=pptl(5,j)
4480           if(pt.gt.0)then
4481             phi2=sign(1.,p2)*acos(p1/pt)
4482           else
4483             phi2=0.
4484           endif
4485 c      write(ifch,*)"jetevent phi",pptl(5,inumj(1)),pptl(5,inumj(2))
4486 c     &                       ,phi1,phi2
4487 c test delta phi 
4488           phi0=pi
4489           if(xpara(8,n).lt.0.)phi0=0.
4490           if(abs(abs(phi1-phi2)-phi0).le.abs(xpara(8,n)))then
4491             ypara(1,n)=1        !event selected (dijet)
4492           endif
4493         else                    !event selected (enough jet with Et>Etmin)
4494           ypara(1,n)=1
4495         endif
4496       endif
4497       return
4498       end
4499 
4500 c----------------------------------------------------------------------
4501       subroutine hardevent(n)
4502 c----------------------------------------------------------------------
4503 c   n = histogram
4504 c input(jet event conditions):
4505 c   xpara(2,n) ... pt1
4506 c   xpara(3,n) ... pt2
4507 c   xpara(4,n) ... absetamax
4508 c   xpara(5,n) ... rmax        (r=sqrt(deltaeta**2+deltaphi**2))
4509 c   xpara(6,n) ... Et_min
4510 c output (jet event found or not):
4511 c   ypara(1,n) ... 1 (found) or 0 if not
4512 c----------------------------------------------------------------------
4513 
4514       include 'epos.inc'
4515       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
4516       parameter (mypara=10,mxpara=10)
4517       logical ilog,icnx,itrevt,idmod
4518       double precision bin,bbin,zcbin,zbbin
4519       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
4520      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
4521      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
4522      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
4523      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
4524      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
4525      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
4526      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
4527      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
4528      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
4529      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
4530 
4531       ypara(1,n)=0
4532       do i=maproj+matarg+1,nptl
4533        if(abs(idptl(i)).ge.100.and.abs(idptl(i)).lt.10000.
4534      $  and.istptl(i).eq.0)then
4535         call idchrg(idptl(i),ch)
4536         if(abs(ch).gt.0.1)then
4537           p1=pptl(1,i)
4538           p2=pptl(2,i)
4539           p3=pptl(3,i)
4540           pt=sqrt(p1**2+p2**2)
4541           if(pt.gt.0)then
4542             eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
4543             phi=sign(1.,p2)*acos(p1/pt)
4544           else
4545             eta=10000
4546             phi=0
4547           endif
4548           if(pt.ge.xpara(2,n).and.abs(eta).lt.xpara(4,n))then
4549             et1=pptl(4,i)*pt/sqrt(p3**2+pt**2)
4550             do j=maproj+matarg+1,nptl
4551               if(j.ne.i
4552      $        .and.abs(idptl(j)).ge.100.and.abs(idptl(j)).lt.10000.
4553      $        .and.istptl(j).eq.0)then
4554                 call idchrg(idptl(j),ch)
4555                 if(abs(ch).gt.0.1.and.abs(idptl(j)).ge.100
4556      $          .and.abs(idptl(j)).lt.10000.and.istptl(j).eq.0)then
4557                   p1=pptl(1,j)
4558                   p2=pptl(2,j)
4559                   p3=pptl(3,j)
4560                   pt=sqrt(p1**2+p2**2)
4561                         if(pt.gt.0)then
4562                    etax=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
4563                    phix=sign(1.,p2)*acos(p1/pt)
4564                   else
4565                     etax=-10000
4566                     phix=0
4567                   endif
4568                   if(pt.ge.xpara(3,n).and.abs(etax).lt.xpara(4,n))then
4569                     deta=eta-etax
4570                     dphi=phi-phix
4571                     if(dphi.lt.-pi)dphi=dphi+2*pi
4572                     if(dphi.gt. pi)dphi=dphi-2*pi
4573                     if(deta**2+dphi**2.lt.xpara(5,n)**2)then
4574                     et2=pptl(4,j)*pt/sqrt(p3**2+pt**2)
4575                      if(et1+et2.gt.xpara(6,n))then
4576                       ypara(1,n)=1
4577                       goto1
4578                      endif
4579                     endif
4580                   endif
4581                 endif
4582               endif
4583             enddo
4584           endif
4585         endif
4586        endif
4587       enddo
4588 
4589    1  continue
4590       return
4591       end
4592 
4593 c----------------------------------------------------------------------
4594       subroutine corrtrig(n)
4595 c----------------------------------------------------------------------
4596 c   n = histogram
4597 c input(trigger conditions):
4598 c   xpara(1,n) ... mode (0,1)
4599 c   xpara(2,n) ... ptmin
4600 c   xpara(3,n) ... ptmax
4601 c   xpara(4,n) ... etamin
4602 c   xpara(5,n) ... etamax
4603 c   xpara(6,n) ...
4604 c   xpara(7,n) ...
4605 c output (triggered particle (the one with highest pt if there are several)):
4606 c   ypara(1,n) ... iptl or 0 if no particle found
4607 c   ypara(2,n) ... phi of particle
4608 c   ypara(3,n) ... phi_null
4609 c   ypara(4,n) ... pt lead
4610 c----------------------------------------------------------------------
4611 
4612       include 'epos.inc'
4613       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
4614       parameter (mypara=10,mxpara=10)
4615       logical ilog,icnx,itrevt,idmod
4616       double precision bin,bbin,zcbin,zbbin
4617       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
4618      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
4619      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
4620      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
4621      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
4622      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
4623      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
4624      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
4625      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
4626      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
4627      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
4628 
4629       pt0=xpara(2,n)
4630 
4631       do i=maproj+matarg+1,nptl
4632        if(abs(idptl(i)).ge.100.and.abs(idptl(i)).lt.10000.
4633      $  and.istptl(i).eq.0)then
4634         call idchrg(idptl(i),ch)
4635         if(abs(ch).gt.0.1)then
4636           p1=pptl(1,i)
4637           p2=pptl(2,i)
4638           p3=pptl(3,i)
4639           pt=sqrt(p1**2+p2**2)
4640           pt=max(pt,1e-20)
4641           eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
4642           phi=sign(1.,p2)*acos(p1/pt)
4643           if(pt.ge.pt0.and.pt.le.xpara(3,n).and.eta.gt.xpara(4,n)
4644      $      .and.eta.lt.xpara(5,n))then
4645             pt0=pt
4646             ypara(1,n)=i
4647             ypara(2,n)=phi
4648             ypara(4,n)=pt
4649           endif
4650         endif
4651        endif
4652       enddo
4653       ypara(3,n)=-0.5
4654       if(nint(xpara(1,n)).eq.1)ypara(3,n)=0.0
4655 
4656       return
4657       end
4658 
4659 c----------------------------------------------------------------------
4660       subroutine caltrig(n)
4661 c----------------------------------------------------------------------
4662 c   n = histogram
4663 c input(trigger conditions):
4664 c   xpara(1,n) ... mode (0,1,2) (one eta side, 2 eta side independently or, 
4665 c                                2 eta side simultaneously)
4666 c   xpara(2,n) ... etamin
4667 c   xpara(3,n) ... etamax
4668 c   xpara(4,n) ... 0 all or 1 charged or 2 charged + photons
4669 c   xpara(5,n) ... ptmin
4670 c   xpara(6,n) ... ptmax
4671 c output (triggered energy):
4672 c   ypara(1,n) ... max E for a particle in eta range (per side)
4673 c----------------------------------------------------------------------
4674 
4675       include 'epos.inc'
4676       parameter (mxhis=500,mxcontr=500,mxidcd=60,mxtri=50,mxbin=405)
4677       parameter (mypara=10,mxpara=10)
4678       logical ilog,icnx,itrevt,idmod
4679       double precision bin,bbin,zcbin,zbbin
4680       common/bins/bin(mxbin,2,mxhis),zcbin(mxbin,2,mxhis)
4681      $     ,bbin(mxbin,2,mxcontr),itrevt(mxhis),zbbin(mxbin,2,mxcontr)
4682      $     ,nac(mxhis),ilog(mxhis),icnx(mxhis),xinc(mxhis),ncevt(mxhis)
4683      $     ,sval(2,mxhis),valtri(mxtri,mxhis),ntrc(mxtri,mxhis)
4684      $     ,xmin(mxhis),xmax(mxhis),nhis,noweak(mxhis)
4685      $     ,ivar(2,mxhis),inorm(mxhis),nbin(mxhis),nidcod(mxhis)
4686      $     ,idcod(mxidcd,mxhis),idmod(mxidcd,mxhis),ntri(mxhis)
4687      $     ,itri(mxtri,mxhis),xmitri(mxtri,mxhis),xmatri(mxtri,mxhis)
4688      $     ,xmitrp(mxtri,mxhis),xmatrp(mxtri,mxhis),xpara(mxpara,mxhis)
4689      $     ,ypara(mypara,mxhis),lookcontr(mxhis)
4690      $     ,lookcontrx(mxhis),ncontrall,icontrtyp(mxhis),nccevt(mxcontr)
4691       logical cont
4692       double precision Eforw,Eback
4693 
4694       mode=nint(xpara(1,n))
4695       etamin=xpara(2,n)
4696       etamax=xpara(3,n)
4697       ptmin=xpara(5,n)
4698       ptmax=xpara(6,n)
4699       ichrd=nint(xpara(4,n))
4700       Eforw=0.d0
4701       Eback=0.d0
4702       etaf=1000.
4703       etab=1000.
4704 
4705       do i=maproj+matarg+1,nptl
4706        if(istptl(i).eq.0)then
4707         call idchrg(idptl(i),ch)
4708         if(ichrd.ge.1)then
4709           cont=abs(ch).gt.0.1.or.(idptl(i).eq.10.and.ichrd.eq.2)
4710         else
4711           cont=.true.
4712         endif
4713         if(cont)then
4714           p1=pptl(1,i)
4715           p2=pptl(2,i)
4716           p3=pptl(3,i)
4717           pt=sqrt(p1**2+p2**2)
4718           pt=max(pt,1e-20)
4719           if(pt.gt.ptmin.and.pt.lt.ptmax)then
4720             eta=sign(1.,p3)*alog((sqrt(p3**2+pt**2)+abs(p3))/pt)
4721 c          phi=sign(1.,p2)*acos(p1/pt)
4722             if(eta.gt.etamin.and.eta.lt.etamax)
4723      *           Eforw=max(Eforw,dble(pptl(4,i)))
4724             if(eta.lt.-etamin.and.eta.gt.-etamax)
4725      *           Eback=max(Eback,dble(pptl(4,i)))
4726             if(eta.le.etamax)etaf=min(etaf,etamax-eta)
4727             if(eta.ge.etamin)etab=min(etab,eta-etamin)
4728           endif
4729         endif
4730        endif
4731       enddo
4732       if(mode.eq.0)then
4733         ypara(1,n)=sngl(Eforw)
4734       elseif(mode.eq.1)then
4735         ypara(1,n)=sngl(max(Eforw,Eback))
4736       elseif(mode.eq.2)then
4737         ypara(1,n)=sngl(min(Eforw,Eback))
4738       endif
4739       ypara(2,n)=etaf
4740       ypara(3,n)=etab
4741       ypara(4,n)=max(etab,etaf)
4742 c      if(typevt.eq.3)print *,mode, ypara(1,n),Eforw,Eback,etaf,etab
4743 c      print *,typevt,mode, ypara(1,n),Eforw,Eback,etaf,etab
4744 
4745 
4746       return
4747       end
4748 
4749 c----------------------------------------------------------------------