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