/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/functionspfa.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/functionspfa.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by pam-fi, Tue May 30 16:30:37 2006 UTC revision 1.20 by pam-fi, Mon Aug 27 12:57:15 2007 UTC
# Line 1  Line 1 
1  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
3  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms:
4  *      *          
5    *     subroutine idtoc(ipfa,cpfa)
6    *
7    *     integer function npfastrips(ic,angle)
8    *
9    *     real function pfaeta(ic,angle)
10    *     real function pfaetal(ic,angle)
11    *     real function pfaeta2(ic,angle)
12    *     real function pfaeta3(ic,angle)
13    *     real function pfaeta4(ic,angle)
14    *     real function cog(ncog,ic)
15    *
16    *     real function fbad_cog(ncog,ic)
17    *     real function fbad_eta(ic,angle)
18    *
19    *     real function riseta(iview,angle)
20    *     FUNCTION risxeta2(x)
21    *     FUNCTION risxeta3(x)
22    *     FUNCTION risxeta4(x)
23    *     FUNCTION risyeta2(x)
24    *     FUNCTION risy_cog(x)
25    *     FUNCTION risx_cog(x)
26    *
27    *     real function pfacorr(ic,angle)
28    *
29    *     NB - The angle is the "effective angle", which is relative
30    *          to the sensor and it takes into account the magnetic field
31  *  *
32  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
33    
34          subroutine idtoc(ipfa,cpfa)
35          
36          integer ipfa
37          character*10 cpfa
38    
39          CPFA='COG4'
40          if(ipfa.eq.0)CPFA='ETA'
41          if(ipfa.eq.2)CPFA='ETA2'
42          if(ipfa.eq.3)CPFA='ETA3'
43          if(ipfa.eq.4)CPFA='ETA4'
44          if(ipfa.eq.5)CPFA='ETAL'
45          if(ipfa.eq.10)CPFA='COG'
46          if(ipfa.eq.11)CPFA='COG1'
47          if(ipfa.eq.12)CPFA='COG2'
48          if(ipfa.eq.13)CPFA='COG3'
49          if(ipfa.eq.14)CPFA='COG4'
50          
51          end
52    
53    
54    
55          integer function npfastrips(ic,angle)
56    *--------------------------------------------------------------
57    *     thid function returns the number of strips used
58    *     to evaluate the position of a cluster, according to the p.f.a.
59    *--------------------------------------------------------------
60          include 'commontracker.f'
61          include 'level1.f'
62          include 'calib.f'
63    
64          character*4 usedPFA
65          
66    
67    
68          call idtoc(pfaid,usedPFA)
69    
70          npfastrips=-1
71    
72          if(usedPFA.eq.'COG1')npfastrips=1
73          if(usedPFA.eq.'COG2')npfastrips=2
74          if(usedPFA.eq.'COG3')npfastrips=3
75          if(usedPFA.eq.'COG4')npfastrips=4
76          if(usedPFA.eq.'ETA2')npfastrips=2
77          if(usedPFA.eq.'ETA3')npfastrips=3
78          if(usedPFA.eq.'ETA4')npfastrips=4
79    *     ----------------------------------------------------------------
80          if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
81    c         print*,VIEW(ic),angle
82             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
83                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
84                   npfastrips=2
85                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
86                   npfastrips=3
87                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
88                   npfastrips=4
89                else
90                   npfastrips=4     !COG4
91                endif                        
92             else                   !X-view
93                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
94                   npfastrips=2
95                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
96                   npfastrips=3
97                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
98                   npfastrips=4
99                else
100                   npfastrips=4     !COG4
101                endif                        
102             endif
103          endif
104    *     ----------------------------------------------------------------
105          if(usedPFA.eq.'COG')then
106    
107             npfastrips=0
108    
109    c$$$         iv=VIEW(ic)
110    c$$$         if(mod(iv,2).eq.1)incut=incuty
111    c$$$         if(mod(iv,2).eq.0)incut=incutx
112    c$$$         istart = INDSTART(IC)
113    c$$$         istop  = TOTCLLENGTH
114    c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
115    c$$$         mu  = 0
116    c$$$         do i = INDMAX(IC),istart,-1
117    c$$$            ipos = i-INDMAX(ic)
118    c$$$            cut  = incut*CLSIGMA(i)
119    c$$$            if(CLSIGNAL(i).ge.cut)then
120    c$$$               mu = mu + 1
121    c$$$               print*,i,mu
122    c$$$            else
123    c$$$               goto 10
124    c$$$            endif
125    c$$$         enddo
126    c$$$ 10      continue
127    c$$$         do i = INDMAX(IC)+1,istop
128    c$$$            ipos = i-INDMAX(ic)
129    c$$$            cut  = incut*CLSIGMA(i)
130    c$$$            if(CLSIGNAL(i).ge.cut)then
131    c$$$               mu = mu + 1
132    c$$$               print*,i,mu
133    c$$$            else
134    c$$$               goto 20
135    c$$$            endif
136    c$$$         enddo
137    c$$$ 20      continue
138    c$$$         npfastrips=mu
139    
140          endif
141    *     ----------------------------------------------------------------
142    
143    c      print*,pfaid,usedPFA,angle,npfastrips
144    
145          return
146          end
147    
148  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
149        real function pfa_eta(ic,angle)        real function pfaeta(ic,angle)
150  *--------------------------------------------------------------  *--------------------------------------------------------------
151  *     this function returns the position (in strip units)  *     this function returns the position (in strip units)
152  *     it calls:  *     it calls:
153  *     - pfa_eta2(ic,angle)  *     - pfaeta2(ic,angle)
154  *     - pfa_eta3(ic,angle)  *     - pfaeta3(ic,angle)
155  *     - pfa_eta4(ic,angle)  *     - pfaeta4(ic,angle)
156  *     according to the angle  *     according to the angle
157  *--------------------------------------------------------------  *--------------------------------------------------------------
158        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
159        include 'level1.f'        include 'level1.f'
160          include 'calib.f'
161                
162        pfa_eta = 0        pfaeta = 0
163    
164        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
165                
166           pfa_eta = pfa_eta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
167                  pfaeta = pfaeta2(ic,angle)
168    cc            print*,pfaeta2(ic,angle)
169             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
170                pfaeta = pfaeta3(ic,angle)
171             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
172                pfaeta = pfaeta4(ic,angle)
173             else
174                pfaeta = cog(4,ic)
175             endif            
176    
177        else                      !X-view        else                      !X-view
178    
179           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
180              pfa_eta = pfa_eta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
181           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
182              pfa_eta = pfa_eta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
183           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
184              pfa_eta = pfa_eta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
185           endif           else
186                pfaeta = cog(4,ic)
187             endif            
188                            
189        endif        endif
190          
   
