/[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.3 by pam-fi, Thu Sep 28 14:04:40 2006 UTC revision 1.12 by pam-fi, Wed May 23 15:31:19 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    
166  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
167        real function ris_eta(ic,angle)  c      real function riseta(ic,angle)
168          real function riseta(iview,angle)
169  *--------------------------------------------------------------  *--------------------------------------------------------------
170  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
171  *     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
172  *     it calls:  *     it calls:
173  *     - risx_eta2(angle)  *     - risx_eta2(angle)
174  *     - risy_eta2(angle)  *     - risy_eta2(angle)
# Line 55  c      include 'calib.f' Line 177  c      include 'calib.f'
177  *     according to the angle  *     according to the angle
178  *--------------------------------------------------------------  *--------------------------------------------------------------
179        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
180        include 'level1.f'        include 'level1.f'
181          include 'calib.f'
182    
183  c$$$      logical DEBUG        riseta = 0
 c$$$      common/dbg/DEBUG  
         
       ris_eta = 0  
184    
185        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
186          if(mod(iview,2).eq.1)then !Y-view
187                
188           ris_eta = risy_eta2(angle)  
189           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
190                riseta = risy_eta2(angle)
191             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
192                riseta = risy_cog(angle) !ATTENZIONE!!
193             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
194                riseta = risy_cog(angle) !ATTENZIONE!!
195             else
196                riseta = risy_cog(angle)
197             endif            
198    
199        else                      !X-view        else                      !X-view
200    
201           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
202              ris_eta = risx_eta2(angle)              riseta = risx_eta2(angle)
203           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
204              ris_eta = risx_eta3(angle)              riseta = risx_eta3(angle)
205           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
206              ris_eta = risx_eta4(angle)              riseta = risx_eta4(angle)
207           elseif(abs(angle).gt.21.)then           else
208              ris_eta = risx_eta4(21.)              riseta = risx_cog(angle)
209           endif           endif            
210                            
211        endif        endif
212    
213  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        print*,'---- ',riseta,iview,angle
 c$$$     $     ,' -->',ris_eta  
   
214    
215   100  return   100  return
216        end        end
# Line 103  c$$$     $     ,' -->',ris_eta Line 229  c$$$     $     ,' -->',ris_eta
229    
230        include 'commontracker.f'        include 'commontracker.f'
231        include 'level1.f'        include 'level1.f'
232  *      include 'calib.f'        include 'calib.f'
233        fbad_eta = 0        fbad_eta = 0
234    
235        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
236                
237           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
238                fbad_eta = fbad_cog(2,ic)
239             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
240                fbad_eta = fbad_cog(3,ic)
241             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
242                fbad_eta = fbad_cog(4,ic)
243             else
244                fbad_eta = fbad_cog(4,ic)
245             endif            
246    
247        else                      !X-view        else                      !X-view
248    
249           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
250              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
251           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
252              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
253           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
254              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
255           endif           else
256                fbad_eta = fbad_cog(4,ic)
257             endif            
258                            
259        endif        endif
260    
# Line 126  c$$$     $     ,' -->',ris_eta Line 262  c$$$     $     ,' -->',ris_eta
262        end        end
263    
264  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
265  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)  
266  *--------------------------------------------------------------  *--------------------------------------------------------------
267  *     this function returns  *     this function returns
268  *  *
# Line 149  c      real function pfa_eta2(cog2,view, Line 281  c      real function pfa_eta2(cog2,view,
281        real cog2,angle        real cog2,angle
282        integer iview,lad        integer iview,lad
283    
284  c      logical DEBUG        iview = VIEW(ic)            
285  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
286          cog2 = cog(2,ic)          
287        iview = VIEW(ic)              !(1)        pfaeta2=cog2
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfa_eta2=cog2  
288    
289  *     find angular bin  *     find angular bin
290  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
291        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
292           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
293              iangle=iang              iangle=iang
294              goto 98              goto 98
295           endif           endif
296        enddo        enddo
297        if(DEBUG)        if(DEBUG)
298       $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
299        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
300        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
301   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 177  c$$$*     find extremes of interpolation Line 305  c$$$*     find extremes of interpolation
305  c$$$      iflag=0  c$$$      iflag=0
306  c$$$*     --------------------------------  c$$$*     --------------------------------
307  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
308  c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2  c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2
309  c$$$*         goto 100  c$$$*         goto 100
310  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
311  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 358  c            print*,'-----',x1,x2,y1,y2
358        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
359        BB=y1-AA*x1        BB=y1-AA*x1
360    
361        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
362        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
363    
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  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
369  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
370  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
371  c$$$      endif  c$$$      endif
372    
373        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'
374       $     ,cog2-iadd,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
375    
376    
377   100  return   100  return
378        end        end
379    
380  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
381  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)  
382  *--------------------------------------------------------------  *--------------------------------------------------------------
383  *     this function returns  *     this function returns
384  *  *
# Line 273  c      real function pfa_eta3(cog3,view, Line 397  c      real function pfa_eta3(cog3,view,
397        real cog3,angle        real cog3,angle
398        integer iview,lad        integer iview,lad
399    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
400    
401        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
402        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
403        cog3 = cog(3,ic)             !(1)        cog3 = cog(3,ic)            
404        pfa_eta3=cog3        pfaeta3=cog3
405    
406  *     find angular bin  *     find angular bin
407  *     (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 413  c         print*,'~~~~~~~~~~~~ ',iang,an
413           endif           endif
414        enddo        enddo
415        if(DEBUG)        if(DEBUG)
416       $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
417        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
418        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
419   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 354  c            print*,'-----',x1,x2,y1,y2 Line 475  c            print*,'-----',x1,x2,y1,y2
475        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
476        BB=y1-AA*x1        BB=y1-AA*x1
477    
478        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
479        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
480    
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  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
486  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
487  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
488  c$$$      endif  c$$$      endif
489    
490        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'
491       $     ,cog3-iadd,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
492    
493   100  return   100  return
494        end        end
495    
496  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
497  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)  
498  *--------------------------------------------------------------  *--------------------------------------------------------------
499  *     this function returns  *     this function returns
500  *  *
# Line 396  c      real function pfa_eta4(cog4,view, Line 513  c      real function pfa_eta4(cog4,view,
513        real cog4,angle        real cog4,angle
514        integer iview,lad        integer iview,lad
515    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
516    
517        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
518        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
519        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
520        pfa_eta4=cog4        pfaeta4=cog4
521    
522  *     find angular bin  *     find angular bin
523  *     (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 529  c         print*,'~~~~~~~~~~~~ ',iang,an
529           endif           endif
530        enddo        enddo
531        if(DEBUG)        if(DEBUG)
532       $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
533        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
534        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
535   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 476  c            print*,'-----',x1,x2,y1,y2 Line 591  c            print*,'-----',x1,x2,y1,y2
591        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
592        BB=y1-AA*x1        BB=y1-AA*x1
593    
594        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
595        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
596    
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  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
602  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
603  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
604  c$$$      endif  c$$$      endif
605    
606        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'
607       $     ,cog4-iadd,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
608    
609   100  return   100  return
610        end        end
# Line 613  c      print *,ncog,ic,cog,'//////////// Line 728  c      print *,ncog,ic,cog,'////////////
728        include 'calib.f'        include 'calib.f'
729        include 'level1.f'        include 'level1.f'
730                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
731    
732    
733        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 655  c      common/dbg/DEBUG Line 768  c      common/dbg/DEBUG
768                    
769           COG = 0.           COG = 0.
770                    
771           if(ncog.eq.2)then  c         print*,'## ',sl2,sl1,sc,sr1,sr2
772    
773             if(ncog.eq.1)then
774                COG = 0.
775             elseif(ncog.eq.2)then
776              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
777                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
778              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
779                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)            
780              endif              endif
781           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
782              COG = (sr1-sl1)/(sl1+sc+sr1)               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)
783           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
784              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
785                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sl1+sc+sr1).ne.0)
786         $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
787              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
788                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)                  if((sr2+sl1+sc+sr1).ne.0)
789         $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
790              endif              endif
791           else           else
792              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
793              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
794              COG = 0.              COG = 0.
795           endif           endif
796    
# Line 685  ' Line 804  '
804           iv=VIEW(ic)           iv=VIEW(ic)
805           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
806           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
807           istart = INDSTART(IC)           istart = INDSTART(IC)
808           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
809           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
810           COG = 0             COG = 0  
811             SGN = 0.
812           mu  = 0           mu  = 0
813           do i = istart,istop  c         print*,'-------'
814             do i = INDMAX(IC),istart,-1
815              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
816              ivk  = nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
817              is   = nst(MAXS(ic)+ipos)              if(CLSIGNAL(i).ge.cut)then              
818  *            print*,'******************',istart,istop,ipos                 COG = COG + ipos*CLSIGNAL(i)
819  *     $           ,MAXS(ic)+ipos,iv,ivk,is                 SGN = SGN + CLSIGNAL(i)
820              cut  = incut*SIGMA(iv,ivk,is)                 mu = mu + 1
821              if(SIGMA(iv,ivk,is).ne.CLSIGMA(i))  c               print*,ipos,CLSIGNAL(i)
822       $           print*,'cog(0,ic) --> hai fatto qualche cazzata'              else
823                   goto 10
824                endif
825             enddo
826     10      continue
827             do i = INDMAX(IC)+1,istop
828                ipos = i-INDMAX(ic)
829                cut  = incut*CLSIGMA(i)
830              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
831                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
832                   SGN = SGN + CLSIGNAL(i)
833                 mu = mu + 1                 mu = mu + 1
834  c               print*,ipos,CLSIGNAL(i),incut,cut  c               print*,ipos,CLSIGNAL(i)
835                else
836                   goto 20
837              endif              endif
838           enddo           enddo
839           if(DEDX(ic).le.0)then   20      continue
840              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)           if(SGN.le.0)then
841                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
842              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
843              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
844              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
845           else           else
846              COG=COG/DEDX(ic)              COG=COG/SGN
847           endif           endif
848  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'  c         print*,'-------'
 c     $        ,cog  
