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