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