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