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