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