191   100  return   100  return
192        end        end
193    
194  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
195        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
196  *--------------------------------------------------------------  *--------------------------------------------------------------
197  *     this function returns the average spatial resolution  *     this function returns the position (in strip units)
 *     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))  
198  *     it calls:  *     it calls:
199  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
200  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
201  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
202  *     according to the angle  *     according to the angle
203  *--------------------------------------------------------------  *--------------------------------------------------------------
204        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
205        include 'level1.f'        include 'level1.f'
206          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
207                
208        ris_eta = 0        pfaetal = 0
209    
210        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
211                
212           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
213           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
214    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
215             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
216                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
217             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
218                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
219             else
220                pfaetal = cog(4,ic)
221             endif            
222    
223        else                      !X-view        else                      !X-view
224    
225           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
226              ris_eta = risx_eta2(angle)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
227           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then  cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
228              ris_eta = risx_eta3(angle)           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
229           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then              pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
230              ris_eta = risx_eta4(angle)           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
231           elseif(abs(angle).gt.21.)then              pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
232              ris_eta = risx_eta4(21.)           else
233           endif              pfaetal = cog(4,ic)
234             endif            
235                            
236        endif        endif
237          
238     100  return
239          end
240    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
241    c      real function riseta(ic,angle)
242          real function riseta(iview,angle)
243    *--------------------------------------------------------------
244    *     this function returns the average spatial resolution
245    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
246    *     it calls:
247    *     - risxeta2(angle)
248    *     - risyeta2(angle)
249    *     - risxeta3(angle)
250    *     - risxeta4(angle)
251    *     according to the angle
252    *--------------------------------------------------------------
253          include 'commontracker.f'
254          include 'level1.f'
255          include 'calib.f'
256    
257          riseta = 0
258    
259    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
260          if(mod(iview,2).eq.1)then !Y-view
261          
262    
263  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
264  c$$$     $     ,' -->',ris_eta              riseta = risyeta2(angle)
265             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
266                riseta = risy_cog(angle) !ATTENZIONE!!
267             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
268                riseta = risy_cog(angle) !ATTENZIONE!!
269             else
270                riseta = risy_cog(angle)
271             endif            
272    
273          else                      !X-view
274    
275             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
276                riseta = risxeta2(angle)
277             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
278                riseta = risxeta3(angle)
279             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
280                riseta = risxeta4(angle)
281             else
282                riseta = risx_cog(angle)
283             endif            
284                
285          endif
286    
287    
288   100  return   100  return
# Line 103  c$$$     $     ,' -->',ris_eta Line 302  c$$$     $     ,' -->',ris_eta
302    
303        include 'commontracker.f'        include 'commontracker.f'
304        include 'level1.f'        include 'level1.f'
305  *      include 'calib.f'        include 'calib.f'
306        fbad_eta = 0        fbad_eta = 0
307    
308        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
309                
310           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
311                fbad_eta = fbad_cog(2,ic)
312             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
313                fbad_eta = fbad_cog(3,ic)
314             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
315                fbad_eta = fbad_cog(4,ic)
316             else
317                fbad_eta = fbad_cog(4,ic)
318             endif            
319    
320        else                      !X-view        else                      !X-view
321    
322           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
323              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
324           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
325              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
326           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
327              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
328           endif           else
329                fbad_eta = fbad_cog(4,ic)
330             endif            
331                            
332        endif        endif
333    
# Line 126  c$$$     $     ,' -->',ris_eta Line 335  c$$$     $     ,' -->',ris_eta
335        end        end
336    
337  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
338  c*****************************************************        real function pfaeta2(ic,angle) !(1)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta2(cog2,view,lad,angle)  
       real function pfa_eta2(ic,angle) !(1)  
