File indexing completed on 2024-04-06 12:14:06
0001
0002 subroutine xiniall
0003
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
0060 subroutine xini
0061
0062
0063
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
0120 indfra=1
0121
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
0189 ivfra(1,nhis)=indfra
0190 ivfra(2,nhis)=indfra
0191 else
0192 inpfra=inl
0193
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
0611 subroutine xana
0612
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
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
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
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
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
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
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
1092
1093 do n=1,nhis
1094 if(ivar(1,n).eq.-1.or.ivar(2,n).eq.-1)goto 99
1095
1096
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
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
1141
1142
1143
1144
1145
1146
1147
1148
1149
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
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
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
1218 subroutine triggercondition(i,n,x,go)
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
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
1293 subroutine fillhistoconex(n,x,y,lf,j) !for conex
1294
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
1369 subroutine xhis(n)
1370
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
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
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
1586 integer function nsdiff(insdif,now)
1587
1588
1589
1590
1591
1592
1593
1594
1595 include 'epos.inc'
1596 integer ncevt,nsdi(0:20)
1597 logical cont
1598 data ncevt/1/
1599 save nsdi,ncevt
1600
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
1628 nowdum=now
1629
1630
1631
1632
1633
1634
1635
1636
1637
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
1730 integer function isdiff(isdif)
1731
1732
1733
1734
1735
1736
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
1763
1764
1765
1766
1767 if(isdif.le.2)yyy=-sign(1.,float(ilprtg))*yyy !trigger on antiproton (target : ilprt=-1)
1768
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
1775
1776
1777 endif
1778
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
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
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
1847 subroutine xtrans(cvar,inom,ifr,n)
1848
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
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
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
2501 stop
2502 endif
2503 end
2504
2505
2506 subroutine xval(n,inom,lf,j,x)
2507
2508
2509
2510
2511
2512
2513
2514
2515
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
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
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
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
2699
2700
2701
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
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
2728
2729
2730
2731 elseif(inom.eq.21)then
2732
2733
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
2767
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
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
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
2914 call utlob5(yhaha,p1,p2,p3,p4,p5)
2915
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
2931 call utlob5(-yhaha,p1,p2,p3,p4,p5)
2932
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
3053
3054
3055
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
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
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
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
3137
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
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
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
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
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
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
3470
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
3479 x=pptl(5,i)**2*(1.-x)**2/x+2*x*pmax*pmax*(1.-cos(theta))
3480
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
3509 else !routine work
3510 do l=1,4
3511 aimuni(l,n)=aimuni(l,n)+p(l,lf)
3512 enddo
3513
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
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
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
3691
3692
3693
3694 x=xsi
3695
3696
3697
3698
3699
3700
3701
3702
3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
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
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
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
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
3812 subroutine mux(n)
3813
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
3874 subroutine PairVariables
3875
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
3984 subroutine StandardVariables
3985
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
4069 multyi=multyi+1
4070
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
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
4132 subroutine jetfind(m,n)
4133
4134
4135
4136
4137
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148
4149
4150
4151
4152
4153
4154
4155
4156
4157
4158
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
4196
4197
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
4305 subroutine fastjet(n)
4306
4307
4308
4309
4310
4311
4312
4313
4314
4315
4316
4317
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
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
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
4398
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
4410 subroutine jetevent(n)
4411
4412
4413
4414
4415
4416
4417
4418
4419
4420
4421
4422
4423
4424
4425
4426
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
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
4455 inumj(numj)=i
4456 else
4457 write(ifmt,*)"Too many jets in jetevent, last are skipped!"
4458 endif
4459 endif
4460 enddo
4461
4462
4463 if(numj.ge.njet)then
4464
4465 if(abs(xpara(8,n)).gt.0)then
4466
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
4486
4487
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
4501 subroutine hardevent(n)
4502
4503
4504
4505
4506
4507
4508
4509
4510
4511
4512
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
4594 subroutine corrtrig(n)
4595
4596
4597
4598
4599
4600
4601
4602
4603
4604
4605
4606
4607
4608
4609
4610
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
4660 subroutine caltrig(n)
4661
4662
4663
4664
4665
4666
4667
4668
4669
4670
4671
4672
4673
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
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
4743
4744
4745
4746 return
4747 end
4748
4749
