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