/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/functionspfa.f
ViewVC logotype

Annotation of /DarthVader/TrackerLevel2/src/F77/functionspfa.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.19 - (hide annotations) (download)
Tue Aug 21 15:21:22 2007 UTC (17 years, 3 months ago) by pam-fi
Branch: MAIN
Changes since 1.18: +13 -10 lines
*** empty log message ***

1 pam-fi 1.9
2 pam-fi 1.10
3     subroutine idtoc(ipfa,cpfa)
4    
5     integer ipfa
6 pam-fi 1.18 character*10 cpfa
7 pam-fi 1.10
8     CPFA='COG4'
9     if(ipfa.eq.0)CPFA='ETA'
10     if(ipfa.eq.2)CPFA='ETA2'
11     if(ipfa.eq.3)CPFA='ETA3'
12     if(ipfa.eq.4)CPFA='ETA4'
13 pam-fi 1.16 if(ipfa.eq.5)CPFA='ETAL'
14 pam-fi 1.10 if(ipfa.eq.10)CPFA='COG'
15     if(ipfa.eq.11)CPFA='COG1'
16     if(ipfa.eq.12)CPFA='COG2'
17     if(ipfa.eq.13)CPFA='COG3'
18     if(ipfa.eq.14)CPFA='COG4'
19    
20     end
21    
22 mocchiut 1.1 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
23     * this file contains all subroutines and functions
24     * that are needed for position finding algorithms
25     *
26     *
27     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
28    
29    
30 pam-fi 1.9 integer function npfastrips(ic,PFA,angle)
31     *--------------------------------------------------------------
32     * thid function returns the number of strips used
33     * to evaluate the position of a cluster, according to the p.f.a.
34     *--------------------------------------------------------------
35     include 'commontracker.f'
36     include 'level1.f'
37     include 'calib.f'
38    
39     character*4 usedPFA,PFA
40    
41    
42     usedPFA=PFA
43    
44     npfastrips=0
45    
46     if(usedPFA.eq.'COG1')npfastrips=1
47     if(usedPFA.eq.'COG2')npfastrips=2
48     if(usedPFA.eq.'COG3')npfastrips=3
49     if(usedPFA.eq.'COG4')npfastrips=4
50     if(usedPFA.eq.'ETA2')npfastrips=2
51     if(usedPFA.eq.'ETA3')npfastrips=3
52     if(usedPFA.eq.'ETA4')npfastrips=4
53     * ----------------------------------------------------------------
54     if(usedPFA.eq.'ETA')then
55     c print*,VIEW(ic),angle
56     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
57     if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
58     npfastrips=2
59     elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
60     npfastrips=3
61     elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
62     npfastrips=4
63     else
64     npfastrips=4
65     c usedPFA='COG'
66     endif
67     else !X-view
68     if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
69     npfastrips=2
70     elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
71     npfastrips=3
72     elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
73     npfastrips=4
74     else
75     npfastrips=4
76     c usedPFA='COG'
77     endif
78     endif
79     endif
80     * ----------------------------------------------------------------
81     if(usedPFA.eq.'COG')then
82    
83     iv=VIEW(ic)
84     if(mod(iv,2).eq.1)incut=incuty
85     if(mod(iv,2).eq.0)incut=incutx
86     istart = INDSTART(IC)
87     istop = TOTCLLENGTH
88     if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
89     mu = 0
90     do i = INDMAX(IC),istart,-1
91     ipos = i-INDMAX(ic)
92     cut = incut*CLSIGMA(i)
93     if(CLSIGNAL(i).ge.cut)then
94     mu = mu + 1
95     print*,i,mu
96     else
97     goto 10
98     endif
99     enddo
100     10 continue
101     do i = INDMAX(IC)+1,istop
102     ipos = i-INDMAX(ic)
103     cut = incut*CLSIGMA(i)
104     if(CLSIGNAL(i).ge.cut)then
105     mu = mu + 1
106     print*,i,mu
107     else
108     goto 20
109     endif
110     enddo
111     20 continue
112     npfastrips=mu
113    
114     endif
115     * ----------------------------------------------------------------
116    
117     c print*,pfastrips
118    
119     return
120     end
121    
122 mocchiut 1.1 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
123 pam-fi 1.4 real function pfaeta(ic,angle)
124 mocchiut 1.1 *--------------------------------------------------------------
125     * this function returns the position (in strip units)
126     * it calls:
127 pam-fi 1.4 * - pfaeta2(ic,angle)
128     * - pfaeta3(ic,angle)
129     * - pfaeta4(ic,angle)
130 mocchiut 1.1 * according to the angle
131     *--------------------------------------------------------------
132     include 'commontracker.f'
133     include 'level1.f'
134 pam-fi 1.9 include 'calib.f'
135 mocchiut 1.1
136 pam-fi 1.4 pfaeta = 0
137 mocchiut 1.1
138     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
139    
140 pam-fi 1.9 if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
141     pfaeta = pfaeta2(ic,angle)
142 pam-fi 1.19 cc print*,pfaeta2(ic,angle)
143 pam-fi 1.9 elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
144     pfaeta = pfaeta3(ic,angle)
145     elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
146     pfaeta = pfaeta4(ic,angle)
147     else
148     pfaeta = cog(4,ic)
149     endif
150    
151 mocchiut 1.1 else !X-view
152    
153 pam-fi 1.9 if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
154 pam-fi 1.4 pfaeta = pfaeta2(ic,angle)
155 pam-fi 1.9 elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
156 pam-fi 1.4 pfaeta = pfaeta3(ic,angle)
157 pam-fi 1.9 elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
158 pam-fi 1.4 pfaeta = pfaeta4(ic,angle)
159 pam-fi 1.9 else
160     pfaeta = cog(4,ic)
161     endif
162 mocchiut 1.1
163     endif
164 pam-fi 1.5
165 mocchiut 1.1 100 return
166     end
167    
168     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
169 pam-fi 1.16 real function pfaetal(ic,angle)
170     *--------------------------------------------------------------
171     * this function returns the position (in strip units)
172     * it calls:
173     * - pfaeta2(ic,angle)+pfcorr(ic,angle)
174     * - pfaeta3(ic,angle)+pfcorr(ic,angle)
175     * - pfaeta4(ic,angle)+pfcorr(ic,angle)
176     * according to the angle
177     *--------------------------------------------------------------
178     include 'commontracker.f'
179     include 'level1.f'
180     include 'calib.f'
181    
182 pam-fi 1.19 pfaetal = 0
183 pam-fi 1.16
184     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
185    
186     if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
187 pam-fi 1.19 pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
188     cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
189 pam-fi 1.16 elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
190 pam-fi 1.19 pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
191 pam-fi 1.16 elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
192 pam-fi 1.19 pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
193 pam-fi 1.16 else
194 pam-fi 1.19 pfaetal = cog(4,ic)
195 pam-fi 1.16 endif
196    
197     else !X-view
198    
199     if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
200 pam-fi 1.19 pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
201     cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
202 pam-fi 1.16 elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
203 pam-fi 1.19 pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
204 pam-fi 1.16 elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
205 pam-fi 1.19 pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
206 pam-fi 1.16 else
207 pam-fi 1.19 pfaetal = cog(4,ic)
208 pam-fi 1.16 endif
209    
210     endif
211    
212     100 return
213     end
214     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
215 pam-fi 1.11 c real function riseta(ic,angle)
216     real function riseta(iview,angle)
217 mocchiut 1.1 *--------------------------------------------------------------
218     * this function returns the average spatial resolution
219 pam-fi 1.4 * (in cm) for the ETA algorithm (function pfaeta(ic,angle))
220 mocchiut 1.1 * it calls:
221 pam-fi 1.13 * - risxeta2(angle)
222     * - risyeta2(angle)
223     * - risxeta3(angle)
224     * - risxeta4(angle)
225 mocchiut 1.1 * according to the angle
226     *--------------------------------------------------------------
227     include 'commontracker.f'
228     include 'level1.f'
229 pam-fi 1.9 include 'calib.f'
230 mocchiut 1.1
231 pam-fi 1.11 riseta = 0
232 mocchiut 1.1
233 pam-fi 1.11 c if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
234     if(mod(iview,2).eq.1)then !Y-view
235 mocchiut 1.1
236 pam-fi 1.9
237     if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
238 pam-fi 1.13 riseta = risyeta2(angle)
239 pam-fi 1.9 elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
240 pam-fi 1.11 riseta = risy_cog(angle) !ATTENZIONE!!
241 pam-fi 1.9 elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
242 pam-fi 1.11 riseta = risy_cog(angle) !ATTENZIONE!!
243 pam-fi 1.9 else
244 pam-fi 1.11 riseta = risy_cog(angle)
245 pam-fi 1.9 endif
246 mocchiut 1.1
247     else !X-view
248    
249 pam-fi 1.9 if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
250 pam-fi 1.13 riseta = risxeta2(angle)
251 pam-fi 1.9 elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
252 pam-fi 1.13 riseta = risxeta3(angle)
253 pam-fi 1.9 elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
254 pam-fi 1.13 riseta = risxeta4(angle)
255 pam-fi 1.9 else
256 pam-fi 1.11 riseta = risx_cog(angle)
257 pam-fi 1.9 endif
258 mocchiut 1.1
259     endif
260    
261 pam-fi 1.11
262 mocchiut 1.1 100 return
263     end
264    
265     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
266     real function fbad_eta(ic,angle)
267     *-------------------------------------------------------
268     * this function returns a factor that takes into
269     * account deterioration of the spatial resolution
270     * in the case BAD strips are included in the cluster.
271     * This factor should multiply the nominal spatial
272     * resolution.
273     * It calls the function FBAD_COG(NCOG,IC),
274     * accordingto the angle
275     *-------------------------------------------------------
276    
277     include 'commontracker.f'
278     include 'level1.f'
279 pam-fi 1.9 include 'calib.f'
280 mocchiut 1.1 fbad_eta = 0
281    
282     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
283    
284 pam-fi 1.9 if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
285     fbad_eta = fbad_cog(2,ic)
286     elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
287     fbad_eta = fbad_cog(3,ic)
288     elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
289     fbad_eta = fbad_cog(4,ic)
290     else
291     fbad_eta = fbad_cog(4,ic)
292     endif
293 mocchiut 1.1
294     else !X-view
295    
296 pam-fi 1.9 if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
297 mocchiut 1.1 fbad_eta = fbad_cog(2,ic)
298 pam-fi 1.9 elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
299 mocchiut 1.1 fbad_eta = fbad_cog(3,ic)
300 pam-fi 1.9 elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
301     fbad_eta = fbad_cog(4,ic)
302     else
303 mocchiut 1.1 fbad_eta = fbad_cog(4,ic)
304 pam-fi 1.9 endif
305 mocchiut 1.1
306     endif
307    
308     return
309     end
310    
311     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
312 pam-fi 1.4 real function pfaeta2(ic,angle) !(1)
313 mocchiut 1.1 *--------------------------------------------------------------
314     * this function returns
315     *
316     * - the position (in strip units)
317     * corrected according to the ETA2 Position Finding Algorithm.
318     * The function performs an interpolation of FETA2%ETA2
319     *
320     * - if the angle is out of range, the calibration parameters
321     * of the lowest or higher bin are used
322     *
323     *--------------------------------------------------------------
324     include 'commontracker.f'
325     include 'calib.f'
326     include 'level1.f'
327    
328     real cog2,angle
329     integer iview,lad
330    
331 pam-fi 1.9 iview = VIEW(ic)
332     lad = nld(MAXS(ic),VIEW(ic))
333     cog2 = cog(2,ic)
334 pam-fi 1.4 pfaeta2=cog2
335 mocchiut 1.1
336     * find angular bin
337     * (in futuro possiamo pensare di interpolare anche sull'angolo)
338     do iang=1,nangbin
339     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
340     iangle=iang
341     goto 98
342     endif
343     enddo
344 pam-fi 1.18 if(DEBUG.EQ.1)
345 pam-fi 1.4 $ print*,'pfaeta2 *** warning *** angle out of range: ',angle
346 mocchiut 1.1 if(angle.lt.angL(1))iang=1
347     if(angle.gt.angR(nangbin))iang=nangbin
348     98 continue !jump here if ok
349    
350    
351     c$$$* find extremes of interpolation
352     c$$$ iflag=0
353     c$$$* --------------------------------
354     c$$$ if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
355 pam-fi 1.4 c$$$c print*,'pfaeta2 *** warning *** argument out of range: ',cog2
356 mocchiut 1.1 c$$$* goto 100
357     c$$$* ----------------------------------------------
358     c$$$* non salto piu`, ma scalo di 1 o -1
359     c$$$* nel caso si tratti di un cluster
360     c$$$* in cui la strip con il segnale massimo non coincide
361     c$$$* con la strip con il rapposto s/n massimo!!!
362     c$$$* ----------------------------------------------
363     c$$$ if(cog2.lt.eta2(1,iang))then !temp
364     c$$$ cog2=cog2+1. !temp
365     c$$$ iflag=1 !temp
366     c$$$ else !temp
367     c$$$ cog2=cog2-1. !temp
368     c$$$ iflag=-1 !temp
369     c$$$ endif !temp
370     c$$$c print*,'shifted >>> ',cog2
371     c$$$ endif
372    
373     iadd=0
374     10 continue
375     if(cog2.lt.eta2(1,iang))then
376     cog2 = cog2 + 1
377     iadd = iadd + 1
378     goto 10
379     endif
380     20 continue
381     if(cog2.gt.eta2(netaval,iang))then
382     cog2 = cog2 - 1
383     iadd = iadd - 1
384     goto 20
385     endif
386    
387     * --------------------------------
388     c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
389     do i=2,netaval
390     if(eta2(i,iang).gt.cog2)then
391    
392     x1 = eta2(i-1,iang)
393     x2 = eta2(i,iang)
394     y1 = feta2(i-1,iview,lad,iang)
395     y2 = feta2(i,iview,lad,iang)
396    
397     c print*,'*****',i,view,lad,iang
398     c print*,'-----',x1,x2,y1,y2
399     goto 99
400     endif
401     enddo
402     99 continue
403    
404    
405     AA=(y2-y1)/(x2-x1)
406     BB=y1-AA*x1
407    
408 pam-fi 1.4 pfaeta2 = AA*cog2+BB
409     pfaeta2 = pfaeta2 - iadd
410 mocchiut 1.1
411     c$$$ if(iflag.eq.1)then
412 pam-fi 1.4 c$$$ pfaeta2=pfaeta2-1. !temp
413 mocchiut 1.1 c$$$ cog2=cog2-1. !temp
414     c$$$ endif
415     c$$$ if(iflag.eq.-1)then
416 pam-fi 1.4 c$$$ pfaeta2=pfaeta2+1. !temp
417 mocchiut 1.1 c$$$ cog2=cog2+1. !temp
418     c$$$ endif
419    
420 pam-fi 1.18 if(DEBUG.EQ.1)print*,'ETA2 (ic ',ic,' ang',angle,')'
421 pam-fi 1.4 $ ,cog2-iadd,' -->',pfaeta2
422 mocchiut 1.1
423    
424     100 return
425     end
426    
427     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
428 pam-fi 1.4 real function pfaeta3(ic,angle) !(1)
429 mocchiut 1.1 *--------------------------------------------------------------
430     * this function returns
431     *
432     * - the position (in strip units)
433     * corrected according to the ETA3 Position Finding Algorithm.
434     * The function performs an interpolation of FETA3%ETA3
435     *
436     * - if the angle is out of range, the calibration parameters
437     * of the lowest or higher bin are used
438     *
439     *--------------------------------------------------------------
440     include 'commontracker.f'
441     include 'calib.f'
442     include 'level1.f'
443    
444     real cog3,angle
445     integer iview,lad
446    
447    
448 pam-fi 1.9 iview = VIEW(ic)
449     lad = nld(MAXS(ic),VIEW(ic))
450     cog3 = cog(3,ic)
451 pam-fi 1.4 pfaeta3=cog3
452 mocchiut 1.1
453     * find angular bin
454     * (in futuro possiamo pensare di interpolare anche sull'angolo)
455     do iang=1,nangbin
456     c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
457     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
458     iangle=iang
459     goto 98
460     endif
461     enddo
462 pam-fi 1.18 if(DEBUG.EQ.1)
463 pam-fi 1.4 $ print*,'pfaeta3 *** warning *** angle out of range: ',angle
464 mocchiut 1.1 if(angle.lt.angL(1))iang=1
465     if(angle.gt.angR(nangbin))iang=nangbin
466     98 continue !jump here if ok
467    
468    
469     c$$$* find extremes of interpolation
470     c$$$ iflag=0
471     c$$$* --------------------------------
472     c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then
473     c$$$* ----------------------------------------------
474     c$$$* non salto piu`, ma scalo di 1 o -1
475     c$$$* nel caso si tratti di un cluster
476     c$$$* in cui la strip con il segnale massimo non coincide
477     c$$$* con la strip con il rapposto s/n massimo!!!
478     c$$$* ----------------------------------------------
479     c$$$ if(cog2.lt.eta2(1,iang))then !temp
480     c$$$ cog2=cog2+1. !temp
481     c$$$ iflag=1 !temp
482     c$$$ else !temp
483     c$$$ cog2=cog2-1. !temp
484     c$$$ iflag=-1 !temp
485     c$$$ endif !temp
486     c$$$c print*,'shifted >>> ',cog2
487     c$$$ endif
488    
489    
490     iadd=0
491     10 continue
492     if(cog3.lt.eta3(1,iang))then
493     cog3 = cog3 + 1
494     iadd = iadd + 1
495     goto 10
496     endif
497     20 continue
498     if(cog3.gt.eta3(netaval,iang))then
499     cog3 = cog3 - 1
500     iadd = iadd - 1
501     goto 20
502     endif
503    
504     * --------------------------------
505     c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
506     do i=2,netaval
507     if(eta3(i,iang).gt.cog3)then
508    
509     x1 = eta3(i-1,iang)
510     x2 = eta3(i,iang)
511     y1 = feta3(i-1,iview,lad,iang)
512     y2 = feta3(i,iview,lad,iang)
513    
514     c print*,'*****',i,view,lad,iang
515     c print*,'-----',x1,x2,y1,y2
516     goto 99
517     endif
518     enddo
519     99 continue
520    
521    
522     AA=(y2-y1)/(x2-x1)
523     BB=y1-AA*x1
524    
525 pam-fi 1.4 pfaeta3 = AA*cog3+BB
526     pfaeta3 = pfaeta3 - iadd
527 mocchiut 1.1
528     c$$$ if(iflag.eq.1)then
529 pam-fi 1.4 c$$$ pfaeta2=pfaeta2-1. !temp
530 mocchiut 1.1 c$$$ cog2=cog2-1. !temp
531     c$$$ endif
532     c$$$ if(iflag.eq.-1)then
533 pam-fi 1.4 c$$$ pfaeta2=pfaeta2+1. !temp
534 mocchiut 1.1 c$$$ cog2=cog2+1. !temp
535     c$$$ endif
536    
537 pam-fi 1.18 if(DEBUG.EQ.1)print*,'ETA3 (ic ',ic,' ang',angle,')'
538 pam-fi 1.4 $ ,cog3-iadd,' -->',pfaeta3
539 mocchiut 1.1
540     100 return
541     end
542    
543     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
544 pam-fi 1.9 real function pfaeta4(ic,angle)
545 mocchiut 1.1 *--------------------------------------------------------------
546     * this function returns
547     *
548     * - the position (in strip units)
549     * corrected according to the ETA4 Position Finding Algorithm.
550     * The function performs an interpolation of FETA3%ETA3
551     *
552     * - if the angle is out of range, the calibration parameters
553     * of the lowest or higher bin are used
554     *
555     *--------------------------------------------------------------
556     include 'commontracker.f'
557     include 'calib.f'
558     include 'level1.f'
559    
560     real cog4,angle
561     integer iview,lad
562    
563 pam-fi 1.6
564 pam-fi 1.9 iview = VIEW(ic)
565     lad = nld(MAXS(ic),VIEW(ic))
566     cog4=cog(4,ic)
567 pam-fi 1.4 pfaeta4=cog4
568 mocchiut 1.1
569     * find angular bin
570     * (in futuro possiamo pensare di interpolare anche sull'angolo)
571     do iang=1,nangbin
572     c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
573     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
574     iangle=iang
575     goto 98
576     endif
577     enddo
578 pam-fi 1.18 if(DEBUG.EQ.1)
579 pam-fi 1.4 $ print*,'pfaeta4 *** warning *** angle out of range: ',angle
580 mocchiut 1.1 if(angle.lt.angL(1))iang=1
581     if(angle.gt.angR(nangbin))iang=nangbin
582     98 continue !jump here if ok
583    
584    
585     c$$$* find extremes of interpolation
586     c$$$ iflag=0
587     c$$$* --------------------------------
588     c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then
589     c$$$* ----------------------------------------------
590     c$$$* non salto piu`, ma scalo di 1 o -1
591     c$$$* nel caso si tratti di un cluster
592     c$$$* in cui la strip con il segnale massimo non coincide
593     c$$$* con la strip con il rapposto s/n massimo!!!
594     c$$$* ----------------------------------------------
595     c$$$ if(cog2.lt.eta2(1,iang))then !temp
596     c$$$ cog2=cog2+1. !temp
597     c$$$ iflag=1 !temp
598     c$$$ else !temp
599     c$$$ cog2=cog2-1. !temp
600     c$$$ iflag=-1 !temp
601     c$$$ endif !temp
602     c$$$c print*,'shifted >>> ',cog2
603     c$$$ endif
604    
605    
606     iadd=0
607     10 continue
608     if(cog4.lt.eta4(1,iang))then
609     cog4 = cog4 + 1
610     iadd = iadd + 1
611     goto 10
612     endif
613     20 continue
614     if(cog4.gt.eta4(netaval,iang))then
615     cog4 = cog4 - 1
616     iadd = iadd - 1
617     goto 20
618     endif
619    
620     * --------------------------------
621     c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
622     do i=2,netaval
623     if(eta4(i,iang).gt.cog4)then
624    
625     x1 = eta4(i-1,iang)
626     x2 = eta4(i,iang)
627     y1 = feta4(i-1,iview,lad,iang)
628     y2 = feta4(i,iview,lad,iang)
629    
630     c print*,'*****',i,view,lad,iang
631     c print*,'-----',x1,x2,y1,y2
632     goto 99
633     endif
634     enddo
635     99 continue
636    
637    
638     AA=(y2-y1)/(x2-x1)
639     BB=y1-AA*x1
640    
641 pam-fi 1.4 pfaeta4 = AA*cog4+BB
642     pfaeta4 = pfaeta4 - iadd
643 mocchiut 1.1
644     c$$$ if(iflag.eq.1)then
645 pam-fi 1.4 c$$$ pfaeta2=pfaeta2-1. !temp
646 mocchiut 1.1 c$$$ cog2=cog2-1. !temp
647     c$$$ endif
648     c$$$ if(iflag.eq.-1)then
649 pam-fi 1.4 c$$$ pfaeta2=pfaeta2+1. !temp
650 mocchiut 1.1 c$$$ cog2=cog2+1. !temp
651     c$$$ endif
652    
653 pam-fi 1.18 if(DEBUG.EQ.1)print*,'ETA4 (ic ',ic,' ang',angle,')'
654 pam-fi 1.4 $ ,cog4-iadd,' -->',pfaeta4
655 mocchiut 1.1
656     100 return
657     end
658    
659    
660    
661     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
662 pam-fi 1.14 c$$$ real function cog0(ncog,ic)
663     c$$$*-------------------------------------------------
664     c$$$* this function returns
665     c$$$*
666     c$$$* - the Center-Of-Gravity of the cluster IC
667     c$$$* evaluated using NCOG strips,
668     c$$$* calculated relative to MAXS(IC)
669     c$$$*
670     c$$$* - zero in case that not enough strips
671     c$$$* have a positive signal
672     c$$$*
673     c$$$* NOTE:
674     c$$$* This is the old definition, used by Straulino.
675     c$$$* The new routine, according to Landi,
676     c$$$* is COG(NCOG,IC)
677     c$$$*-------------------------------------------------
678     c$$$
679     c$$$
680     c$$$ include 'commontracker.f'
681     c$$$ include 'level1.f'
682     c$$$
683     c$$$* --> signal of the central strip
684     c$$$ sc = CLSIGNAL(INDMAX(ic)) !center
685     c$$$
686     c$$$* signal of adjacent strips
687     c$$$* --> left
688     c$$$ sl1 = 0 !left 1
689     c$$$ if(
690     c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic)
691     c$$$ $ )
692     c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
693     c$$$
694     c$$$ sl2 = 0 !left 2
695     c$$$ if(
696     c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic)
697     c$$$ $ )
698     c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
699     c$$$
700     c$$$* --> right
701     c$$$ sr1 = 0 !right 1
702     c$$$ if(
703     c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
704     c$$$ $ .or.
705     c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
706     c$$$ $ )
707     c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
708     c$$$
709     c$$$ sr2 = 0 !right 2
710     c$$$ if(
711     c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
712     c$$$ $ .or.
713     c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
714     c$$$ $ )
715     c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
716     c$$$
717     c$$$************************************************************
718     c$$$* COG computation
719     c$$$************************************************************
720     c$$$
721     c$$$c print*,sl2,sl1,sc,sr1,sr2
722     c$$$
723     c$$$ COG = 0.
724     c$$$
725     c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then
726     c$$$
727     c$$$ if(ncog.eq.2.and.sl1.ne.0)then
728     c$$$ COG = -sl1/(sl1+sc)
729     c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
730     c$$$ COG = (sr1-sl1)/(sl1+sc+sr1)
731     c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
732     c$$$ COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
733     c$$$ else
734     c$$$ COG = 0.
735     c$$$ endif
736     c$$$
737     c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then
738     c$$$
739     c$$$ if(ncog.eq.2.and.sr1.ne.0)then
740     c$$$ COG = sr1/(sc+sr1)
741     c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
742     c$$$ COG = (sr1-sl1)/(sl1+sc+sr1)
743     c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
744     c$$$ COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
745     c$$$ else
746     c$$$ COG = 0.
747     c$$$ endif
748     c$$$
749     c$$$ endif
750     c$$$
751     c$$$ COG0 = COG
752     c$$$
753     c$$$c print *,ncog,ic,cog,'/////////////'
754     c$$$
755     c$$$ return
756     c$$$ end
757 mocchiut 1.1
758     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
759     real function cog(ncog,ic)
760     *-------------------------------------------------
761     * this function returns
762     *
763     * - if NCOG=0, the Center-Of-Gravity of the
764     * cluster IC, relative to MAXS(IC), according to
765     * the cluster multiplicity
766     *
767     * - if NCOG>0, the Center-Of-Gravity of the cluster IC
768     * evaluated using NCOG strips, even if they have a
769     * negative signal (according to Landi)
770     *
771     *-------------------------------------------------
772    
773    
774     include 'commontracker.f'
775     include 'calib.f'
776     include 'level1.f'
777    
778    
779    
780     if (ncog.gt.0) then
781     * ===========================
782     * ETA2 ETA3 ETA4 computation
783     * ===========================
784    
785     * --> signal of the central strip
786     sc = CLSIGNAL(INDMAX(ic)) !center
787     * signal of adjacent strips
788 pam-fi 1.14 sl1 = -9999. !left 1
789 mocchiut 1.1 if(
790     $ (INDMAX(ic)-1).ge.INDSTART(ic)
791     $ )
792     $ sl1 = CLSIGNAL(INDMAX(ic)-1)
793    
794 pam-fi 1.14 sl2 = -9999. !left 2
795 mocchiut 1.1 if(
796     $ (INDMAX(ic)-2).ge.INDSTART(ic)
797     $ )
798     $ sl2 = CLSIGNAL(INDMAX(ic)-2)
799    
800 pam-fi 1.14 sr1 = -9999. !right 1
801 mocchiut 1.1 if(
802     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
803     $ .or.
804     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
805     $ )
806     $ sr1 = CLSIGNAL(INDMAX(ic)+1)
807    
808 pam-fi 1.14 sr2 = -9999. !right 2
809 mocchiut 1.1 if(
810     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
811     $ .or.
812     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
813     $ )
814     $ sr2 = CLSIGNAL(INDMAX(ic)+2)
815    
816     COG = 0.
817    
818 pam-fi 1.6 c print*,'## ',sl2,sl1,sc,sr1,sr2
819    
820 pam-fi 1.14 c ==============================================================
821 pam-fi 1.4 if(ncog.eq.1)then
822     COG = 0.
823 pam-fi 1.14 if(sr1.gt.sc)cog=1. !NEW
824     if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
825     c ==============================================================
826 pam-fi 1.4 elseif(ncog.eq.2)then
827 mocchiut 1.1 if(sl1.gt.sr1)then
828 pam-fi 1.6 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)
829 pam-fi 1.14 elseif(sl1.lt.sr1)then
830     if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
831     elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
832     if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
833     $ .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
834     if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
835     $ .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
836     endif
837 pam-fi 1.15 c if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
838     c $ ,' : ',sl2,sl1,sc,sr1,sr2
839 pam-fi 1.14 c ==============================================================
840 mocchiut 1.1 elseif(ncog.eq.3)then
841 pam-fi 1.14 if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
842 pam-fi 1.15 c if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
843     c $ ,' : ',sl2,sl1,sc,sr1,sr2
844 pam-fi 1.14 c ==============================================================
845 mocchiut 1.1 elseif(ncog.eq.4)then
846     if(sl2.gt.sr2)then
847 pam-fi 1.6 if((sl2+sl1+sc+sr1).ne.0)
848     $ COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
849 pam-fi 1.14 elseif(sl2.lt.sr2)then
850     if((sr2+sl1+sc+sr1).ne.0)
851 pam-fi 1.12 $ COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
852 pam-fi 1.14 elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
853     if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
854     $ .and.(sl2+sl1+sc+sr1).ne.0 )
855     $ cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
856     if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
857     $ .and.(sr2+sl1+sc+sr1).ne.0 )
858     $ cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1) !NEW
859 mocchiut 1.1 endif
860     else
861     print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
862 pam-fi 1.7 $ ,' not implemented'
863 mocchiut 1.1 COG = 0.
864     endif
865    
866     c print*,'NCOG ',ncog,ic,' @@@ ',sl1,sc,sr1,' @@@ ',cog
867    
868     elseif(ncog.eq.0)then
869     * =========================
870     * COG computation
871     * =========================
872    
873     iv=VIEW(ic)
874     if(mod(iv,2).eq.1)incut=incuty
875     if(mod(iv,2).eq.0)incut=incutx
876 pam-fi 1.3 istart = INDSTART(IC)
877     istop = TOTCLLENGTH
878 mocchiut 1.1 if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
879 pam-fi 1.3 COG = 0
880 pam-fi 1.9 SGN = 0.
881 pam-fi 1.3 mu = 0
882 pam-fi 1.9 c print*,'-------'
883     do i = INDMAX(IC),istart,-1
884     ipos = i-INDMAX(ic)
885     cut = incut*CLSIGMA(i)
886     if(CLSIGNAL(i).ge.cut)then
887     COG = COG + ipos*CLSIGNAL(i)
888     SGN = SGN + CLSIGNAL(i)
889     mu = mu + 1
890 pam-fi 1.12 c print*,ipos,CLSIGNAL(i)
891 pam-fi 1.9 else
892     goto 10
893     endif
894     enddo
895     10 continue
896     do i = INDMAX(IC)+1,istop
897 pam-fi 1.3 ipos = i-INDMAX(ic)
898 pam-fi 1.4 cut = incut*CLSIGMA(i)
899 mocchiut 1.1 if(CLSIGNAL(i).ge.cut)then
900     COG = COG + ipos*CLSIGNAL(i)
901 pam-fi 1.9 SGN = SGN + CLSIGNAL(i)
902 mocchiut 1.1 mu = mu + 1
903 pam-fi 1.12 c print*,ipos,CLSIGNAL(i)
904 pam-fi 1.9 else
905     goto 20
906 mocchiut 1.1 endif
907     enddo
908 pam-fi 1.9 20 continue
909     if(SGN.le.0)then
910 pam-fi 1.12 print*,'cog(0,ic) --> ic, dedx ',ic,SGN
911 pam-fi 1.3 print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
912     print*,(CLSIGNAL(i),i=istart,istop)
913 pam-fi 1.9 c print*,'cog(0,ic) --> NOT EVALUATED '
914 pam-fi 1.3 else
915 pam-fi 1.9 COG=COG/SGN
916 pam-fi 1.3 endif
917 pam-fi 1.9 c print*,'-------'
918 mocchiut 1.1
919     else
920    
921     COG=0
922     print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
923     print*,' (NCOG must be >= 0)'
924    
925    
926     endif
927    
928 pam-fi 1.6 c print *,'## cog ',ncog,ic,cog,'/////////////'
929 mocchiut 1.1
930     return
931     end
932    
933     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
934 pam-fi 1.9
935 mocchiut 1.1 real function fbad_cog(ncog,ic)
936     *-------------------------------------------------------
937     * this function returns a factor that takes into
938     * account deterioration of the spatial resolution
939     * in the case BAD strips are included in the cluster.
940     * This factor should multiply the nominal spatial
941     * resolution.
942     *
943     *-------------------------------------------------------
944    
945     include 'commontracker.f'
946     include 'level1.f'
947     include 'calib.f'
948    
949     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
950 pam-fi 1.9 si = 8.4 !average good-strip noise
951     f = 4. !average bad-strip noise: f*si
952 mocchiut 1.1 incut=incuty
953     else !X-view
954 pam-fi 1.9 si = 3.9 !average good-strip noise
955     f = 6. !average bad-strip noise: f*si
956 mocchiut 1.1 incut=incutx
957     endif
958    
959     fbad_cog = 1.
960    
961     if (ncog.gt.0) then
962    
963     * --> signal of the central strip
964     sc = CLSIGNAL(INDMAX(ic)) !center
965     fsc = 1
966 pam-fi 1.9 c if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
967     fsc = clsigma(INDMAX(ic))/si
968 mocchiut 1.1 * --> signal of adjacent strips
969     sl1 = 0 !left 1
970     fsl1 = 1 !left 1
971     if(
972     $ (INDMAX(ic)-1).ge.INDSTART(ic)
973     $ )then
974     sl1 = CLSIGNAL(INDMAX(ic)-1)
975 pam-fi 1.9 c if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f
976     fsl1 = clsigma(INDMAX(ic)-1)/si
977 mocchiut 1.1 endif
978    
979     sl2 = 0 !left 2
980     fsl2 = 1 !left 2
981     if(
982     $ (INDMAX(ic)-2).ge.INDSTART(ic)
983     $ )then
984     sl2 = CLSIGNAL(INDMAX(ic)-2)
985 pam-fi 1.9 c if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f
986     fsl2 = clsigma(INDMAX(ic)-2)/si
987 mocchiut 1.1 endif
988     sr1 = 0 !right 1
989     fsr1 = 1 !right 1
990     if(
991     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
992     $ .or.
993     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
994     $ )then
995     sr1 = CLSIGNAL(INDMAX(ic)+1)
996 pam-fi 1.9 c if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f
997     fsr1 = clsigma(INDMAX(ic)+1)/si
998 mocchiut 1.1 endif
999     sr2 = 0 !right 2
1000     fsr2 = 1 !right 2
1001     if(
1002     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1003     $ .or.
1004     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1005     $ )then
1006     sr2 = CLSIGNAL(INDMAX(ic)+2)
1007 pam-fi 1.9 c if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f
1008     fsr2 = clsigma(INDMAX(ic)+2)/si
1009 mocchiut 1.1 endif
1010    
1011    
1012    
1013     ************************************************************
1014 pam-fi 1.9 * COG2-3-4 computation
1015 mocchiut 1.1 ************************************************************
1016    
1017     c print*,sl2,sl1,sc,sr1,sr2
1018    
1019 pam-fi 1.9 vCOG = cog(ncog,ic)!0.
1020 mocchiut 1.1
1021     if(ncog.eq.2)then
1022     if(sl1.gt.sr1)then
1023 pam-fi 1.9 c COG = -sl1/(sl1+sc)
1024     fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1025     fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1026 mocchiut 1.1 elseif(sl1.le.sr1)then
1027 pam-fi 1.9 c COG = sr1/(sc+sr1)
1028     fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1029     fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1030 mocchiut 1.1 endif
1031     elseif(ncog.eq.3)then
1032 pam-fi 1.9 c COG = (sr1-sl1)/(sl1+sc+sr1)
1033 mocchiut 1.1 fbad_cog =
1034 pam-fi 1.9 $ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1035 mocchiut 1.1 fbad_cog =
1036 pam-fi 1.9 $ fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1037 mocchiut 1.1 elseif(ncog.eq.4)then
1038     if(sl2.gt.sr2)then
1039 pam-fi 1.9 c COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1040 mocchiut 1.1 fbad_cog =
1041 pam-fi 1.9 $ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1042     $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1043 mocchiut 1.1 fbad_cog =
1044 pam-fi 1.9 $ fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1045     $ +(-vCOG)**2+(1-vCOG)**2)
1046 mocchiut 1.1 elseif(sl2.le.sr2)then
1047 pam-fi 1.9 c COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1048 mocchiut 1.1 fbad_cog =
1049 pam-fi 1.9 $ (fsl1*(-1-vCOG)**2
1050     $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1051 mocchiut 1.1 fbad_cog =
1052 pam-fi 1.9 $ fbad_cog / ((-1-vCOG)**2
1053     $ +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1054 mocchiut 1.1 endif
1055     else
1056     print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1057     print*,' (NCOG must be <= 4)'
1058 pam-fi 1.9 c COG = 0.
1059 mocchiut 1.1 endif
1060    
1061     elseif(ncog.eq.0)then
1062 pam-fi 1.9 * =========================
1063     * COG computation
1064     * =========================
1065    
1066     vCOG = cog(0,ic)
1067 mocchiut 1.1
1068 pam-fi 1.9 iv = VIEW(ic)
1069     istart = INDSTART(IC)
1070     istop = TOTCLLENGTH
1071     if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1072     SGN = 0.
1073     SNU = 0.
1074     SDE = 0.
1075     c$$$ do i=INDMAX(IC),istart,-1
1076     c$$$ ipos = i-INDMAX(ic)
1077     c$$$ cut = incut*CLSIGMA(i)
1078     c$$$ if(CLSIGNAL(i).gt.cut)then
1079     c$$$ COG = COG + ipos*CLSIGNAL(i)
1080     c$$$ SGN = SGN + CLSIGNAL(i)
1081     c$$$ else
1082     c$$$ goto 10
1083     c$$$ endif
1084     c$$$ enddo
1085     c$$$ 10 continue
1086     c$$$ do i=INDMAX(IC)+1,istop
1087     c$$$ ipos = i-INDMAX(ic)
1088     c$$$ cut = incut*CLSIGMA(i)
1089     c$$$ if(CLSIGNAL(i).gt.cut)then
1090     c$$$ COG = COG + ipos*CLSIGNAL(i)
1091     c$$$ SGN = SGN + CLSIGNAL(i)
1092     c$$$ else
1093     c$$$ goto 20
1094     c$$$ endif
1095     c$$$ enddo
1096     c$$$ 20 continue
1097     c$$$ if(SGN.le.0)then
1098     c$$$ print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1099     c$$$ print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1100     c$$$ print*,(CLSIGNAL(i),i=istart,istop)
1101     c$$$ print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1102     c$$$ else
1103     c$$$ COG=COG/SGN
1104     c$$$ endif
1105    
1106     do i=INDMAX(IC),istart,-1
1107     ipos = i-INDMAX(ic)
1108     cut = incut*CLSIGMA(i)
1109 mocchiut 1.1 if(CLSIGNAL(i).gt.cut)then
1110 pam-fi 1.9 fs = clsigma(i)/si
1111     SNU = SNU + fs*(ipos-vCOG)**2
1112     SDE = SDE + (ipos-vCOG)**2
1113     else
1114     goto 10
1115     endif
1116 mocchiut 1.1 enddo
1117 pam-fi 1.9 10 continue
1118     do i=INDMAX(IC)+1,istop
1119     ipos = i-INDMAX(ic)
1120     cut = incut*CLSIGMA(i)
1121 mocchiut 1.1 if(CLSIGNAL(i).gt.cut)then
1122 pam-fi 1.9 fs = clsigma(i)/si
1123     SNU = SNU + fs*(ipos-vCOG)**2
1124     SDE = SDE + (ipos-vCOG)**2
1125     else
1126     goto 20
1127 mocchiut 1.1 endif
1128     enddo
1129 pam-fi 1.9 20 continue
1130     if(SDE.ne.0)then
1131     FBAD_COG=SNU/SDE
1132     else
1133    
1134     endif
1135 mocchiut 1.1
1136     else
1137    
1138     FBAD_COG=0
1139     print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1140     print*,' (NCOG must be >= 0)'
1141    
1142    
1143     endif
1144    
1145    
1146     fbad_cog = sqrt(fbad_cog)
1147    
1148     return
1149     end
1150    
1151    
1152     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1153     real function fbad_cog0(ncog,ic)
1154     *-------------------------------------------------------
1155     * this function returns a factor that takes into
1156     * account deterioration of the spatial resolution
1157     * in the case BAD strips are included in the cluster.
1158     * This factor should multiply the nominal spatial
1159     * resolution.
1160     *
1161     * NB!!!
1162     * (this is the old version. It consider only the two
1163     * strips with the greatest signal. The new one is
1164     * fbad_cog(ncog,ic) )
1165     *
1166     *-------------------------------------------------------
1167    
1168     include 'commontracker.f'
1169     include 'level1.f'
1170     include 'calib.f'
1171    
1172     * --> signal of the central strip
1173     sc = CLSIGNAL(INDMAX(ic)) !center
1174    
1175     * signal of adjacent strips
1176     * --> left
1177     sl1 = 0 !left 1
1178     if(
1179     $ (INDMAX(ic)-1).ge.INDSTART(ic)
1180     $ )
1181     $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
1182    
1183     sl2 = 0 !left 2
1184     if(
1185     $ (INDMAX(ic)-2).ge.INDSTART(ic)
1186     $ )
1187     $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
1188    
1189     * --> right
1190     sr1 = 0 !right 1
1191     if(
1192     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1193     $ .or.
1194     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1195     $ )
1196     $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
1197    
1198     sr2 = 0 !right 2
1199     if(
1200     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1201     $ .or.
1202     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1203     $ )
1204     $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
1205    
1206    
1207     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1208     f = 4.
1209     si = 8.4
1210     else !X-view
1211     f = 6.
1212     si = 3.9
1213     endif
1214    
1215     fbad_cog = 1.
1216     f0 = 1
1217     f1 = 1
1218     f2 = 1
1219     f3 = 1
1220     if(sl1.gt.sr1.and.sl1.gt.0.)then
1221    
1222     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
1223     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
1224     c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
1225    
1226     if(ncog.eq.2.and.sl1.ne.0)then
1227     fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
1228     elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
1229     fbad_cog = 1.
1230     elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
1231     fbad_cog = 1.
1232     else
1233     fbad_cog = 1.
1234     endif
1235    
1236     elseif(sl1.le.sr1.and.sr1.gt.0.)then
1237    
1238    
1239     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
1240     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
1241     c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
1242    
1243     if(ncog.eq.2.and.sr1.ne.0)then
1244     fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
1245     elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
1246     fbad_cog = 1.
1247     elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
1248     fbad_cog = 1.
1249     else
1250     fbad_cog = 1.
1251     endif
1252    
1253     endif
1254    
1255     fbad_cog0 = sqrt(fbad_cog)
1256    
1257     return
1258     end
1259    
1260    
1261    
1262    
1263     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1264    
1265 pam-fi 1.13 FUNCTION risxeta2(x)
1266 mocchiut 1.1
1267     DOUBLE PRECISION V( 1)
1268     INTEGER NPAR, NDIM, IMQFUN, I, J
1269     DOUBLE PRECISION HQDJ, VV, VCONST
1270     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1271     DOUBLE PRECISION SIGV( 18, 1)
1272     DOUBLE PRECISION SIGDEL( 18)
1273     DOUBLE PRECISION SIGA( 18)
1274     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1275     DATA VCONST / 0.000000000000 /
1276     DATA SIGVMI / -20.50000000000 /
1277     DATA SIGVT / 41.00000000000 /
1278     DATA SIGV / 0.6097560748458E-01
1279     +, 0.1097560971975
1280     +, 0.1341463327408
1281     +, 0.1829268187284
1282     +, 0.2317073047161
1283     +, 0.4268292486668
1284     +, 0.4756097495556
1285     +, 0.4999999701977
1286     +, 0.5243902206421
1287     +, 0.5731707215309
1288     +, 0.7682926654816
1289     +, 0.8170731663704
1290     +, 0.8658536076546
1291     +, 0.8902438879013
1292     +, 0.9390243291855
1293     +, 0.000000000000
1294     +, 1.000000000000
1295     +, 0.3658536374569
1296     +/
1297     DATA SIGDEL / 0.4878048598766E-01
1298     +, 0.4878048598766E-01
1299     +, 0.4878048598766E-01
1300     +, 0.4878048598766E-01
1301     +, 0.4878048598766E-01
1302     +, 0.4878048598766E-01
1303     +, 0.4878048598766E-01
1304     +, 0.4878048598766E-01
1305     +, 0.4878048598766E-01
1306     +, 0.4878048598766E-01
1307     +, 0.4878048598766E-01
1308     +, 0.4878048598766E-01
1309     +, 0.4878048598766E-01
1310     +, 0.4878048598766E-01
1311     +, 0.4878048598766E-01
1312     +, 0.1999999994950E-05
1313     +, 0.1999999994950E-05
1314     +, 0.9756097197533E-01
1315     +/
1316     DATA SIGA / 51.65899502118
1317     +, -150.4733247841
1318     +, 143.0468613786
1319     +, -16.56096738997
1320     +, 5.149319798083
1321     +, 21.57149712673
1322     +, -39.46652322782
1323     +, 47.13181632948
1324     +, -32.93197883680
1325     +, 16.38645317092
1326     +, 1.453688482992
1327     +, -10.00547244421
1328     +, 131.3517670587
1329     +, -140.6351538257
1330     +, 49.05515749582
1331     +, -23.00028974788
1332     +, -22.58470403729
1333     +, -3.824682486418
1334     +/
1335    
1336     V(1)= abs(x)
1337 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1338 mocchiut 1.1
1339     HQUADF = 0.
1340     DO 20 J = 1, NPAR
1341     HQDJ = 0.
1342     DO 10 I = 1, NDIM
1343     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1344     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1345     10 CONTINUE
1346     HQDJ = HQDJ + SIGDEL (J) ** 2
1347     HQDJ = SQRT (HQDJ)
1348     HQUADF = HQUADF + SIGA (J) * HQDJ
1349     20 CONTINUE
1350     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1351    
1352 pam-fi 1.13 risxeta2=HQUADF* 1e-4
1353 mocchiut 1.1
1354     END
1355    
1356     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1357 pam-fi 1.13 FUNCTION risxeta3(x)
1358 mocchiut 1.1 DOUBLE PRECISION V( 1)
1359     INTEGER NPAR, NDIM, IMQFUN, I, J
1360     DOUBLE PRECISION HQDJ, VV, VCONST
1361     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1362     DOUBLE PRECISION SIGV( 18, 1)
1363     DOUBLE PRECISION SIGDEL( 18)
1364     DOUBLE PRECISION SIGA( 18)
1365     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1366     DATA VCONST / 0.000000000000 /
1367     DATA SIGVMI / -20.50000000000 /
1368     DATA SIGVT / 41.00000000000 /
1369     DATA SIGV / 0.6097560748458E-01
1370     +, 0.1097560971975
1371     +, 0.1341463327408
1372     +, 0.1829268187284
1373     +, 0.2317073047161
1374     +, 0.4756097495556
1375     +, 0.4999999701977
1376     +, 0.5243902206421
1377     +, 0.7682926654816
1378     +, 0.8170731663704
1379     +, 0.8658536076546
1380     +, 0.8902438879013
1381     +, 0.9390243291855
1382     +, 0.000000000000
1383     +, 1.000000000000
1384     +, 0.3658536374569
1385     +, 0.4146341383457
1386     +, 0.6097560524940
1387     +/
1388     DATA SIGDEL / 0.4878048598766E-01
1389     +, 0.4878048598766E-01
1390     +, 0.4878048598766E-01
1391     +, 0.4878048598766E-01
1392     +, 0.4878048598766E-01
1393     +, 0.4878048598766E-01
1394     +, 0.4878048598766E-01
1395     +, 0.4878048598766E-01
1396     +, 0.4878048598766E-01
1397     +, 0.4878048598766E-01
1398     +, 0.4878048598766E-01
1399     +, 0.4878048598766E-01
1400     +, 0.4878048598766E-01
1401     +, 0.1999999994950E-05
1402     +, 0.1999999994950E-05
1403     +, 0.9756097197533E-01
1404     +, 0.9756097197533E-01
1405     +, 0.9756097197533E-01
1406     +/
1407     DATA SIGA / 55.18284054458
1408     +, -160.3358431242
1409     +, 144.6939185763
1410     +, -20.45200854118
1411     +, 5.223570087108
1412     +,-0.4171476953945
1413     +, -27.67911907462
1414     +, 17.70327157495
1415     +, -1.867165491707
1416     +, -8.884458169181
1417     +, 124.3526608791
1418     +, -143.3309398345
1419     +, 50.80345027122
1420     +, -16.44454904415
1421     +, -15.73785568450
1422     +, -22.71810502561
1423     +, 36.86170101430
1424     +, 2.437918198452
1425     +/
1426    
1427     V(1) = abs(x)
1428 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1429    
1430 mocchiut 1.1 HQUADF = 0.
1431     DO 20 J = 1, NPAR
1432     HQDJ = 0.
1433     DO 10 I = 1, NDIM
1434     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1435     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1436     10 CONTINUE
1437     HQDJ = HQDJ + SIGDEL (J) ** 2
1438     HQDJ = SQRT (HQDJ)
1439     HQUADF = HQUADF + SIGA (J) * HQDJ
1440     20 CONTINUE
1441     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1442    
1443 pam-fi 1.13 risxeta3 = HQUADF* 1e-4
1444 mocchiut 1.1
1445     END
1446     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1447 pam-fi 1.13 FUNCTION risxeta4(x)
1448 mocchiut 1.1 DOUBLE PRECISION V( 1)
1449     INTEGER NPAR, NDIM, IMQFUN, I, J
1450     DOUBLE PRECISION HQDJ, VV, VCONST
1451     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1452     DOUBLE PRECISION SIGV( 18, 1)
1453     DOUBLE PRECISION SIGDEL( 18)
1454     DOUBLE PRECISION SIGA( 18)
1455     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1456     DATA VCONST / 0.000000000000 /
1457     DATA SIGVMI / -20.50000000000 /
1458     DATA SIGVT / 41.00000000000 /
1459     DATA SIGV / 0.3658536449075E-01
1460     +, 0.6097560748458E-01
1461     +, 0.1097560971975
1462     +, 0.1341463327408
1463     +, 0.4756097495556
1464     +, 0.5243902206421
1465     +, 0.8658536076546
1466     +, 0.8902438879013
1467     +, 0.9390243291855
1468     +, 0.9634146094322
1469     +, 0.000000000000
1470     +, 1.000000000000
1471     +, 0.3658536374569
1472     +, 0.4146341383457
1473     +, 0.6097560524940
1474     +, 0.6585365533829
1475     +, 0.7560975551605
1476     +, 0.2439024299383
1477     +/
1478     DATA SIGDEL / 0.4878048598766E-01
1479     +, 0.4878048598766E-01
1480     +, 0.4878048598766E-01
1481     +, 0.4878048598766E-01
1482     +, 0.4878048598766E-01
1483     +, 0.4878048598766E-01
1484     +, 0.4878048598766E-01
1485     +, 0.4878048598766E-01
1486     +, 0.4878048598766E-01
1487     +, 0.4878048598766E-01
1488     +, 0.1999999994950E-05
1489     +, 0.1999999994950E-05
1490     +, 0.9756097197533E-01
1491     +, 0.9756097197533E-01
1492     +, 0.9756097197533E-01
1493     +, 0.9756097197533E-01
1494     +, 0.9756097197533E-01
1495     +, 0.1951219439507
1496     +/
1497     DATA SIGA / -43.61551887895
1498     +, 57.88466995373
1499     +, -92.04113299504
1500     +, 74.08166649890
1501     +, -9.768686062558
1502     +, -4.304496875334
1503     +, 72.62237333937
1504     +, -91.21920840618
1505     +, 56.75519978630
1506     +, -43.21115751243
1507     +, 12.79984505413
1508     +, 12.10074868595
1509     +, -6.238587250860
1510     +, 23.43447356326
1511     +, 17.98221401495
1512     +, -7.980332610975
1513     +, -3.426733307051
1514     +, -8.683439558751
1515     +/
1516    
1517     V(1)=abs(x)
1518 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1519 mocchiut 1.1
1520     HQUADF = 0.
1521     DO 20 J = 1, NPAR
1522     HQDJ = 0.
1523     DO 10 I = 1, NDIM
1524     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1525     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1526     10 CONTINUE
1527     HQDJ = HQDJ + SIGDEL (J) ** 2
1528     HQDJ = SQRT (HQDJ)
1529     HQUADF = HQUADF + SIGA (J) * HQDJ
1530     20 CONTINUE
1531     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1532    
1533 pam-fi 1.13 risxeta4=HQUADF* 1e-4
1534 mocchiut 1.1
1535     END
1536     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1537 pam-fi 1.13 FUNCTION risyeta2(x)
1538 mocchiut 1.1 DOUBLE PRECISION V( 1)
1539     INTEGER NPAR, NDIM, IMQFUN, I, J
1540     DOUBLE PRECISION HQDJ, VV, VCONST
1541     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1542     DOUBLE PRECISION SIGV( 12, 1)
1543     DOUBLE PRECISION SIGDEL( 12)
1544     DOUBLE PRECISION SIGA( 12)
1545     DATA NPAR, NDIM, IMQFUN / 12, 1, 1/
1546     DATA VCONST / 0.000000000000 /
1547     DATA SIGVMI / -20.50000000000 /
1548     DATA SIGVT / 41.00000000000 /
1549     DATA SIGV / 0.1585365831852
1550     +, 0.4024389982224
1551     +, 0.4756097495556
1552     +, 0.5243902206421
1553     +, 0.5975609421730
1554     +, 0.8414633870125
1555     +, 0.000000000000
1556     +, 1.000000000000
1557     +, 0.2682926654816
1558     +, 0.3170731663704
1559     +, 0.7073170542717
1560     +, 0.7560975551605
1561     +/
1562     DATA SIGDEL / 0.4878048598766E-01
1563     +, 0.4878048598766E-01
1564     +, 0.4878048598766E-01
1565     +, 0.4878048598766E-01
1566     +, 0.4878048598766E-01
1567     +, 0.4878048598766E-01
1568     +, 0.1999999994950E-05
1569     +, 0.1999999994950E-05
1570     +, 0.9756097197533E-01
1571     +, 0.9756097197533E-01
1572     +, 0.9756097197533E-01
1573     +, 0.9756097197533E-01
1574     +/
1575     DATA SIGA / 14.57433603529
1576     +, -15.93532436156
1577     +, -13.24628335221
1578     +, -14.31193855410
1579     +, -12.67339684488
1580     +, 18.19876051780
1581     +, -5.270493486725
1582     +, -5.107670990828
1583     +, -9.553262933901
1584     +, 43.34150727448
1585     +, 55.91366786432
1586     +, -29.38037318563
1587     +/
1588    
1589     v(1)= abs(x)
1590 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1591 mocchiut 1.1
1592     HQUADF = 0.
1593     DO 20 J = 1, NPAR
1594     HQDJ = 0.
1595     DO 10 I = 1, NDIM
1596     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1597     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1598     10 CONTINUE
1599     HQDJ = HQDJ + SIGDEL (J) ** 2
1600     HQDJ = SQRT (HQDJ)
1601     HQUADF = HQUADF + SIGA (J) * HQDJ
1602     20 CONTINUE
1603     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1604    
1605 pam-fi 1.13 risyeta2=HQUADF* 1e-4
1606 mocchiut 1.1
1607     END
1608     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1609    
1610     FUNCTION risy_cog(x)
1611     DOUBLE PRECISION V( 1)
1612     INTEGER NPAR, NDIM, IMQFUN, I, J
1613     DOUBLE PRECISION HQDJ, VV, VCONST
1614     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1615     DOUBLE PRECISION SIGV( 10, 1)
1616     DOUBLE PRECISION SIGDEL( 10)
1617     DOUBLE PRECISION SIGA( 10)
1618     DATA NPAR, NDIM, IMQFUN / 10, 1, 1/
1619     DATA VCONST / 0.000000000000 /
1620     DATA SIGVMI / -20.50000000000 /
1621     DATA SIGVT / 41.00000000000 /
1622     DATA SIGV / 0.1585365831852
1623     +, 0.8414633870125
1624     +, 0.000000000000
1625     +, 1.000000000000
1626     +, 0.4634146094322
1627     +, 0.5121951103210
1628     +, 0.5609756112099
1629     +, 0.6585365533829
1630     +, 0.7073170542717
1631     +, 0.3414633870125
1632     +/
1633     DATA SIGDEL / 0.4878048598766E-01
1634     +, 0.4878048598766E-01
1635     +, 0.1999999994950E-05
1636     +, 0.1999999994950E-05
1637     +, 0.9756097197533E-01
1638     +, 0.9756097197533E-01
1639     +, 0.9756097197533E-01
1640     +, 0.9756097197533E-01
1641     +, 0.9756097197533E-01
1642     +, 0.1951219439507
1643     +/
1644     DATA SIGA / 23.73833445988
1645     +, 24.10182100013
1646     +, 1.865894323190
1647     +, 1.706006262931
1648     +, -1.075607857202
1649     +, -22.11489493403
1650     +, 1.663100707801
1651     +, 4.089852595440
1652     +, -4.314993873697
1653     +, -2.174479487744
1654     +/
1655    
1656     V(1)=abs(x)
1657 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1658 mocchiut 1.1
1659     HQUADF = 0.
1660     DO 20 J = 1, NPAR
1661     HQDJ = 0.
1662     DO 10 I = 1, NDIM
1663     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1664     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1665     10 CONTINUE
1666     HQDJ = HQDJ + SIGDEL (J) ** 2
1667     HQDJ = SQRT (HQDJ)
1668     HQUADF = HQUADF + SIGA (J) * HQDJ
1669     20 CONTINUE
1670     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1671    
1672     risy_cog=HQUADF* 1e-4
1673    
1674     END
1675     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1676     FUNCTION risx_cog(x)
1677     DOUBLE PRECISION V( 1)
1678     INTEGER NPAR, NDIM, IMQFUN, I, J
1679     DOUBLE PRECISION HQDJ, VV, VCONST
1680     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1681     DOUBLE PRECISION SIGV( 15, 1)
1682     DOUBLE PRECISION SIGDEL( 15)
1683     DOUBLE PRECISION SIGA( 15)
1684     DATA NPAR, NDIM, IMQFUN / 15, 1, 1/
1685     DATA VCONST / 0.000000000000 /
1686     DATA SIGVMI / -20.50000000000 /
1687     DATA SIGVT / 41.00000000000 /
1688     DATA SIGV / 0.6097560748458E-01
1689     +, 0.8536584675312E-01
1690     +, 0.1341463327408
1691     +, 0.2317073047161
1692     +, 0.2804878056049
1693     +, 0.3780487775803
1694     +, 0.6219512224197
1695     +, 0.7195121645927
1696     +, 0.7682926654816
1697     +, 0.8658536076546
1698     +, 0.9146341085434
1699     +, 0.9390243291855
1700     +, 0.000000000000
1701     +, 1.000000000000
1702     +, 0.5121951103210
1703     +/
1704     DATA SIGDEL / 0.4878048598766E-01
1705     +, 0.4878048598766E-01
1706     +, 0.4878048598766E-01
1707     +, 0.4878048598766E-01
1708     +, 0.4878048598766E-01
1709     +, 0.4878048598766E-01
1710     +, 0.4878048598766E-01
1711     +, 0.4878048598766E-01
1712     +, 0.4878048598766E-01
1713     +, 0.4878048598766E-01
1714     +, 0.4878048598766E-01
1715     +, 0.4878048598766E-01
1716     +, 0.1999999994950E-05
1717     +, 0.1999999994950E-05
1718     +, 0.9756097197533E-01
1719     +/
1720     DATA SIGA / 31.95672945139
1721     +, -34.23286209245
1722     +, -6.298459168211
1723     +, 10.98847700545
1724     +,-0.3052213535054
1725     +, 13.10517991464
1726     +, 15.60290821679
1727     +, -1.956118448507
1728     +, 12.41453816720
1729     +, -7.354056408553
1730     +, -32.32512668778
1731     +, 30.61116178966
1732     +, 1.418505329236
1733     +, 1.583492573619
1734     +, -18.48799977042
1735     +/
1736    
1737     V(1)=abs(x)
1738 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1739 mocchiut 1.1
1740     HQUADF = 0.
1741     DO 20 J = 1, NPAR
1742     HQDJ = 0.
1743     DO 10 I = 1, NDIM
1744     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1745     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1746     10 CONTINUE
1747     HQDJ = HQDJ + SIGDEL (J) ** 2
1748     HQDJ = SQRT (HQDJ)
1749     HQUADF = HQUADF + SIGA (J) * HQDJ
1750     20 CONTINUE
1751     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1752    
1753     risx_cog = HQUADF * 1e-4
1754    
1755     END
1756 pam-fi 1.16
1757    
1758     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1759 pam-fi 1.19 real function pfacorr(ic,angle)
1760 pam-fi 1.16 *--------------------------------------------------------------
1761     * this function returns the landi correction for this cluster
1762     *--------------------------------------------------------------
1763     include 'commontracker.f'
1764     include 'calib.f'
1765     include 'level1.f'
1766    
1767     real angle
1768     integer iview,lad
1769    
1770     iview = VIEW(ic)
1771     lad = nld(MAXS(ic),VIEW(ic))
1772    
1773     * find angular bin
1774     * (in futuro possiamo pensare di interpolare anche sull'angolo)
1775     do iang=1,nangbin
1776     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1777     iangle=iang
1778     goto 98
1779     endif
1780     enddo
1781 pam-fi 1.18 if(DEBUG.eq.1)
1782 pam-fi 1.16 $ print*,'pfacorr *** warning *** angle out of range: ',angle
1783     if(angle.lt.angL(1))iang=1
1784     if(angle.gt.angR(nangbin))iang=nangbin
1785     98 continue !jump here if ok
1786    
1787     pfacorr = fcorr(iview,lad,iang)
1788    
1789 pam-fi 1.18 if(DEBUG.eq.1)print*,'CORR (ic ',ic,' ang',angle,') -->',pfacorr
1790 pam-fi 1.16
1791    
1792     100 return
1793     end

  ViewVC Help
Powered by ViewVC 1.1.23