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