/[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.27 - (hide annotations) (download)
Tue Mar 12 11:02:02 2013 UTC (12 years, 7 months ago) by pam-fi
Branch: MAIN
Changes since 1.26: +52 -3 lines 
Add digital position finding for saturated clusters.

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