/[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.10 by pam-fi, Mon May 14 11:03:06 2007 UTC
# Line 1  Line 1 
1    
2    
3          subroutine idtoc(ipfa,cpfa)
4          
5          integer ipfa
6          character*4 cpfa
7    
8          CPFA='COG4'
9          if(ipfa.eq.0)CPFA='ETA'
10          if(ipfa.eq.2)CPFA='ETA2'
11          if(ipfa.eq.3)CPFA='ETA3'
12          if(ipfa.eq.4)CPFA='ETA4'
13          if(ipfa.eq.10)CPFA='COG'
14          if(ipfa.eq.11)CPFA='COG1'
15          if(ipfa.eq.12)CPFA='COG2'
16          if(ipfa.eq.13)CPFA='COG3'
17          if(ipfa.eq.14)CPFA='COG4'
18          
19          end
20    
21  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
22  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
23  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms
# Line 6  Line 26 
26  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
27    
28    
29          integer function npfastrips(ic,PFA,angle)
30    *--------------------------------------------------------------
31    *     thid function returns the number of strips used
32    *     to evaluate the position of a cluster, according to the p.f.a.
33    *--------------------------------------------------------------
34          include 'commontracker.f'
35          include 'level1.f'
36          include 'calib.f'
37    
38          character*4 usedPFA,PFA
39    
40    
41          usedPFA=PFA
42    
43          npfastrips=0
44    
45          if(usedPFA.eq.'COG1')npfastrips=1
46          if(usedPFA.eq.'COG2')npfastrips=2
47          if(usedPFA.eq.'COG3')npfastrips=3
48          if(usedPFA.eq.'COG4')npfastrips=4
49          if(usedPFA.eq.'ETA2')npfastrips=2
50          if(usedPFA.eq.'ETA3')npfastrips=3
51          if(usedPFA.eq.'ETA4')npfastrips=4
52    *     ----------------------------------------------------------------
53          if(usedPFA.eq.'ETA')then
54    c         print*,VIEW(ic),angle
55             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
56                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
57                   npfastrips=2
58                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
59                   npfastrips=3
60                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
61                   npfastrips=4
62                else
63                   npfastrips=4
64    c               usedPFA='COG'
65                endif                        
66             else                   !X-view
67                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
68                   npfastrips=2
69                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
70                   npfastrips=3
71                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
72                   npfastrips=4
73                else
74                   npfastrips=4
75    c               usedPFA='COG'
76                endif                        
77             endif
78          endif
79    *     ----------------------------------------------------------------
80          if(usedPFA.eq.'COG')then
81    
82             iv=VIEW(ic)
83             if(mod(iv,2).eq.1)incut=incuty
84             if(mod(iv,2).eq.0)incut=incutx
85             istart = INDSTART(IC)
86             istop  = TOTCLLENGTH
87             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
88             mu  = 0
89             do i = INDMAX(IC),istart,-1
90                ipos = i-INDMAX(ic)
91                cut  = incut*CLSIGMA(i)
92                if(CLSIGNAL(i).ge.cut)then
93                   mu = mu + 1
94                   print*,i,mu
95                else
96                   goto 10
97                endif
98             enddo
99     10      continue
100             do i = INDMAX(IC)+1,istop
101                ipos = i-INDMAX(ic)
102                cut  = incut*CLSIGMA(i)
103                if(CLSIGNAL(i).ge.cut)then
104                   mu = mu + 1
105                   print*,i,mu
106                else
107                   goto 20
108                endif
109             enddo
110     20      continue
111             npfastrips=mu
112    
113          endif
114    *     ----------------------------------------------------------------
115    
116    c      print*,pfastrips
117    
118          return
119          end
120    
121  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
122        real function pfa_eta(ic,angle)        real function pfaeta(ic,angle)
123  *--------------------------------------------------------------  *--------------------------------------------------------------
124  *     this function returns the position (in strip units)  *     this function returns the position (in strip units)
125  *     it calls:  *     it calls:
126  *     - pfa_eta2(ic,angle)  *     - pfaeta2(ic,angle)
127  *     - pfa_eta3(ic,angle)  *     - pfaeta3(ic,angle)
128  *     - pfa_eta4(ic,angle)  *     - pfaeta4(ic,angle)
129  *     according to the angle  *     according to the angle
130  *--------------------------------------------------------------  *--------------------------------------------------------------
131        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
132        include 'level1.f'        include 'level1.f'
133          include 'calib.f'
134                
135        pfa_eta = 0        pfaeta = 0
136    
137        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
138                
139           pfa_eta = pfa_eta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
140                  pfaeta = pfaeta2(ic,angle)
141             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
142                pfaeta = pfaeta3(ic,angle)
143             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
144                pfaeta = pfaeta4(ic,angle)
145             else
146                pfaeta = cog(4,ic)
147             endif            
148    
149        else                      !X-view        else                      !X-view
150    
151           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
152              pfa_eta = pfa_eta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
153           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
154              pfa_eta = pfa_eta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
155           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
156              pfa_eta = pfa_eta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
157           endif           else
158                pfaeta = cog(4,ic)
159             endif            
160                            
161        endif        endif
162          
   
163   100  return   100  return
164        end        end
165    
# Line 46  c      include 'calib.f' Line 167  c      include 'calib.f'
167        real function ris_eta(ic,angle)        real function ris_eta(ic,angle)
168  *--------------------------------------------------------------  *--------------------------------------------------------------
169  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
170  *     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
171  *     it calls:  *     it calls:
172  *     - risx_eta2(angle)  *     - risx_eta2(angle)
173  *     - risy_eta2(angle)  *     - risy_eta2(angle)
# Line 55  c      include 'calib.f' Line 176  c      include 'calib.f'
176  *     according to the angle  *     according to the angle
177  *--------------------------------------------------------------  *--------------------------------------------------------------
178        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
179        include 'level1.f'        include 'level1.f'
180          include 'calib.f'
181    
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
         
182        ris_eta = 0        ris_eta = 0
183    
184        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
185                
186           ris_eta = risy_eta2(angle)  
187           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
188                ris_eta = risy_eta2(angle)
189             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
190                ris_eta = risy_cog(angle) !ATTENZIONE!!
191             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
192                ris_eta = risy_cog(angle) !ATTENZIONE!!
193             else
194                ris_eta = risy_cog(angle)
195             endif            
196    
197        else                      !X-view        else                      !X-view
198    
199           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
200              ris_eta = risx_eta2(angle)              ris_eta = risx_eta2(angle)
201           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
202              ris_eta = risx_eta3(angle)              ris_eta = risx_eta3(angle)
203           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
204              ris_eta = risx_eta4(angle)              ris_eta = risx_eta4(angle)
205           elseif(abs(angle).gt.21.)then           else
206              ris_eta = risx_eta4(21.)              ris_eta = risx_cog(angle)
207           endif           endif            
208                            
209        endif        endif
210    
 c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'  
 c$$$     $     ,' -->',ris_eta  
   
   
211   100  return   100  return
212        end        end
213    
# Line 103  c$$$     $     ,' -->',ris_eta Line 225  c$$$     $     ,' -->',ris_eta
225    
226        include 'commontracker.f'        include 'commontracker.f'
227        include 'level1.f'        include 'level1.f'
228  *      include 'calib.f'        include 'calib.f'
229        fbad_eta = 0        fbad_eta = 0
230    
231        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
232                
233           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
234                fbad_eta = fbad_cog(2,ic)
235             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
236                fbad_eta = fbad_cog(3,ic)
237             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
238                fbad_eta = fbad_cog(4,ic)
239             else
240                fbad_eta = fbad_cog(4,ic)
241             endif            
242    
243        else                      !X-view        else                      !X-view
244    
245           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
246              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
247           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
248              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
249           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
250              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
251           endif           else
252                fbad_eta = fbad_cog(4,ic)
253             endif            
254                            
255        endif        endif
256    
# Line 126  c$$$     $     ,' -->',ris_eta Line 258  c$$$     $     ,' -->',ris_eta
258        end        end
259    
260  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
261  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)  
262  *--------------------------------------------------------------  *--------------------------------------------------------------
263  *     this function returns  *     this function returns
264  *  *
# Line 149  c      real function pfa_eta2(cog2,view, Line 277  c      real function pfa_eta2(cog2,view,
277        real cog2,angle        real cog2,angle
278        integer iview,lad        integer iview,lad
279    
280  c      logical DEBUG        iview = VIEW(ic)            
281  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
282          cog2 = cog(2,ic)          
283        iview = VIEW(ic)              !(1)        pfaeta2=cog2
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfa_eta2=cog2  
284    
285  *     find angular bin  *     find angular bin
286  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
287        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
288           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
289              iangle=iang              iangle=iang
290              goto 98              goto 98
291           endif           endif
292        enddo        enddo
293        if(DEBUG)        if(DEBUG)
294       $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
295        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
296        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
297   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 177  c$$$*     find extremes of interpolation Line 301  c$$$*     find extremes of interpolation
301  c$$$      iflag=0  c$$$      iflag=0
302  c$$$*     --------------------------------  c$$$*     --------------------------------
303  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
304  c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2  c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2
305  c$$$*         goto 100  c$$$*         goto 100
306  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
307  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 354  c            print*,'-----',x1,x2,y1,y2
354        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
355        BB=y1-AA*x1        BB=y1-AA*x1
356    
357        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
358        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
359    
360  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
361  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
362  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
363  c$$$      endif  c$$$      endif
364  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
365  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
366  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
367  c$$$      endif  c$$$      endif
368    
369        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'
370       $     ,cog2-iadd,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
371    
372    
373   100  return   100  return
374        end        end
375    
376  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
377  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)  
378  *--------------------------------------------------------------  *--------------------------------------------------------------
379  *     this function returns  *     this function returns
380  *  *
# Line 273  c      real function pfa_eta3(cog3,view, Line 393  c      real function pfa_eta3(cog3,view,
393        real cog3,angle        real cog3,angle
394        integer iview,lad        integer iview,lad
395    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
396    
397        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
398        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
399        cog3 = cog(3,ic)             !(1)        cog3 = cog(3,ic)            
400        pfa_eta3=cog3        pfaeta3=cog3
401    
402  *     find angular bin  *     find angular bin
403  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 292  c         print*,'~~~~~~~~~~~~ ',iang,an Line 409  c         print*,'~~~~~~~~~~~~ ',iang,an
409           endif           endif
410        enddo        enddo
411        if(DEBUG)        if(DEBUG)
412       $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
413        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
414        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
415   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 354  c            print*,'-----',x1,x2,y1,y2 Line 471  c            print*,'-----',x1,x2,y1,y2
471        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
472        BB=y1-AA*x1        BB=y1-AA*x1
473    
474        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
475        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
476    
477  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
478  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
479  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
480  c$$$      endif  c$$$      endif
481  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
482  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
483  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
484  c$$$      endif  c$$$      endif
485    
486        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'
487       $     ,cog3-iadd,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
488    
489   100  return   100  return
490        end        end
491    
492  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
493  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)  
494  *--------------------------------------------------------------  *--------------------------------------------------------------
495  *     this function returns  *     this function returns
496  *  *
# Line 396  c      real function pfa_eta4(cog4,view, Line 509  c      real function pfa_eta4(cog4,view,
509        real cog4,angle        real cog4,angle
510        integer iview,lad        integer iview,lad
511    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
512    
513        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
514        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
515        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
516        pfa_eta4=cog4        pfaeta4=cog4
517    
518  *     find angular bin  *     find angular bin
519  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 414  c         print*,'~~~~~~~~~~~~ ',iang,an Line 525  c         print*,'~~~~~~~~~~~~ ',iang,an
525           endif           endif
526        enddo        enddo
527        if(DEBUG)        if(DEBUG)
528       $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
529        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
530        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
531   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 476  c            print*,'-----',x1,x2,y1,y2 Line 587  c            print*,'-----',x1,x2,y1,y2
587        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
588        BB=y1-AA*x1        BB=y1-AA*x1
589    
590        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
591        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
592    
593  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
594  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
595  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
596  c$$$      endif  c$$$      endif
597  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
598  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
599  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
600  c$$$      endif  c$$$      endif
601    
602        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'
603       $     ,cog4-iadd,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
604    
605   100  return   100  return
606        end        end
# Line 613  c      print *,ncog,ic,cog,'//////////// Line 724  c      print *,ncog,ic,cog,'////////////
724        include 'calib.f'        include 'calib.f'
725        include 'level1.f'        include 'level1.f'
726                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
727    
728    
729        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 655  c      common/dbg/DEBUG Line 764  c      common/dbg/DEBUG
764                    
765           COG = 0.           COG = 0.
766                    
767           if(ncog.eq.2)then  c         print*,'## ',sl2,sl1,sc,sr1,sr2
768    
769             if(ncog.eq.1)then
770                COG = 0.
771             elseif(ncog.eq.2)then
772              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
773                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
774              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
775                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)            
776              endif              endif
777           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
778              COG = (sr1-sl1)/(sl1+sc+sr1)               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)
779           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
780              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
781                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sl1+sc+sr1).ne.0)
782         $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
783              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
784                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)                  if((sl2+sl1+sc+sr1).ne.0)
785         $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
786              endif              endif
787           else           else
788              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
789              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
790              COG = 0.              COG = 0.
791           endif           endif
792    
# Line 685  ' Line 800  '
800           iv=VIEW(ic)           iv=VIEW(ic)
801           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
802           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
803             istart = INDSTART(IC)
804           istart=INDSTART(IC)           istop  = TOTCLLENGTH
          istop=TOTCLLENGTH  