339  *--------------------------------------------------------------  *--------------------------------------------------------------
340  *     this function returns  *     this function returns
341  *  *
# Line 149  c      real function pfa_eta2(cog2,view, Line 354  c      real function pfa_eta2(cog2,view,
354        real cog2,angle        real cog2,angle
355        integer iview,lad        integer iview,lad
356    
357  c      logical DEBUG        iview = VIEW(ic)            
358  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
359          cog2 = cog(2,ic)          
360        iview = VIEW(ic)              !(1)        pfaeta2=cog2
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfa_eta2=cog2  
361    
362  *     find angular bin  *     find angular bin
363  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
364        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
365           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
366              iangle=iang              iangle=iang
367              goto 98              goto 98
368           endif           endif
369        enddo        enddo
370        if(DEBUG)        if(DEBUG.EQ.1)
371       $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
372        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
373        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
374   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 177  c$$$*     find extremes of interpolation Line 378  c$$$*     find extremes of interpolation
378  c$$$      iflag=0  c$$$      iflag=0
379  c$$$*     --------------------------------  c$$$*     --------------------------------
380  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
381  c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2  c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2
382  c$$$*         goto 100  c$$$*         goto 100
383  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
384  c$$$*     non salto piu`, ma scalo di 1 o -1  c$$$*     non salto piu`, ma scalo di 1 o -1
# Line 230  c            print*,'-----',x1,x2,y1,y2 Line 431  c            print*,'-----',x1,x2,y1,y2
431        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
432        BB=y1-AA*x1        BB=y1-AA*x1
433    
434        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
435        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
436    
437  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
438  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
439  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
440  c$$$      endif  c$$$      endif
441  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
442  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
443  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
444  c$$$      endif  c$$$      endif
445    
446        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
447       $     ,cog2-iadd,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
448    
449    
450   100  return   100  return
451        end        end
452    
453  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
454  c*****************************************************        real function pfaeta3(ic,angle) !(1)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta3(cog3,view,lad,angle)  
       real function pfa_eta3(ic,angle) !(1)  
455  *--------------------------------------------------------------  *--------------------------------------------------------------
456  *     this function returns  *     this function returns
457  *  *
# Line 273  c      real function pfa_eta3(cog3,view, Line 470  c      real function pfa_eta3(cog3,view,
470        real cog3,angle        real cog3,angle
471        integer iview,lad        integer iview,lad
472    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
473    
474          iview = VIEW(ic)            
475        iview = VIEW(ic)              !(1)        lad = nld(MAXS(ic),VIEW(ic))
476        lad = nld(MAXS(ic),VIEW(ic)) !(1)        cog3 = cog(3,ic)            
477        cog3 = cog(3,ic)             !(1)        pfaeta3=cog3
       pfa_eta3=cog3  
