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