849                    
850        else        else
851                    
# Line 726  c     $        ,cog Line 856  c     $        ,cog
856    
857        endif        endif
858    
859  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
860    
861        return        return
862        end        end
863    
864  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
865    
866        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
867  *-------------------------------------------------------  *-------------------------------------------------------
868  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 747  c      print *,ncog,ic,cog,'//////////// Line 878  c      print *,ncog,ic,cog,'////////////
878        include 'calib.f'        include 'calib.f'
879    
880        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
881           f  = 4.           si = 8.4  !average good-strip noise
882           si = 8.4           f  = 4.   !average bad-strip noise: f*si
883           incut=incuty           incut=incuty
884        else                      !X-view        else                      !X-view
885           f  = 6.           si = 3.9  !average good-strip noise
886           si = 3.9           f  = 6.   !average bad-strip noise: f*si
887           incut=incutx           incut=incutx
888        endif        endif
889                
# Line 763  c      print *,ncog,ic,cog,'//////////// Line 894  c      print *,ncog,ic,cog,'////////////
894  *     --> signal of the central strip  *     --> signal of the central strip
895           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
896           fsc = 1           fsc = 1
897           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
898             fsc = clsigma(INDMAX(ic))/si
899  *     --> signal of adjacent strips  *     --> signal of adjacent strips
900           sl1 = 0                !left 1           sl1 = 0                !left 1
901           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 771  c      print *,ncog,ic,cog,'//////////// Line 903  c      print *,ncog,ic,cog,'////////////
903       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
904       $        )then       $        )then
905              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
906              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
907  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
908           endif           endif
909    
910           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 782  c            fsl1 = 0 Line 913  c            fsl1 = 0
913       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
914       $        )then       $        )then
915              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
916              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
917  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
918           endif           endif
919           sr1 = 0                !right 1           sr1 = 0                !right 1
920           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 794  c            fsl2 = 0 Line 924  c            fsl2 = 0
924       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
925       $        )then       $        )then
926              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
927              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
928  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
929           endif               endif    
930           sr2 = 0                !right 2           sr2 = 0                !right 2
931           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 806  c            fsr1 = 0 Line 935  c            fsr1 = 0
935       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
936       $        )then       $        )then
937              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
938              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
939  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
940           endif           endif
941    
942    
943    
944  ************************************************************  ************************************************************
945  *     COG computation  *     COG2-3-4 computation
946  ************************************************************  ************************************************************
947    
948  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
949                    
950           COG = 0.           vCOG = cog(ncog,ic)!0.
951                    
952           if(ncog.eq.2)then           if(ncog.eq.2)then
953              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
954                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
955                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
956                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
957              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
958                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
959                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
960                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
961              endif              endif
962           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
963              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
964              fbad_cog =              fbad_cog =
965       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
966              fbad_cog =              fbad_cog =
967       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
968           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
969              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
970                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
971                 fbad_cog =                 fbad_cog =
972       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
973       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
974                 fbad_cog =                 fbad_cog =
975       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
976       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
977              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
978                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
979                 fbad_cog =                 fbad_cog =
980       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
981       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
982                 fbad_cog =                 fbad_cog =
983       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
984       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
985              endif              endif
986           else           else
987              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
988              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
989              COG = 0.  c            COG = 0.
990           endif           endif
991                    
992        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
993    *     =========================
994    *     COG computation
995    *     =========================
996                    
997           iv=VIEW(ic)           vCOG = cog(0,ic)
998           istart=INDSTART(IC)  
999           istop=TOTCLLENGTH           iv     = VIEW(ic)
1000           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1001           COG=0.           istop  = TOTCLLENGTH
1002           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1003           SDE=0.           SGN = 0.
1004           do i=istart,istop           SNU = 0.
1005              ipos=i-INDMAX(ic)           SDE = 0.
1006              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1007              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1008              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1009    c$$$            if(CLSIGNAL(i).gt.cut)then
1010    c$$$               COG = COG + ipos*CLSIGNAL(i)
1011    c$$$               SGN = SGN + CLSIGNAL(i)
1012    c$$$            else
1013    c$$$               goto 10
1014    c$$$            endif
1015    c$$$         enddo
1016    c$$$ 10      continue
1017    c$$$         do i=INDMAX(IC)+1,istop
1018    c$$$            ipos = i-INDMAX(ic)
1019    c$$$            cut  = incut*CLSIGMA(i)
1020    c$$$            if(CLSIGNAL(i).gt.cut)then
1021    c$$$               COG = COG + ipos*CLSIGNAL(i)
1022    c$$$               SGN = SGN + CLSIGNAL(i)
1023    c$$$            else
1024    c$$$               goto 20
1025    c$$$            endif
1026    c$$$         enddo
1027    c$$$ 20      continue
1028    c$$$         if(SGN.le.0)then
1029    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1030    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1031    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1032    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1033    c$$$         else
1034    c$$$            COG=COG/SGN
1035    c$$$         endif
1036    
1037             do i=INDMAX(IC),istart,-1
1038                ipos = i-INDMAX(ic)
1039                cut  = incut*CLSIGMA(i)
1040              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1041                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1042              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1043                   SDE = SDE + (ipos-vCOG)**2
1044                else
1045                   goto 10
1046                endif            
1047           enddo           enddo
1048           COG=COG/DEDX(ic)   10      continue
1049           do i=istart,istop           do i=INDMAX(IC)+1,istop
1050              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1051              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1052              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1053                 fs=1                 fs = clsigma(i)/si
1054                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1055                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1056                 SDE = SDE + (ipos-COG)**2              else
1057                   goto 20
1058              endif                          endif            
1059           enddo           enddo
1060           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1061             if(SDE.ne.0)then
1062                FBAD_COG=SNU/SDE
1063             else
1064                
1065             endif
1066    
1067        else        else
1068                                        
# Line 1096  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1265  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1265       +/       +/
1266    
1267        V(1)= abs(x)        V(1)= abs(x)
1268          if(V(1).gt.20.)V(1)=20.
1269    
1270        HQUADF = 0.        HQUADF = 0.
1271        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1186  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1356  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1356       +/       +/
1357    
1358        V(1) =  abs(x)        V(1) =  abs(x)
1359          if(V(1).gt.20.)V(1)=20.
1360    
1361        HQUADF = 0.        HQUADF = 0.
1362        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1363           HQDJ = 0.           HQDJ = 0.
# Line 1274  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1446  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1446       +/       +/
1447    
1448        V(1)=abs(x)        V(1)=abs(x)
1449          if(V(1).gt.20.)V(1)=20.
1450    
1451        HQUADF = 0.        HQUADF = 0.
1452        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1345  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1518  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1518       +/       +/
1519    
1520        v(1)= abs(x)        v(1)= abs(x)
1521          if(V(1).gt.20.)V(1)=20.
1522                
1523        HQUADF = 0.        HQUADF = 0.
1524        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1411  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1585  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1585       +/       +/
1586                
1587        V(1)=abs(x)        V(1)=abs(x)
1588          if(V(1).gt.20.)V(1)=20.
1589    
1590        HQUADF = 0.        HQUADF = 0.
1591        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1491  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1666  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1666       +/       +/
1667    
1668        V(1)=abs(x)        V(1)=abs(x)
1669          if(V(1).gt.20.)V(1)=20.
1670    
1671        HQUADF = 0.        HQUADF = 0.
1672        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1508  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1684  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1684        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1685    
1686        END        END
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23