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