1 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
2 |
* this file contains all subroutines and functions |
* this file contains all subroutines and functions |
3 |
* that are needed for position finding algorithms |
* that are needed for position finding algorithms: |
4 |
* |
* |
5 |
|
* subroutine idtoc(ipfa,cpfa) |
6 |
|
* |
7 |
|
* subroutine applypfa(PFAtt,ic,ang,corr,res) |
8 |
|
* |
9 |
|
* integer function npfastrips(ic,angle) |
10 |
|
* |
11 |
|
* ----------------------------------------------------------------- |
12 |
|
* p.f.a. |
13 |
|
* ----------------------------------------------------------------- |
14 |
|
* real function pfaeta(ic,angle) |
15 |
|
* real function pfaetal(ic,angle) |
16 |
|
* real function pfaeta2(ic,angle) |
17 |
|
* real function pfaeta3(ic,angle) |
18 |
|
* real function pfaeta4(ic,angle) |
19 |
|
* real function cog(ncog,ic) |
20 |
|
* |
21 |
|
* ----------------------------------------------------------------- |
22 |
|
* risoluzione spaziale media, stimata dalla simulazione (samuele) |
23 |
|
* ----------------------------------------------------------------- |
24 |
|
* FUNCTION risxeta2(angle) |
25 |
|
* FUNCTION risxeta3(angle) |
26 |
|
* FUNCTION risxeta4(angle) |
27 |
|
* FUNCTION risyeta2(angle) |
28 |
|
* FUNCTION risy_cog(angle) |
29 |
|
* FUNCTION risx_cog(angle) |
30 |
|
* real function riseta(iview,angle) |
31 |
|
* ----------------------------------------------------------------- |
32 |
|
* fattore moltiplicativo per tenere conto della dipendenza della |
33 |
|
* risoluzione dal rumore delle strip |
34 |
|
* ----------------------------------------------------------------- |
35 |
|
* real function fbad_cog(ncog,ic) |
36 |
|
* real function fbad_eta(ic,angle) |
37 |
|
* |
38 |
|
* ----------------------------------------------------------------- |
39 |
|
* NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE |
40 |
|
* ----------------------------------------------------------------- |
41 |
|
* real function riscogtheor(ncog,ic) |
42 |
|
* real function risetatheor(ncog,ic,angle) |
43 |
|
* |
44 |
|
* ----------------------------------------------------------------- |
45 |
|
* correzione landi |
46 |
|
* ----------------------------------------------------------------- |
47 |
|
* real function pfacorr(ic,angle) |
48 |
* |
* |
49 |
|
* real function effectiveangle(ang,iview,bbb) |
50 |
|
* real function fieldcorr(iview,bbb) |
51 |
|
* |
52 |
|
* NB - The angle is the "effective angle", which is relative |
53 |
|
* to the sensor and it takes into account the magnetic field |
54 |
|
* |
55 |
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
56 |
|
|
57 |
|
subroutine idtoc(ipfa,cpfa) |
58 |
|
|
59 |
|
integer ipfa |
60 |
|
character*10 cpfa |
61 |
|
|
62 |
|
CPFA='COG4' |
63 |
|
if(ipfa.eq.0)CPFA='ETA' |
64 |
|
if(ipfa.eq.2)CPFA='ETA2' |
65 |
|
if(ipfa.eq.3)CPFA='ETA3' |
66 |
|
if(ipfa.eq.4)CPFA='ETA4' |
67 |
|
if(ipfa.eq.5)CPFA='ETAL' |
68 |
|
if(ipfa.eq.10)CPFA='COG' |
69 |
|
if(ipfa.eq.11)CPFA='COG1' |
70 |
|
if(ipfa.eq.12)CPFA='COG2' |
71 |
|
if(ipfa.eq.13)CPFA='COG3' |
72 |
|
if(ipfa.eq.14)CPFA='COG4' |
73 |
|
|
74 |
|
end |
75 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
76 |
|
real function effectiveangle(ang,iview,bbb) |
77 |
|
|
78 |
|
include 'commontracker.f' |
79 |
|
|
80 |
|
effectiveangle = 0. |
81 |
|
|
82 |
|
if(mod(iview,2).eq.0)then |
83 |
|
c ================================================= |
84 |
|
c X view |
85 |
|
c ================================================= |
86 |
|
c here bbb is the y component of the m.field |
87 |
|
angx = ang |
88 |
|
by = bbb |
89 |
|
if(iview.eq.12) angx = -1. * ang |
90 |
|
if(iview.eq.12) by = -1. * bbb |
91 |
|
cc tgtemp = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!! |
92 |
|
tgtemp = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001 |
93 |
|
|
94 |
|
elseif(mod(iview,2).eq.1)then |
95 |
|
c ================================================= |
96 |
|
c Y view |
97 |
|
c ================================================= |
98 |
|
c here bbb is the x component of the m.filed |
99 |
|
angy = ang |
100 |
|
bx = bbb |
101 |
|
tgtemp = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001 |
102 |
|
|
103 |
|
endif |
104 |
|
effectiveangle = 180.*atan(tgtemp)/acos(-1.) |
105 |
|
|
106 |
|
return |
107 |
|
end |
108 |
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
109 |
|
real function fieldcorr(iview,bbb) |
110 |
|
|
111 |
|
include 'commontracker.f' |
112 |
|
|
113 |
|
fieldcorr = 0. |
114 |
|
|
115 |
|
if(mod(iview,2).eq.0)then |
116 |
|
|
117 |
|
c ================================================= |
118 |
|
c X view |
119 |
|
c ================================================= |
120 |
|
c here bbb is the y component of the m.field |
121 |
|
by = bbb |
122 |
|
if(iview.eq.12) by = -1. * bbb |
123 |
|
fieldcorr = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX |
124 |
|
|
125 |
|
elseif(mod(iview,2).eq.1)then |
126 |
|
c ================================================= |
127 |
|
c Y view |
128 |
|
c ================================================= |
129 |
|
c here bbb is the x component of the m.filed |
130 |
|
bx = bbb |
131 |
|
fieldcorr = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY |
132 |
|
|
133 |
|
endif |
134 |
|
|
135 |
|
return |
136 |
|
end |
137 |
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
138 |
|
|
139 |
|
subroutine applypfa(PFAtt,ic,ang,corr,res) |
140 |
|
*--------------------------------------------------------------- |
141 |
|
* this subroutine calculate the coordinate of cluster ic (in |
142 |
|
* strip units), relative to the strip with the maximum signal, |
143 |
|
* and its spatial resolution (in cm), applying PFAtt. |
144 |
|
* ang is the effective angle, relative to the sensor |
145 |
|
*--------------------------------------------------------------- |
146 |
|
|
147 |
|
character*4 PFAtt |
148 |
|
include 'commontracker.f' |
149 |
|
include 'level1.f' |
150 |
|
|
151 |
|
corr = 0 |
152 |
|
res = 0 |
153 |
|
|
154 |
|
if(ic.le.0)return |
155 |
|
|
156 |
|
iview = VIEW(ic) |
157 |
|
|
158 |
|
if(mod(iview,2).eq.0)then |
159 |
|
c ================================================= |
160 |
|
c X view |
161 |
|
c ================================================= |
162 |
|
|
163 |
|
res = RESXAV |
164 |
|
|
165 |
|
if(PFAtt.eq.'COG1')then |
166 |
|
|
167 |
|
corr = 0 |
168 |
|
res = 1e-4*pitchX/sqrt(12.)!!res |
169 |
|
|
170 |
|
elseif(PFAtt.eq.'COG2')then |
171 |
|
|
172 |
|
corr = cog(2,ic) |
173 |
|
res = risx_cog(abs(ang))!TEMPORANEO |
174 |
|
res = res*fbad_cog(2,ic) |
175 |
|
|
176 |
|
elseif(PFAtt.eq.'COG3')then |
177 |
|
|
178 |
|
corr = cog(3,ic) |
179 |
|
res = risx_cog(abs(ang))!TEMPORANEO |
180 |
|
res = res*fbad_cog(3,ic) |
181 |
|
|
182 |
|
elseif(PFAtt.eq.'COG4')then |
183 |
|
|
184 |
|
corr = cog(4,ic) |
185 |
|
res = risx_cog(abs(ang))!TEMPORANEO |
186 |
|
res = res*fbad_cog(4,ic) |
187 |
|
|
188 |
|
elseif(PFAtt.eq.'ETA2')then |
189 |
|
|
190 |
|
corr = pfaeta2(ic,ang) |
191 |
|
res = risxeta2(abs(ang)) |
192 |
|
res = res*fbad_cog(2,ic) |
193 |
|
|
194 |
|
elseif(PFAtt.eq.'ETA3')then |
195 |
|
|
196 |
|
corr = pfaeta3(ic,ang) |
197 |
|
res = risxeta3(abs(ang)) |
198 |
|
res = res*fbad_cog(3,ic) |
199 |
|
|
200 |
|
elseif(PFAtt.eq.'ETA4')then |
201 |
|
|
202 |
|
corr = pfaeta4(ic,ang) |
203 |
|
res = risxeta4(abs(ang)) |
204 |
|
res = res*fbad_cog(4,ic) |
205 |
|
|
206 |
|
elseif(PFAtt.eq.'ETA')then |
207 |
|
|
208 |
|
corr = pfaeta(ic,ang) |
209 |
|
c res = riseta(ic,ang) |
210 |
|
res = riseta(iview,ang) |
211 |
|
res = res*fbad_eta(ic,ang) |
212 |
|
|
213 |
|
elseif(PFAtt.eq.'ETAL')then |
214 |
|
|
215 |
|
corr = pfaetal(ic,ang) |
216 |
|
res = riseta(iview,ang) |
217 |
|
res = res*fbad_eta(ic,ang) |
218 |
|
|
219 |
|
elseif(PFAtt.eq.'COG')then |
220 |
|
|
221 |
|
corr = cog(0,ic) |
222 |
|
res = risx_cog(abs(ang)) |
223 |
|
res = res*fbad_cog(0,ic) |
224 |
|
|
225 |
|
else |
226 |
|
if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt |
227 |
|
endif |
228 |
|
|
229 |
|
|
230 |
|
* ====================================== |
231 |
|
* temporary patch for saturated clusters |
232 |
|
* ====================================== |
233 |
|
if( nsatstrips(ic).gt.0 )then |
234 |
|
corr = cog(4,ic) |
235 |
|
res = pitchX*1e-4/sqrt(12.) |
236 |
|
cc cc=cog(4,ic) |
237 |
|
c$$$ print*,ic,' *** ',cc |
238 |
|
c$$$ print*,ic,' *** ',res |
239 |
|
endif |
240 |
|
|
241 |
|
|
242 |
|
elseif(mod(iview,2).eq.1)then |
243 |
|
c ================================================= |
244 |
|
c Y view |
245 |
|
c ================================================= |
246 |
|
|
247 |
|
res = RESYAV |
248 |
|
|
249 |
|
if(PFAtt.eq.'COG1')then |
250 |
|
|
251 |
|
corr = 0 |
252 |
|
res = 1e-4*pitchY/sqrt(12.)!res |
253 |
|
|
254 |
|
elseif(PFAtt.eq.'COG2')then |
255 |
|
|
256 |
|
corr = cog(2,ic) |
257 |
|
res = risy_cog(abs(ang))!TEMPORANEO |
258 |
|
res = res*fbad_cog(2,ic) |
259 |
|
|
260 |
|
elseif(PFAtt.eq.'COG3')then |
261 |
|
|
262 |
|
corr = cog(3,ic) |
263 |
|
res = risy_cog(abs(ang))!TEMPORANEO |
264 |
|
res = res*fbad_cog(3,ic) |
265 |
|
|
266 |
|
elseif(PFAtt.eq.'COG4')then |
267 |
|
|
268 |
|
corr = cog(4,ic) |
269 |
|
res = risy_cog(abs(ang))!TEMPORANEO |
270 |
|
res = res*fbad_cog(4,ic) |
271 |
|
|
272 |
|
elseif(PFAtt.eq.'ETA2')then |
273 |
|
|
274 |
|
corr = pfaeta2(ic,ang) |
275 |
|
res = risyeta2(abs(ang)) |
276 |
|
res = res*fbad_cog(2,ic) |
277 |
|
|
278 |
|
elseif(PFAtt.eq.'ETA3')then |
279 |
|
|
280 |
|
corr = pfaeta3(ic,ang) |
281 |
|
res = res*fbad_cog(3,ic) |
282 |
|
|
283 |
|
elseif(PFAtt.eq.'ETA4')then |
284 |
|
|
285 |
|
corr = pfaeta4(ic,ang) |
286 |
|
res = res*fbad_cog(4,ic) |
287 |
|
|
288 |
|
elseif(PFAtt.eq.'ETA')then |
289 |
|
|
290 |
|
corr = pfaeta(ic,ang) |
291 |
|
c res = riseta(ic,ang) |
292 |
|
res = riseta(iview,ang) |
293 |
|
res = res*fbad_eta(ic,ang) |
294 |
|
|
295 |
|
elseif(PFAtt.eq.'ETAL')then |
296 |
|
|
297 |
|
corr = pfaetal(ic,ang) |
298 |
|
res = riseta(iview,ang) |
299 |
|
res = res*fbad_eta(ic,ang) |
300 |
|
|
301 |
|
elseif(PFAtt.eq.'COG')then |
302 |
|
|
303 |
|
corr = cog(0,ic) |
304 |
|
res = risy_cog(abs(ang)) |
305 |
|
res = res*fbad_cog(0,ic) |
306 |
|
|
307 |
|
else |
308 |
|
if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt |
309 |
|
endif |
310 |
|
|
311 |
|
|
312 |
|
* ====================================== |
313 |
|
* temporary patch for saturated clusters |
314 |
|
* ====================================== |
315 |
|
if( nsatstrips(ic).gt.0 )then |
316 |
|
corr = cog(4,ic) |
317 |
|
res = pitchY*1e-4/sqrt(12.) |
318 |
|
cc cc=cog(4,ic) |
319 |
|
c$$$ print*,ic,' *** ',cc |
320 |
|
c$$$ print*,ic,' *** ',res |
321 |
|
endif |
322 |
|
|
323 |
|
endif |
324 |
|
end |
325 |
|
|
326 |
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
327 |
|
integer function npfastrips(ic,angle) |
328 |
|
*-------------------------------------------------------------- |
329 |
|
* thid function returns the number of strips used |
330 |
|
* to evaluate the position of a cluster, according to the p.f.a. |
331 |
|
*-------------------------------------------------------------- |
332 |
|
include 'commontracker.f' |
333 |
|
include 'level1.f' |
334 |
|
include 'calib.f' |
335 |
|
|
336 |
|
character*4 usedPFA |
337 |
|
|
338 |
|
|
339 |
|
|
340 |
|
call idtoc(pfaid,usedPFA) |
341 |
|
|
342 |
|
npfastrips=-1 |
343 |
|
|
344 |
|
if(usedPFA.eq.'COG1')npfastrips=1 |
345 |
|
if(usedPFA.eq.'COG2')npfastrips=2 |
346 |
|
if(usedPFA.eq.'COG3')npfastrips=3 |
347 |
|
if(usedPFA.eq.'COG4')npfastrips=4 |
348 |
|
if(usedPFA.eq.'ETA2')npfastrips=2 |
349 |
|
if(usedPFA.eq.'ETA3')npfastrips=3 |
350 |
|
if(usedPFA.eq.'ETA4')npfastrips=4 |
351 |
|
* ---------------------------------------------------------------- |
352 |
|
if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then |
353 |
|
c print*,VIEW(ic),angle |
354 |
|
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
355 |
|
if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then |
356 |
|
npfastrips=2 |
357 |
|
elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then |
358 |
|
npfastrips=3 |
359 |
|
elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then |
360 |
|
npfastrips=4 |
361 |
|
else |
362 |
|
npfastrips=4 !COG4 |
363 |
|
endif |
364 |
|
else !X-view |
365 |
|
if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then |
366 |
|
npfastrips=2 |
367 |
|
elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then |
368 |
|
npfastrips=3 |
369 |
|
elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then |
370 |
|
npfastrips=4 |
371 |
|
else |
372 |
|
npfastrips=4 !COG4 |
373 |
|
endif |
374 |
|
endif |
375 |
|
endif |
376 |
|
* ---------------------------------------------------------------- |
377 |
|
if(usedPFA.eq.'COG')then |
378 |
|
|
379 |
|
npfastrips=0 |
380 |
|
|
381 |
|
c$$$ iv=VIEW(ic) |
382 |
|
c$$$ if(mod(iv,2).eq.1)incut=incuty |
383 |
|
c$$$ if(mod(iv,2).eq.0)incut=incutx |
384 |
|
c$$$ istart = INDSTART(IC) |
385 |
|
c$$$ istop = TOTCLLENGTH |
386 |
|
c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 |
387 |
|
c$$$ mu = 0 |
388 |
|
c$$$ do i = INDMAX(IC),istart,-1 |
389 |
|
c$$$ ipos = i-INDMAX(ic) |
390 |
|
c$$$ cut = incut*CLSIGMA(i) |
391 |
|
c$$$ if(CLSIGNAL(i).ge.cut)then |
392 |
|
c$$$ mu = mu + 1 |
393 |
|
c$$$ print*,i,mu |
394 |
|
c$$$ else |
395 |
|
c$$$ goto 10 |
396 |
|
c$$$ endif |
397 |
|
c$$$ enddo |
398 |
|
c$$$ 10 continue |
399 |
|
c$$$ do i = INDMAX(IC)+1,istop |
400 |
|
c$$$ ipos = i-INDMAX(ic) |
401 |
|
c$$$ cut = incut*CLSIGMA(i) |
402 |
|
c$$$ if(CLSIGNAL(i).ge.cut)then |
403 |
|
c$$$ mu = mu + 1 |
404 |
|
c$$$ print*,i,mu |
405 |
|
c$$$ else |
406 |
|
c$$$ goto 20 |
407 |
|
c$$$ endif |
408 |
|
c$$$ enddo |
409 |
|
c$$$ 20 continue |
410 |
|
c$$$ npfastrips=mu |
411 |
|
|
412 |
|
endif |
413 |
|
* ---------------------------------------------------------------- |
414 |
|
|
415 |
|
c print*,pfaid,usedPFA,angle,npfastrips |
416 |
|
|
417 |
|
return |
418 |
|
end |
419 |
|
|
420 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
421 |
real function pfaeta(ic,angle) |
real function pfaeta(ic,angle) |
428 |
* according to the angle |
* according to the angle |
429 |
*-------------------------------------------------------------- |
*-------------------------------------------------------------- |
430 |
include 'commontracker.f' |
include 'commontracker.f' |
|
c include 'calib.f' |
|
431 |
include 'level1.f' |
include 'level1.f' |
432 |
|
include 'calib.f' |
433 |
|
|
434 |
pfaeta = 0 |
pfaeta = 0 |
435 |
|
|
436 |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
437 |
|
|
438 |
pfaeta = pfaeta2(ic,angle) |
if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then |
439 |
|
pfaeta = pfaeta2(ic,angle) |
440 |
|
cc print*,pfaeta2(ic,angle) |
441 |
|
elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then |
442 |
|
pfaeta = pfaeta3(ic,angle) |
443 |
|
elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then |
444 |
|
pfaeta = pfaeta4(ic,angle) |
445 |
|
else |
446 |
|
pfaeta = cog(4,ic) |
447 |
|
endif |
448 |
|
|
449 |
else !X-view |
else !X-view |
450 |
|
|
451 |
if(abs(angle).le.10.)then |
if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then |
452 |
pfaeta = pfaeta2(ic,angle) |
pfaeta = pfaeta2(ic,angle) |
453 |
elseif(abs(angle).gt.10..and.abs(angle).le.15.)then |
elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then |
454 |
pfaeta = pfaeta3(ic,angle) |
pfaeta = pfaeta3(ic,angle) |
455 |
elseif(abs(angle).gt.15.)then |
elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then |
456 |
pfaeta = pfaeta4(ic,angle) |
pfaeta = pfaeta4(ic,angle) |
457 |
endif |
else |
458 |
|
pfaeta = cog(4,ic) |
459 |
|
endif |
460 |
|
|
461 |
endif |
endif |
462 |
|
|
463 |
c print*,'pfaeta ',pfaeta, angle |
c 100 return |
464 |
|
return |
|
100 return |
|
465 |
end |
end |
466 |
|
|
467 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
468 |
real function ris_eta(ic,angle) |
real function pfaetal(ic,angle) |
469 |
*-------------------------------------------------------------- |
*-------------------------------------------------------------- |
470 |
* this function returns the average spatial resolution |
* this function returns the position (in strip units) |
|
* (in cm) for the ETA algorithm (function pfaeta(ic,angle)) |
|
471 |
* it calls: |
* it calls: |
472 |
* - risx_eta2(angle) |
* - pfaeta2(ic,angle)+pfcorr(ic,angle) |
473 |
* - risy_eta2(angle) |
* - pfaeta3(ic,angle)+pfcorr(ic,angle) |
474 |
* - risx_eta3(angle) |
* - pfaeta4(ic,angle)+pfcorr(ic,angle) |
|
* - risx_eta4(angle) |
|
475 |
* according to the angle |
* according to the angle |
476 |
*-------------------------------------------------------------- |
*-------------------------------------------------------------- |
477 |
include 'commontracker.f' |
include 'commontracker.f' |
|
c include 'calib.f' |
|
478 |
include 'level1.f' |
include 'level1.f' |
479 |
|
include 'calib.f' |
|
c$$$ logical DEBUG |
|
|
c$$$ common/dbg/DEBUG |
|
480 |
|
|
481 |
ris_eta = 0 |
pfaetal = 0 |
482 |
|
|
483 |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
484 |
|
|
485 |
ris_eta = risy_eta2(angle) |
if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then |
486 |
if(abs(angle).gt.21.)ris_eta = risy_eta2(21.) |
pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle) |
487 |
|
cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle) |
488 |
|
elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then |
489 |
|
pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle) |
490 |
|
elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then |
491 |
|
pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle) |
492 |
|
else |
493 |
|
pfaetal = cog(4,ic) |
494 |
|
endif |
495 |
|
|
496 |
else !X-view |
else !X-view |
497 |
|
|
498 |
if(abs(angle).le.10.)then |
if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then |
499 |
ris_eta = risx_eta2(angle) |
pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle) |
500 |
elseif(abs(angle).gt.10..and.abs(angle).le.15.)then |
cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle) |
501 |
ris_eta = risx_eta3(angle) |
elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then |
502 |
elseif(abs(angle).gt.15..and.abs(angle).le.21.)then |
pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle) |
503 |
ris_eta = risx_eta4(angle) |
elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then |
504 |
elseif(abs(angle).gt.21.)then |
pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle) |
505 |
ris_eta = risx_eta4(21.) |
else |
506 |
endif |
pfaetal = cog(4,ic) |
507 |
|
endif |
508 |
|
|
509 |
endif |
endif |
510 |
|
|
511 |
|
c 100 return |
512 |
|
return |
513 |
|
end |
514 |
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
515 |
|
c real function riseta(ic,angle) |
516 |
|
real function riseta(iview,angle) |
517 |
|
*-------------------------------------------------------------- |
518 |
|
* this function returns the average spatial resolution |
519 |
|
* (in cm) for the ETA algorithm (function pfaeta(ic,angle)) |
520 |
|
* it calls: |
521 |
|
* - risxeta2(angle) |
522 |
|
* - risyeta2(angle) |
523 |
|
* - risxeta3(angle) |
524 |
|
* - risxeta4(angle) |
525 |
|
* according to the angle |
526 |
|
*-------------------------------------------------------------- |
527 |
|
include 'commontracker.f' |
528 |
|
include 'level1.f' |
529 |
|
include 'calib.f' |
530 |
|
|
531 |
c$$$ if(DEBUG)print*,'ris (ic ',ic,' ang',angle,')' |
riseta = 0 |
|
c$$$ $ ,' -->',ris_eta |
|
532 |
|
|
533 |
|
c if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
534 |
|
if(mod(iview,2).eq.1)then !Y-view |
535 |
|
|
536 |
|
|
537 |
|
if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then |
538 |
|
riseta = risyeta2(angle) |
539 |
|
elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then |
540 |
|
riseta = risy_cog(angle) !ATTENZIONE!! |
541 |
|
elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then |
542 |
|
riseta = risy_cog(angle) !ATTENZIONE!! |
543 |
|
else |
544 |
|
riseta = risy_cog(angle) |
545 |
|
endif |
546 |
|
|
547 |
100 return |
else !X-view |
548 |
|
|
549 |
|
if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then |
550 |
|
riseta = risxeta2(angle) |
551 |
|
elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then |
552 |
|
riseta = risxeta3(angle) |
553 |
|
elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then |
554 |
|
riseta = risxeta4(angle) |
555 |
|
else |
556 |
|
riseta = risx_cog(angle) |
557 |
|
endif |
558 |
|
|
559 |
|
endif |
560 |
|
|
561 |
|
|
562 |
|
c 100 return |
563 |
|
return |
564 |
end |
end |
565 |
|
|
566 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
573 |
* resolution. |
* resolution. |
574 |
* It calls the function FBAD_COG(NCOG,IC), |
* It calls the function FBAD_COG(NCOG,IC), |
575 |
* accordingto the angle |
* accordingto the angle |
576 |
|
* |
577 |
|
* >>> cosi` non e` corretto!! |
578 |
|
* >>> l'errore sulla coordinata eta si ottiene moltiplicando |
579 |
|
* >>> l'errore sulla coordinata cog per la derivata della |
580 |
|
* >>> distribuzione eta... pur sapendolo l'ho sempre ignorato... |
581 |
|
* >>> deve essere modificato!!!! |
582 |
|
* |
583 |
*------------------------------------------------------- |
*------------------------------------------------------- |
584 |
|
|
585 |
include 'commontracker.f' |
include 'commontracker.f' |
586 |
include 'level1.f' |
include 'level1.f' |
587 |
* include 'calib.f' |
include 'calib.f' |
588 |
fbad_eta = 0 |
fbad_eta = 0 |
589 |
|
|
590 |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
591 |
|
|
592 |
fbad_eta = fbad_cog(2,ic) |
if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then |
593 |
|
fbad_eta = fbad_cog(2,ic) |
594 |
|
elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then |
595 |
|
fbad_eta = fbad_cog(3,ic) |
596 |
|
elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then |
597 |
|
fbad_eta = fbad_cog(4,ic) |
598 |
|
else |
599 |
|
fbad_eta = fbad_cog(4,ic) |
600 |
|
endif |
601 |
|
|
602 |
else !X-view |
else !X-view |
603 |
|
|
604 |
if(abs(angle).le.10.)then |
if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then |
605 |
fbad_eta = fbad_cog(2,ic) |
fbad_eta = fbad_cog(2,ic) |
606 |
elseif(abs(angle).gt.10..and.abs(angle).le.15.)then |
elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then |
607 |
fbad_eta = fbad_cog(3,ic) |
fbad_eta = fbad_cog(3,ic) |
608 |
elseif(abs(angle).gt.15.)then |
elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then |
609 |
fbad_eta = fbad_cog(4,ic) |
fbad_eta = fbad_cog(4,ic) |
610 |
endif |
else |
611 |
|
fbad_eta = fbad_cog(4,ic) |
612 |
|
endif |
613 |
|
|
614 |
endif |
endif |
615 |
|
|
617 |
end |
end |
618 |
|
|
619 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
620 |
c***************************************************** |
real function pfaeta2(ic,angle) |
|
cccccc 02/02/2006 modified by Elena Vannuccini --> (1) |
|
|
c***************************************************** |
|
|
c real function pfaeta2(cog2,view,lad,angle) |
|
|
real function pfaeta2(ic,angle) !(1) |
|
621 |
*-------------------------------------------------------------- |
*-------------------------------------------------------------- |
622 |
* this function returns |
* this function returns |
623 |
* |
* |
636 |
real cog2,angle |
real cog2,angle |
637 |
integer iview,lad |
integer iview,lad |
638 |
|
|
639 |
c logical DEBUG |
iview = VIEW(ic) |
640 |
c common/dbg/DEBUG |
lad = nld(MAXS(ic),VIEW(ic)) |
641 |
|
cog2 = cog(2,ic) |
642 |
c print*,'## pfaeta2 ',ic,angle |
pfaeta2 = cog2 |
|
iview = VIEW(ic) !(1) |
|
|
lad = nld(MAXS(ic),VIEW(ic)) !(1) |
|
|
cog2 = cog(2,ic) !(1) |
|
|
pfaeta2=cog2 |
|
643 |
|
|
644 |
|
* ---------------- |
645 |
* find angular bin |
* find angular bin |
646 |
|
* ---------------- |
647 |
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
648 |
do iang=1,nangbin |
do iang=1,nangbin |
|
c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle |
|
649 |
if(angL(iang).lt.angle.and.angR(iang).ge.angle)then |
if(angL(iang).lt.angle.and.angR(iang).ge.angle)then |
650 |
iangle=iang |
iangle=iang |
651 |
goto 98 |
goto 98 |
652 |
endif |
endif |
653 |
enddo |
enddo |
654 |
if(DEBUG) |
if(DEBUG.EQ.1) |
655 |
$ print*,'pfaeta2 *** warning *** angle out of range: ',angle |
$ print*,'pfaeta2 *** warning *** angle out of range: ',angle |
656 |
if(angle.lt.angL(1))iang=1 |
if(angle.le.angL(1))iang=1 |
657 |
if(angle.gt.angR(nangbin))iang=nangbin |
if(angle.ge.angR(nangbin))iang=nangbin |
658 |
98 continue !jump here if ok |
98 continue !jump here if ok |
659 |
|
|
660 |
|
* ------------- |
661 |
|
* within +/-0.5 |
662 |
|
* ------------- |
663 |
|
|
664 |
c$$$* find extremes of interpolation |
iaddmax=10 |
|
c$$$ iflag=0 |
|
|
c$$$* -------------------------------- |
|
|
c$$$ if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then |
|
|
c$$$c print*,'pfaeta2 *** warning *** argument out of range: ',cog2 |
|
|
c$$$* goto 100 |
|
|
c$$$* ---------------------------------------------- |
|
|
c$$$* non salto piu`, ma scalo di 1 o -1 |
|
|
c$$$* nel caso si tratti di un cluster |
|
|
c$$$* in cui la strip con il segnale massimo non coincide |
|
|
c$$$* con la strip con il rapposto s/n massimo!!! |
|
|
c$$$* ---------------------------------------------- |
|
|
c$$$ if(cog2.lt.eta2(1,iang))then !temp |
|
|
c$$$ cog2=cog2+1. !temp |
|
|
c$$$ iflag=1 !temp |
|
|
c$$$ else !temp |
|
|
c$$$ cog2=cog2-1. !temp |
|
|
c$$$ iflag=-1 !temp |
|
|
c$$$ endif !temp |
|
|
c$$$c print*,'shifted >>> ',cog2 |
|
|
c$$$ endif |
|
|
|
|
665 |
iadd=0 |
iadd=0 |
666 |
10 continue |
10 continue |
667 |
if(cog2.lt.eta2(1,iang))then |
if(cog2.lt.eta2(1,iang))then |
668 |
cog2 = cog2 + 1 |
cog2 = cog2 + 1 |
669 |
iadd = iadd + 1 |
iadd = iadd + 1 |
670 |
|
if(iadd>iaddmax)goto 111 |
671 |
goto 10 |
goto 10 |
672 |
endif |
endif |
673 |
20 continue |
20 continue |
674 |
if(cog2.gt.eta2(netaval,iang))then |
if(cog2.gt.eta2(netaval,iang))then |
675 |
cog2 = cog2 - 1 |
cog2 = cog2 - 1 |
676 |
iadd = iadd - 1 |
iadd = iadd - 1 |
677 |
|
if(iadd<-1*iaddmax)goto 111 |
678 |
goto 20 |
goto 20 |
679 |
endif |
endif |
680 |
|
goto 1111 |
681 |
|
111 continue |
682 |
|
if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster' |
683 |
|
if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)' |
684 |
|
cog2=0 |
685 |
|
1111 continue |
686 |
|
|
687 |
* -------------------------------- |
* -------------------------------- |
688 |
c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 |
c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 |
717 |
c$$$ cog2=cog2+1. !temp |
c$$$ cog2=cog2+1. !temp |
718 |
c$$$ endif |
c$$$ endif |
719 |
|
|
720 |
if(DEBUG)print*,'ETA2 (ic ',ic,' ang',angle,')' |
if(DEBUG.EQ.1)print*,'ETA2 (ic ',ic,' ang',angle,')' |
721 |
$ ,cog2-iadd,' -->',pfaeta2 |
$ ,cog2-iadd,' -->',pfaeta2 |
722 |
|
|
723 |
|
|
724 |
100 return |
c 100 return |
725 |
|
return |
726 |
end |
end |
727 |
|
|
728 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
|
c***************************************************** |
|
|
cccccc 02/02/2006 modified by Elena Vannuccini --> (1) |
|
|
c***************************************************** |
|
|
c real function pfaeta3(cog3,view,lad,angle) |
|
729 |
real function pfaeta3(ic,angle) !(1) |
real function pfaeta3(ic,angle) !(1) |
730 |
*-------------------------------------------------------------- |
*-------------------------------------------------------------- |
731 |
* this function returns |
* this function returns |
745 |
real cog3,angle |
real cog3,angle |
746 |
integer iview,lad |
integer iview,lad |
747 |
|
|
|
c logical DEBUG |
|
|
c common/dbg/DEBUG |
|
748 |
|
|
749 |
c print*,'## pfaeta3 ',ic,angle |
iview = VIEW(ic) |
750 |
|
lad = nld(MAXS(ic),VIEW(ic)) |
751 |
iview = VIEW(ic) !(1) |
cog3 = cog(3,ic) |
752 |
lad = nld(MAXS(ic),VIEW(ic)) !(1) |
cc = cog3 |
753 |
cog3 = cog(3,ic) !(1) |
cog3 = cc |
754 |
pfaeta3=cog3 |
pfaeta3=cog3 |
755 |
|
|
756 |
|
* ---------------- |
757 |
* find angular bin |
* find angular bin |
758 |
|
* ---------------- |
759 |
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
760 |
do iang=1,nangbin |
do iang=1,nangbin |
761 |
c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle |
c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle |
764 |
goto 98 |
goto 98 |
765 |
endif |
endif |
766 |
enddo |
enddo |
767 |
if(DEBUG) |
if(DEBUG.EQ.1) |
768 |
$ print*,'pfaeta3 *** warning *** angle out of range: ',angle |
$ print*,'pfaeta3 *** warning *** angle out of range: ',angle |
769 |
if(angle.lt.angL(1))iang=1 |
if(angle.le.angL(1))iang=1 |
770 |
if(angle.gt.angR(nangbin))iang=nangbin |
if(angle.ge.angR(nangbin))iang=nangbin |
771 |
98 continue !jump here if ok |
98 continue !jump here if ok |
772 |
|
|
773 |
|
* ------------- |
774 |
|
* within +/-0.5 |
775 |
|
* ------------- |
776 |
|
|
777 |
c$$$* find extremes of interpolation |
iaddmax=10 |
|
c$$$ iflag=0 |
|
|
c$$$* -------------------------------- |
|
|
c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then |
|
|
c$$$* ---------------------------------------------- |
|
|
c$$$* non salto piu`, ma scalo di 1 o -1 |
|
|
c$$$* nel caso si tratti di un cluster |
|
|
c$$$* in cui la strip con il segnale massimo non coincide |
|
|
c$$$* con la strip con il rapposto s/n massimo!!! |
|
|
c$$$* ---------------------------------------------- |
|
|
c$$$ if(cog2.lt.eta2(1,iang))then !temp |
|
|
c$$$ cog2=cog2+1. !temp |
|
|
c$$$ iflag=1 !temp |
|
|
c$$$ else !temp |
|
|
c$$$ cog2=cog2-1. !temp |
|
|
c$$$ iflag=-1 !temp |
|
|
c$$$ endif !temp |
|
|
c$$$c print*,'shifted >>> ',cog2 |
|
|
c$$$ endif |
|
|
|
|
|
|
|
778 |
iadd=0 |
iadd=0 |
779 |
10 continue |
10 continue |
780 |
if(cog3.lt.eta3(1,iang))then |
if(cog3.lt.eta3(1,iang))then |
781 |
cog3 = cog3 + 1 |
cog3 = cog3 + 1. |
782 |
iadd = iadd + 1 |
iadd = iadd + 1 |
783 |
|
if(iadd>iaddmax) goto 111 |
784 |
goto 10 |
goto 10 |
785 |
endif |
endif |
786 |
20 continue |
20 continue |
787 |
if(cog3.gt.eta3(netaval,iang))then |
if(cog3.gt.eta3(netaval,iang))then |
788 |
cog3 = cog3 - 1 |
cog3 = cog3 - 1. |
789 |
iadd = iadd - 1 |
iadd = iadd - 1 |
790 |
|
if(iadd<-1*iaddmax) goto 111 |
791 |
goto 20 |
goto 20 |
792 |
endif |
endif |
793 |
|
goto 1111 |
794 |
|
111 continue |
795 |
|
if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster' |
796 |
|
if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)' |
797 |
|
cog3=0 |
798 |
|
1111 continue |
799 |
|
|
800 |
* -------------------------------- |
* -------------------------------- |
801 |
c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 |
c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 |
821 |
pfaeta3 = AA*cog3+BB |
pfaeta3 = AA*cog3+BB |
822 |
pfaeta3 = pfaeta3 - iadd |
pfaeta3 = pfaeta3 - iadd |
823 |
|
|
|
c$$$ if(iflag.eq.1)then |
|
|
c$$$ pfaeta2=pfaeta2-1. !temp |
|
|
c$$$ cog2=cog2-1. !temp |
|
|
c$$$ endif |
|
|
c$$$ if(iflag.eq.-1)then |
|
|
c$$$ pfaeta2=pfaeta2+1. !temp |
|
|
c$$$ cog2=cog2+1. !temp |
|
|
c$$$ endif |
|
824 |
|
|
825 |
if(DEBUG)print*,'ETA3 (ic ',ic,' ang',angle,')' |
if(DEBUG.EQ.1)print*,'ETA3 (ic ',ic,' ang',angle,')' |
826 |
$ ,cog3-iadd,' -->',pfaeta3 |
$ ,cog3-iadd,' -->',pfaeta3 |
827 |
|
|
828 |
100 return |
c 100 return |
829 |
|
return |
830 |
end |
end |
831 |
|
|
832 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
833 |
c***************************************************** |
real function pfaeta4(ic,angle) |
|
cccccc 02/02/2006 modified by Elena Vannuccini --> (1) |
|
|
c***************************************************** |
|
|
c real function pfaeta4(cog4,view,lad,angle) |
|
|
real function pfaeta4(ic,angle) !(1) |
|
834 |
*-------------------------------------------------------------- |
*-------------------------------------------------------------- |
835 |
* this function returns |
* this function returns |
836 |
* |
* |
849 |
real cog4,angle |
real cog4,angle |
850 |
integer iview,lad |
integer iview,lad |
851 |
|
|
|
c logical DEBUG |
|
|
c common/dbg/DEBUG |
|
852 |
|
|
853 |
c print*,'## pfaeta4 ',ic,angle |
iview = VIEW(ic) |
854 |
|
lad = nld(MAXS(ic),VIEW(ic)) |
855 |
iview = VIEW(ic) !(1) |
cog4=cog(4,ic) |
|
lad = nld(MAXS(ic),VIEW(ic)) !(1) |
|
|
cog4=cog(4,ic) !(1) |
|
856 |
pfaeta4=cog4 |
pfaeta4=cog4 |
857 |
|
|
858 |
|
* ---------------- |
859 |
* find angular bin |
* find angular bin |
860 |
|
* ---------------- |
861 |
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
862 |
do iang=1,nangbin |
do iang=1,nangbin |
863 |
c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle |
c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle |
866 |
goto 98 |
goto 98 |
867 |
endif |
endif |
868 |
enddo |
enddo |
869 |
if(DEBUG) |
if(DEBUG.EQ.1) |
870 |
$ print*,'pfaeta4 *** warning *** angle out of range: ',angle |
$ print*,'pfaeta4 *** warning *** angle out of range: ',angle |
871 |
if(angle.lt.angL(1))iang=1 |
if(angle.le.angL(1))iang=1 |
872 |
if(angle.gt.angR(nangbin))iang=nangbin |
if(angle.ge.angR(nangbin))iang=nangbin |
873 |
98 continue !jump here if ok |
98 continue !jump here if ok |
874 |
|
|
875 |
|
* ------------- |
876 |
|
* within +/-0.5 |
877 |
|
* ------------- |
878 |
|
|
879 |
c$$$* find extremes of interpolation |
iaddmax=10 |
|
c$$$ iflag=0 |
|
|
c$$$* -------------------------------- |
|
|
c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then |
|
|
c$$$* ---------------------------------------------- |
|
|
c$$$* non salto piu`, ma scalo di 1 o -1 |
|
|
c$$$* nel caso si tratti di un cluster |
|
|
c$$$* in cui la strip con il segnale massimo non coincide |
|
|
c$$$* con la strip con il rapposto s/n massimo!!! |
|
|
c$$$* ---------------------------------------------- |
|
|
c$$$ if(cog2.lt.eta2(1,iang))then !temp |
|
|
c$$$ cog2=cog2+1. !temp |
|
|
c$$$ iflag=1 !temp |
|
|
c$$$ else !temp |
|
|
c$$$ cog2=cog2-1. !temp |
|
|
c$$$ iflag=-1 !temp |
|
|
c$$$ endif !temp |
|
|
c$$$c print*,'shifted >>> ',cog2 |
|
|
c$$$ endif |
|
|
|
|
|
|
|
880 |
iadd=0 |
iadd=0 |
881 |
10 continue |
10 continue |
882 |
if(cog4.lt.eta4(1,iang))then |
if(cog4.lt.eta4(1,iang))then |
883 |
cog4 = cog4 + 1 |
cog4 = cog4 + 1 |
884 |
iadd = iadd + 1 |
iadd = iadd + 1 |
885 |
|
if(iadd>iaddmax)goto 111 |
886 |
goto 10 |
goto 10 |
887 |
endif |
endif |
888 |
20 continue |
20 continue |
889 |
if(cog4.gt.eta4(netaval,iang))then |
if(cog4.gt.eta4(netaval,iang))then |
890 |
cog4 = cog4 - 1 |
cog4 = cog4 - 1 |
891 |
iadd = iadd - 1 |
iadd = iadd - 1 |
892 |
|
if(iadd<-1*iaddmax)goto 111 |
893 |
goto 20 |
goto 20 |
894 |
endif |
endif |
895 |
|
goto 1111 |
896 |
|
111 continue |
897 |
|
if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster' |
898 |
|
if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)' |
899 |
|
cog4=0 |
900 |
|
1111 continue |
901 |
|
|
902 |
* -------------------------------- |
* -------------------------------- |
903 |
c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 |
c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2 |
932 |
c$$$ cog2=cog2+1. !temp |
c$$$ cog2=cog2+1. !temp |
933 |
c$$$ endif |
c$$$ endif |
934 |
|
|
935 |
if(DEBUG)print*,'ETA4 (ic ',ic,' ang',angle,')' |
if(DEBUG.EQ.1)print*,'ETA4 (ic ',ic,' ang',angle,')' |
936 |
$ ,cog4-iadd,' -->',pfaeta4 |
$ ,cog4-iadd,' -->',pfaeta4 |
937 |
|
|
938 |
100 return |
c 100 return |
939 |
|
return |
940 |
end |
end |
941 |
|
|
942 |
|
|
943 |
|
|
944 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
|
real function cog0(ncog,ic) |
|
|
*------------------------------------------------- |
|
|
* this function returns |
|
|
* |
|
|
* - the Center-Of-Gravity of the cluster IC |
|
|
* evaluated using NCOG strips, |
|
|
* calculated relative to MAXS(IC) |
|
|
* |
|
|
* - zero in case that not enough strips |
|
|
* have a positive signal |
|
|
* |
|
|
* NOTE: |
|
|
* This is the old definition, used by Straulino. |
|
|
* The new routine, according to Landi, |
|
|
* is COG(NCOG,IC) |
|
|
*------------------------------------------------- |
|
|
|
|
|
|
|
|
include 'commontracker.f' |
|
|
include 'level1.f' |
|
|
|
|
|
* --> signal of the central strip |
|
|
sc = CLSIGNAL(INDMAX(ic)) !center |
|
|
|
|
|
* signal of adjacent strips |
|
|
* --> left |
|
|
sl1 = 0 !left 1 |
|
|
if( |
|
|
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
|
|
$ ) |
|
|
$ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) |
|
|
|
|
|
sl2 = 0 !left 2 |
|
|
if( |
|
|
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
|
|
$ ) |
|
|
$ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) |
|
|
|
|
|
* --> right |
|
|
sr1 = 0 !right 1 |
|
|
if( |
|
|
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) |
|
|
$ .or. |
|
|
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) |
|
|
$ ) |
|
|
$ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) |
|
|
|
|
|
sr2 = 0 !right 2 |
|
|
if( |
|
|
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) |
|
|
$ .or. |
|
|
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) |
|
|
$ ) |
|
|
$ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) |
|
|
|
|
|
************************************************************ |
|
|
* COG computation |
|
|
************************************************************ |
|
|
|
|
|
c print*,sl2,sl1,sc,sr1,sr2 |
|
|
|
|
|
COG = 0. |
|
|
|
|
|
if(sl1.gt.sr1.and.sl1.gt.0.)then |
|
|
|
|
|
if(ncog.eq.2.and.sl1.ne.0)then |
|
|
COG = -sl1/(sl1+sc) |
|
|
elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then |
|
|
COG = (sr1-sl1)/(sl1+sc+sr1) |
|
|
elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then |
|
|
COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) |
|
|
else |
|
|
COG = 0. |
|
|
endif |
|
|
|
|
|
elseif(sl1.le.sr1.and.sr1.gt.0.)then |
|
|
|
|
|
if(ncog.eq.2.and.sr1.ne.0)then |
|
|
COG = sr1/(sc+sr1) |
|
|
elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then |
|
|
COG = (sr1-sl1)/(sl1+sc+sr1) |
|
|
elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then |
|
|
COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) |
|
|
else |
|
|
COG = 0. |
|
|
endif |
|
|
|
|
|
endif |
|
|
|
|
|
COG0 = COG |
|
|
|
|
|
c print *,ncog,ic,cog,'/////////////' |
|
|
|
|
|
return |
|
|
end |
|
|
|
|
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
|
945 |
real function cog(ncog,ic) |
real function cog(ncog,ic) |
946 |
*------------------------------------------------- |
*------------------------------------------------- |
947 |
* this function returns |
* this function returns |
961 |
include 'calib.f' |
include 'calib.f' |
962 |
include 'level1.f' |
include 'level1.f' |
963 |
|
|
|
c logical DEBUG |
|
|
c common/dbg/DEBUG |
|
964 |
|
|
965 |
|
|
966 |
if (ncog.gt.0) then |
if (ncog.gt.0) then |
971 |
* --> signal of the central strip |
* --> signal of the central strip |
972 |
sc = CLSIGNAL(INDMAX(ic)) !center |
sc = CLSIGNAL(INDMAX(ic)) !center |
973 |
* signal of adjacent strips |
* signal of adjacent strips |
974 |
sl1 = 0 !left 1 |
sl1 = -9999. !left 1 |
975 |
if( |
if( |
976 |
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
977 |
$ ) |
$ ) |
978 |
$ sl1 = CLSIGNAL(INDMAX(ic)-1) |
$ sl1 = CLSIGNAL(INDMAX(ic)-1) |
979 |
|
|
980 |
sl2 = 0 !left 2 |
sl2 = -9999. !left 2 |
981 |
if( |
if( |
982 |
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
983 |
$ ) |
$ ) |
984 |
$ sl2 = CLSIGNAL(INDMAX(ic)-2) |
$ sl2 = CLSIGNAL(INDMAX(ic)-2) |
985 |
|
|
986 |
sr1 = 0 !right 1 |
sr1 = -9999. !right 1 |
987 |
if( |
if( |
988 |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) |
989 |
$ .or. |
$ .or. |
991 |
$ ) |
$ ) |
992 |
$ sr1 = CLSIGNAL(INDMAX(ic)+1) |
$ sr1 = CLSIGNAL(INDMAX(ic)+1) |
993 |
|
|
994 |
sr2 = 0 !right 2 |
sr2 = -9999. !right 2 |
995 |
if( |
if( |
996 |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) |
997 |
$ .or. |
$ .or. |
1001 |
|
|
1002 |
COG = 0. |
COG = 0. |
1003 |
|
|
1004 |
c print*,'## ',sl2,sl1,sc,sr1,sr2 |
c print *,'## ',sl2,sl1,sc,sr1,sr2 |
1005 |
|
|
1006 |
|
c ============================================================== |
1007 |
if(ncog.eq.1)then |
if(ncog.eq.1)then |
1008 |
COG = 0. |
COG = 0. |
1009 |
|
if(sr1.gt.sc)cog=1. |
1010 |
|
if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. |
1011 |
|
c ============================================================== |
1012 |
elseif(ncog.eq.2)then |
elseif(ncog.eq.2)then |
1013 |
|
COG = 0. |
1014 |
if(sl1.gt.sr1)then |
if(sl1.gt.sr1)then |
1015 |
if((sl1+sc).ne.0)COG = -sl1/(sl1+sc) |
if((sl1+sc).ne.0)COG = -sl1/(sl1+sc) |
1016 |
elseif(sl1.le.sr1)then |
elseif(sl1.lt.sr1)then |
1017 |
if((sc+sr1).ne.0)COG = sr1/(sc+sr1) |
if((sc+sr1).ne.0)COG = sr1/(sc+sr1) |
1018 |
endif |
elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then |
1019 |
|
if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1) |
1020 |
|
$ .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) |
1021 |
|
if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1) |
1022 |
|
$ .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) |
1023 |
|
endif |
1024 |
|
c if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic) |
1025 |
|
c $ ,' : ',sl2,sl1,sc,sr1,sr2 |
1026 |
|
c ============================================================== |
1027 |
elseif(ncog.eq.3)then |
elseif(ncog.eq.3)then |
1028 |
if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1) |
COG = 0 |
1029 |
|
sss = sc |
1030 |
|
if( sl1.ne.-9999. )COG = COG-sl1 |
1031 |
|
if( sl1.ne.-9999. )sss = sss+sl1 |
1032 |
|
if( sr1.ne.-9999. )COG = COG+sr1 |
1033 |
|
if( sr1.ne.-9999. )sss = sss+sr1 |
1034 |
|
if(sss.ne.0)COG=COG/sss |
1035 |
|
|
1036 |
|
c if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1) |
1037 |
|
c if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic) |
1038 |
|
c $ ,' : ',sl2,sl1,sc,sr1,sr2 |
1039 |
|
c ============================================================== |
1040 |
elseif(ncog.eq.4)then |
elseif(ncog.eq.4)then |
1041 |
|
|
1042 |
|
COG = 0 |
1043 |
|
sss = sc |
1044 |
|
if( sl1.ne.-9999. )COG = COG-sl1 |
1045 |
|
if( sl1.ne.-9999. )sss = sss+sl1 |
1046 |
|
if( sr1.ne.-9999. )COG = COG+sr1 |
1047 |
|
if( sr1.ne.-9999. )sss = sss+sr1 |
1048 |
if(sl2.gt.sr2)then |
if(sl2.gt.sr2)then |
1049 |
if((sl2+sl1+sc+sr1).ne.0) |
if((sl2+sss).ne.0) |
1050 |
$ COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) |
$ COG = (COG-2*sl2)/(sl2+sss) |
1051 |
elseif(sl2.le.sr2)then |
elseif(sl2.lt.sr2)then |
1052 |
if((sl2+sl1+sc+sr1).ne.0) |
if((sr2+sss).ne.0) |
1053 |
$ COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) |
$ COG = (2*sr2+COG)/(sr2+sss) |
1054 |
|
elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then |
1055 |
|
if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2) |
1056 |
|
$ .and.(sl2+sss).ne.0 ) |
1057 |
|
$ cog = (cog-2*sl2)/(sl2+sss) |
1058 |
|
if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2) |
1059 |
|
$ .and.(sr2+sss).ne.0 ) |
1060 |
|
$ cog = (2*sr2+cog)/(sr2+sss) |
1061 |
endif |
endif |
1062 |
|
c ============================================================== |
1063 |
|
elseif(ncog.eq.5)then |
1064 |
|
COG = 0 |
1065 |
|
sss = sc |
1066 |
|
if( sl1.ne.-9999. )COG = COG-sl1 |
1067 |
|
if( sl1.ne.-9999. )sss = sss+sl1 |
1068 |
|
if( sr1.ne.-9999. )COG = COG+sr1 |
1069 |
|
if( sr1.ne.-9999. )sss = sss+sr1 |
1070 |
|
if( sl2.ne.-9999. )COG = COG-2*sl2 |
1071 |
|
if( sl2.ne.-9999. )sss = sss+sl2 |
1072 |
|
if( sr2.ne.-9999. )COG = COG+2*sr2 |
1073 |
|
if( sr2.ne.-9999. )sss = sss+sr2 |
1074 |
|
if(sss.ne.0)COG=COG/sss |
1075 |
else |
else |
1076 |
print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG |
print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG |
1077 |
print*,' (NCOG must be 0-4)' |
$ ,' not implemented' |
1078 |
COG = 0. |
COG = 0. |
1079 |
endif |
endif |
1080 |
|
|
1088 |
iv=VIEW(ic) |
iv=VIEW(ic) |
1089 |
if(mod(iv,2).eq.1)incut=incuty |
if(mod(iv,2).eq.1)incut=incuty |
1090 |
if(mod(iv,2).eq.0)incut=incutx |
if(mod(iv,2).eq.0)incut=incutx |
|
|
|
1091 |
istart = INDSTART(IC) |
istart = INDSTART(IC) |
1092 |
istop = TOTCLLENGTH |
istop = TOTCLLENGTH |
1093 |
if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 |
if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 |
1094 |
COG = 0 |
COG = 0 |
1095 |
|
SGN = 0. |
1096 |
mu = 0 |
mu = 0 |
1097 |
do i = istart,istop |
c print*,'-------' |
1098 |
|
do i = INDMAX(IC),istart,-1 |
1099 |
|
ipos = i-INDMAX(ic) |
1100 |
|
cut = incut*CLSIGMA(i) |
1101 |
|
if(CLSIGNAL(i).ge.cut)then |
1102 |
|
COG = COG + ipos*CLSIGNAL(i) |
1103 |
|
SGN = SGN + CLSIGNAL(i) |
1104 |
|
mu = mu + 1 |
1105 |
|
c print*,ipos,CLSIGNAL(i) |
1106 |
|
else |
1107 |
|
goto 10 |
1108 |
|
endif |
1109 |
|
enddo |
1110 |
|
10 continue |
1111 |
|
do i = INDMAX(IC)+1,istop |
1112 |
ipos = i-INDMAX(ic) |
ipos = i-INDMAX(ic) |
|
cc ivk = nvk(MAXS(ic)+ipos) |
|
|
cc is = nst(MAXS(ic)+ipos) |
|
|
* print*,'******************',istart,istop,ipos |
|
|
* $ ,MAXS(ic)+ipos,iv,ivk,is |
|
|
cc cut = incut*SIGMA(iv,ivk,is) |
|
1113 |
cut = incut*CLSIGMA(i) |
cut = incut*CLSIGMA(i) |
|
c if(SIGMA(iv,ivk,is).ne.CLSIGMA(i)) |
|
|
c $ print*,'cog(0,ic) --> hai fatto qualche cazzata' |
|
1114 |
if(CLSIGNAL(i).ge.cut)then |
if(CLSIGNAL(i).ge.cut)then |
1115 |
COG = COG + ipos*CLSIGNAL(i) |
COG = COG + ipos*CLSIGNAL(i) |
1116 |
|
SGN = SGN + CLSIGNAL(i) |
1117 |
mu = mu + 1 |
mu = mu + 1 |
1118 |
c print*,ipos,CLSIGNAL(i),incut,cut |
c print*,ipos,CLSIGNAL(i) |
1119 |
|
else |
1120 |
|
goto 20 |
1121 |
endif |
endif |
1122 |
enddo |
enddo |
1123 |
if(DEDX(ic).le.0)then |
20 continue |
1124 |
print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic) |
if(SGN.le.0)then |
1125 |
|
print*,'cog(0,ic) --> ic, dedx ',ic,SGN |
1126 |
print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) |
print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) |
1127 |
print*,(CLSIGNAL(i),i=istart,istop) |
print*,(CLSIGNAL(i),i=istart,istop) |
1128 |
print*,'cog(0,ic) --> NOT EVALUATED ' |
c print*,'cog(0,ic) --> NOT EVALUATED ' |
1129 |
else |
else |
1130 |
COG=COG/DEDX(ic) |
COG=COG/SGN |
1131 |
endif |
endif |
1132 |
c if(DEBUG)print*,'COG (ic ',ic,' m',mu,')' |
c print*,'-------' |
|
c $ ,cog |
|
1133 |
|
|
1134 |
else |
else |
1135 |
|
|
1142 |
|
|
1143 |
c print *,'## cog ',ncog,ic,cog,'/////////////' |
c print *,'## cog ',ncog,ic,cog,'/////////////' |
1144 |
|
|
1145 |
|
if(COG.lt.-0.75.or.COG.gt.+0.75)then |
1146 |
|
if(DEBUG.eq.1) |
1147 |
|
$ print*,'cog *** warning *** anomalous cluster ??? --> ' |
1148 |
|
if(DEBUG.eq.1) |
1149 |
|
$ print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG |
1150 |
|
endif |
1151 |
|
|
1152 |
|
|
1153 |
return |
return |
1154 |
end |
end |
1155 |
|
|
1156 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1157 |
|
|
1158 |
real function fbad_cog(ncog,ic) |
real function fbad_cog(ncog,ic) |
1159 |
*------------------------------------------------------- |
*------------------------------------------------------- |
1160 |
* this function returns a factor that takes into |
* this function returns a factor that takes into |
1170 |
include 'calib.f' |
include 'calib.f' |
1171 |
|
|
1172 |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
1173 |
f = 4. |
si = 8.4 !average good-strip noise |
1174 |
si = 8.4 |
f = 4. !average bad-strip noise: f*si |
1175 |
incut=incuty |
incut=incuty |
1176 |
else !X-view |
else !X-view |
1177 |
f = 6. |
si = 3.9 !average good-strip noise |
1178 |
si = 3.9 |
f = 6. !average bad-strip noise: f*si |
1179 |
incut=incutx |
incut=incutx |
1180 |
endif |
endif |
1181 |
|
|
1186 |
* --> signal of the central strip |
* --> signal of the central strip |
1187 |
sc = CLSIGNAL(INDMAX(ic)) !center |
sc = CLSIGNAL(INDMAX(ic)) !center |
1188 |
fsc = 1 |
fsc = 1 |
1189 |
if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)fsc=f |
c if( CLBAD(INDMAX(ic)).eq.0 )fsc=f |
1190 |
|
fsc = clsigma(INDMAX(ic))/si |
1191 |
* --> signal of adjacent strips |
* --> signal of adjacent strips |
1192 |
sl1 = 0 !left 1 |
sl1 = 0 !left 1 |
1193 |
fsl1 = 1 !left 1 |
fsl1 = 1 !left 1 |
1195 |
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
1196 |
$ )then |
$ )then |
1197 |
sl1 = CLSIGNAL(INDMAX(ic)-1) |
sl1 = CLSIGNAL(INDMAX(ic)-1) |
1198 |
if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f |
c if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f |
1199 |
c else |
fsl1 = clsigma(INDMAX(ic)-1)/si |
|
c fsl1 = 0 |
|
1200 |
endif |
endif |
1201 |
|
|
1202 |
sl2 = 0 !left 2 |
sl2 = 0 !left 2 |
1205 |
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
1206 |
$ )then |
$ )then |
1207 |
sl2 = CLSIGNAL(INDMAX(ic)-2) |
sl2 = CLSIGNAL(INDMAX(ic)-2) |
1208 |
if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f |
c if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f |
1209 |
c else |
fsl2 = clsigma(INDMAX(ic)-2)/si |
|
c fsl2 = 0 |
|
1210 |
endif |
endif |
1211 |
sr1 = 0 !right 1 |
sr1 = 0 !right 1 |
1212 |
fsr1 = 1 !right 1 |
fsr1 = 1 !right 1 |
1216 |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) |
1217 |
$ )then |
$ )then |
1218 |
sr1 = CLSIGNAL(INDMAX(ic)+1) |
sr1 = CLSIGNAL(INDMAX(ic)+1) |
1219 |
if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f |
c if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f |
1220 |
c else |
fsr1 = clsigma(INDMAX(ic)+1)/si |
|
c fsr1 = 0 |
|
1221 |
endif |
endif |
1222 |
sr2 = 0 !right 2 |
sr2 = 0 !right 2 |
1223 |
fsr2 = 1 !right 2 |
fsr2 = 1 !right 2 |
1227 |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) |
1228 |
$ )then |
$ )then |
1229 |
sr2 = CLSIGNAL(INDMAX(ic)+2) |
sr2 = CLSIGNAL(INDMAX(ic)+2) |
1230 |
if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f |
c if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f |
1231 |
c else |
fsr2 = clsigma(INDMAX(ic)+2)/si |
|
c fsr2 = 0 |
|
1232 |
endif |
endif |
1233 |
|
|
1234 |
|
|
1235 |
|
|
1236 |
************************************************************ |
************************************************************ |
1237 |
* COG computation |
* COG2-3-4 computation |
1238 |
************************************************************ |
************************************************************ |
1239 |
|
|
1240 |
c print*,sl2,sl1,sc,sr1,sr2 |
c print*,sl2,sl1,sc,sr1,sr2 |
1241 |
|
|
1242 |
COG = 0. |
vCOG = cog(ncog,ic)!0. |
1243 |
|
|
1244 |
if(ncog.eq.2)then |
if(ncog.eq.2)then |
1245 |
if(sl1.gt.sr1)then |
if(sl1.gt.sr1)then |
1246 |
COG = -sl1/(sl1+sc) |
c COG = -sl1/(sl1+sc) |
1247 |
fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2) |
fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2) |
1248 |
fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2) |
fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2) |
1249 |
elseif(sl1.le.sr1)then |
elseif(sl1.le.sr1)then |
1250 |
COG = sr1/(sc+sr1) |
c COG = sr1/(sc+sr1) |
1251 |
fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2) |
fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) |
1252 |
fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2) |
fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2) |
1253 |
endif |
endif |
1254 |
elseif(ncog.eq.3)then |
elseif(ncog.eq.3)then |
1255 |
COG = (sr1-sl1)/(sl1+sc+sr1) |
c COG = (sr1-sl1)/(sl1+sc+sr1) |
1256 |
fbad_cog = |
fbad_cog = |
1257 |
$ (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2) |
$ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) |
1258 |
fbad_cog = |
fbad_cog = |
1259 |
$ fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2) |
$ fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2) |
1260 |
elseif(ncog.eq.4)then |
elseif(ncog.eq.4)then |
1261 |
if(sl2.gt.sr2)then |
if(sl2.gt.sr2)then |
1262 |
COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) |
c COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) |
1263 |
fbad_cog = |
fbad_cog = |
1264 |
$ (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2 |
$ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2 |
1265 |
$ +fsc*(-COG)**2+fsr1*(1-COG)**2) |
$ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) |
1266 |
fbad_cog = |
fbad_cog = |
1267 |
$ fbad_cog / ((-2-COG)**2+(-1-COG)**2 |
$ fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2 |
1268 |
$ +(-COG)**2+(1-COG)**2) |
$ +(-vCOG)**2+(1-vCOG)**2) |
1269 |
elseif(sl2.le.sr2)then |
elseif(sl2.le.sr2)then |
1270 |
COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) |
c COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1) |
1271 |
fbad_cog = |
fbad_cog = |
1272 |
$ (fsl1*(-1-COG)**2 |
$ (fsl1*(-1-vCOG)**2 |
1273 |
$ +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2) |
$ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2) |
1274 |
fbad_cog = |
fbad_cog = |
1275 |
$ fbad_cog / ((-1-COG)**2 |
$ fbad_cog / ((-1-vCOG)**2 |
1276 |
$ +(-COG)**2+(1-COG)**2+(2-COG)**2) |
$ +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2) |
1277 |
endif |
endif |
1278 |
else |
else |
1279 |
print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG |
print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG |
1280 |
print*,' (NCOG must be <= 4)' |
print*,' (NCOG must be <= 4)' |
1281 |
COG = 0. |
c COG = 0. |
1282 |
endif |
endif |
1283 |
|
|
1284 |
elseif(ncog.eq.0)then |
elseif(ncog.eq.0)then |
1285 |
|
* ========================= |
1286 |
|
* COG computation |
1287 |
|
* ========================= |
1288 |
|
|
1289 |
iv=VIEW(ic) |
vCOG = cog(0,ic) |
1290 |
istart=INDSTART(IC) |
|
1291 |
istop=TOTCLLENGTH |
iv = VIEW(ic) |
1292 |
if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 |
istart = INDSTART(IC) |
1293 |
COG=0. |
istop = TOTCLLENGTH |
1294 |
SNU=0. |
if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1 |
1295 |
SDE=0. |
SGN = 0. |
1296 |
do i=istart,istop |
SNU = 0. |
1297 |
ipos=i-INDMAX(ic) |
SDE = 0. |
1298 |
il=nvk(MAXS(ic)+ipos) |
|
1299 |
is=nst(MAXS(ic)+ipos) |
do i=INDMAX(IC),istart,-1 |
1300 |
cut=incut*SIGMA(iv,il,is) |
ipos = i-INDMAX(ic) |
1301 |
|
cut = incut*CLSIGMA(i) |
1302 |
if(CLSIGNAL(i).gt.cut)then |
if(CLSIGNAL(i).gt.cut)then |
1303 |
COG = COG + ipos*CLSIGNAL(i) |
fs = clsigma(i)/si |
1304 |
endif |
SNU = SNU + fs*(ipos-vCOG)**2 |
1305 |
|
SDE = SDE + (ipos-vCOG)**2 |
1306 |
|
else |
1307 |
|
goto 10 |
1308 |
|
endif |
1309 |
enddo |
enddo |
1310 |
COG=COG/DEDX(ic) |
10 continue |
1311 |
do i=istart,istop |
do i=INDMAX(IC)+1,istop |
1312 |
ipos=i-INDMAX(ic) |
ipos = i-INDMAX(ic) |
1313 |
il=nvk(MAXS(ic)+ipos) |
cut = incut*CLSIGMA(i) |
|
is=nst(MAXS(ic)+ipos) |
|
|
cut=incut*SIGMA(iv,il,is) |
|
1314 |
if(CLSIGNAL(i).gt.cut)then |
if(CLSIGNAL(i).gt.cut)then |
1315 |
fs=1 |
fs = clsigma(i)/si |
1316 |
if(BAD(iv,il,is).eq.0)fs=f |
SNU = SNU + fs*(ipos-vCOG)**2 |
1317 |
SNU = SNU + fs*(ipos-COG)**2 |
SDE = SDE + (ipos-vCOG)**2 |
1318 |
SDE = SDE + (ipos-COG)**2 |
else |
1319 |
|
goto 20 |
1320 |
endif |
endif |
1321 |
enddo |
enddo |
1322 |
if(SDE.ne.0)FBAD_COG=SNU/SDE |
20 continue |
1323 |
|
if(SDE.ne.0)then |
1324 |
|
FBAD_COG=SNU/SDE |
1325 |
|
else |
1326 |
|
|
1327 |
|
endif |
1328 |
|
|
1329 |
else |
else |
1330 |
|
|
1343 |
|
|
1344 |
|
|
1345 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1346 |
real function fbad_cog0(ncog,ic) |
|
1347 |
|
real function riscogtheor(ncog,ic) |
1348 |
*------------------------------------------------------- |
*------------------------------------------------------- |
|
* this function returns a factor that takes into |
|
|
* account deterioration of the spatial resolution |
|
|
* in the case BAD strips are included in the cluster. |
|
|
* This factor should multiply the nominal spatial |
|
|
* resolution. |
|
1349 |
* |
* |
1350 |
* NB!!! |
* this function returns the expected resolution |
1351 |
* (this is the old version. It consider only the two |
* obtained by propagating the strip noise |
1352 |
* strips with the greatest signal. The new one is |
* to the center-of-gravity coordinate |
1353 |
* fbad_cog(ncog,ic) ) |
* |
1354 |
* |
* ncog = n.strip used in the coordinate evaluation |
1355 |
|
* (ncog=0 => all strips above threshold) |
1356 |
|
* |
1357 |
*------------------------------------------------------- |
*------------------------------------------------------- |
1358 |
|
|
1359 |
include 'commontracker.f' |
include 'commontracker.f' |
1360 |
include 'level1.f' |
include 'level1.f' |
1361 |
include 'calib.f' |
include 'calib.f' |
1362 |
|
|
1363 |
|
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
1364 |
|
incut = incuty |
1365 |
|
pitch = pitchY / 1.e4 |
1366 |
|
else !X-view |
1367 |
|
incut = incutx |
1368 |
|
pitch = pitchX / 1.e4 |
1369 |
|
endif |
1370 |
|
|
1371 |
|
func = 100000. |
1372 |
|
stot = 0. |
1373 |
|
|
1374 |
|
if (ncog.gt.0) then |
1375 |
|
|
1376 |
* --> signal of the central strip |
* --> signal of the central strip |
1377 |
sc = CLSIGNAL(INDMAX(ic)) !center |
sc = CLSIGNAL(INDMAX(ic)) !center |
1378 |
|
fsc = clsigma(INDMAX(ic)) |
1379 |
|
* --> signal of adjacent strips |
1380 |
|
sl1 = 0 !left 1 |
1381 |
|
fsl1 = 1 !left 1 |
1382 |
|
if( |
1383 |
|
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
1384 |
|
$ )then |
1385 |
|
sl1 = CLSIGNAL(INDMAX(ic)-1) |
1386 |
|
fsl1 = clsigma(INDMAX(ic)-1) |
1387 |
|
endif |
1388 |
|
|
1389 |
* signal of adjacent strips |
sl2 = 0 !left 2 |
1390 |
* --> left |
fsl2 = 1 !left 2 |
1391 |
sl1 = 0 !left 1 |
if( |
1392 |
if( |
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
1393 |
$ (INDMAX(ic)-1).ge.INDSTART(ic) |
$ )then |
1394 |
$ ) |
sl2 = CLSIGNAL(INDMAX(ic)-2) |
1395 |
$ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) |
fsl2 = clsigma(INDMAX(ic)-2) |
1396 |
|
endif |
1397 |
sl2 = 0 !left 2 |
sr1 = 0 !right 1 |
1398 |
if( |
fsr1 = 1 !right 1 |
1399 |
$ (INDMAX(ic)-2).ge.INDSTART(ic) |
if( |
1400 |
$ ) |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) |
1401 |
$ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) |
$ .or. |
1402 |
|
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) |
1403 |
* --> right |
$ )then |
1404 |
sr1 = 0 !right 1 |
sr1 = CLSIGNAL(INDMAX(ic)+1) |
1405 |
if( |
fsr1 = clsigma(INDMAX(ic)+1) |
1406 |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) |
endif |
1407 |
$ .or. |
sr2 = 0 !right 2 |
1408 |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) |
fsr2 = 1 !right 2 |
1409 |
$ ) |
if( |
1410 |
$ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) |
1411 |
|
$ .or. |
1412 |
sr2 = 0 !right 2 |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) |
1413 |
if( |
$ )then |
1414 |
$ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) |
sr2 = CLSIGNAL(INDMAX(ic)+2) |
1415 |
$ .or. |
fsr2 = clsigma(INDMAX(ic)+2) |
1416 |
$ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) |
endif |
|
$ ) |
|
|
$ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) |
|
1417 |
|
|
1418 |
|
|
|
if(mod(int(VIEW(ic)),2).eq.1)then !Y-view |
|
|
f = 4. |
|
|
si = 8.4 |
|
|
else !X-view |
|
|
f = 6. |
|
|
si = 3.9 |
|
|
endif |
|
1419 |
|
|
1420 |
fbad_cog = 1. |
************************************************************ |
1421 |
f0 = 1 |
* COG2-3-4 computation |
1422 |
f1 = 1 |
************************************************************ |
1423 |
f2 = 1 |
|
1424 |
f3 = 1 |
c print*,sl2,sl1,sc,sr1,sr2 |
|
if(sl1.gt.sr1.and.sl1.gt.0.)then |
|
1425 |
|
|
1426 |
if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f |
vCOG = cog(ncog,ic)!0. |
1427 |
if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f |
|
1428 |
c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f |
if(ncog.eq.1)then |
1429 |
|
func = 1./12. |
1430 |
if(ncog.eq.2.and.sl1.ne.0)then |
stot = 1. |
1431 |
fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.) |
elseif(ncog.eq.2)then |
1432 |
elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then |
if(sl1.gt.sr1)then |
1433 |
fbad_cog = 1. |
func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2) |
1434 |
elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then |
stot = sl1+sc |
1435 |
fbad_cog = 1. |
elseif(sl1.le.sr1)then |
1436 |
|
func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) |
1437 |
|
stot = sc+sr1 |
1438 |
|
endif |
1439 |
|
elseif(ncog.eq.3)then |
1440 |
|
func = |
1441 |
|
$ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) |
1442 |
|
stot = sl1+sc+sr1 |
1443 |
|
elseif(ncog.eq.4)then |
1444 |
|
if(sl2.gt.sr2)then |
1445 |
|
func = |
1446 |
|
$ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2 |
1447 |
|
$ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2) |
1448 |
|
stot = sl2+sl1+sc+sr1 |
1449 |
|
elseif(sl2.le.sr2)then |
1450 |
|
func = |
1451 |
|
$ (fsl1*(-1-vCOG)**2 |
1452 |
|
$ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2) |
1453 |
|
stot = sl2+sl1+sc+sr1 |
1454 |
|
endif |
1455 |
else |
else |
1456 |
fbad_cog = 1. |
print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG |
1457 |
|
$ ,' not implemented' |
1458 |
endif |
endif |
1459 |
|
|
1460 |
elseif(sl1.le.sr1.and.sr1.gt.0.)then |
elseif(ncog.eq.0)then |
1461 |
|
* ========================= |
1462 |
|
* COG computation |
1463 |
|
* ========================= |
1464 |
|
|
1465 |
|
vCOG = cog(0,ic) |
1466 |
|
|
1467 |
|
iv = VIEW(ic) |
1468 |
|
istart = INDSTART(IC) |
1469 |
|
istop = TOTCLLENGTH |
1470 |
|
if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1 |
1471 |
|
ccc SGN = 0. |
1472 |
|
SNU = 0. |
1473 |
|
ccc SDE = 0. |
1474 |
|
|
1475 |
if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f |
do i=INDMAX(IC),istart,-1 |
1476 |
if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f |
ipos = i-INDMAX(ic) |
1477 |
c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f |
cut = incut*CLSIGMA(i) |
1478 |
|
if(CLSIGNAL(i).gt.cut)then |
1479 |
if(ncog.eq.2.and.sr1.ne.0)then |
fs = clsigma(i) |
1480 |
fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.) |
SNU = SNU + fs*(ipos-vCOG)**2 |
1481 |
elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then |
stot = stot + CLSIGNAL(i) |
1482 |
fbad_cog = 1. |
else |
1483 |
elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then |
goto 10 |
1484 |
fbad_cog = 1. |
endif |
1485 |
|
enddo |
1486 |
|
10 continue |
1487 |
|
do i=INDMAX(IC)+1,istop |
1488 |
|
ipos = i-INDMAX(ic) |
1489 |
|
cut = incut*CLSIGMA(i) |
1490 |
|
if(CLSIGNAL(i).gt.cut)then |
1491 |
|
fs = clsigma(i) |
1492 |
|
SNU = SNU + fs*(ipos-vCOG)**2 |
1493 |
|
stot = stot + CLSIGNAL(i) |
1494 |
|
else |
1495 |
|
goto 20 |
1496 |
|
endif |
1497 |
|
enddo |
1498 |
|
20 continue |
1499 |
|
if(SDE.ne.0)then |
1500 |
|
FUNC=SNU |
1501 |
else |
else |
1502 |
fbad_cog = 1. |
|
1503 |
endif |
endif |
1504 |
|
|
1505 |
|
else |
1506 |
|
|
1507 |
|
FUNC=0 |
1508 |
|
print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG |
1509 |
|
print*,' (NCOG must be >= 0)' |
1510 |
|
|
1511 |
|
|
1512 |
endif |
endif |
1513 |
|
|
1514 |
fbad_cog0 = sqrt(fbad_cog) |
|
1515 |
|
if(stot.gt.0..and.func.gt.0.)then |
1516 |
|
func = sqrt(func) |
1517 |
|
func = pitch * func/stot |
1518 |
|
endif |
1519 |
|
|
1520 |
|
riscogtheor = func |
1521 |
|
|
1522 |
return |
return |
1523 |
end |
end |
1524 |
|
|
1525 |
|
|
1526 |
|
|
1527 |
|
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1528 |
|
|
1529 |
|
real function risetatheor(ncog,ic,angle) |
1530 |
|
*------------------------------------------------------- |
1531 |
|
* |
1532 |
|
* this function returns the expected resolution |
1533 |
|
* obtained by propagating the strip noise |
1534 |
|
* to the coordinate evaluated with non-linear eta-function |
1535 |
|
* |
1536 |
|
* ncog = n.strip used in the coordinate evaluation |
1537 |
|
* (ncog=0 => ncog=2,3,4 according to angle) |
1538 |
|
* |
1539 |
|
*------------------------------------------------------- |
1540 |
|
|
1541 |
|
include 'commontracker.f' |
1542 |
|
include 'level1.f' |
1543 |
|
include 'calib.f' |
1544 |
|
|
1545 |
|
|
1546 |
|
func = 1. |
1547 |
|
|
1548 |
|
iview = VIEW(ic) |
1549 |
|
lad = nld(MAXS(ic),VIEW(ic)) |
1550 |
|
|
1551 |
|
* ------------------------------------------------ |
1552 |
|
* number of strip to be used (in case of ncog = 0) |
1553 |
|
* ------------------------------------------------ |
1554 |
|
|
1555 |
|
inoeta = 0 |
1556 |
|
|
1557 |
|
if(mod(int(iview),2).eq.1)then !Y-view |
1558 |
|
|
1559 |
|
pitch = pitchY / 1.e4 |
1560 |
|
|
1561 |
|
if(ncog.eq.0)then |
1562 |
|
if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then |
1563 |
|
ncog = 2 |
1564 |
|
elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then |
1565 |
|
ncog = 3 |
1566 |
|
elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then |
1567 |
|
ncog = 4 |
1568 |
|
else |
1569 |
|
ncog = 4 |
1570 |
|
inoeta = 1 |
1571 |
|
endif |
1572 |
|
endif |
1573 |
|
|
1574 |
|
else !X-view |
1575 |
|
|
1576 |
|
pitch = pitchX / 1.e4 |
1577 |
|
|
1578 |
|
if(ncog.eq.0)then |
1579 |
|
if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then |
1580 |
|
ncog = 2 |
1581 |
|
elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then |
1582 |
|
ncog = 3 |
1583 |
|
elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then |
1584 |
|
ncog = 4 |
1585 |
|
else |
1586 |
|
ncog = 4 |
1587 |
|
inoeta = 1 |
1588 |
|
endif |
1589 |
|
endif |
1590 |
|
|
1591 |
|
endif |
1592 |
|
|
1593 |
|
func = riscogtheor(ncog,ic) |
1594 |
|
|
1595 |
|
risetatheor = func |
1596 |
|
|
1597 |
|
if(inoeta.eq.1)return ! no eta correction is applied --> exit |
1598 |
|
if(ncog.lt.1.or.ncog.gt.4)return |
1599 |
|
|
1600 |
|
* ---------------- |
1601 |
|
* find angular bin |
1602 |
|
* ---------------- |
1603 |
|
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
1604 |
|
do iang=1,nangbin |
1605 |
|
if(angL(iang).lt.angle.and.angR(iang).ge.angle)then |
1606 |
|
iangle=iang |
1607 |
|
goto 98 |
1608 |
|
endif |
1609 |
|
enddo |
1610 |
|
if(DEBUG.EQ.1)print* |
1611 |
|
$ ,'risetatheor *** warning *** angle out of range: ',angle |
1612 |
|
if(angle.le.angL(1))iang=1 |
1613 |
|
if(angle.ge.angR(nangbin))iang=nangbin |
1614 |
|
98 continue !jump here if ok |
1615 |
|
|
1616 |
|
* ------------- |
1617 |
|
* within +/-0.5 |
1618 |
|
* ------------- |
1619 |
|
|
1620 |
|
vcog = cog(ncog,ic) |
1621 |
|
|
1622 |
|
etamin = eta2(1,iang) |
1623 |
|
etamax = eta2(netaval,iang) |
1624 |
|
|
1625 |
|
iaddmax=10 |
1626 |
|
iadd=0 |
1627 |
|
10 continue |
1628 |
|
if(vcog.lt.etamin)then |
1629 |
|
vcog = vcog + 1 |
1630 |
|
iadd = iadd + 1 |
1631 |
|
if(iadd>iaddmax)goto 111 |
1632 |
|
goto 10 |
1633 |
|
endif |
1634 |
|
20 continue |
1635 |
|
if(vcog.gt.etamax)then |
1636 |
|
vcog = vcog - 1 |
1637 |
|
iadd = iadd - 1 |
1638 |
|
if(iadd<-1*iaddmax)goto 111 |
1639 |
|
goto 20 |
1640 |
|
endif |
1641 |
|
goto 1111 |
1642 |
|
111 continue |
1643 |
|
if(DEBUG.eq.1) |
1644 |
|
$ print*,'risetatheor *** warning *** anomalous cluster' |
1645 |
|
if(DEBUG.eq.1) |
1646 |
|
$ print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)' |
1647 |
|
vcog=0 |
1648 |
|
1111 continue |
1649 |
|
|
1650 |
|
* ------------------------------------------------ |
1651 |
|
* interpolation |
1652 |
|
* ------------------------------------------------ |
1653 |
|
|
1654 |
|
|
1655 |
|
ibin = netaval |
1656 |
|
do i=2,netaval |
1657 |
|
if(ncog.eq.2)eta=eta2(i,iang) |
1658 |
|
if(ncog.eq.3)eta=eta3(i,iang) |
1659 |
|
if(ncog.eq.4)eta=eta4(i,iang) |
1660 |
|
if(eta.ge.vcog)then |
1661 |
|
ibin = i |
1662 |
|
goto 99 |
1663 |
|
endif |
1664 |
|
enddo |
1665 |
|
99 continue |
1666 |
|
|
1667 |
|
if(ncog.eq.2)then |
1668 |
|
x1 = eta2(ibin-1,iang) |
1669 |
|
x2 = eta2(ibin,iang) |
1670 |
|
y1 = feta2(ibin-1,iview,lad,iang) |
1671 |
|
y2 = feta2(ibin,iview,lad,iang) |
1672 |
|
elseif(ncog.eq.3)then |
1673 |
|
x1 = eta3(ibin-1,iang) |
1674 |
|
x2 = eta3(ibin,iang) |
1675 |
|
y1 = feta3(ibin-1,iview,lad,iang) |
1676 |
|
y2 = feta3(ibin,iview,lad,iang) |
1677 |
|
elseif(ncog.eq.4)then |
1678 |
|
x1 = eta4(ibin-1,iang) |
1679 |
|
x2 = eta4(ibin,iang) |
1680 |
|
y1 = feta4(ibin-1,iview,lad,iang) |
1681 |
|
y2 = feta4(ibin,iview,lad,iang) |
1682 |
|
endif |
1683 |
|
|
1684 |
|
func = func * (y2-y1)/(x2-x1) |
1685 |
|
|
1686 |
|
risetatheor = func |
1687 |
|
|
1688 |
|
return |
1689 |
|
end |
1690 |
|
|
1691 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1692 |
|
|
1693 |
FUNCTION risx_eta2(x) |
FUNCTION risxeta2(x) |
1694 |
|
|
1695 |
DOUBLE PRECISION V( 1) |
DOUBLE PRECISION V( 1) |
1696 |
INTEGER NPAR, NDIM, IMQFUN, I, J |
INTEGER NPAR, NDIM, IMQFUN, I, J |
1762 |
+/ |
+/ |
1763 |
|
|
1764 |
V(1)= abs(x) |
V(1)= abs(x) |
1765 |
|
if(V(1).gt.20.)V(1)=20. |
1766 |
|
|
1767 |
HQUADF = 0. |
HQUADF = 0. |
1768 |
DO 20 J = 1, NPAR |
DO 20 J = 1, NPAR |
1777 |
20 CONTINUE |
20 CONTINUE |
1778 |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
1779 |
|
|
1780 |
risx_eta2=HQUADF* 1e-4 |
risxeta2=HQUADF* 1e-4 |
1781 |
|
|
1782 |
END |
END |
1783 |
|
|
1784 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1785 |
FUNCTION risx_eta3(x) |
FUNCTION risxeta3(x) |
1786 |
DOUBLE PRECISION V( 1) |
DOUBLE PRECISION V( 1) |
1787 |
INTEGER NPAR, NDIM, IMQFUN, I, J |
INTEGER NPAR, NDIM, IMQFUN, I, J |
1788 |
DOUBLE PRECISION HQDJ, VV, VCONST |
DOUBLE PRECISION HQDJ, VV, VCONST |
1853 |
+/ |
+/ |
1854 |
|
|
1855 |
V(1) = abs(x) |
V(1) = abs(x) |
1856 |
|
if(V(1).gt.20.)V(1)=20. |
1857 |
|
|
1858 |
HQUADF = 0. |
HQUADF = 0. |
1859 |
DO 20 J = 1, NPAR |
DO 20 J = 1, NPAR |
1860 |
HQDJ = 0. |
HQDJ = 0. |
1868 |
20 CONTINUE |
20 CONTINUE |
1869 |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
1870 |
|
|
1871 |
risx_eta3 = HQUADF* 1e-4 |
risxeta3 = HQUADF* 1e-4 |
1872 |
|
|
1873 |
END |
END |
1874 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1875 |
FUNCTION risx_eta4(x) |
FUNCTION risxeta4(x) |
1876 |
DOUBLE PRECISION V( 1) |
DOUBLE PRECISION V( 1) |
1877 |
INTEGER NPAR, NDIM, IMQFUN, I, J |
INTEGER NPAR, NDIM, IMQFUN, I, J |
1878 |
DOUBLE PRECISION HQDJ, VV, VCONST |
DOUBLE PRECISION HQDJ, VV, VCONST |
1943 |
+/ |
+/ |
1944 |
|
|
1945 |
V(1)=abs(x) |
V(1)=abs(x) |
1946 |
|
if(V(1).gt.20.)V(1)=20. |
1947 |
|
|
1948 |
HQUADF = 0. |
HQUADF = 0. |
1949 |
DO 20 J = 1, NPAR |
DO 20 J = 1, NPAR |
1958 |
20 CONTINUE |
20 CONTINUE |
1959 |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
1960 |
|
|
1961 |
risx_eta4=HQUADF* 1e-4 |
risxeta4=HQUADF* 1e-4 |
1962 |
|
|
1963 |
END |
END |
1964 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
1965 |
FUNCTION risy_eta2(x) |
FUNCTION risyeta2(x) |
1966 |
DOUBLE PRECISION V( 1) |
DOUBLE PRECISION V( 1) |
1967 |
INTEGER NPAR, NDIM, IMQFUN, I, J |
INTEGER NPAR, NDIM, IMQFUN, I, J |
1968 |
DOUBLE PRECISION HQDJ, VV, VCONST |
DOUBLE PRECISION HQDJ, VV, VCONST |
2015 |
+/ |
+/ |
2016 |
|
|
2017 |
v(1)= abs(x) |
v(1)= abs(x) |
2018 |
|
if(V(1).gt.20.)V(1)=20. |
2019 |
|
|
2020 |
HQUADF = 0. |
HQUADF = 0. |
2021 |
DO 20 J = 1, NPAR |
DO 20 J = 1, NPAR |
2030 |
20 CONTINUE |
20 CONTINUE |
2031 |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF) |
2032 |
|
|
2033 |
risy_eta2=HQUADF* 1e-4 |
risyeta2=HQUADF* 1e-4 |
2034 |
|
|
2035 |
END |
END |
2036 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
2082 |
+/ |
+/ |
2083 |
|
|
2084 |
V(1)=abs(x) |
V(1)=abs(x) |
2085 |
|
if(V(1).gt.20.)V(1)=20. |
2086 |
|
|
2087 |
HQUADF = 0. |
HQUADF = 0. |
2088 |
DO 20 J = 1, NPAR |
DO 20 J = 1, NPAR |
2163 |
+/ |
+/ |
2164 |
|
|
2165 |
V(1)=abs(x) |
V(1)=abs(x) |
2166 |
|
if(V(1).gt.20.)V(1)=20. |
2167 |
|
|
2168 |
HQUADF = 0. |
HQUADF = 0. |
2169 |
DO 20 J = 1, NPAR |
DO 20 J = 1, NPAR |
2181 |
risx_cog = HQUADF * 1e-4 |
risx_cog = HQUADF * 1e-4 |
2182 |
|
|
2183 |
END |
END |
2184 |
|
|
2185 |
|
|
2186 |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** |
2187 |
|
real function pfacorr(ic,angle) |
2188 |
|
*-------------------------------------------------------------- |
2189 |
|
* this function returns the landi correction for this cluster |
2190 |
|
*-------------------------------------------------------------- |
2191 |
|
include 'commontracker.f' |
2192 |
|
include 'calib.f' |
2193 |
|
include 'level1.f' |
2194 |
|
|
2195 |
|
real angle |
2196 |
|
integer iview,lad |
2197 |
|
|
2198 |
|
iview = VIEW(ic) |
2199 |
|
lad = nld(MAXS(ic),VIEW(ic)) |
2200 |
|
|
2201 |
|
* find angular bin |
2202 |
|
* (in futuro possiamo pensare di interpolare anche sull'angolo) |
2203 |
|
do iang=1,nangbin |
2204 |
|
if(angL(iang).lt.angle.and.angR(iang).ge.angle)then |
2205 |
|
iangle=iang |
2206 |
|
goto 98 |
2207 |
|
endif |
2208 |
|
enddo |
2209 |
|
if(DEBUG.eq.1) |
2210 |
|
$ print*,'pfacorr *** warning *** angle out of range: ',angle |
2211 |
|
if(angle.le.angL(1))iang=1 |
2212 |
|
if(angle.ge.angR(nangbin))iang=nangbin |
2213 |
|
98 continue !jump here if ok |
2214 |
|
|
2215 |
|
pfacorr = fcorr(iview,lad,iang) |
2216 |
|
|
2217 |
|
if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr |
2218 |
|
|
2219 |
|
|
2220 |
|
c 100 return |
2221 |
|
return |
2222 |
|
end |