/[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.21 - (hide annotations) (download)
Fri Aug 31 14:56:52 2007 UTC (17 years, 3 months ago) by pam-fi
Branch: MAIN
Changes since 1.20: +298 -16 lines
new variables added to TrkTrack + other changes

1 pam-fi 1.20 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2     * this file contains all subroutines and functions
3     * that are needed for position finding algorithms:
4     *
5     * subroutine idtoc(ipfa,cpfa)
6     *
7 pam-fi 1.21 * subroutine applypfa(PFAtt,ic,ang,corr,res)
8     *
9 pam-fi 1.20 * integer function npfastrips(ic,angle)
10     *
11     * real function pfaeta(ic,angle)
12     * real function pfaetal(ic,angle)
13     * real function pfaeta2(ic,angle)
14     * real function pfaeta3(ic,angle)
15     * real function pfaeta4(ic,angle)
16     * real function cog(ncog,ic)
17     *
18     * real function fbad_cog(ncog,ic)
19     * real function fbad_eta(ic,angle)
20     *
21     * real function riseta(iview,angle)
22     * FUNCTION risxeta2(x)
23     * FUNCTION risxeta3(x)
24     * FUNCTION risxeta4(x)
25     * FUNCTION risyeta2(x)
26     * FUNCTION risy_cog(x)
27     * FUNCTION risx_cog(x)
28     *
29     * real function pfacorr(ic,angle)
30     *
31 pam-fi 1.21 * real function effectiveangle(ang,iview,bbb)
32     * real function fieldcorr(iview,bbb)
33     *
34 pam-fi 1.20 * NB - The angle is the "effective angle", which is relative
35     * to the sensor and it takes into account the magnetic field
36     *
37     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
38 pam-fi 1.10
39     subroutine idtoc(ipfa,cpfa)
40    
41     integer ipfa
42 pam-fi 1.18 character*10 cpfa
43 pam-fi 1.10
44     CPFA='COG4'
45     if(ipfa.eq.0)CPFA='ETA'
46     if(ipfa.eq.2)CPFA='ETA2'
47     if(ipfa.eq.3)CPFA='ETA3'
48     if(ipfa.eq.4)CPFA='ETA4'
49 pam-fi 1.16 if(ipfa.eq.5)CPFA='ETAL'
50 pam-fi 1.10 if(ipfa.eq.10)CPFA='COG'
51     if(ipfa.eq.11)CPFA='COG1'
52     if(ipfa.eq.12)CPFA='COG2'
53     if(ipfa.eq.13)CPFA='COG3'
54     if(ipfa.eq.14)CPFA='COG4'
55    
56     end
57 pam-fi 1.21 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
58     real function effectiveangle(ang,iview,bbb)
59    
60     include 'commontracker.f'
61    
62     effectiveangle = 0.
63    
64     if(mod(iview,2).eq.0)then
65     c =================================================
66     c X view
67     c =================================================
68     c here bbb is the y component of the m.field
69     angx = ang
70     by = bbb
71     if(iview.eq.12) angx = -1. * ang
72     if(iview.eq.12) by = -1. * bbb
73     tgtemp = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001
74    
75     elseif(mod(iview,2).eq.1)then
76     c =================================================
77     c Y view
78     c =================================================
79     c here bbb is the x component of the m.filed
80     angy = ang
81     bx = bbb
82     tgtemp = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001
83    
84     endif
85     effectiveangle = 180.*atan(tgtemp)/acos(-1.)
86    
87     return
88     end
89     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
90     real function fieldcorr(iview,bbb)
91    
92     include 'commontracker.f'
93    
94     fieldcorr = 0.
95    
96     if(mod(iview,2).eq.0)then
97    
98     c =================================================
99     c X view
100     c =================================================
101     c here bbb is the y component of the m.field
102     by = bbb
103     if(iview.eq.12) by = -1. * bbb
104     fieldcorr = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
105    
106     elseif(mod(iview,2).eq.1)then
107     c =================================================
108     c Y view
109     c =================================================
110     c here bbb is the x component of the m.filed
111     bx = bbb
112     fieldcorr = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
113    
114     endif
115    
116     return
117     end
118     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
119    
120     subroutine applypfa(PFAtt,ic,ang,corr,res)
121     *---------------------------------------------------------------
122     * this subroutine calculate the coordinate of cluster ic (in
123     * strip units), relative to the strip with the maximum signal,
124     * and its spatial resolution (in cm), applying PFAtt.
125     * ang is the effective angle, relative to the sensor
126     *---------------------------------------------------------------
127    
128     character*4 PFAtt
129     include 'commontracker.f'
130     include 'level1.f'
131    
132     corr = 0
133     res = 0
134    
135     if(ic.le.0)return
136    
137     iview = VIEW(ic)
138    
139     if(mod(iview,2).eq.0)then
140     c =================================================
141     c X view
142     c =================================================
143    
144     res = RESXAV
145    
146     if(PFAtt.eq.'COG1')then
147    
148     corr = 0
149     res = 1e-4*pitchX/sqrt(12.)!!res
150    
151     elseif(PFAtt.eq.'COG2')then
152    
153     corr = cog(2,ic)
154     res = risx_cog(abs(ang))!TEMPORANEO
155     res = res*fbad_cog(2,ic)
156    
157     elseif(PFAtt.eq.'COG3')then
158    
159     corr = cog(3,ic)
160     res = risx_cog(abs(ang))!TEMPORANEO
161     res = res*fbad_cog(3,ic)
162    
163     elseif(PFAtt.eq.'COG4')then
164    
165     corr = cog(4,ic)
166     res = risx_cog(abs(ang))!TEMPORANEO
167     res = res*fbad_cog(4,ic)
168    
169     elseif(PFAtt.eq.'ETA2')then
170    
171     corr = pfaeta2(ic,ang)
172     res = risxeta2(abs(ang))
173     res = res*fbad_cog(2,ic)
174    
175     elseif(PFAtt.eq.'ETA3')then
176    
177     corr = pfaeta3(ic,ang)
178     res = risxeta3(abs(ang))
179     res = res*fbad_cog(3,ic)
180    
181     elseif(PFAtt.eq.'ETA4')then
182    
183     corr = pfaeta4(ic,ang)
184     res = risxeta4(abs(ang))
185     res = res*fbad_cog(4,ic)
186    
187     elseif(PFAtt.eq.'ETA')then
188 pam-fi 1.10
189 pam-fi 1.21 corr = pfaeta(ic,ang)
190     c res = riseta(ic,ang)
191     res = riseta(iview,ang)
192     res = res*fbad_eta(ic,ang)
193 mocchiut 1.1
194 pam-fi 1.21 elseif(PFAtt.eq.'ETAL')then
195 mocchiut 1.1
196 pam-fi 1.21 corr = pfaetal(ic,ang)
197     res = riseta(iview,ang)
198     res = res*fbad_eta(ic,ang)
199    
200     elseif(PFAtt.eq.'COG')then
201    
202     corr = cog(0,ic)
203     res = risx_cog(abs(ang))
204     res = res*fbad_cog(0,ic)
205    
206     else
207     if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
208     endif
209    
210    
211     * ======================================
212     * temporary patch for saturated clusters
213     * ======================================
214     if( nsatstrips(ic).gt.0 )then
215     corr = cog(4,ic)
216     res = pitchX*1e-4/sqrt(12.)
217     cc cc=cog(4,ic)
218     c$$$ print*,ic,' *** ',cc
219     c$$$ print*,ic,' *** ',res
220     endif
221    
222    
223     elseif(mod(iview,2).eq.1)then
224     c =================================================
225     c Y view
226     c =================================================
227    
228     res = RESYAV
229    
230     if(PFAtt.eq.'COG1')then
231    
232     corr = 0
233     res = 1e-4*pitchY/sqrt(12.)!res
234    
235     elseif(PFAtt.eq.'COG2')then
236    
237     corr = cog(2,ic)
238     res = risy_cog(abs(ang))!TEMPORANEO
239     res = res*fbad_cog(2,ic)
240    
241     elseif(PFAtt.eq.'COG3')then
242    
243     corr = cog(3,ic)
244     res = risy_cog(abs(ang))!TEMPORANEO
245     res = res*fbad_cog(3,ic)
246    
247     elseif(PFAtt.eq.'COG4')then
248    
249     corr = cog(4,ic)
250     res = risy_cog(abs(ang))!TEMPORANEO
251     res = res*fbad_cog(4,ic)
252    
253     elseif(PFAtt.eq.'ETA2')then
254    
255     corr = pfaeta2(ic,ang)
256     res = risyeta2(abs(ang))
257     res = res*fbad_cog(2,ic)
258    
259     elseif(PFAtt.eq.'ETA3')then
260    
261     corr = pfaeta3(ic,ang)
262     res = res*fbad_cog(3,ic)
263    
264     elseif(PFAtt.eq.'ETA4')then
265    
266     corr = pfaeta4(ic,ang)
267     res = res*fbad_cog(4,ic)
268    
269     elseif(PFAtt.eq.'ETA')then
270    
271     corr = pfaeta(ic,ang)
272     c res = riseta(ic,ang)
273     res = riseta(iview,ang)
274     res = res*fbad_eta(ic,ang)
275    
276     elseif(PFAtt.eq.'ETAL')then
277    
278     corr = pfaetal(ic,ang)
279     res = riseta(iview,ang)
280     res = res*fbad_eta(ic,ang)
281    
282     elseif(PFAtt.eq.'COG')then
283    
284     corr = cog(0,ic)
285     res = risy_cog(abs(ang))
286     res = res*fbad_cog(0,ic)
287    
288     else
289     if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
290     endif
291    
292    
293     * ======================================
294     * temporary patch for saturated clusters
295     * ======================================
296     if( nsatstrips(ic).gt.0 )then
297     corr = cog(4,ic)
298     res = pitchY*1e-4/sqrt(12.)
299     cc cc=cog(4,ic)
300     c$$$ print*,ic,' *** ',cc
301     c$$$ print*,ic,' *** ',res
302     endif
303    
304     endif
305     end
306    
307     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
308 pam-fi 1.20 integer function npfastrips(ic,angle)
309 pam-fi 1.9 *--------------------------------------------------------------
310     * thid function returns the number of strips used
311     * to evaluate the position of a cluster, according to the p.f.a.
312     *--------------------------------------------------------------
313     include 'commontracker.f'
314     include 'level1.f'
315     include 'calib.f'
316    
317 pam-fi 1.20 character*4 usedPFA
318    
319 pam-fi 1.9
320    
321 pam-fi 1.20 call idtoc(pfaid,usedPFA)
322 pam-fi 1.9
323 pam-fi 1.20 npfastrips=-1
324 pam-fi 1.9
325     if(usedPFA.eq.'COG1')npfastrips=1
326     if(usedPFA.eq.'COG2')npfastrips=2
327     if(usedPFA.eq.'COG3')npfastrips=3
328     if(usedPFA.eq.'COG4')npfastrips=4
329     if(usedPFA.eq.'ETA2')npfastrips=2
330     if(usedPFA.eq.'ETA3')npfastrips=3
331     if(usedPFA.eq.'ETA4')npfastrips=4
332     * ----------------------------------------------------------------
333 pam-fi 1.20 if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
334 pam-fi 1.9 c print*,VIEW(ic),angle
335     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
336     if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
337     npfastrips=2
338     elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
339     npfastrips=3
340     elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
341     npfastrips=4
342     else
343 pam-fi 1.20 npfastrips=4 !COG4
344 pam-fi 1.9 endif
345     else !X-view
346     if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
347     npfastrips=2
348     elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
349     npfastrips=3
350     elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
351     npfastrips=4
352     else
353 pam-fi 1.20 npfastrips=4 !COG4
354 pam-fi 1.9 endif
355     endif
356     endif
357     * ----------------------------------------------------------------
358     if(usedPFA.eq.'COG')then
359    
360 pam-fi 1.20 npfastrips=0
361    
362     c$$$ iv=VIEW(ic)
363     c$$$ if(mod(iv,2).eq.1)incut=incuty
364     c$$$ if(mod(iv,2).eq.0)incut=incutx
365     c$$$ istart = INDSTART(IC)
366     c$$$ istop = TOTCLLENGTH
367     c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
368     c$$$ mu = 0
369     c$$$ do i = INDMAX(IC),istart,-1
370     c$$$ ipos = i-INDMAX(ic)
371     c$$$ cut = incut*CLSIGMA(i)
372     c$$$ if(CLSIGNAL(i).ge.cut)then
373     c$$$ mu = mu + 1
374     c$$$ print*,i,mu
375     c$$$ else
376     c$$$ goto 10
377     c$$$ endif
378     c$$$ enddo
379     c$$$ 10 continue
380     c$$$ do i = INDMAX(IC)+1,istop
381     c$$$ ipos = i-INDMAX(ic)
382     c$$$ cut = incut*CLSIGMA(i)
383     c$$$ if(CLSIGNAL(i).ge.cut)then
384     c$$$ mu = mu + 1
385     c$$$ print*,i,mu
386     c$$$ else
387     c$$$ goto 20
388     c$$$ endif
389     c$$$ enddo
390     c$$$ 20 continue
391     c$$$ npfastrips=mu
392 pam-fi 1.9
393     endif
394     * ----------------------------------------------------------------
395    
396 pam-fi 1.20 c print*,pfaid,usedPFA,angle,npfastrips
397 pam-fi 1.9
398     return
399     end
400    
401 mocchiut 1.1 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
402 pam-fi 1.4 real function pfaeta(ic,angle)
403 mocchiut 1.1 *--------------------------------------------------------------
404     * this function returns the position (in strip units)
405     * it calls:
406 pam-fi 1.4 * - pfaeta2(ic,angle)
407     * - pfaeta3(ic,angle)
408     * - pfaeta4(ic,angle)
409 mocchiut 1.1 * according to the angle
410     *--------------------------------------------------------------
411     include 'commontracker.f'
412     include 'level1.f'
413 pam-fi 1.9 include 'calib.f'
414 mocchiut 1.1
415 pam-fi 1.4 pfaeta = 0
416 mocchiut 1.1
417     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
418    
419 pam-fi 1.20 if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
420 pam-fi 1.9 pfaeta = pfaeta2(ic,angle)
421 pam-fi 1.19 cc print*,pfaeta2(ic,angle)
422 pam-fi 1.20 elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
423 pam-fi 1.9 pfaeta = pfaeta3(ic,angle)
424 pam-fi 1.20 elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
425 pam-fi 1.9 pfaeta = pfaeta4(ic,angle)
426     else
427     pfaeta = cog(4,ic)
428     endif
429    
430 mocchiut 1.1 else !X-view
431    
432 pam-fi 1.20 if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
433 pam-fi 1.4 pfaeta = pfaeta2(ic,angle)
434 pam-fi 1.20 elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
435 pam-fi 1.4 pfaeta = pfaeta3(ic,angle)
436 pam-fi 1.20 elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
437 pam-fi 1.4 pfaeta = pfaeta4(ic,angle)
438 pam-fi 1.9 else
439     pfaeta = cog(4,ic)
440     endif
441 mocchiut 1.1
442     endif
443 pam-fi 1.5
444 mocchiut 1.1 100 return
445     end
446    
447     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
448 pam-fi 1.16 real function pfaetal(ic,angle)
449     *--------------------------------------------------------------
450     * this function returns the position (in strip units)
451     * it calls:
452     * - pfaeta2(ic,angle)+pfcorr(ic,angle)
453     * - pfaeta3(ic,angle)+pfcorr(ic,angle)
454     * - pfaeta4(ic,angle)+pfcorr(ic,angle)
455     * according to the angle
456     *--------------------------------------------------------------
457     include 'commontracker.f'
458     include 'level1.f'
459     include 'calib.f'
460    
461 pam-fi 1.19 pfaetal = 0
462 pam-fi 1.16
463     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
464    
465 pam-fi 1.20 if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
466 pam-fi 1.19 pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
467     cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
468 pam-fi 1.20 elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
469 pam-fi 1.19 pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
470 pam-fi 1.20 elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
471 pam-fi 1.19 pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
472 pam-fi 1.16 else
473 pam-fi 1.19 pfaetal = cog(4,ic)
474 pam-fi 1.16 endif
475    
476     else !X-view
477    
478 pam-fi 1.20 if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
479 pam-fi 1.19 pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
480     cc print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
481 pam-fi 1.20 elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
482 pam-fi 1.19 pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
483 pam-fi 1.20 elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
484 pam-fi 1.19 pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
485 pam-fi 1.16 else
486 pam-fi 1.19 pfaetal = cog(4,ic)
487 pam-fi 1.16 endif
488    
489     endif
490    
491     100 return
492     end
493     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
494 pam-fi 1.11 c real function riseta(ic,angle)
495     real function riseta(iview,angle)
496 mocchiut 1.1 *--------------------------------------------------------------
497     * this function returns the average spatial resolution
498 pam-fi 1.4 * (in cm) for the ETA algorithm (function pfaeta(ic,angle))
499 mocchiut 1.1 * it calls:
500 pam-fi 1.13 * - risxeta2(angle)
501     * - risyeta2(angle)
502     * - risxeta3(angle)
503     * - risxeta4(angle)
504 mocchiut 1.1 * according to the angle
505     *--------------------------------------------------------------
506     include 'commontracker.f'
507     include 'level1.f'
508 pam-fi 1.9 include 'calib.f'
509 mocchiut 1.1
510 pam-fi 1.11 riseta = 0
511 mocchiut 1.1
512 pam-fi 1.11 c if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
513     if(mod(iview,2).eq.1)then !Y-view
514 mocchiut 1.1
515 pam-fi 1.9
516     if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
517 pam-fi 1.13 riseta = risyeta2(angle)
518 pam-fi 1.9 elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
519 pam-fi 1.11 riseta = risy_cog(angle) !ATTENZIONE!!
520 pam-fi 1.9 elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
521 pam-fi 1.11 riseta = risy_cog(angle) !ATTENZIONE!!
522 pam-fi 1.9 else
523 pam-fi 1.11 riseta = risy_cog(angle)
524 pam-fi 1.9 endif
525 mocchiut 1.1
526     else !X-view
527    
528 pam-fi 1.9 if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
529 pam-fi 1.13 riseta = risxeta2(angle)
530 pam-fi 1.9 elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
531 pam-fi 1.13 riseta = risxeta3(angle)
532 pam-fi 1.9 elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
533 pam-fi 1.13 riseta = risxeta4(angle)
534 pam-fi 1.9 else
535 pam-fi 1.11 riseta = risx_cog(angle)
536 pam-fi 1.9 endif
537 mocchiut 1.1
538     endif
539    
540 pam-fi 1.11
541 mocchiut 1.1 100 return
542     end
543    
544     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
545     real function fbad_eta(ic,angle)
546     *-------------------------------------------------------
547     * this function returns a factor that takes into
548     * account deterioration of the spatial resolution
549     * in the case BAD strips are included in the cluster.
550     * This factor should multiply the nominal spatial
551     * resolution.
552     * It calls the function FBAD_COG(NCOG,IC),
553     * accordingto the angle
554     *-------------------------------------------------------
555    
556     include 'commontracker.f'
557     include 'level1.f'
558 pam-fi 1.9 include 'calib.f'
559 mocchiut 1.1 fbad_eta = 0
560    
561     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
562    
563 pam-fi 1.9 if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
564     fbad_eta = fbad_cog(2,ic)
565     elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
566     fbad_eta = fbad_cog(3,ic)
567     elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
568     fbad_eta = fbad_cog(4,ic)
569     else
570     fbad_eta = fbad_cog(4,ic)
571     endif
572 mocchiut 1.1
573     else !X-view
574    
575 pam-fi 1.9 if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
576 mocchiut 1.1 fbad_eta = fbad_cog(2,ic)
577 pam-fi 1.9 elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
578 mocchiut 1.1 fbad_eta = fbad_cog(3,ic)
579 pam-fi 1.9 elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
580     fbad_eta = fbad_cog(4,ic)
581     else
582 mocchiut 1.1 fbad_eta = fbad_cog(4,ic)
583 pam-fi 1.9 endif
584 mocchiut 1.1
585     endif
586    
587     return
588     end
589    
590     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
591 pam-fi 1.4 real function pfaeta2(ic,angle) !(1)
592 mocchiut 1.1 *--------------------------------------------------------------
593     * this function returns
594     *
595     * - the position (in strip units)
596     * corrected according to the ETA2 Position Finding Algorithm.
597     * The function performs an interpolation of FETA2%ETA2
598     *
599     * - if the angle is out of range, the calibration parameters
600     * of the lowest or higher bin are used
601     *
602     *--------------------------------------------------------------
603     include 'commontracker.f'
604     include 'calib.f'
605     include 'level1.f'
606    
607     real cog2,angle
608     integer iview,lad
609    
610 pam-fi 1.9 iview = VIEW(ic)
611     lad = nld(MAXS(ic),VIEW(ic))
612     cog2 = cog(2,ic)
613 pam-fi 1.4 pfaeta2=cog2
614 mocchiut 1.1
615     * find angular bin
616     * (in futuro possiamo pensare di interpolare anche sull'angolo)
617     do iang=1,nangbin
618     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
619     iangle=iang
620     goto 98
621     endif
622     enddo
623 pam-fi 1.18 if(DEBUG.EQ.1)
624 pam-fi 1.4 $ print*,'pfaeta2 *** warning *** angle out of range: ',angle
625 mocchiut 1.1 if(angle.lt.angL(1))iang=1
626     if(angle.gt.angR(nangbin))iang=nangbin
627     98 continue !jump here if ok
628    
629    
630     c$$$* find extremes of interpolation
631     c$$$ iflag=0
632     c$$$* --------------------------------
633     c$$$ if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
634 pam-fi 1.4 c$$$c print*,'pfaeta2 *** warning *** argument out of range: ',cog2
635 mocchiut 1.1 c$$$* goto 100
636     c$$$* ----------------------------------------------
637     c$$$* non salto piu`, ma scalo di 1 o -1
638     c$$$* nel caso si tratti di un cluster
639     c$$$* in cui la strip con il segnale massimo non coincide
640     c$$$* con la strip con il rapposto s/n massimo!!!
641     c$$$* ----------------------------------------------
642     c$$$ if(cog2.lt.eta2(1,iang))then !temp
643     c$$$ cog2=cog2+1. !temp
644     c$$$ iflag=1 !temp
645     c$$$ else !temp
646     c$$$ cog2=cog2-1. !temp
647     c$$$ iflag=-1 !temp
648     c$$$ endif !temp
649     c$$$c print*,'shifted >>> ',cog2
650     c$$$ endif
651    
652     iadd=0
653     10 continue
654     if(cog2.lt.eta2(1,iang))then
655     cog2 = cog2 + 1
656     iadd = iadd + 1
657     goto 10
658     endif
659     20 continue
660     if(cog2.gt.eta2(netaval,iang))then
661     cog2 = cog2 - 1
662     iadd = iadd - 1
663     goto 20
664     endif
665    
666     * --------------------------------
667     c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
668     do i=2,netaval
669     if(eta2(i,iang).gt.cog2)then
670    
671     x1 = eta2(i-1,iang)
672     x2 = eta2(i,iang)
673     y1 = feta2(i-1,iview,lad,iang)
674     y2 = feta2(i,iview,lad,iang)
675    
676     c print*,'*****',i,view,lad,iang
677     c print*,'-----',x1,x2,y1,y2
678     goto 99
679     endif
680     enddo
681     99 continue
682    
683    
684     AA=(y2-y1)/(x2-x1)
685     BB=y1-AA*x1
686    
687 pam-fi 1.4 pfaeta2 = AA*cog2+BB
688     pfaeta2 = pfaeta2 - iadd
689 mocchiut 1.1
690     c$$$ if(iflag.eq.1)then
691 pam-fi 1.4 c$$$ pfaeta2=pfaeta2-1. !temp
692 mocchiut 1.1 c$$$ cog2=cog2-1. !temp
693     c$$$ endif
694     c$$$ if(iflag.eq.-1)then
695 pam-fi 1.4 c$$$ pfaeta2=pfaeta2+1. !temp
696 mocchiut 1.1 c$$$ cog2=cog2+1. !temp
697     c$$$ endif
698    
699 pam-fi 1.18 if(DEBUG.EQ.1)print*,'ETA2 (ic ',ic,' ang',angle,')'
700 pam-fi 1.4 $ ,cog2-iadd,' -->',pfaeta2
701 mocchiut 1.1
702    
703     100 return
704     end
705    
706     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
707 pam-fi 1.4 real function pfaeta3(ic,angle) !(1)
708 mocchiut 1.1 *--------------------------------------------------------------
709     * this function returns
710     *
711     * - the position (in strip units)
712     * corrected according to the ETA3 Position Finding Algorithm.
713     * The function performs an interpolation of FETA3%ETA3
714     *
715     * - if the angle is out of range, the calibration parameters
716     * of the lowest or higher bin are used
717     *
718     *--------------------------------------------------------------
719     include 'commontracker.f'
720     include 'calib.f'
721     include 'level1.f'
722    
723     real cog3,angle
724     integer iview,lad
725    
726    
727 pam-fi 1.9 iview = VIEW(ic)
728     lad = nld(MAXS(ic),VIEW(ic))
729     cog3 = cog(3,ic)
730 pam-fi 1.4 pfaeta3=cog3
731 mocchiut 1.1
732     * find angular bin
733     * (in futuro possiamo pensare di interpolare anche sull'angolo)
734     do iang=1,nangbin
735     c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
736     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
737     iangle=iang
738     goto 98
739     endif
740     enddo
741 pam-fi 1.18 if(DEBUG.EQ.1)
742 pam-fi 1.4 $ print*,'pfaeta3 *** warning *** angle out of range: ',angle
743 mocchiut 1.1 if(angle.lt.angL(1))iang=1
744     if(angle.gt.angR(nangbin))iang=nangbin
745     98 continue !jump here if ok
746    
747    
748     c$$$* find extremes of interpolation
749     c$$$ iflag=0
750     c$$$* --------------------------------
751     c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then
752     c$$$* ----------------------------------------------
753     c$$$* non salto piu`, ma scalo di 1 o -1
754     c$$$* nel caso si tratti di un cluster
755     c$$$* in cui la strip con il segnale massimo non coincide
756     c$$$* con la strip con il rapposto s/n massimo!!!
757     c$$$* ----------------------------------------------
758     c$$$ if(cog2.lt.eta2(1,iang))then !temp
759     c$$$ cog2=cog2+1. !temp
760     c$$$ iflag=1 !temp
761     c$$$ else !temp
762     c$$$ cog2=cog2-1. !temp
763     c$$$ iflag=-1 !temp
764     c$$$ endif !temp
765     c$$$c print*,'shifted >>> ',cog2
766     c$$$ endif
767    
768    
769     iadd=0
770     10 continue
771     if(cog3.lt.eta3(1,iang))then
772     cog3 = cog3 + 1
773     iadd = iadd + 1
774     goto 10
775     endif
776     20 continue
777     if(cog3.gt.eta3(netaval,iang))then
778     cog3 = cog3 - 1
779     iadd = iadd - 1
780     goto 20
781     endif
782    
783     * --------------------------------
784     c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
785     do i=2,netaval
786     if(eta3(i,iang).gt.cog3)then
787    
788     x1 = eta3(i-1,iang)
789     x2 = eta3(i,iang)
790     y1 = feta3(i-1,iview,lad,iang)
791     y2 = feta3(i,iview,lad,iang)
792    
793     c print*,'*****',i,view,lad,iang
794     c print*,'-----',x1,x2,y1,y2
795     goto 99
796     endif
797     enddo
798     99 continue
799    
800    
801     AA=(y2-y1)/(x2-x1)
802     BB=y1-AA*x1
803    
804 pam-fi 1.4 pfaeta3 = AA*cog3+BB
805     pfaeta3 = pfaeta3 - iadd
806 mocchiut 1.1
807     c$$$ if(iflag.eq.1)then
808 pam-fi 1.4 c$$$ pfaeta2=pfaeta2-1. !temp
809 mocchiut 1.1 c$$$ cog2=cog2-1. !temp
810     c$$$ endif
811     c$$$ if(iflag.eq.-1)then
812 pam-fi 1.4 c$$$ pfaeta2=pfaeta2+1. !temp
813 mocchiut 1.1 c$$$ cog2=cog2+1. !temp
814     c$$$ endif
815    
816 pam-fi 1.18 if(DEBUG.EQ.1)print*,'ETA3 (ic ',ic,' ang',angle,')'
817 pam-fi 1.4 $ ,cog3-iadd,' -->',pfaeta3
818 mocchiut 1.1
819     100 return
820     end
821    
822     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
823 pam-fi 1.9 real function pfaeta4(ic,angle)
824 mocchiut 1.1 *--------------------------------------------------------------
825     * this function returns
826     *
827     * - the position (in strip units)
828     * corrected according to the ETA4 Position Finding Algorithm.
829     * The function performs an interpolation of FETA3%ETA3
830     *
831     * - if the angle is out of range, the calibration parameters
832     * of the lowest or higher bin are used
833     *
834     *--------------------------------------------------------------
835     include 'commontracker.f'
836     include 'calib.f'
837     include 'level1.f'
838    
839     real cog4,angle
840     integer iview,lad
841    
842 pam-fi 1.6
843 pam-fi 1.9 iview = VIEW(ic)
844     lad = nld(MAXS(ic),VIEW(ic))
845     cog4=cog(4,ic)
846 pam-fi 1.4 pfaeta4=cog4
847 mocchiut 1.1
848     * find angular bin
849     * (in futuro possiamo pensare di interpolare anche sull'angolo)
850     do iang=1,nangbin
851     c print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
852     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
853     iangle=iang
854     goto 98
855     endif
856     enddo
857 pam-fi 1.18 if(DEBUG.EQ.1)
858 pam-fi 1.4 $ print*,'pfaeta4 *** warning *** angle out of range: ',angle
859 mocchiut 1.1 if(angle.lt.angL(1))iang=1
860     if(angle.gt.angR(nangbin))iang=nangbin
861     98 continue !jump here if ok
862    
863    
864     c$$$* find extremes of interpolation
865     c$$$ iflag=0
866     c$$$* --------------------------------
867     c$$$ if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then
868     c$$$* ----------------------------------------------
869     c$$$* non salto piu`, ma scalo di 1 o -1
870     c$$$* nel caso si tratti di un cluster
871     c$$$* in cui la strip con il segnale massimo non coincide
872     c$$$* con la strip con il rapposto s/n massimo!!!
873     c$$$* ----------------------------------------------
874     c$$$ if(cog2.lt.eta2(1,iang))then !temp
875     c$$$ cog2=cog2+1. !temp
876     c$$$ iflag=1 !temp
877     c$$$ else !temp
878     c$$$ cog2=cog2-1. !temp
879     c$$$ iflag=-1 !temp
880     c$$$ endif !temp
881     c$$$c print*,'shifted >>> ',cog2
882     c$$$ endif
883    
884    
885     iadd=0
886     10 continue
887     if(cog4.lt.eta4(1,iang))then
888     cog4 = cog4 + 1
889     iadd = iadd + 1
890     goto 10
891     endif
892     20 continue
893     if(cog4.gt.eta4(netaval,iang))then
894     cog4 = cog4 - 1
895     iadd = iadd - 1
896     goto 20
897     endif
898    
899     * --------------------------------
900     c print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
901     do i=2,netaval
902     if(eta4(i,iang).gt.cog4)then
903    
904     x1 = eta4(i-1,iang)
905     x2 = eta4(i,iang)
906     y1 = feta4(i-1,iview,lad,iang)
907     y2 = feta4(i,iview,lad,iang)
908    
909     c print*,'*****',i,view,lad,iang
910     c print*,'-----',x1,x2,y1,y2
911     goto 99
912     endif
913     enddo
914     99 continue
915    
916    
917     AA=(y2-y1)/(x2-x1)
918     BB=y1-AA*x1
919    
920 pam-fi 1.4 pfaeta4 = AA*cog4+BB
921     pfaeta4 = pfaeta4 - iadd
922 mocchiut 1.1
923     c$$$ if(iflag.eq.1)then
924 pam-fi 1.4 c$$$ pfaeta2=pfaeta2-1. !temp
925 mocchiut 1.1 c$$$ cog2=cog2-1. !temp
926     c$$$ endif
927     c$$$ if(iflag.eq.-1)then
928 pam-fi 1.4 c$$$ pfaeta2=pfaeta2+1. !temp
929 mocchiut 1.1 c$$$ cog2=cog2+1. !temp
930     c$$$ endif
931    
932 pam-fi 1.18 if(DEBUG.EQ.1)print*,'ETA4 (ic ',ic,' ang',angle,')'
933 pam-fi 1.4 $ ,cog4-iadd,' -->',pfaeta4
934 mocchiut 1.1
935     100 return
936     end
937    
938    
939    
940     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
941     real function cog(ncog,ic)
942     *-------------------------------------------------
943     * this function returns
944     *
945     * - if NCOG=0, the Center-Of-Gravity of the
946     * cluster IC, relative to MAXS(IC), according to
947     * the cluster multiplicity
948     *
949     * - if NCOG>0, the Center-Of-Gravity of the cluster IC
950     * evaluated using NCOG strips, even if they have a
951     * negative signal (according to Landi)
952     *
953     *-------------------------------------------------
954    
955    
956     include 'commontracker.f'
957     include 'calib.f'
958     include 'level1.f'
959    
960    
961    
962     if (ncog.gt.0) then
963     * ===========================
964     * ETA2 ETA3 ETA4 computation
965     * ===========================
966    
967     * --> signal of the central strip
968     sc = CLSIGNAL(INDMAX(ic)) !center
969     * signal of adjacent strips
970 pam-fi 1.14 sl1 = -9999. !left 1
971 mocchiut 1.1 if(
972     $ (INDMAX(ic)-1).ge.INDSTART(ic)
973     $ )
974     $ sl1 = CLSIGNAL(INDMAX(ic)-1)
975    
976 pam-fi 1.14 sl2 = -9999. !left 2
977 mocchiut 1.1 if(
978     $ (INDMAX(ic)-2).ge.INDSTART(ic)
979     $ )
980     $ sl2 = CLSIGNAL(INDMAX(ic)-2)
981    
982 pam-fi 1.14 sr1 = -9999. !right 1
983 mocchiut 1.1 if(
984     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
985     $ .or.
986     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
987     $ )
988     $ sr1 = CLSIGNAL(INDMAX(ic)+1)
989    
990 pam-fi 1.14 sr2 = -9999. !right 2
991 mocchiut 1.1 if(
992     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
993     $ .or.
994     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
995     $ )
996     $ sr2 = CLSIGNAL(INDMAX(ic)+2)
997    
998     COG = 0.
999    
1000 pam-fi 1.6 c print*,'## ',sl2,sl1,sc,sr1,sr2
1001    
1002 pam-fi 1.14 c ==============================================================
1003 pam-fi 1.4 if(ncog.eq.1)then
1004     COG = 0.
1005 pam-fi 1.21 if(sr1.gt.sc)cog=1.
1006     if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1007 pam-fi 1.14 c ==============================================================
1008 pam-fi 1.4 elseif(ncog.eq.2)then
1009 pam-fi 1.21 COG = 0.
1010 mocchiut 1.1 if(sl1.gt.sr1)then
1011 pam-fi 1.6 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)
1012 pam-fi 1.14 elseif(sl1.lt.sr1)then
1013     if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
1014 pam-fi 1.21 elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1015 pam-fi 1.14 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1016 pam-fi 1.21 $ .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1017 pam-fi 1.14 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1018 pam-fi 1.21 $ .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1019 pam-fi 1.14 endif
1020 pam-fi 1.15 c if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1021     c $ ,' : ',sl2,sl1,sc,sr1,sr2
1022 pam-fi 1.14 c ==============================================================
1023 mocchiut 1.1 elseif(ncog.eq.3)then
1024 pam-fi 1.21 COG = 0
1025     sss = sc
1026     if( sl1.ne.-9999. )COG = COG-sl1
1027     if( sl1.ne.-9999. )sss = sss+sl1
1028     if( sr1.ne.-9999. )COG = COG+sr1
1029     if( sr1.ne.-9999. )sss = sss+sr1
1030     if(sss.ne.0)COG=COG/sss
1031    
1032     c if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1033     c if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1034 pam-fi 1.15 c $ ,' : ',sl2,sl1,sc,sr1,sr2
1035 pam-fi 1.14 c ==============================================================
1036 mocchiut 1.1 elseif(ncog.eq.4)then
1037 pam-fi 1.21
1038     COG = 0
1039     sss = sc
1040     if( sl1.ne.-9999. )COG = COG-sl1
1041     if( sl1.ne.-9999. )sss = sss+sl1
1042     if( sr1.ne.-9999. )COG = COG+sr1
1043     if( sr1.ne.-9999. )sss = sss+sr1
1044 mocchiut 1.1 if(sl2.gt.sr2)then
1045 pam-fi 1.21 if((sl2+sss).ne.0)
1046     $ COG = (COG-2*sl2)/(sl2+sss)
1047 pam-fi 1.14 elseif(sl2.lt.sr2)then
1048 pam-fi 1.21 if((sr2+sss).ne.0)
1049     $ COG = (2*sr2+COG)/(sr2+sss)
1050     elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1051 pam-fi 1.14 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1052 pam-fi 1.21 $ .and.(sl2+sss).ne.0 )
1053     $ cog = (cog-2*sl2)/(sl2+sss)
1054 pam-fi 1.14 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1055 pam-fi 1.21 $ .and.(sr2+sss).ne.0 )
1056     $ cog = (2*sr2+cog)/(sr2+sss)
1057 mocchiut 1.1 endif
1058 pam-fi 1.21 c ==============================================================
1059     elseif(ncog.eq.5)then
1060     COG = 0
1061     sss = sc
1062     if( sl1.ne.-9999. )COG = COG-sl1
1063     if( sl1.ne.-9999. )sss = sss+sl1
1064     if( sr1.ne.-9999. )COG = COG+sr1
1065     if( sr1.ne.-9999. )sss = sss+sr1
1066     if( sl2.ne.-9999. )COG = COG-2*sl2
1067     if( sl2.ne.-9999. )sss = sss+sl2
1068     if( sr2.ne.-9999. )COG = COG+2*sr2
1069     if( sr2.ne.-9999. )sss = sss+sr2
1070     if(sss.ne.0)COG=COG/sss
1071 mocchiut 1.1 else
1072     print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1073 pam-fi 1.7 $ ,' not implemented'
1074 mocchiut 1.1 COG = 0.
1075     endif
1076    
1077     c print*,'NCOG ',ncog,ic,' @@@ ',sl1,sc,sr1,' @@@ ',cog
1078    
1079     elseif(ncog.eq.0)then
1080     * =========================
1081     * COG computation
1082     * =========================
1083    
1084     iv=VIEW(ic)
1085     if(mod(iv,2).eq.1)incut=incuty
1086     if(mod(iv,2).eq.0)incut=incutx
1087 pam-fi 1.3 istart = INDSTART(IC)
1088     istop = TOTCLLENGTH
1089 mocchiut 1.1 if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
1090 pam-fi 1.3 COG = 0
1091 pam-fi 1.9 SGN = 0.
1092 pam-fi 1.3 mu = 0
1093 pam-fi 1.9 c print*,'-------'
1094     do i = INDMAX(IC),istart,-1
1095     ipos = i-INDMAX(ic)
1096     cut = incut*CLSIGMA(i)
1097     if(CLSIGNAL(i).ge.cut)then
1098     COG = COG + ipos*CLSIGNAL(i)
1099     SGN = SGN + CLSIGNAL(i)
1100     mu = mu + 1
1101 pam-fi 1.12 c print*,ipos,CLSIGNAL(i)
1102 pam-fi 1.9 else
1103     goto 10
1104     endif
1105     enddo
1106     10 continue
1107     do i = INDMAX(IC)+1,istop
1108 pam-fi 1.3 ipos = i-INDMAX(ic)
1109 pam-fi 1.4 cut = incut*CLSIGMA(i)
1110 mocchiut 1.1 if(CLSIGNAL(i).ge.cut)then
1111     COG = COG + ipos*CLSIGNAL(i)
1112 pam-fi 1.9 SGN = SGN + CLSIGNAL(i)
1113 mocchiut 1.1 mu = mu + 1
1114 pam-fi 1.12 c print*,ipos,CLSIGNAL(i)
1115 pam-fi 1.9 else
1116     goto 20
1117 mocchiut 1.1 endif
1118     enddo
1119 pam-fi 1.9 20 continue
1120     if(SGN.le.0)then
1121 pam-fi 1.12 print*,'cog(0,ic) --> ic, dedx ',ic,SGN
1122 pam-fi 1.3 print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1123     print*,(CLSIGNAL(i),i=istart,istop)
1124 pam-fi 1.9 c print*,'cog(0,ic) --> NOT EVALUATED '
1125 pam-fi 1.3 else
1126 pam-fi 1.9 COG=COG/SGN
1127 pam-fi 1.3 endif
1128 pam-fi 1.9 c print*,'-------'
1129 mocchiut 1.1
1130     else
1131    
1132     COG=0
1133     print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1134     print*,' (NCOG must be >= 0)'
1135    
1136    
1137     endif
1138    
1139 pam-fi 1.6 c print *,'## cog ',ncog,ic,cog,'/////////////'
1140 mocchiut 1.1
1141     return
1142     end
1143    
1144     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1145 pam-fi 1.9
1146 mocchiut 1.1 real function fbad_cog(ncog,ic)
1147     *-------------------------------------------------------
1148     * this function returns a factor that takes into
1149     * account deterioration of the spatial resolution
1150     * in the case BAD strips are included in the cluster.
1151     * This factor should multiply the nominal spatial
1152     * resolution.
1153     *
1154     *-------------------------------------------------------
1155    
1156     include 'commontracker.f'
1157     include 'level1.f'
1158     include 'calib.f'
1159    
1160     if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1161 pam-fi 1.9 si = 8.4 !average good-strip noise
1162     f = 4. !average bad-strip noise: f*si
1163 mocchiut 1.1 incut=incuty
1164     else !X-view
1165 pam-fi 1.9 si = 3.9 !average good-strip noise
1166     f = 6. !average bad-strip noise: f*si
1167 mocchiut 1.1 incut=incutx
1168     endif
1169    
1170     fbad_cog = 1.
1171    
1172     if (ncog.gt.0) then
1173    
1174     * --> signal of the central strip
1175     sc = CLSIGNAL(INDMAX(ic)) !center
1176     fsc = 1
1177 pam-fi 1.9 c if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
1178     fsc = clsigma(INDMAX(ic))/si
1179 mocchiut 1.1 * --> signal of adjacent strips
1180     sl1 = 0 !left 1
1181     fsl1 = 1 !left 1
1182     if(
1183     $ (INDMAX(ic)-1).ge.INDSTART(ic)
1184     $ )then
1185     sl1 = CLSIGNAL(INDMAX(ic)-1)
1186 pam-fi 1.9 c if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f
1187     fsl1 = clsigma(INDMAX(ic)-1)/si
1188 mocchiut 1.1 endif
1189    
1190     sl2 = 0 !left 2
1191     fsl2 = 1 !left 2
1192     if(
1193     $ (INDMAX(ic)-2).ge.INDSTART(ic)
1194     $ )then
1195     sl2 = CLSIGNAL(INDMAX(ic)-2)
1196 pam-fi 1.9 c if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f
1197     fsl2 = clsigma(INDMAX(ic)-2)/si
1198 mocchiut 1.1 endif
1199     sr1 = 0 !right 1
1200     fsr1 = 1 !right 1
1201     if(
1202     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1203     $ .or.
1204     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1205     $ )then
1206     sr1 = CLSIGNAL(INDMAX(ic)+1)
1207 pam-fi 1.9 c if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f
1208     fsr1 = clsigma(INDMAX(ic)+1)/si
1209 mocchiut 1.1 endif
1210     sr2 = 0 !right 2
1211     fsr2 = 1 !right 2
1212     if(
1213     $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1214     $ .or.
1215     $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1216     $ )then
1217     sr2 = CLSIGNAL(INDMAX(ic)+2)
1218 pam-fi 1.9 c if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f
1219     fsr2 = clsigma(INDMAX(ic)+2)/si
1220 mocchiut 1.1 endif
1221    
1222    
1223    
1224     ************************************************************
1225 pam-fi 1.9 * COG2-3-4 computation
1226 mocchiut 1.1 ************************************************************
1227    
1228     c print*,sl2,sl1,sc,sr1,sr2
1229    
1230 pam-fi 1.9 vCOG = cog(ncog,ic)!0.
1231 mocchiut 1.1
1232     if(ncog.eq.2)then
1233     if(sl1.gt.sr1)then
1234 pam-fi 1.9 c COG = -sl1/(sl1+sc)
1235     fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1236     fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1237 mocchiut 1.1 elseif(sl1.le.sr1)then
1238 pam-fi 1.9 c COG = sr1/(sc+sr1)
1239     fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1240     fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1241 mocchiut 1.1 endif
1242     elseif(ncog.eq.3)then
1243 pam-fi 1.9 c COG = (sr1-sl1)/(sl1+sc+sr1)
1244 mocchiut 1.1 fbad_cog =
1245 pam-fi 1.9 $ (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1246 mocchiut 1.1 fbad_cog =
1247 pam-fi 1.9 $ fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1248 mocchiut 1.1 elseif(ncog.eq.4)then
1249     if(sl2.gt.sr2)then
1250 pam-fi 1.9 c COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1251 mocchiut 1.1 fbad_cog =
1252 pam-fi 1.9 $ (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1253     $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1254 mocchiut 1.1 fbad_cog =
1255 pam-fi 1.9 $ fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1256     $ +(-vCOG)**2+(1-vCOG)**2)
1257 mocchiut 1.1 elseif(sl2.le.sr2)then
1258 pam-fi 1.9 c COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1259 mocchiut 1.1 fbad_cog =
1260 pam-fi 1.9 $ (fsl1*(-1-vCOG)**2
1261     $ +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1262 mocchiut 1.1 fbad_cog =
1263 pam-fi 1.9 $ fbad_cog / ((-1-vCOG)**2
1264     $ +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1265 mocchiut 1.1 endif
1266     else
1267     print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1268     print*,' (NCOG must be <= 4)'
1269 pam-fi 1.9 c COG = 0.
1270 mocchiut 1.1 endif
1271    
1272     elseif(ncog.eq.0)then
1273 pam-fi 1.9 * =========================
1274     * COG computation
1275     * =========================
1276    
1277     vCOG = cog(0,ic)
1278 mocchiut 1.1
1279 pam-fi 1.9 iv = VIEW(ic)
1280     istart = INDSTART(IC)
1281     istop = TOTCLLENGTH
1282     if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1283     SGN = 0.
1284     SNU = 0.
1285     SDE = 0.
1286    
1287     do i=INDMAX(IC),istart,-1
1288     ipos = i-INDMAX(ic)
1289     cut = incut*CLSIGMA(i)
1290 mocchiut 1.1 if(CLSIGNAL(i).gt.cut)then
1291 pam-fi 1.9 fs = clsigma(i)/si
1292     SNU = SNU + fs*(ipos-vCOG)**2
1293     SDE = SDE + (ipos-vCOG)**2
1294     else
1295     goto 10
1296     endif
1297 mocchiut 1.1 enddo
1298 pam-fi 1.9 10 continue
1299     do i=INDMAX(IC)+1,istop
1300     ipos = i-INDMAX(ic)
1301     cut = incut*CLSIGMA(i)
1302 mocchiut 1.1 if(CLSIGNAL(i).gt.cut)then
1303 pam-fi 1.9 fs = clsigma(i)/si
1304     SNU = SNU + fs*(ipos-vCOG)**2
1305     SDE = SDE + (ipos-vCOG)**2
1306     else
1307     goto 20
1308 mocchiut 1.1 endif
1309     enddo
1310 pam-fi 1.9 20 continue
1311     if(SDE.ne.0)then
1312     FBAD_COG=SNU/SDE
1313     else
1314    
1315     endif
1316 mocchiut 1.1
1317     else
1318    
1319     FBAD_COG=0
1320     print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1321     print*,' (NCOG must be >= 0)'
1322    
1323    
1324     endif
1325    
1326    
1327     fbad_cog = sqrt(fbad_cog)
1328    
1329     return
1330     end
1331    
1332    
1333 pam-fi 1.20 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1334     c$$$ real function fbad_cog0(ncog,ic)
1335     c$$$*-------------------------------------------------------
1336     c$$$* this function returns a factor that takes into
1337     c$$$* account deterioration of the spatial resolution
1338     c$$$* in the case BAD strips are included in the cluster.
1339     c$$$* This factor should multiply the nominal spatial
1340     c$$$* resolution.
1341     c$$$*
1342     c$$$* NB!!!
1343     c$$$* (this is the old version. It consider only the two
1344     c$$$* strips with the greatest signal. The new one is
1345     c$$$* fbad_cog(ncog,ic) )
1346     c$$$*
1347     c$$$*-------------------------------------------------------
1348     c$$$
1349     c$$$ include 'commontracker.f'
1350     c$$$ include 'level1.f'
1351     c$$$ include 'calib.f'
1352     c$$$
1353     c$$$* --> signal of the central strip
1354     c$$$ sc = CLSIGNAL(INDMAX(ic)) !center
1355     c$$$
1356     c$$$* signal of adjacent strips
1357     c$$$* --> left
1358     c$$$ sl1 = 0 !left 1
1359     c$$$ if(
1360     c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic)
1361     c$$$ $ )
1362     c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
1363     c$$$
1364     c$$$ sl2 = 0 !left 2
1365     c$$$ if(
1366     c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic)
1367     c$$$ $ )
1368     c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
1369     c$$$
1370     c$$$* --> right
1371     c$$$ sr1 = 0 !right 1
1372     c$$$ if(
1373     c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1374     c$$$ $ .or.
1375     c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1376     c$$$ $ )
1377     c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
1378     c$$$
1379     c$$$ sr2 = 0 !right 2
1380     c$$$ if(
1381     c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1382     c$$$ $ .or.
1383     c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1384     c$$$ $ )
1385     c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
1386     c$$$
1387     c$$$
1388     c$$$ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1389     c$$$ f = 4.
1390     c$$$ si = 8.4
1391     c$$$ else !X-view
1392     c$$$ f = 6.
1393     c$$$ si = 3.9
1394     c$$$ endif
1395     c$$$
1396     c$$$ fbad_cog = 1.
1397     c$$$ f0 = 1
1398     c$$$ f1 = 1
1399     c$$$ f2 = 1
1400     c$$$ f3 = 1
1401     c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then
1402     c$$$
1403     c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
1404     c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
1405     c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
1406     c$$$
1407     c$$$ if(ncog.eq.2.and.sl1.ne.0)then
1408     c$$$ fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
1409     c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
1410     c$$$ fbad_cog = 1.
1411     c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
1412     c$$$ fbad_cog = 1.
1413     c$$$ else
1414     c$$$ fbad_cog = 1.
1415     c$$$ endif
1416     c$$$
1417     c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then
1418     c$$$
1419     c$$$
1420     c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
1421     c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
1422     c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
1423     c$$$
1424     c$$$ if(ncog.eq.2.and.sr1.ne.0)then
1425     c$$$ fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
1426     c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
1427     c$$$ fbad_cog = 1.
1428     c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
1429     c$$$ fbad_cog = 1.
1430     c$$$ else
1431     c$$$ fbad_cog = 1.
1432     c$$$ endif
1433     c$$$
1434     c$$$ endif
1435     c$$$
1436     c$$$ fbad_cog0 = sqrt(fbad_cog)
1437     c$$$
1438     c$$$ return
1439     c$$$ end
1440     c$$$
1441     c$$$
1442     c$$$
1443 mocchiut 1.1
1444     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1445    
1446 pam-fi 1.13 FUNCTION risxeta2(x)
1447 mocchiut 1.1
1448     DOUBLE PRECISION V( 1)
1449     INTEGER NPAR, NDIM, IMQFUN, I, J
1450     DOUBLE PRECISION HQDJ, VV, VCONST
1451     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1452     DOUBLE PRECISION SIGV( 18, 1)
1453     DOUBLE PRECISION SIGDEL( 18)
1454     DOUBLE PRECISION SIGA( 18)
1455     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1456     DATA VCONST / 0.000000000000 /
1457     DATA SIGVMI / -20.50000000000 /
1458     DATA SIGVT / 41.00000000000 /
1459     DATA SIGV / 0.6097560748458E-01
1460     +, 0.1097560971975
1461     +, 0.1341463327408
1462     +, 0.1829268187284
1463     +, 0.2317073047161
1464     +, 0.4268292486668
1465     +, 0.4756097495556
1466     +, 0.4999999701977
1467     +, 0.5243902206421
1468     +, 0.5731707215309
1469     +, 0.7682926654816
1470     +, 0.8170731663704
1471     +, 0.8658536076546
1472     +, 0.8902438879013
1473     +, 0.9390243291855
1474     +, 0.000000000000
1475     +, 1.000000000000
1476     +, 0.3658536374569
1477     +/
1478     DATA SIGDEL / 0.4878048598766E-01
1479     +, 0.4878048598766E-01
1480     +, 0.4878048598766E-01
1481     +, 0.4878048598766E-01
1482     +, 0.4878048598766E-01
1483     +, 0.4878048598766E-01
1484     +, 0.4878048598766E-01
1485     +, 0.4878048598766E-01
1486     +, 0.4878048598766E-01
1487     +, 0.4878048598766E-01
1488     +, 0.4878048598766E-01
1489     +, 0.4878048598766E-01
1490     +, 0.4878048598766E-01
1491     +, 0.4878048598766E-01
1492     +, 0.4878048598766E-01
1493     +, 0.1999999994950E-05
1494     +, 0.1999999994950E-05
1495     +, 0.9756097197533E-01
1496     +/
1497     DATA SIGA / 51.65899502118
1498     +, -150.4733247841
1499     +, 143.0468613786
1500     +, -16.56096738997
1501     +, 5.149319798083
1502     +, 21.57149712673
1503     +, -39.46652322782
1504     +, 47.13181632948
1505     +, -32.93197883680
1506     +, 16.38645317092
1507     +, 1.453688482992
1508     +, -10.00547244421
1509     +, 131.3517670587
1510     +, -140.6351538257
1511     +, 49.05515749582
1512     +, -23.00028974788
1513     +, -22.58470403729
1514     +, -3.824682486418
1515     +/
1516    
1517     V(1)= abs(x)
1518 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1519 mocchiut 1.1
1520     HQUADF = 0.
1521     DO 20 J = 1, NPAR
1522     HQDJ = 0.
1523     DO 10 I = 1, NDIM
1524     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1525     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1526     10 CONTINUE
1527     HQDJ = HQDJ + SIGDEL (J) ** 2
1528     HQDJ = SQRT (HQDJ)
1529     HQUADF = HQUADF + SIGA (J) * HQDJ
1530     20 CONTINUE
1531     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1532    
1533 pam-fi 1.13 risxeta2=HQUADF* 1e-4
1534 mocchiut 1.1
1535     END
1536    
1537     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1538 pam-fi 1.13 FUNCTION risxeta3(x)
1539 mocchiut 1.1 DOUBLE PRECISION V( 1)
1540     INTEGER NPAR, NDIM, IMQFUN, I, J
1541     DOUBLE PRECISION HQDJ, VV, VCONST
1542     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1543     DOUBLE PRECISION SIGV( 18, 1)
1544     DOUBLE PRECISION SIGDEL( 18)
1545     DOUBLE PRECISION SIGA( 18)
1546     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1547     DATA VCONST / 0.000000000000 /
1548     DATA SIGVMI / -20.50000000000 /
1549     DATA SIGVT / 41.00000000000 /
1550     DATA SIGV / 0.6097560748458E-01
1551     +, 0.1097560971975
1552     +, 0.1341463327408
1553     +, 0.1829268187284
1554     +, 0.2317073047161
1555     +, 0.4756097495556
1556     +, 0.4999999701977
1557     +, 0.5243902206421
1558     +, 0.7682926654816
1559     +, 0.8170731663704
1560     +, 0.8658536076546
1561     +, 0.8902438879013
1562     +, 0.9390243291855
1563     +, 0.000000000000
1564     +, 1.000000000000
1565     +, 0.3658536374569
1566     +, 0.4146341383457
1567     +, 0.6097560524940
1568     +/
1569     DATA SIGDEL / 0.4878048598766E-01
1570     +, 0.4878048598766E-01
1571     +, 0.4878048598766E-01
1572     +, 0.4878048598766E-01
1573     +, 0.4878048598766E-01
1574     +, 0.4878048598766E-01
1575     +, 0.4878048598766E-01
1576     +, 0.4878048598766E-01
1577     +, 0.4878048598766E-01
1578     +, 0.4878048598766E-01
1579     +, 0.4878048598766E-01
1580     +, 0.4878048598766E-01
1581     +, 0.4878048598766E-01
1582     +, 0.1999999994950E-05
1583     +, 0.1999999994950E-05
1584     +, 0.9756097197533E-01
1585     +, 0.9756097197533E-01
1586     +, 0.9756097197533E-01
1587     +/
1588     DATA SIGA / 55.18284054458
1589     +, -160.3358431242
1590     +, 144.6939185763
1591     +, -20.45200854118
1592     +, 5.223570087108
1593     +,-0.4171476953945
1594     +, -27.67911907462
1595     +, 17.70327157495
1596     +, -1.867165491707
1597     +, -8.884458169181
1598     +, 124.3526608791
1599     +, -143.3309398345
1600     +, 50.80345027122
1601     +, -16.44454904415
1602     +, -15.73785568450
1603     +, -22.71810502561
1604     +, 36.86170101430
1605     +, 2.437918198452
1606     +/
1607    
1608     V(1) = abs(x)
1609 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1610    
1611 mocchiut 1.1 HQUADF = 0.
1612     DO 20 J = 1, NPAR
1613     HQDJ = 0.
1614     DO 10 I = 1, NDIM
1615     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1616     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1617     10 CONTINUE
1618     HQDJ = HQDJ + SIGDEL (J) ** 2
1619     HQDJ = SQRT (HQDJ)
1620     HQUADF = HQUADF + SIGA (J) * HQDJ
1621     20 CONTINUE
1622     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1623    
1624 pam-fi 1.13 risxeta3 = HQUADF* 1e-4
1625 mocchiut 1.1
1626     END
1627     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1628 pam-fi 1.13 FUNCTION risxeta4(x)
1629 mocchiut 1.1 DOUBLE PRECISION V( 1)
1630     INTEGER NPAR, NDIM, IMQFUN, I, J
1631     DOUBLE PRECISION HQDJ, VV, VCONST
1632     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1633     DOUBLE PRECISION SIGV( 18, 1)
1634     DOUBLE PRECISION SIGDEL( 18)
1635     DOUBLE PRECISION SIGA( 18)
1636     DATA NPAR, NDIM, IMQFUN / 18, 1, 1/
1637     DATA VCONST / 0.000000000000 /
1638     DATA SIGVMI / -20.50000000000 /
1639     DATA SIGVT / 41.00000000000 /
1640     DATA SIGV / 0.3658536449075E-01
1641     +, 0.6097560748458E-01
1642     +, 0.1097560971975
1643     +, 0.1341463327408
1644     +, 0.4756097495556
1645     +, 0.5243902206421
1646     +, 0.8658536076546
1647     +, 0.8902438879013
1648     +, 0.9390243291855
1649     +, 0.9634146094322
1650     +, 0.000000000000
1651     +, 1.000000000000
1652     +, 0.3658536374569
1653     +, 0.4146341383457
1654     +, 0.6097560524940
1655     +, 0.6585365533829
1656     +, 0.7560975551605
1657     +, 0.2439024299383
1658     +/
1659     DATA SIGDEL / 0.4878048598766E-01
1660     +, 0.4878048598766E-01
1661     +, 0.4878048598766E-01
1662     +, 0.4878048598766E-01
1663     +, 0.4878048598766E-01
1664     +, 0.4878048598766E-01
1665     +, 0.4878048598766E-01
1666     +, 0.4878048598766E-01
1667     +, 0.4878048598766E-01
1668     +, 0.4878048598766E-01
1669     +, 0.1999999994950E-05
1670     +, 0.1999999994950E-05
1671     +, 0.9756097197533E-01
1672     +, 0.9756097197533E-01
1673     +, 0.9756097197533E-01
1674     +, 0.9756097197533E-01
1675     +, 0.9756097197533E-01
1676     +, 0.1951219439507
1677     +/
1678     DATA SIGA / -43.61551887895
1679     +, 57.88466995373
1680     +, -92.04113299504
1681     +, 74.08166649890
1682     +, -9.768686062558
1683     +, -4.304496875334
1684     +, 72.62237333937
1685     +, -91.21920840618
1686     +, 56.75519978630
1687     +, -43.21115751243
1688     +, 12.79984505413
1689     +, 12.10074868595
1690     +, -6.238587250860
1691     +, 23.43447356326
1692     +, 17.98221401495
1693     +, -7.980332610975
1694     +, -3.426733307051
1695     +, -8.683439558751
1696     +/
1697    
1698     V(1)=abs(x)
1699 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1700 mocchiut 1.1
1701     HQUADF = 0.
1702     DO 20 J = 1, NPAR
1703     HQDJ = 0.
1704     DO 10 I = 1, NDIM
1705     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1706     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1707     10 CONTINUE
1708     HQDJ = HQDJ + SIGDEL (J) ** 2
1709     HQDJ = SQRT (HQDJ)
1710     HQUADF = HQUADF + SIGA (J) * HQDJ
1711     20 CONTINUE
1712     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1713    
1714 pam-fi 1.13 risxeta4=HQUADF* 1e-4
1715 mocchiut 1.1
1716     END
1717     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1718 pam-fi 1.13 FUNCTION risyeta2(x)
1719 mocchiut 1.1 DOUBLE PRECISION V( 1)
1720     INTEGER NPAR, NDIM, IMQFUN, I, J
1721     DOUBLE PRECISION HQDJ, VV, VCONST
1722     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1723     DOUBLE PRECISION SIGV( 12, 1)
1724     DOUBLE PRECISION SIGDEL( 12)
1725     DOUBLE PRECISION SIGA( 12)
1726     DATA NPAR, NDIM, IMQFUN / 12, 1, 1/
1727     DATA VCONST / 0.000000000000 /
1728     DATA SIGVMI / -20.50000000000 /
1729     DATA SIGVT / 41.00000000000 /
1730     DATA SIGV / 0.1585365831852
1731     +, 0.4024389982224
1732     +, 0.4756097495556
1733     +, 0.5243902206421
1734     +, 0.5975609421730
1735     +, 0.8414633870125
1736     +, 0.000000000000
1737     +, 1.000000000000
1738     +, 0.2682926654816
1739     +, 0.3170731663704
1740     +, 0.7073170542717
1741     +, 0.7560975551605
1742     +/
1743     DATA SIGDEL / 0.4878048598766E-01
1744     +, 0.4878048598766E-01
1745     +, 0.4878048598766E-01
1746     +, 0.4878048598766E-01
1747     +, 0.4878048598766E-01
1748     +, 0.4878048598766E-01
1749     +, 0.1999999994950E-05
1750     +, 0.1999999994950E-05
1751     +, 0.9756097197533E-01
1752     +, 0.9756097197533E-01
1753     +, 0.9756097197533E-01
1754     +, 0.9756097197533E-01
1755     +/
1756     DATA SIGA / 14.57433603529
1757     +, -15.93532436156
1758     +, -13.24628335221
1759     +, -14.31193855410
1760     +, -12.67339684488
1761     +, 18.19876051780
1762     +, -5.270493486725
1763     +, -5.107670990828
1764     +, -9.553262933901
1765     +, 43.34150727448
1766     +, 55.91366786432
1767     +, -29.38037318563
1768     +/
1769    
1770     v(1)= abs(x)
1771 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1772 mocchiut 1.1
1773     HQUADF = 0.
1774     DO 20 J = 1, NPAR
1775     HQDJ = 0.
1776     DO 10 I = 1, NDIM
1777     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1778     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1779     10 CONTINUE
1780     HQDJ = HQDJ + SIGDEL (J) ** 2
1781     HQDJ = SQRT (HQDJ)
1782     HQUADF = HQUADF + SIGA (J) * HQDJ
1783     20 CONTINUE
1784     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1785    
1786 pam-fi 1.13 risyeta2=HQUADF* 1e-4
1787 mocchiut 1.1
1788     END
1789     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1790    
1791     FUNCTION risy_cog(x)
1792     DOUBLE PRECISION V( 1)
1793     INTEGER NPAR, NDIM, IMQFUN, I, J
1794     DOUBLE PRECISION HQDJ, VV, VCONST
1795     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1796     DOUBLE PRECISION SIGV( 10, 1)
1797     DOUBLE PRECISION SIGDEL( 10)
1798     DOUBLE PRECISION SIGA( 10)
1799     DATA NPAR, NDIM, IMQFUN / 10, 1, 1/
1800     DATA VCONST / 0.000000000000 /
1801     DATA SIGVMI / -20.50000000000 /
1802     DATA SIGVT / 41.00000000000 /
1803     DATA SIGV / 0.1585365831852
1804     +, 0.8414633870125
1805     +, 0.000000000000
1806     +, 1.000000000000
1807     +, 0.4634146094322
1808     +, 0.5121951103210
1809     +, 0.5609756112099
1810     +, 0.6585365533829
1811     +, 0.7073170542717
1812     +, 0.3414633870125
1813     +/
1814     DATA SIGDEL / 0.4878048598766E-01
1815     +, 0.4878048598766E-01
1816     +, 0.1999999994950E-05
1817     +, 0.1999999994950E-05
1818     +, 0.9756097197533E-01
1819     +, 0.9756097197533E-01
1820     +, 0.9756097197533E-01
1821     +, 0.9756097197533E-01
1822     +, 0.9756097197533E-01
1823     +, 0.1951219439507
1824     +/
1825     DATA SIGA / 23.73833445988
1826     +, 24.10182100013
1827     +, 1.865894323190
1828     +, 1.706006262931
1829     +, -1.075607857202
1830     +, -22.11489493403
1831     +, 1.663100707801
1832     +, 4.089852595440
1833     +, -4.314993873697
1834     +, -2.174479487744
1835     +/
1836    
1837     V(1)=abs(x)
1838 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1839 mocchiut 1.1
1840     HQUADF = 0.
1841     DO 20 J = 1, NPAR
1842     HQDJ = 0.
1843     DO 10 I = 1, NDIM
1844     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1845     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1846     10 CONTINUE
1847     HQDJ = HQDJ + SIGDEL (J) ** 2
1848     HQDJ = SQRT (HQDJ)
1849     HQUADF = HQUADF + SIGA (J) * HQDJ
1850     20 CONTINUE
1851     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1852    
1853     risy_cog=HQUADF* 1e-4
1854    
1855     END
1856     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1857     FUNCTION risx_cog(x)
1858     DOUBLE PRECISION V( 1)
1859     INTEGER NPAR, NDIM, IMQFUN, I, J
1860     DOUBLE PRECISION HQDJ, VV, VCONST
1861     DOUBLE PRECISION SIGVMI( 1), SIGVT( 1)
1862     DOUBLE PRECISION SIGV( 15, 1)
1863     DOUBLE PRECISION SIGDEL( 15)
1864     DOUBLE PRECISION SIGA( 15)
1865     DATA NPAR, NDIM, IMQFUN / 15, 1, 1/
1866     DATA VCONST / 0.000000000000 /
1867     DATA SIGVMI / -20.50000000000 /
1868     DATA SIGVT / 41.00000000000 /
1869     DATA SIGV / 0.6097560748458E-01
1870     +, 0.8536584675312E-01
1871     +, 0.1341463327408
1872     +, 0.2317073047161
1873     +, 0.2804878056049
1874     +, 0.3780487775803
1875     +, 0.6219512224197
1876     +, 0.7195121645927
1877     +, 0.7682926654816
1878     +, 0.8658536076546
1879     +, 0.9146341085434
1880     +, 0.9390243291855
1881     +, 0.000000000000
1882     +, 1.000000000000
1883     +, 0.5121951103210
1884     +/
1885     DATA SIGDEL / 0.4878048598766E-01
1886     +, 0.4878048598766E-01
1887     +, 0.4878048598766E-01
1888     +, 0.4878048598766E-01
1889     +, 0.4878048598766E-01
1890     +, 0.4878048598766E-01
1891     +, 0.4878048598766E-01
1892     +, 0.4878048598766E-01
1893     +, 0.4878048598766E-01
1894     +, 0.4878048598766E-01
1895     +, 0.4878048598766E-01
1896     +, 0.4878048598766E-01
1897     +, 0.1999999994950E-05
1898     +, 0.1999999994950E-05
1899     +, 0.9756097197533E-01
1900     +/
1901     DATA SIGA / 31.95672945139
1902     +, -34.23286209245
1903     +, -6.298459168211
1904     +, 10.98847700545
1905     +,-0.3052213535054
1906     +, 13.10517991464
1907     +, 15.60290821679
1908     +, -1.956118448507
1909     +, 12.41453816720
1910     +, -7.354056408553
1911     +, -32.32512668778
1912     +, 30.61116178966
1913     +, 1.418505329236
1914     +, 1.583492573619
1915     +, -18.48799977042
1916     +/
1917    
1918     V(1)=abs(x)
1919 pam-fi 1.9 if(V(1).gt.20.)V(1)=20.
1920 mocchiut 1.1
1921     HQUADF = 0.
1922     DO 20 J = 1, NPAR
1923     HQDJ = 0.
1924     DO 10 I = 1, NDIM
1925     VV = (V (I) - SIGVMI (I)) / SIGVT (I)
1926     HQDJ = HQDJ + (VV - SIGV (J, I)) ** 2
1927     10 CONTINUE
1928     HQDJ = HQDJ + SIGDEL (J) ** 2
1929     HQDJ = SQRT (HQDJ)
1930     HQUADF = HQUADF + SIGA (J) * HQDJ
1931     20 CONTINUE
1932     IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1933    
1934     risx_cog = HQUADF * 1e-4
1935    
1936     END
1937 pam-fi 1.16
1938    
1939     *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1940 pam-fi 1.19 real function pfacorr(ic,angle)
1941 pam-fi 1.16 *--------------------------------------------------------------
1942     * this function returns the landi correction for this cluster
1943     *--------------------------------------------------------------
1944     include 'commontracker.f'
1945     include 'calib.f'
1946     include 'level1.f'
1947    
1948     real angle
1949     integer iview,lad
1950    
1951     iview = VIEW(ic)
1952     lad = nld(MAXS(ic),VIEW(ic))
1953    
1954     * find angular bin
1955     * (in futuro possiamo pensare di interpolare anche sull'angolo)
1956     do iang=1,nangbin
1957     if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1958     iangle=iang
1959     goto 98
1960     endif
1961     enddo
1962 pam-fi 1.18 if(DEBUG.eq.1)
1963 pam-fi 1.16 $ print*,'pfacorr *** warning *** angle out of range: ',angle
1964     if(angle.lt.angL(1))iang=1
1965     if(angle.gt.angR(nangbin))iang=nangbin
1966     98 continue !jump here if ok
1967    
1968     pfacorr = fcorr(iview,lad,iang)
1969    
1970 pam-fi 1.18 if(DEBUG.eq.1)print*,'CORR (ic ',ic,' ang',angle,') -->',pfacorr
1971 pam-fi 1.16
1972    
1973     100 return
1974     end

  ViewVC Help
Powered by ViewVC 1.1.23