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