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