478    
479  *     find angular bin  *     find angular bin
480  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 291  c         print*,'~~~~~~~~~~~~ ',iang,an Line 485  c         print*,'~~~~~~~~~~~~ ',iang,an
485              goto 98              goto 98
486           endif           endif
487        enddo        enddo
488        if(DEBUG)        if(DEBUG.EQ.1)
489       $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
490        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
491        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
492   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 354  c            print*,'-----',x1,x2,y1,y2 Line 548  c            print*,'-----',x1,x2,y1,y2
548        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
549        BB=y1-AA*x1        BB=y1-AA*x1
550    
551        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
552        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
553    
554  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
555  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
556  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
557  c$$$      endif  c$$$      endif
558  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
559  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
560  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
561  c$$$      endif  c$$$      endif
562    
563        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
564       $     ,cog3-iadd,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
565    
566   100  return   100  return
567        end        end
568    
569  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
570  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta4(cog4,view,lad,angle)  
       real function pfa_eta4(ic,angle) !(1)  
571  *--------------------------------------------------------------  *--------------------------------------------------------------
572  *     this function returns  *     this function returns
573  *  *
# Line 396  c      real function pfa_eta4(cog4,view, Line 586  c      real function pfa_eta4(cog4,view,
586        real cog4,angle        real cog4,angle
587        integer iview,lad        integer iview,lad
588    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
589    
590        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
591        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
592        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
593        pfa_eta4=cog4        pfaeta4=cog4
594    
595  *     find angular bin  *     find angular bin
596  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 413  c         print*,'~~~~~~~~~~~~ ',iang,an Line 601  c         print*,'~~~~~~~~~~~~ ',iang,an
601              goto 98              goto 98
602           endif           endif
603        enddo        enddo
604        if(DEBUG)        if(DEBUG.EQ.1)
605       $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
606        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
607        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
608   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 476  c            print*,'-----',x1,x2,y1,y2 Line 664  c            print*,'-----',x1,x2,y1,y2
664        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
665        BB=y1-AA*x1        BB=y1-AA*x1
666    
667        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
668        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
669    
670  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
671  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
672  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
673  c$$$      endif  c$$$      endif
674  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
675  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
676  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
677  c$$$      endif  c$$$      endif
678    
679        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
680       $     ,cog4-iadd,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
681    
682   100  return   100  return
683        end        end
# Line 497  c$$$      endif Line 685  c$$$      endif
685    
686    
687  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
       real function cog0(ncog,ic)  
 *-------------------------------------------------  
 *     this function returns  
 *  
 *     - the Center-Of-Gravity of the cluster IC  
 *     evaluated using NCOG strips,  
 *     calculated relative to MAXS(IC)  
 *      
 *     - zero in case that not  enough strips  
 *     have a positive signal  
 *      
 *     NOTE:  
 *     This is the old definition, used by Straulino.  
 *     The new routine, according to Landi,  
 *     is COG(NCOG,IC)  
 *-------------------------------------------------  
   
   
       include 'commontracker.f'  
       include 'level1.f'  
         
 *     --> signal of the central strip  
       sc = CLSIGNAL(INDMAX(ic)) !center  
   
 *     signal of adjacent strips  
 *     --> left  
       sl1 = 0                  !left 1  
       if(  
      $     (INDMAX(ic)-1).ge.INDSTART(ic)  
      $     )  
      $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
   
       sl2 = 0                  !left 2  
       if(  
      $     (INDMAX(ic)-2).ge.INDSTART(ic)  
      $     )  
      $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
   
 *     --> right  
       sr1 = 0                  !right 1  
       if(  
      $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
      $     .or.  
      $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
      $     )  
      $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
   
       sr2 = 0                  !right 2  
       if(  
      $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
      $     .or.  
      $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
      $     )  
      $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
         
 ************************************************************  
 *     COG computation  
 ************************************************************  
   
 c      print*,sl2,sl1,sc,sr1,sr2  
   
       COG = 0.  
           
       if(sl1.gt.sr1.and.sl1.gt.0.)then  
           
          if(ncog.eq.2.and.sl1.ne.0)then  
             COG = -sl1/(sl1+sc)          
          elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
             COG = (sr1-sl1)/(sl1+sc+sr1)  
          elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
             COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  
          else  
             COG = 0.  
          endif  
           
       elseif(sl1.le.sr1.and.sr1.gt.0.)then  
   
          if(ncog.eq.2.and.sr1.ne.0)then  
             COG = sr1/(sc+sr1)              
          elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
             COG = (sr1-sl1)/(sl1+sc+sr1)  
          elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
             COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  
          else  
             COG = 0.  
          endif  
   
       endif  
   
       COG0 = COG  
   
 c      print *,ncog,ic,cog,'/////////////'  
   
       return  
       end  
   
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
688        real function cog(ncog,ic)        real function cog(ncog,ic)
689  *-------------------------------------------------  *-------------------------------------------------
690  *     this function returns  *     this function returns
# Line 613  c      print *,ncog,ic,cog,'//////////// Line 704  c      print *,ncog,ic,cog,'////////////
704        include 'calib.f'        include 'calib.f'
705        include 'level1.f'        include 'level1.f'
706                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
707    
708    
709        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 625  c      common/dbg/DEBUG Line 714  c      common/dbg/DEBUG
714  *     --> signal of the central strip  *     --> signal of the central strip
715           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
716  *     signal of adjacent strips  *     signal of adjacent strips
717           sl1 = 0                !left 1           sl1 = -9999.           !left 1
718           if(           if(
719       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
720       $        )       $        )
721       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
722                    
723           sl2 = 0                !left 2           sl2 = -9999.           !left 2
724           if(           if(
725       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
726       $        )       $        )
727       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
728                    
729           sr1 = 0                !right 1           sr1 = -9999.           !right 1
730           if(           if(
731       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
732       $        .or.       $        .or.
# Line 645  c      common/dbg/DEBUG Line 734  c      common/dbg/DEBUG
734       $        )       $        )
735       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
736                    
737           sr2 = 0                !right 2           sr2 = -9999.           !right 2
738           if(           if(
739       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
740       $        .or.       $        .or.
# Line 655  c      common/dbg/DEBUG Line 744  c      common/dbg/DEBUG
744                    
745           COG = 0.           COG = 0.
746                    
747           if(ncog.eq.2)then  c         print*,'## ',sl2,sl1,sc,sr1,sr2
748    
749    c     ==============================================================
750             if(ncog.eq.1)then
751                COG = 0.
752                if(sr1.gt.sc)cog=1. !NEW
753                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
754    c     ==============================================================
755             elseif(ncog.eq.2)then
756              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
757                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
758              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
759                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
760              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
761                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
762         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
763                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
764         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
765                endif
766    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
767    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
768    c     ==============================================================
769           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
770              COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
771    c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
772    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
773    c     ==============================================================
774           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
775              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
776                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sl1+sc+sr1).ne.0)
777              elseif(sl2.le.sr2)then       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
778                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)              elseif(sl2.lt.sr2)then
779                   if((sr2+sl1+sc+sr1).ne.0)
780         $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
781                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
782                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
783         $              .and.(sl2+sl1+sc+sr1).ne.0 )
784         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
785                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
786         $              .and.(sr2+sl1+sc+sr1).ne.0 )
787         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
788              endif              endif
789           else           else
790              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
791              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
792              COG = 0.              COG = 0.
793           endif           endif
794    
# Line 685  ' Line 802  '
802           iv=VIEW(ic)           iv=VIEW(ic)
803           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
804           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
805             istart = INDSTART(IC)
806           istart=INDSTART(IC)           istop  = TOTCLLENGTH
          istop=TOTCLLENGTH  
