/[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.10 by pam-fi, Mon May 14 11:03:06 2007 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
3    *     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)        subroutine idtoc(ipfa,cpfa)
35                
36        integer ipfa        integer ipfa
37        character*4 cpfa        character*10 cpfa
38    
39        CPFA='COG4'        CPFA='COG4'
40        if(ipfa.eq.0)CPFA='ETA'        if(ipfa.eq.0)CPFA='ETA'
41        if(ipfa.eq.2)CPFA='ETA2'        if(ipfa.eq.2)CPFA='ETA2'
42        if(ipfa.eq.3)CPFA='ETA3'        if(ipfa.eq.3)CPFA='ETA3'
43        if(ipfa.eq.4)CPFA='ETA4'        if(ipfa.eq.4)CPFA='ETA4'
44          if(ipfa.eq.5)CPFA='ETAL'
45        if(ipfa.eq.10)CPFA='COG'        if(ipfa.eq.10)CPFA='COG'
46        if(ipfa.eq.11)CPFA='COG1'        if(ipfa.eq.11)CPFA='COG1'
47        if(ipfa.eq.12)CPFA='COG2'        if(ipfa.eq.12)CPFA='COG2'
# Line 18  Line 50 
50                
51        end        end
52    
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 *     this file contains all subroutines and functions  
 *     that are needed for position finding algorithms  
 *      
 *  
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
53    
54    
55        integer function npfastrips(ic,PFA,angle)        integer function npfastrips(ic,angle)
56  *--------------------------------------------------------------  *--------------------------------------------------------------
57  *     thid function returns the number of strips used  *     thid function returns the number of strips used
58  *     to evaluate the position of a cluster, according to the p.f.a.  *     to evaluate the position of a cluster, according to the p.f.a.
# Line 35  Line 61 
61        include 'level1.f'        include 'level1.f'
62        include 'calib.f'        include 'calib.f'
63    
64        character*4 usedPFA,PFA        character*4 usedPFA
65          
66    
67    
68        usedPFA=PFA        call idtoc(pfaid,usedPFA)
69    
70        npfastrips=0        npfastrips=-1
71    
72        if(usedPFA.eq.'COG1')npfastrips=1        if(usedPFA.eq.'COG1')npfastrips=1
73        if(usedPFA.eq.'COG2')npfastrips=2        if(usedPFA.eq.'COG2')npfastrips=2
# Line 50  Line 77 
77        if(usedPFA.eq.'ETA3')npfastrips=3        if(usedPFA.eq.'ETA3')npfastrips=3
78        if(usedPFA.eq.'ETA4')npfastrips=4        if(usedPFA.eq.'ETA4')npfastrips=4
79  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
80        if(usedPFA.eq.'ETA')then        if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
81  c         print*,VIEW(ic),angle  c         print*,VIEW(ic),angle
82           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
83              if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then              if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
# Line 60  c         print*,VIEW(ic),angle Line 87  c         print*,VIEW(ic),angle
87              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
88                 npfastrips=4                 npfastrips=4
89              else              else
90                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
91              endif                                      endif                        
92           else                   !X-view           else                   !X-view
93              if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then              if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
# Line 71  c               usedPFA='COG' Line 97  c               usedPFA='COG'
97              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
98                 npfastrips=4                 npfastrips=4
99              else              else
100                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
101              endif                                      endif                        
102           endif           endif
103        endif        endif
104  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
105        if(usedPFA.eq.'COG')then        if(usedPFA.eq.'COG')then
106    
107           iv=VIEW(ic)           npfastrips=0
108           if(mod(iv,2).eq.1)incut=incuty  
109           if(mod(iv,2).eq.0)incut=incutx  c$$$         iv=VIEW(ic)
110           istart = INDSTART(IC)  c$$$         if(mod(iv,2).eq.1)incut=incuty
111           istop  = TOTCLLENGTH  c$$$         if(mod(iv,2).eq.0)incut=incutx
112           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1  c$$$         istart = INDSTART(IC)
113           mu  = 0  c$$$         istop  = TOTCLLENGTH
114           do i = INDMAX(IC),istart,-1  c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
115              ipos = i-INDMAX(ic)  c$$$         mu  = 0
116              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC),istart,-1
117              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
118                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
119                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
120              else  c$$$               mu = mu + 1
121                 goto 10  c$$$               print*,i,mu
122              endif  c$$$            else
123           enddo  c$$$               goto 10
124   10      continue  c$$$            endif
125           do i = INDMAX(IC)+1,istop  c$$$         enddo
126              ipos = i-INDMAX(ic)  c$$$ 10      continue
127              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC)+1,istop
128              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
129                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
130                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
131              else  c$$$               mu = mu + 1
132                 goto 20  c$$$               print*,i,mu
133              endif  c$$$            else
134           enddo  c$$$               goto 20
135   20      continue  c$$$            endif
136           npfastrips=mu  c$$$         enddo
137    c$$$ 20      continue
138    c$$$         npfastrips=mu
139    
140        endif        endif
141  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
142    
143  c      print*,pfastrips  c      print*,pfaid,usedPFA,angle,npfastrips
144    
145        return        return
146        end        end
# Line 136  c      print*,pfastrips Line 163  c      print*,pfastrips
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           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
167              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
168           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then  cc            print*,pfaeta2(ic,angle)
169             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
170              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
171           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
172              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
173           else           else
174              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 148  c      print*,pfastrips Line 176  c      print*,pfastrips
176    
177        else                      !X-view        else                      !X-view
178    
179           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
180              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
181           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
182              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
183           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
184              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
185           else           else
186              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 164  c      print*,pfastrips Line 192  c      print*,pfastrips
192        end        end
193    
194  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
195        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
196    *--------------------------------------------------------------
197    *     this function returns the position (in strip units)
198    *     it calls:
199    *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
200    *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
201    *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
202    *     according to the angle
203    *--------------------------------------------------------------
204          include 'commontracker.f'
205          include 'level1.f'
206          include 'calib.f'
207          
208          pfaetal = 0
209    
210          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
211          
212             if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
213                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
224    
225             if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
226                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
227    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
228             elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
229                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
230             elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
231                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
232             else
233                pfaetal = cog(4,ic)
234             endif            
235                
236          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  *     this function returns the average spatial resolution
245  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
246  *     it calls:  *     it calls:
247  *     - risx_eta2(angle)  *     - risxeta2(angle)
248  *     - risy_eta2(angle)  *     - risyeta2(angle)
249  *     - risx_eta3(angle)  *     - risxeta3(angle)
250  *     - risx_eta4(angle)  *     - risxeta4(angle)
251  *     according to the angle  *     according to the angle
252  *--------------------------------------------------------------  *--------------------------------------------------------------
253        include 'commontracker.f'        include 'commontracker.f'
254        include 'level1.f'        include 'level1.f'
255        include 'calib.f'        include 'calib.f'
256    
257        ris_eta = 0        riseta = 0
258    
259        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  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           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
264              ris_eta = risy_eta2(angle)              riseta = risyeta2(angle)
265           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
266              ris_eta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
267           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
268              ris_eta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
269           else           else
270              ris_eta = risy_cog(angle)              riseta = risy_cog(angle)
271           endif                       endif            
272    
273        else                      !X-view        else                      !X-view
274    
275           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
276              ris_eta = risx_eta2(angle)              riseta = risxeta2(angle)
277           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
278              ris_eta = risx_eta3(angle)              riseta = risxeta3(angle)
279           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
280              ris_eta = risx_eta4(angle)              riseta = risxeta4(angle)
281           else           else
282              ris_eta = risx_cog(angle)              riseta = risx_cog(angle)
283           endif                       endif            
284                            
285        endif        endif
286    
287    
288   100  return   100  return
289        end        end
290    
# Line 290  c      print*,pfastrips Line 367  c      print*,pfastrips
367              goto 98              goto 98
368           endif           endif
369        enddo        enddo
370        if(DEBUG)        if(DEBUG.EQ.1)
371       $     print*,'pfaeta2 *** 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
# Line 366  c$$$         pfaeta2=pfaeta2+1.   !temp Line 443  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,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
448    
449    
# Line 408  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*,'pfaeta3 *** 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
# Line 483  c$$$         pfaeta2=pfaeta2+1.   !temp Line 560  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,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
565    
566   100  return   100  return
# Line 524  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*,'pfaeta4 *** 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
# Line 599  c$$$         pfaeta2=pfaeta2+1.   !temp Line 676  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,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
681    
682   100  return   100  return
# Line 608  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 734  c      print *,ncog,ic,cog,'//////////// Line 714  c      print *,ncog,ic,cog,'////////////
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 754  c      print *,ncog,ic,cog,'//////////// Line 734  c      print *,ncog,ic,cog,'////////////
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 766  c      print *,ncog,ic,cog,'//////////// Line 746  c      print *,ncog,ic,cog,'////////////
746                    
747  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c         print*,'## ',sl2,sl1,sc,sr1,sr2
748    
749    c     ==============================================================
750           if(ncog.eq.1)then           if(ncog.eq.1)then
751              COG = 0.              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           elseif(ncog.eq.2)then
756              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
757                 if((sl1+sc).ne.0)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                 if((sc+sr1).ne.0)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               if((sl1+sc+sr1).ne.0)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                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sl1+sc+sr1).ne.0)
777       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
778              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
779                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sl1+sc+sr1).ne.0)
780       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              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
# Line 814  c         print*,'-------' Line 816  c         print*,'-------'
816                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
817                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
818                 mu = mu + 1                 mu = mu + 1
819                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
820              else              else
821                 goto 10                 goto 10
822              endif              endif
# Line 827  c         print*,'-------' Line 829  c         print*,'-------'
829                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
830                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
831                 mu = mu + 1                 mu = mu + 1
832                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
833              else              else
834                 goto 20                 goto 20
835              endif              endif
836           enddo           enddo
837   20      continue   20      continue
838           if(SGN.le.0)then           if(SGN.le.0)then
839  c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
840              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
841              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
842  c            print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
# Line 999  c            COG = 0. Line 1001  c            COG = 0.
1001           SGN = 0.           SGN = 0.
1002           SNU = 0.           SNU = 0.
1003           SDE = 0.           SDE = 0.
 c$$$         do i=INDMAX(IC),istart,-1  
 c$$$            ipos = i-INDMAX(ic)  
 c$$$            cut  = incut*CLSIGMA(i)  
 c$$$            if(CLSIGNAL(i).gt.cut)then  
 c$$$               COG = COG + ipos*CLSIGNAL(i)  
 c$$$               SGN = SGN + CLSIGNAL(i)  
 c$$$            else  
 c$$$               goto 10  
 c$$$            endif  
 c$$$         enddo  
 c$$$ 10      continue  
 c$$$         do i=INDMAX(IC)+1,istop  
 c$$$            ipos = i-INDMAX(ic)  
 c$$$            cut  = incut*CLSIGMA(i)  
 c$$$            if(CLSIGNAL(i).gt.cut)then  
 c$$$               COG = COG + ipos*CLSIGNAL(i)  
 c$$$               SGN = SGN + CLSIGNAL(i)  
 c$$$            else  
 c$$$               goto 20  
 c$$$            endif  
 c$$$         enddo  
 c$$$ 20      continue  
 c$$$         if(SGN.le.0)then  
 c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN  
 c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$            print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '  
 c$$$         else  
 c$$$            COG=COG/SGN  
 c$$$         endif  
1004    
1005           do i=INDMAX(IC),istart,-1           do i=INDMAX(IC),istart,-1
1006              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
# Line 1076  c$$$         endif Line 1048  c$$$         endif
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 1276  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 1367  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 1457  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 1529  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 1680  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.10  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23