/[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.3 - (hide annotations) (download)
Thu Sep 28 14:04:40 2006 UTC (18 years, 2 months ago) by pam-fi
Branch: MAIN
Changes since 1.2: +19 -10 lines 
some bugs fixed, some changings in the classes:
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 pam-fi 1.3 istart = INDSTART(IC)
690     istop = TOTCLLENGTH
691 mocchiut 1.1 if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
692 pam-fi 1.3 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 mocchiut 1.1 * print*,'******************',istart,istop,ipos
699     * $ ,MAXS(ic)+ipos,iv,ivk,is
700 pam-fi 1.3 cut = incut*SIGMA(iv,ivk,is)
701     if(SIGMA(iv,ivk,is).ne.CLSIGMA(i))
702     $ print*,'cog(0,ic) --> hai fatto qualche cazzata'
703 mocchiut 1.1 if(CLSIGNAL(i).ge.cut)then
704     COG = COG + ipos*CLSIGNAL(i)
705     mu = mu + 1
706     c print*,ipos,CLSIGNAL(i),incut,cut
707     endif
708     enddo
709 pam-fi 1.3 if(DEDX(ic).le.0)then
710     print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)
711     print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
712     print*,(CLSIGNAL(i),i=istart,istop)
713     print*,'cog(0,ic) --> NOT EVALUATED '
714     else
715     COG=COG/DEDX(ic)
716     endif
717 mocchiut 1.1 c if(DEBUG)print*,'COG (ic ',ic,' m',mu,')'
718     c $ ,cog
719    
720     else
721    
722     COG=0
723     print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
724     print*,' (NCOG must be >= 0)'
725    
726    
727     endif
728    
729     c print *,ncog,ic,cog,'/////////////'
730    
731     return
732     end
733    
734     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
735     real function fbad_cog(ncog,ic)
736     *-------------------------------------------------------
737     * this function returns a factor that takes into
738     * account deterioration of the spatial resolution
739     * in the case BAD strips are included in the cluster.
740     * This factor should multiply the nominal spatial
741     * resolution.
742     *
743     *-------------------------------------------------------
744    
745     include 'commontracker.f'
746     include 'level1.f'
747     include 'calib.f'
748    
749     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
750     f = 4.
751     si = 8.4
752     incut=incuty
753     else !X-view
754     f = 6.
755     si = 3.9
756     incut=incutx
757     endif
758    
759     fbad_cog = 1.
760    
761     if (ncog.gt.0) then
762    
763     * --> signal of the central strip
764     sc = CLSIGNAL(INDMAX(ic)) !center
765     fsc = 1
766     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)fsc=f
767     * --> signal of adjacent strips
768     sl1 = 0 !left 1
769     fsl1 = 1 !left 1
770     if(
771     $ (INDMAX(ic)-1).ge.INDSTART(ic)
772     $ )then
773     sl1 = CLSIGNAL(INDMAX(ic)-1)
774     if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f
775     c else
776     c fsl1 = 0
777     endif
778    
779     sl2 = 0 !left 2
780     fsl2 = 1 !left 2
781     if(
782     $ (INDMAX(ic)-2).ge.INDSTART(ic)
783     $ )then
784     sl2 = CLSIGNAL(INDMAX(ic)-2)
785     if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f
786     c else
787     c fsl2 = 0
788     endif
789     sr1 = 0 !right 1
790     fsr1 = 1 !right 1
791     if(
792     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
793     $ .or.
794     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
795     $ )then
796     sr1 = CLSIGNAL(INDMAX(ic)+1)
797     if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f
798     c else
799     c fsr1 = 0
800     endif
801     sr2 = 0 !right 2
802     fsr2 = 1 !right 2
803     if(
804     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
805     $ .or.
806     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
807     $ )then
808     sr2 = CLSIGNAL(INDMAX(ic)+2)
809     if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f
810     c else
811     c fsr2 = 0
812     endif
813    
814    
815    
816     ************************************************************
817     * COG computation
818     ************************************************************
819    
820     c print*,sl2,sl1,sc,sr1,sr2
821    
822     COG = 0.
823    
824     if(ncog.eq.2)then
825     if(sl1.gt.sr1)then
826     COG = -sl1/(sl1+sc)
827     fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)
828     fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)
829     elseif(sl1.le.sr1)then
830     COG = sr1/(sc+sr1)
831     fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)
832     fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)
833     endif
834     elseif(ncog.eq.3)then
835     COG = (sr1-sl1)/(sl1+sc+sr1)
836     fbad_cog =
837     $ (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)
838     fbad_cog =
839     $ fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)
840     elseif(ncog.eq.4)then
841     if(sl2.gt.sr2)then
842     COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
843     fbad_cog =
844     $ (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2
845     $ +fsc*(-COG)**2+fsr1*(1-COG)**2)
846     fbad_cog =
847     $ fbad_cog / ((-2-COG)**2+(-1-COG)**2
848     $ +(-COG)**2+(1-COG)**2)
849     elseif(sl2.le.sr2)then
850     COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
851     fbad_cog =
852     $ (fsl1*(-1-COG)**2
853     $ +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)
854     fbad_cog =
855     $ fbad_cog / ((-1-COG)**2
856     $ +(-COG)**2+(1-COG)**2+(2-COG)**2)
857     endif
858     else
859     print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
860     print*,' (NCOG must be <= 4)'
861     COG = 0.
862     endif
863    
864     elseif(ncog.eq.0)then
865    
866    
867     iv=VIEW(ic)
868     istart=INDSTART(IC)
869     istop=TOTCLLENGTH
870     if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
871     COG=0.
872     SNU=0.
873     SDE=0.
874     do i=istart,istop
875     ipos=i-INDMAX(ic)
876     il=nvk(MAXS(ic)+ipos)
877     is=nst(MAXS(ic)+ipos)
878     cut=incut*SIGMA(iv,il,is)
879     if(CLSIGNAL(i).gt.cut)then
880     COG = COG + ipos*CLSIGNAL(i)
881     endif
882     enddo
883     COG=COG/DEDX(ic)
884     do i=istart,istop
885     ipos=i-INDMAX(ic)
886     il=nvk(MAXS(ic)+ipos)
887     is=nst(MAXS(ic)+ipos)
888     cut=incut*SIGMA(iv,il,is)
889     if(CLSIGNAL(i).gt.cut)then
890     fs=1
891     if(BAD(iv,il,is).eq.0)fs=f
892     SNU = SNU + fs*(ipos-COG)**2
893     SDE = SDE + (ipos-COG)**2
894     endif
895     enddo
896     if(SDE.ne.0)FBAD_COG=SNU/SDE
897    
898     else
899    
900     FBAD_COG=0
901     print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
902     print*,' (NCOG must be >= 0)'
903    
904    
905     endif
906    
907    
908     fbad_cog = sqrt(fbad_cog)
909    
910     return
911     end
912    
913    
914     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
915     real function fbad_cog0(ncog,ic)
916     *-------------------------------------------------------
917     * this function returns a factor that takes into
918     * account deterioration of the spatial resolution
919     * in the case BAD strips are included in the cluster.
920     * This factor should multiply the nominal spatial
921     * resolution.
922     *
923     * NB!!!
924     * (this is the old version. It consider only the two
925     * strips with the greatest signal. The new one is
926     * fbad_cog(ncog,ic) )
927     *
928     *-------------------------------------------------------
929    
930     include 'commontracker.f'
931     include 'level1.f'
932     include 'calib.f'
933    
934     * --> signal of the central strip
935     sc = CLSIGNAL(INDMAX(ic)) !center
936    
937     * signal of adjacent strips
938     * --> left
939     sl1 = 0 !left 1
940     if(
941     $ (INDMAX(ic)-1).ge.INDSTART(ic)
942     $ )
943     $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
944    
945     sl2 = 0 !left 2
946     if(
947     $ (INDMAX(ic)-2).ge.INDSTART(ic)
948     $ )
949     $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
950    
951     * --> right
952     sr1 = 0 !right 1
953     if(
954     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
955     $ .or.
956     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
957     $ )
958     $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
959    
960     sr2 = 0 !right 2
961     if(
962     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
963     $ .or.
964     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
965     $ )
966     $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
967    
968    
969     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
970     f = 4.
971     si = 8.4
972     else !X-view
973     f = 6.
974     si = 3.9
975     endif
976    
977     fbad_cog = 1.
978     f0 = 1
979     f1 = 1
980     f2 = 1
981     f3 = 1
982     if(sl1.gt.sr1.and.sl1.gt.0.)then
983    
984     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
985     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
986     c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
987    
988     if(ncog.eq.2.and.sl1.ne.0)then
989     fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
990     elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
991     fbad_cog = 1.
992     elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
993     fbad_cog = 1.
994     else
995     fbad_cog = 1.
996     endif
997    
998     elseif(sl1.le.sr1.and.sr1.gt.0.)then
999    
1000    
1001     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
1002     if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
1003     c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
1004    
1005     if(ncog.eq.2.and.sr1.ne.0)then
1006     fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
1007     elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
1008     fbad_cog = 1.
1009     elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
1010     fbad_cog = 1.
1011     else
1012     fbad_cog = 1.
1013     endif
1014    
1015     endif
1016    
1017     fbad_cog0 = sqrt(fbad_cog)
1018    
1019     return
1020     end
1021    
1022    
1023    
1024    
1025     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1026    
1027     FUNCTION risx_eta2(x)
1028    
1029     DOUBLE PRECISION V( 1)
1030     INTEGER NPAR, NDIM, IMQFUN, I, J
1031     DOUBLE PRECISION HQDJ, VV, VCONST
1032     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1033     DOUBLE PRECISION SIGV( 18, 1)
1034     DOUBLE PRECISION SIGDEL( 18)
1035     DOUBLE PRECISION SIGA( 18)
1036     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1037     DATA VCONST / 0.000000000000 /
1038     DATA SIGVMI / -20.50000000000 /
1039     DATA SIGVT / 41.00000000000 /
1040     DATA SIGV / 0.6097560748458E-01
1041     +, 0.1097560971975
1042     +, 0.1341463327408
1043     +, 0.1829268187284
1044     +, 0.2317073047161
1045     +, 0.4268292486668
1046     +, 0.4756097495556
1047     +, 0.4999999701977
1048     +, 0.5243902206421
1049     +, 0.5731707215309
1050     +, 0.7682926654816
1051     +, 0.8170731663704
1052     +, 0.8658536076546
1053     +, 0.8902438879013
1054     +, 0.9390243291855
1055     +, 0.000000000000
1056     +, 1.000000000000
1057     +, 0.3658536374569
1058     +/
1059     DATA SIGDEL / 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.4878048598766E-01
1066     +, 0.4878048598766E-01
1067     +, 0.4878048598766E-01
1068     +, 0.4878048598766E-01
1069     +, 0.4878048598766E-01
1070     +, 0.4878048598766E-01
1071     +, 0.4878048598766E-01
1072     +, 0.4878048598766E-01
1073     +, 0.4878048598766E-01
1074     +, 0.1999999994950E-05
1075     +, 0.1999999994950E-05
1076     +, 0.9756097197533E-01
1077     +/
1078     DATA SIGA / 51.65899502118
1079     +, -150.4733247841
1080     +, 143.0468613786
1081     +, -16.56096738997
1082     +, 5.149319798083
1083     +, 21.57149712673
1084     +, -39.46652322782
1085     +, 47.13181632948
1086     +, -32.93197883680
1087     +, 16.38645317092
1088     +, 1.453688482992
1089     +, -10.00547244421
1090     +, 131.3517670587
1091     +, -140.6351538257
1092     +, 49.05515749582
1093     +, -23.00028974788
1094     +, -22.58470403729
1095     +, -3.824682486418
1096     +/
1097    
1098     V(1)= abs(x)
1099    
1100     HQUADF = 0.
1101     DO 20 J = 1, NPAR
1102     HQDJ = 0.
1103     DO 10 I = 1, NDIM
1104     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1105     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1106     10 CONTINUE
1107     HQDJ = HQDJ + SIGDEL (J) ** 2
1108     HQDJ = SQRT (HQDJ)
1109     HQUADF = HQUADF + SIGA (J) * HQDJ
1110     20 CONTINUE
1111     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1112    
1113     risx_eta2=HQUADF* 1e-4
1114    
1115     END
1116    
1117     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1118     FUNCTION risx_eta3(x)
1119     DOUBLE PRECISION V( 1)
1120     INTEGER NPAR, NDIM, IMQFUN, I, J
1121     DOUBLE PRECISION HQDJ, VV, VCONST
1122     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1123     DOUBLE PRECISION SIGV( 18, 1)
1124     DOUBLE PRECISION SIGDEL( 18)
1125     DOUBLE PRECISION SIGA( 18)
1126     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1127     DATA VCONST / 0.000000000000 /
1128     DATA SIGVMI / -20.50000000000 /
1129     DATA SIGVT / 41.00000000000 /
1130     DATA SIGV / 0.6097560748458E-01
1131     +, 0.1097560971975
1132     +, 0.1341463327408
1133     +, 0.1829268187284
1134     +, 0.2317073047161
1135     +, 0.4756097495556
1136     +, 0.4999999701977
1137     +, 0.5243902206421
1138     +, 0.7682926654816
1139     +, 0.8170731663704
1140     +, 0.8658536076546
1141     +, 0.8902438879013
1142     +, 0.9390243291855
1143     +, 0.000000000000
1144     +, 1.000000000000
1145     +, 0.3658536374569
1146     +, 0.4146341383457
1147     +, 0.6097560524940
1148     +/
1149     DATA SIGDEL / 0.4878048598766E-01
1150     +, 0.4878048598766E-01
1151     +, 0.4878048598766E-01
1152     +, 0.4878048598766E-01
1153     +, 0.4878048598766E-01
1154     +, 0.4878048598766E-01
1155     +, 0.4878048598766E-01
1156     +, 0.4878048598766E-01
1157     +, 0.4878048598766E-01
1158     +, 0.4878048598766E-01
1159     +, 0.4878048598766E-01
1160     +, 0.4878048598766E-01
1161     +, 0.4878048598766E-01
1162     +, 0.1999999994950E-05
1163     +, 0.1999999994950E-05
1164     +, 0.9756097197533E-01
1165     +, 0.9756097197533E-01
1166     +, 0.9756097197533E-01
1167     +/
1168     DATA SIGA / 55.18284054458
1169     +, -160.3358431242
1170     +, 144.6939185763
1171     +, -20.45200854118
1172     +, 5.223570087108
1173     +,-0.4171476953945
1174     +, -27.67911907462
1175     +, 17.70327157495
1176     +, -1.867165491707
1177     +, -8.884458169181
1178     +, 124.3526608791
1179     +, -143.3309398345
1180     +, 50.80345027122
1181     +, -16.44454904415
1182     +, -15.73785568450
1183     +, -22.71810502561
1184     +, 36.86170101430
1185     +, 2.437918198452
1186     +/
1187    
1188     V(1) = abs(x)
1189     HQUADF = 0.
1190     DO 20 J = 1, NPAR
1191     HQDJ = 0.
1192     DO 10 I = 1, NDIM
1193     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1194     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1195     10 CONTINUE
1196     HQDJ = HQDJ + SIGDEL (J) ** 2
1197     HQDJ = SQRT (HQDJ)
1198     HQUADF = HQUADF + SIGA (J) * HQDJ
1199     20 CONTINUE
1200     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1201    
1202     risx_eta3 = HQUADF* 1e-4
1203    
1204     END
1205     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1206     FUNCTION risx_eta4(x)
1207     DOUBLE PRECISION V( 1)
1208     INTEGER NPAR, NDIM, IMQFUN, I, J
1209     DOUBLE PRECISION HQDJ, VV, VCONST
1210     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1211     DOUBLE PRECISION SIGV( 18, 1)
1212     DOUBLE PRECISION SIGDEL( 18)
1213     DOUBLE PRECISION SIGA( 18)
1214     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1215     DATA VCONST / 0.000000000000 /
1216     DATA SIGVMI / -20.50000000000 /
1217     DATA SIGVT / 41.00000000000 /
1218     DATA SIGV / 0.3658536449075E-01
1219     +, 0.6097560748458E-01
1220     +, 0.1097560971975
1221     +, 0.1341463327408
1222     +, 0.4756097495556
1223     +, 0.5243902206421
1224     +, 0.8658536076546
1225     +, 0.8902438879013
1226     +, 0.9390243291855
1227     +, 0.9634146094322
1228     +, 0.000000000000
1229     +, 1.000000000000
1230     +, 0.3658536374569
1231     +, 0.4146341383457
1232     +, 0.6097560524940
1233     +, 0.6585365533829
1234     +, 0.7560975551605
1235     +, 0.2439024299383
1236     +/
1237     DATA SIGDEL / 0.4878048598766E-01
1238     +, 0.4878048598766E-01
1239     +, 0.4878048598766E-01
1240     +, 0.4878048598766E-01
1241     +, 0.4878048598766E-01
1242     +, 0.4878048598766E-01
1243     +, 0.4878048598766E-01
1244     +, 0.4878048598766E-01
1245     +, 0.4878048598766E-01
1246     +, 0.4878048598766E-01
1247     +, 0.1999999994950E-05
1248     +, 0.1999999994950E-05
1249     +, 0.9756097197533E-01
1250     +, 0.9756097197533E-01
1251     +, 0.9756097197533E-01
1252     +, 0.9756097197533E-01
1253     +, 0.9756097197533E-01
1254     +, 0.1951219439507
1255     +/
1256     DATA SIGA / -43.61551887895
1257     +, 57.88466995373
1258     +, -92.04113299504
1259     +, 74.08166649890
1260     +, -9.768686062558
1261     +, -4.304496875334
1262     +, 72.62237333937
1263     +, -91.21920840618
1264     +, 56.75519978630
1265     +, -43.21115751243
1266     +, 12.79984505413
1267     +, 12.10074868595
1268     +, -6.238587250860
1269     +, 23.43447356326
1270     +, 17.98221401495
1271     +, -7.980332610975
1272     +, -3.426733307051
1273     +, -8.683439558751
1274     +/
1275    
1276     V(1)=abs(x)
1277    
1278     HQUADF = 0.
1279     DO 20 J = 1, NPAR
1280     HQDJ = 0.
1281     DO 10 I = 1, NDIM
1282     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1283     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1284     10 CONTINUE
1285     HQDJ = HQDJ + SIGDEL (J) ** 2
1286     HQDJ = SQRT (HQDJ)
1287     HQUADF = HQUADF + SIGA (J) * HQDJ
1288     20 CONTINUE
1289     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1290    
1291     risx_eta4=HQUADF* 1e-4
1292    
1293     END
1294     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1295     FUNCTION risy_eta2(x)
1296     DOUBLE PRECISION V( 1)
1297     INTEGER NPAR, NDIM, IMQFUN, I, J
1298     DOUBLE PRECISION HQDJ, VV, VCONST
1299     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1300     DOUBLE PRECISION SIGV( 12, 1)
1301     DOUBLE PRECISION SIGDEL( 12)
1302     DOUBLE PRECISION SIGA( 12)
1303     DATA NPAR, NDIM, IMQFUN / 12, 1, 1/
1304     DATA VCONST / 0.000000000000 /
1305     DATA SIGVMI / -20.50000000000 /
1306     DATA SIGVT / 41.00000000000 /
1307     DATA SIGV / 0.1585365831852
1308     +, 0.4024389982224
1309     +, 0.4756097495556
1310     +, 0.5243902206421
1311     +, 0.5975609421730
1312     +, 0.8414633870125
1313     +, 0.000000000000
1314     +, 1.000000000000
1315     +, 0.2682926654816
1316     +, 0.3170731663704
1317     +, 0.7073170542717
1318     +, 0.7560975551605
1319     +/
1320     DATA SIGDEL / 0.4878048598766E-01
1321     +, 0.4878048598766E-01
1322     +, 0.4878048598766E-01
1323     +, 0.4878048598766E-01
1324     +, 0.4878048598766E-01
1325     +, 0.4878048598766E-01
1326     +, 0.1999999994950E-05
1327     +, 0.1999999994950E-05
1328     +, 0.9756097197533E-01
1329     +, 0.9756097197533E-01
1330     +, 0.9756097197533E-01
1331     +, 0.9756097197533E-01
1332     +/
1333     DATA SIGA / 14.57433603529
1334     +, -15.93532436156
1335     +, -13.24628335221
1336     +, -14.31193855410
1337     +, -12.67339684488
1338     +, 18.19876051780
1339     +, -5.270493486725
1340     +, -5.107670990828
1341     +, -9.553262933901
1342     +, 43.34150727448
1343     +, 55.91366786432
1344     +, -29.38037318563
1345     +/
1346    
1347     v(1)= abs(x)
1348    
1349     HQUADF = 0.
1350     DO 20 J = 1, NPAR
1351     HQDJ = 0.
1352     DO 10 I = 1, NDIM
1353     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1354     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1355     10 CONTINUE
1356     HQDJ = HQDJ + SIGDEL (J) ** 2
1357     HQDJ = SQRT (HQDJ)
1358     HQUADF = HQUADF + SIGA (J) * HQDJ
1359     20 CONTINUE
1360     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1361    
1362     risy_eta2=HQUADF* 1e-4
1363    
1364     END
1365     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1366    
1367     FUNCTION risy_cog(x)
1368     DOUBLE PRECISION V( 1)
1369     INTEGER NPAR, NDIM, IMQFUN, I, J
1370     DOUBLE PRECISION HQDJ, VV, VCONST
1371     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1372     DOUBLE PRECISION SIGV( 10, 1)
1373     DOUBLE PRECISION SIGDEL( 10)
1374     DOUBLE PRECISION SIGA( 10)
1375     DATA NPAR, NDIM, IMQFUN / 10, 1, 1/
1376     DATA VCONST / 0.000000000000 /
1377     DATA SIGVMI / -20.50000000000 /
1378     DATA SIGVT / 41.00000000000 /
1379     DATA SIGV / 0.1585365831852
1380     +, 0.8414633870125
1381     +, 0.000000000000
1382     +, 1.000000000000
1383     +, 0.4634146094322
1384     +, 0.5121951103210
1385     +, 0.5609756112099
1386     +, 0.6585365533829
1387     +, 0.7073170542717
1388     +, 0.3414633870125
1389     +/
1390     DATA SIGDEL / 0.4878048598766E-01
1391     +, 0.4878048598766E-01
1392     +, 0.1999999994950E-05
1393     +, 0.1999999994950E-05
1394     +, 0.9756097197533E-01
1395     +, 0.9756097197533E-01
1396     +, 0.9756097197533E-01
1397     +, 0.9756097197533E-01
1398     +, 0.9756097197533E-01
1399     +, 0.1951219439507
1400     +/
1401     DATA SIGA / 23.73833445988
1402     +, 24.10182100013
1403     +, 1.865894323190
1404     +, 1.706006262931
1405     +, -1.075607857202
1406     +, -22.11489493403
1407     +, 1.663100707801
1408     +, 4.089852595440
1409     +, -4.314993873697
1410     +, -2.174479487744
1411     +/
1412    
1413     V(1)=abs(x)
1414    
1415     HQUADF = 0.
1416     DO 20 J = 1, NPAR
1417     HQDJ = 0.
1418     DO 10 I = 1, NDIM
1419     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1420     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1421     10 CONTINUE
1422     HQDJ = HQDJ + SIGDEL (J) ** 2
1423     HQDJ = SQRT (HQDJ)
1424     HQUADF = HQUADF + SIGA (J) * HQDJ
1425     20 CONTINUE
1426     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1427    
1428     risy_cog=HQUADF* 1e-4
1429    
1430     END
1431     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1432     FUNCTION risx_cog(x)
1433     DOUBLE PRECISION V( 1)
1434     INTEGER NPAR, NDIM, IMQFUN, I, J
1435     DOUBLE PRECISION HQDJ, VV, VCONST
1436     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1437     DOUBLE PRECISION SIGV( 15, 1)
1438     DOUBLE PRECISION SIGDEL( 15)
1439     DOUBLE PRECISION SIGA( 15)
1440     DATA NPAR, NDIM, IMQFUN / 15, 1, 1/
1441     DATA VCONST / 0.000000000000 /
1442     DATA SIGVMI / -20.50000000000 /
1443     DATA SIGVT / 41.00000000000 /
1444     DATA SIGV / 0.6097560748458E-01
1445     +, 0.8536584675312E-01
1446     +, 0.1341463327408
1447     +, 0.2317073047161
1448     +, 0.2804878056049
1449     +, 0.3780487775803
1450     +, 0.6219512224197
1451     +, 0.7195121645927
1452     +, 0.7682926654816
1453     +, 0.8658536076546
1454     +, 0.9146341085434
1455     +, 0.9390243291855
1456     +, 0.000000000000
1457     +, 1.000000000000
1458     +, 0.5121951103210
1459     +/
1460     DATA SIGDEL / 0.4878048598766E-01
1461     +, 0.4878048598766E-01
1462     +, 0.4878048598766E-01
1463     +, 0.4878048598766E-01
1464     +, 0.4878048598766E-01
1465     +, 0.4878048598766E-01
1466     +, 0.4878048598766E-01
1467     +, 0.4878048598766E-01
1468     +, 0.4878048598766E-01
1469     +, 0.4878048598766E-01
1470     +, 0.4878048598766E-01
1471     +, 0.4878048598766E-01
1472     +, 0.1999999994950E-05
1473     +, 0.1999999994950E-05
1474     +, 0.9756097197533E-01
1475     +/
1476     DATA SIGA / 31.95672945139
1477     +, -34.23286209245
1478     +, -6.298459168211
1479     +, 10.98847700545
1480     +,-0.3052213535054
1481     +, 13.10517991464
1482     +, 15.60290821679
1483     +, -1.956118448507
1484     +, 12.41453816720
1485     +, -7.354056408553
1486     +, -32.32512668778
1487     +, 30.61116178966
1488     +, 1.418505329236
1489     +, 1.583492573619
1490     +, -18.48799977042
1491     +/
1492    
1493     V(1)=abs(x)
1494    
1495     HQUADF = 0.
1496     DO 20 J = 1, NPAR
1497     HQDJ = 0.
1498     DO 10 I = 1, NDIM
1499     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1500     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1501     10 CONTINUE
1502     HQDJ = HQDJ + SIGDEL (J) ** 2
1503     HQDJ = SQRT (HQDJ)
1504     HQUADF = HQUADF + SIGA (J) * HQDJ
1505     20 CONTINUE
1506     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1507    
1508     risx_cog = HQUADF * 1e-4
1509    
1510     END
1511     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
  ViewVC Help
Powered by ViewVC 1.1.23