807           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
808           COG=0             COG = 0  
809           mu=0           SGN = 0.
810           do i=istart,istop           mu  = 0
811              ipos=i-INDMAX(ic)  c         print*,'-------'
812              ivk=nvk(MAXS(ic)+ipos)           do i = INDMAX(IC),istart,-1
813              is=nst(MAXS(ic)+ipos)              ipos = i-INDMAX(ic)
814  *            print*,'******************',istart,istop,ipos              cut  = incut*CLSIGMA(i)
815  *     $           ,MAXS(ic)+ipos,iv,ivk,is              if(CLSIGNAL(i).ge.cut)then              
816              cut=incut*SIGMA(iv,ivk,is)                 COG = COG + ipos*CLSIGNAL(i)
817                   SGN = SGN + CLSIGNAL(i)
818                   mu = mu + 1
819    c               print*,ipos,CLSIGNAL(i)
820                else
821                   goto 10
822                endif
823             enddo
824     10      continue
825             do i = INDMAX(IC)+1,istop
826                ipos = i-INDMAX(ic)
827                cut  = incut*CLSIGMA(i)
828              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
829                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
830                   SGN = SGN + CLSIGNAL(i)
831                 mu = mu + 1                 mu = mu + 1
832  c               print*,ipos,CLSIGNAL(i),incut,cut  c               print*,ipos,CLSIGNAL(i)
833                else
834                   goto 20
835              endif              endif
836           enddo           enddo
837           COG=COG/DEDX(ic)   20      continue
838  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'           if(SGN.le.0)then
839  c     $        ,cog              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
840                print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
841                print*,(CLSIGNAL(i),i=istart,istop)
842    c            print*,'cog(0,ic) --> NOT EVALUATED '
843             else
844                COG=COG/SGN
845             endif
846    c         print*,'-------'
847                    
848        else        else
849                    
# Line 717  c     $        ,cog Line 854  c     $        ,cog
854    
855        endif        endif
856    
857  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
858    
859        return        return
860        end        end
861    
862  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
863    
864        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
865  *-------------------------------------------------------  *-------------------------------------------------------
866  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 876  c      print *,ncog,ic,cog,'////////////
876        include 'calib.f'        include 'calib.f'
877    
878        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
879           f  = 4.           si = 8.4  !average good-strip noise
880           si = 8.4           f  = 4.   !average bad-strip noise: f*si
881           incut=incuty           incut=incuty
882        else                      !X-view        else                      !X-view
883           f  = 6.           si = 3.9  !average good-strip noise
884           si = 3.9           f  = 6.   !average bad-strip noise: f*si
885           incut=incutx           incut=incutx
886        endif        endif
887                
# Line 754  c      print *,ncog,ic,cog,'//////////// Line 892  c      print *,ncog,ic,cog,'////////////
892  *     --> signal of the central strip  *     --> signal of the central strip
893           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
894           fsc = 1           fsc = 1
895           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
896             fsc = clsigma(INDMAX(ic))/si
897  *     --> signal of adjacent strips  *     --> signal of adjacent strips
898           sl1 = 0                !left 1           sl1 = 0                !left 1
899           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 762  c      print *,ncog,ic,cog,'//////////// Line 901  c      print *,ncog,ic,cog,'////////////
901       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
902       $        )then       $        )then
903              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
904              if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f  c            if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f
905  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
906           endif           endif
907    
908           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 773  c            fsl1 = 0 Line 911  c            fsl1 = 0
911       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
912       $        )then       $        )then
913              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
914              if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f  c            if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f
915  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
916           endif           endif
917           sr1 = 0                !right 1           sr1 = 0                !right 1
918           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 785  c            fsl2 = 0 Line 922  c            fsl2 = 0
922       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
923       $        )then       $        )then
924              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
925              if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f  c            if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f
926  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
927           endif               endif    
928           sr2 = 0                !right 2           sr2 = 0                !right 2
929           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 797  c            fsr1 = 0 Line 933  c            fsr1 = 0
933       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
934       $        )then       $        )then
935              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
936              if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f  c            if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f
937  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
938           endif           endif
939    
940    
941    
942  ************************************************************  ************************************************************
943  *     COG computation  *     COG2-3-4 computation
944  ************************************************************  ************************************************************
945    
946  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
947                    
948           COG = 0.           vCOG = cog(ncog,ic)!0.
949                    
950           if(ncog.eq.2)then           if(ncog.eq.2)then
951              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
952                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
953                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
954                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
955              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
956                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
957                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
958                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
959              endif              endif
960           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
961              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
962              fbad_cog =              fbad_cog =
963       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
964              fbad_cog =              fbad_cog =
965       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
966           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
967              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
968                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
969                 fbad_cog =                 fbad_cog =
970       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
971       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
972                 fbad_cog =                 fbad_cog =
973       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
974       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
975              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
976                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
977                 fbad_cog =                 fbad_cog =
978       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
979       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
980                 fbad_cog =                 fbad_cog =
981       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
982       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
983              endif              endif
984           else           else
985              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
986              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
987              COG = 0.  c            COG = 0.
988           endif           endif
989                    
990        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
991    *     =========================
992    *     COG computation
993    *     =========================
994                    
995           iv=VIEW(ic)           vCOG = cog(0,ic)
996           istart=INDSTART(IC)  
997           istop=TOTCLLENGTH           iv     = VIEW(ic)
998           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
999           COG=0.           istop  = TOTCLLENGTH
1000           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1001           SDE=0.           SGN = 0.
1002           do i=istart,istop           SNU = 0.
1003              ipos=i-INDMAX(ic)           SDE = 0.
1004              il=nvk(MAXS(ic)+ipos)  
1005              is=nst(MAXS(ic)+ipos)           do i=INDMAX(IC),istart,-1
1006              cut=incut*SIGMA(iv,il,is)              ipos = i-INDMAX(ic)
1007                cut  = incut*CLSIGMA(i)
1008              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1009                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1010              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1011                   SDE = SDE + (ipos-vCOG)**2
1012                else
1013                   goto 10
1014                endif            
1015           enddo           enddo
1016           COG=COG/DEDX(ic)   10      continue
1017           do i=istart,istop           do i=INDMAX(IC)+1,istop
1018              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1019              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1020              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1021                 fs=1                 fs = clsigma(i)/si
1022                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1023                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1024                 SDE = SDE + (ipos-COG)**2              else
1025                   goto 20
1026              endif                          endif            
1027           enddo           enddo
1028           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1029             if(SDE.ne.0)then
1030                FBAD_COG=SNU/SDE
1031             else
1032                
1033             endif
1034    
1035        else        else
1036                                        
# Line 902  c      print*,sl2,sl1,sc,sr1,sr2 Line 1048  c      print*,sl2,sl1,sc,sr1,sr2
1048        end        end
1049    
1050    
1051  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1052        real function fbad_cog0(ncog,ic)  c$$$      real function fbad_cog0(ncog,ic)
1053  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1054  *     this function returns a factor that takes into  c$$$*     this function returns a factor that takes into
1055  *     account deterioration of the spatial resolution  c$$$*     account deterioration of the spatial resolution
1056  *     in the case BAD strips are included in the cluster.  c$$$*     in the case BAD strips are included in the cluster.
1057  *     This factor should multiply the nominal spatial  c$$$*     This factor should multiply the nominal spatial
1058  *     resolution.  c$$$*     resolution.
1059  *  c$$$*
1060  *     NB!!!  c$$$*     NB!!!
1061  *     (this is the old version. It consider only the two  c$$$*     (this is the old version. It consider only the two
1062  *     strips with the greatest signal. The new one is  c$$$*     strips with the greatest signal. The new one is
1063  *     fbad_cog(ncog,ic) )  c$$$*     fbad_cog(ncog,ic) )
1064  *      c$$$*    
1065  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1066    c$$$
1067        include 'commontracker.f'  c$$$      include 'commontracker.f'
1068        include 'level1.f'  c$$$      include 'level1.f'
1069        include 'calib.f'  c$$$      include 'calib.f'
1070    c$$$
1071  *     --> signal of the central strip  c$$$*     --> signal of the central strip
1072        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
1073    c$$$
1074  *     signal of adjacent strips  c$$$*     signal of adjacent strips
1075  *     --> left  c$$$*     --> left
1076        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
1077        if(  c$$$      if(
1078       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
1079       $     )  c$$$     $     )
1080       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
1081    c$$$
1082        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
1083        if(  c$$$      if(
1084       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
1085       $     )  c$$$     $     )
1086       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
1087    c$$$
1088  *     --> right  c$$$*     --> right
1089        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
1090        if(  c$$$      if(
1091       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1092       $     .or.  c$$$     $     .or.
1093       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1094       $     )  c$$$     $     )
1095       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
1096    c$$$
1097        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
1098        if(  c$$$      if(
1099       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1100       $     .or.  c$$$     $     .or.
1101       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1102       $     )  c$$$     $     )
1103       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
1104    c$$$
1105    c$$$
1106        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1107           f  = 4.  c$$$         f  = 4.
1108           si = 8.4  c$$$         si = 8.4
1109        else                              !X-view  c$$$      else                              !X-view
1110           f  = 6.  c$$$         f  = 6.
1111           si = 3.9  c$$$         si = 3.9
1112        endif  c$$$      endif
1113    c$$$
1114        fbad_cog = 1.  c$$$      fbad_cog = 1.
1115        f0 = 1  c$$$      f0 = 1
1116        f1 = 1  c$$$      f1 = 1
1117        f2 = 1  c$$$      f2 = 1
1118        f3 = 1    c$$$      f3 = 1  
1119        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
1120            c$$$        
1121           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
1122           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
1123  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
1124    c$$$
1125           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
1126              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
1127           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
1128              fbad_cog = 1.  c$$$            fbad_cog = 1.
1129           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
1130              fbad_cog = 1.  c$$$            fbad_cog = 1.
1131           else  c$$$         else
1132              fbad_cog = 1.  c$$$            fbad_cog = 1.
1133           endif  c$$$         endif
1134            c$$$        
1135        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
1136    c$$$
1137    c$$$
1138           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
1139           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
1140  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
1141    c$$$
1142           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
1143              fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
1144           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
1145              fbad_cog = 1.  c$$$            fbad_cog = 1.
1146           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
1147              fbad_cog = 1.  c$$$            fbad_cog = 1.
1148           else  c$$$         else
1149              fbad_cog = 1.  c$$$            fbad_cog = 1.
1150           endif  c$$$         endif
1151    c$$$
1152        endif  c$$$      endif
1153    c$$$
1154        fbad_cog0 = sqrt(fbad_cog)  c$$$      fbad_cog0 = sqrt(fbad_cog)
1155    c$$$
1156        return  c$$$      return
1157        end  c$$$      end
1158    c$$$
1159    c$$$
1160    c$$$
1161    
1162  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1163    
1164        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1165    
1166        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1167        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1087  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1233  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1233       +/       +/
1234    
1235        V(1)= abs(x)        V(1)= abs(x)
1236          if(V(1).gt.20.)V(1)=20.
1237    
1238        HQUADF = 0.        HQUADF = 0.
1239        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1101  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1248  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1248     20 CONTINUE     20 CONTINUE
1249        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1250                
1251        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1252    
1253        END        END
1254    
1255  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1256        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1257        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1258        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1259        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1177  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1324  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1324       +/       +/
1325    
1326        V(1) =  abs(x)        V(1) =  abs(x)
1327          if(V(1).gt.20.)V(1)=20.
1328    
1329        HQUADF = 0.        HQUADF = 0.
1330        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1331           HQDJ = 0.           HQDJ = 0.
# Line 1190  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1339  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1339     20 CONTINUE     20 CONTINUE
1340        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1341                
1342        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1343    
1344        END        END
1345  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1346        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1347        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1348        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1349        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1265  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1414  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1414       +/       +/
1415    
1416        V(1)=abs(x)        V(1)=abs(x)
1417          if(V(1).gt.20.)V(1)=20.
1418    
1419        HQUADF = 0.        HQUADF = 0.
1420        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1279  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1429  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1429     20 CONTINUE     20 CONTINUE
1430        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1431                
1432        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1433    
1434        END        END
1435  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1436        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1437        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1438        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1439        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1336  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1486  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1486       +/       +/
1487    
1488        v(1)= abs(x)        v(1)= abs(x)
1489          if(V(1).gt.20.)V(1)=20.
1490                
1491        HQUADF = 0.        HQUADF = 0.
1492        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1350  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1501  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1501     20 CONTINUE     20 CONTINUE
1502        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1503    
1504        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1505    
1506        END        END
1507  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1402  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1553  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1553       +/       +/
1554                
1555        V(1)=abs(x)        V(1)=abs(x)
1556          if(V(1).gt.20.)V(1)=20.
1557    
1558        HQUADF = 0.        HQUADF = 0.
1559        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1482  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1634  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1634       +/       +/
1635    
1636        V(1)=abs(x)        V(1)=abs(x)
1637          if(V(1).gt.20.)V(1)=20.
1638    
1639        HQUADF = 0.        HQUADF = 0.
1640        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1499  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1652  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1652        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1653    
1654        END        END
1655    
1656    
1657  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1658          real function pfacorr(ic,angle)
1659    *--------------------------------------------------------------
1660    *     this function returns the landi correction for this cluster
1661    *--------------------------------------------------------------
1662          include 'commontracker.f'
1663          include 'calib.f'
1664          include 'level1.f'
1665    
1666          real angle
1667          integer iview,lad
1668    
1669          iview = VIEW(ic)            
1670          lad = nld(MAXS(ic),VIEW(ic))
1671    
1672    *     find angular bin
1673    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1674          do iang=1,nangbin
1675             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1676                iangle=iang
1677                goto 98
1678             endif
1679          enddo
1680          if(DEBUG.eq.1)
1681         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1682          if(angle.lt.angL(1))iang=1
1683          if(angle.gt.angR(nangbin))iang=nangbin
1684     98   continue                  !jump here if ok
1685    
1686          pfacorr = fcorr(iview,lad,iang)
1687    
1688          if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1689    
1690    
1691     100  return
1692          end
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23