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