805           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
806           COG=0             COG = 0  
807           mu=0           SGN = 0.
808           do i=istart,istop           mu  = 0
809              ipos=i-INDMAX(ic)  c         print*,'-------'
810              ivk=nvk(MAXS(ic)+ipos)           do i = INDMAX(IC),istart,-1
811              is=nst(MAXS(ic)+ipos)              ipos = i-INDMAX(ic)
812  *            print*,'******************',istart,istop,ipos              cut  = incut*CLSIGMA(i)
813  *     $           ,MAXS(ic)+ipos,iv,ivk,is              if(CLSIGNAL(i).ge.cut)then              
814              cut=incut*SIGMA(iv,ivk,is)                 COG = COG + ipos*CLSIGNAL(i)
815                   SGN = SGN + CLSIGNAL(i)
816                   mu = mu + 1
817                   print*,ipos,CLSIGNAL(i)
818                else
819                   goto 10
820                endif
821             enddo
822     10      continue
823             do i = INDMAX(IC)+1,istop
824                ipos = i-INDMAX(ic)
825                cut  = incut*CLSIGMA(i)
826              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
827                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
828                   SGN = SGN + CLSIGNAL(i)
829                 mu = mu + 1                 mu = mu + 1
830  c               print*,ipos,CLSIGNAL(i),incut,cut                 print*,ipos,CLSIGNAL(i)
831                else
832                   goto 20
833              endif              endif
834           enddo           enddo
835           COG=COG/DEDX(ic)   20      continue
836  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'           if(SGN.le.0)then
837  c     $        ,cog  c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN
838                print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
839                print*,(CLSIGNAL(i),i=istart,istop)
840    c            print*,'cog(0,ic) --> NOT EVALUATED '
841             else
842                COG=COG/SGN
843             endif
844    c         print*,'-------'
845                    
846        else        else
847                    
# Line 717  c     $        ,cog Line 852  c     $        ,cog
852    
853        endif        endif
854    
855  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
856    
857        return        return
858        end        end
859    
860  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
861    
862        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
863  *-------------------------------------------------------  *-------------------------------------------------------
864  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 874  c      print *,ncog,ic,cog,'////////////
874        include 'calib.f'        include 'calib.f'
875    
876        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
877           f  = 4.           si = 8.4  !average good-strip noise
878           si = 8.4           f  = 4.   !average bad-strip noise: f*si
879           incut=incuty           incut=incuty
880        else                      !X-view        else                      !X-view
881           f  = 6.           si = 3.9  !average good-strip noise
882           si = 3.9           f  = 6.   !average bad-strip noise: f*si
883           incut=incutx           incut=incutx
884        endif        endif
885                
# Line 754  c      print *,ncog,ic,cog,'//////////// Line 890  c      print *,ncog,ic,cog,'////////////
890  *     --> signal of the central strip  *     --> signal of the central strip
891           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
892           fsc = 1           fsc = 1
893           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
894             fsc = clsigma(INDMAX(ic))/si
895  *     --> signal of adjacent strips  *     --> signal of adjacent strips
896           sl1 = 0                !left 1           sl1 = 0                !left 1
897           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 762  c      print *,ncog,ic,cog,'//////////// Line 899  c      print *,ncog,ic,cog,'////////////
899       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
900       $        )then       $        )then
901              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
902              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
903  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
904           endif           endif
905    
906           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 773  c            fsl1 = 0 Line 909  c            fsl1 = 0
909       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
910       $        )then       $        )then
911              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
912              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
913  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
914           endif           endif
915           sr1 = 0                !right 1           sr1 = 0                !right 1
916           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 785  c            fsl2 = 0 Line 920  c            fsl2 = 0
920       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
921       $        )then       $        )then
922              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
923              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
924  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
925           endif               endif    
926           sr2 = 0                !right 2           sr2 = 0                !right 2
927           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 797  c            fsr1 = 0 Line 931  c            fsr1 = 0
931       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
932       $        )then       $        )then
933              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
934              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
935  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
936           endif           endif
937    
938    
939    
940  ************************************************************  ************************************************************
941  *     COG computation  *     COG2-3-4 computation
942  ************************************************************  ************************************************************
943    
944  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
945                    
946           COG = 0.           vCOG = cog(ncog,ic)!0.
947                    
948           if(ncog.eq.2)then           if(ncog.eq.2)then
949              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
950                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
951                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
952                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
953              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
954                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
955                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
956                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
957              endif              endif
958           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
959              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
960              fbad_cog =              fbad_cog =
961       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
962              fbad_cog =              fbad_cog =
963       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
964           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
965              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
966                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
967                 fbad_cog =                 fbad_cog =
968       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
969       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
970                 fbad_cog =                 fbad_cog =
971       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
972       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
973              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
974                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
975                 fbad_cog =                 fbad_cog =
976       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
977       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
978                 fbad_cog =                 fbad_cog =
979       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
980       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
981              endif              endif
982           else           else
983              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
984              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
985              COG = 0.  c            COG = 0.
986           endif           endif
987                    
988        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
989    *     =========================
990    *     COG computation
991    *     =========================
992                    
993           iv=VIEW(ic)           vCOG = cog(0,ic)
994           istart=INDSTART(IC)  
995           istop=TOTCLLENGTH           iv     = VIEW(ic)
996           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
997           COG=0.           istop  = TOTCLLENGTH
998           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
999           SDE=0.           SGN = 0.
1000           do i=istart,istop           SNU = 0.
1001              ipos=i-INDMAX(ic)           SDE = 0.
1002              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1003              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1004              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1005    c$$$            if(CLSIGNAL(i).gt.cut)then
1006    c$$$               COG = COG + ipos*CLSIGNAL(i)
1007    c$$$               SGN = SGN + CLSIGNAL(i)
1008    c$$$            else
1009    c$$$               goto 10
1010    c$$$            endif
1011    c$$$         enddo
1012    c$$$ 10      continue
1013    c$$$         do i=INDMAX(IC)+1,istop
1014    c$$$            ipos = i-INDMAX(ic)
1015    c$$$            cut  = incut*CLSIGMA(i)
1016    c$$$            if(CLSIGNAL(i).gt.cut)then
1017    c$$$               COG = COG + ipos*CLSIGNAL(i)
1018    c$$$               SGN = SGN + CLSIGNAL(i)
1019    c$$$            else
1020    c$$$               goto 20
1021    c$$$            endif
1022    c$$$         enddo
1023    c$$$ 20      continue
1024    c$$$         if(SGN.le.0)then
1025    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1026    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1027    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1028    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1029    c$$$         else
1030    c$$$            COG=COG/SGN
1031    c$$$         endif
1032    
1033             do i=INDMAX(IC),istart,-1
1034                ipos = i-INDMAX(ic)
1035                cut  = incut*CLSIGMA(i)
1036              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1037                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1038              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1039                   SDE = SDE + (ipos-vCOG)**2
1040                else
1041                   goto 10
1042                endif            
1043           enddo           enddo
1044           COG=COG/DEDX(ic)   10      continue
1045           do i=istart,istop           do i=INDMAX(IC)+1,istop
1046              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1047              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1048              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1049                 fs=1                 fs = clsigma(i)/si
1050                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1051                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1052                 SDE = SDE + (ipos-COG)**2              else
1053                   goto 20
1054              endif                          endif            
1055           enddo           enddo
1056           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1057             if(SDE.ne.0)then
1058                FBAD_COG=SNU/SDE
1059             else
1060                
1061             endif
1062    
1063        else        else
1064                                        
# Line 1087  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1261  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1261       +/       +/
1262    
1263        V(1)= abs(x)        V(1)= abs(x)
1264          if(V(1).gt.20.)V(1)=20.
1265    
1266        HQUADF = 0.        HQUADF = 0.
1267        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1177  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1352  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1352       +/       +/
1353    
1354        V(1) =  abs(x)        V(1) =  abs(x)
1355          if(V(1).gt.20.)V(1)=20.
1356    
1357        HQUADF = 0.        HQUADF = 0.
1358        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1359           HQDJ = 0.           HQDJ = 0.
# Line 1265  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1442  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1442       +/       +/
1443    
1444        V(1)=abs(x)        V(1)=abs(x)
1445          if(V(1).gt.20.)V(1)=20.
1446    
1447        HQUADF = 0.        HQUADF = 0.
1448        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1336  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1514  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1514       +/       +/
1515    
1516        v(1)= abs(x)        v(1)= abs(x)
1517          if(V(1).gt.20.)V(1)=20.
1518                
1519        HQUADF = 0.        HQUADF = 0.
1520        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1402  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1581  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1581       +/       +/
1582                
1583        V(1)=abs(x)        V(1)=abs(x)
1584          if(V(1).gt.20.)V(1)=20.
1585    
1586        HQUADF = 0.        HQUADF = 0.
1587        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1482  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1662  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1662       +/       +/
1663    
1664        V(1)=abs(x)        V(1)=abs(x)
1665          if(V(1).gt.20.)V(1)=20.
1666    
1667        HQUADF = 0.        HQUADF = 0.
1668        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1499  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1680  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1680        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1681    
1682        END        END
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23