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