/[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.20 - (hide annotations) (download)
Mon Aug 27 12:57:15 2007 UTC (17 years, 3 months ago) by pam-fi
Branch: MAIN
Changes since 1.19: +195 -296 lines
fixed bug on pfa + new methods added

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

  ViewVC Help
Powered by ViewVC 1.1.23