/[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.5 by pam-fi, Thu Oct 26 16:22:37 2006 UTC revision 1.17 by pam-fi, Fri Aug 17 16:08:15 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.5)CPFA='ETAL'
14          if(ipfa.eq.10)CPFA='COG'
15          if(ipfa.eq.11)CPFA='COG1'
16          if(ipfa.eq.12)CPFA='COG2'
17          if(ipfa.eq.13)CPFA='COG3'
18          if(ipfa.eq.14)CPFA='COG4'
19          
20          end
21    
22  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
23  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
24  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms
# Line 6  Line 27 
27  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
28    
29    
30          integer function npfastrips(ic,PFA,angle)
31    *--------------------------------------------------------------
32    *     thid function returns the number of strips used
33    *     to evaluate the position of a cluster, according to the p.f.a.
34    *--------------------------------------------------------------
35          include 'commontracker.f'
36          include 'level1.f'
37          include 'calib.f'
38    
39          character*4 usedPFA,PFA
40    
41    
42          usedPFA=PFA
43    
44          npfastrips=0
45    
46          if(usedPFA.eq.'COG1')npfastrips=1
47          if(usedPFA.eq.'COG2')npfastrips=2
48          if(usedPFA.eq.'COG3')npfastrips=3
49          if(usedPFA.eq.'COG4')npfastrips=4
50          if(usedPFA.eq.'ETA2')npfastrips=2
51          if(usedPFA.eq.'ETA3')npfastrips=3
52          if(usedPFA.eq.'ETA4')npfastrips=4
53    *     ----------------------------------------------------------------
54          if(usedPFA.eq.'ETA')then
55    c         print*,VIEW(ic),angle
56             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
57                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
58                   npfastrips=2
59                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
60                   npfastrips=3
61                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
62                   npfastrips=4
63                else
64                   npfastrips=4
65    c               usedPFA='COG'
66                endif                        
67             else                   !X-view
68                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
69                   npfastrips=2
70                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
71                   npfastrips=3
72                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
73                   npfastrips=4
74                else
75                   npfastrips=4
76    c               usedPFA='COG'
77                endif                        
78             endif
79          endif
80    *     ----------------------------------------------------------------
81          if(usedPFA.eq.'COG')then
82    
83             iv=VIEW(ic)
84             if(mod(iv,2).eq.1)incut=incuty
85             if(mod(iv,2).eq.0)incut=incutx
86             istart = INDSTART(IC)
87             istop  = TOTCLLENGTH
88             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
89             mu  = 0
90             do i = INDMAX(IC),istart,-1
91                ipos = i-INDMAX(ic)
92                cut  = incut*CLSIGMA(i)
93                if(CLSIGNAL(i).ge.cut)then
94                   mu = mu + 1
95                   print*,i,mu
96                else
97                   goto 10
98                endif
99             enddo
100     10      continue
101             do i = INDMAX(IC)+1,istop
102                ipos = i-INDMAX(ic)
103                cut  = incut*CLSIGMA(i)
104                if(CLSIGNAL(i).ge.cut)then
105                   mu = mu + 1
106                   print*,i,mu
107                else
108                   goto 20
109                endif
110             enddo
111     20      continue
112             npfastrips=mu
113    
114          endif
115    *     ----------------------------------------------------------------
116    
117    c      print*,pfastrips
118    
119          return
120          end
121    
122  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
123        real function pfaeta(ic,angle)        real function pfaeta(ic,angle)
124  *--------------------------------------------------------------  *--------------------------------------------------------------
# Line 17  Line 130 
130  *     according to the angle  *     according to the angle
131  *--------------------------------------------------------------  *--------------------------------------------------------------
132        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
133        include 'level1.f'        include 'level1.f'
134          include 'calib.f'
135                
136        pfaeta = 0        pfaeta = 0
137    
138        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
139                
140           pfaeta = pfaeta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
141                  pfaeta = pfaeta2(ic,angle)
142             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
143                pfaeta = pfaeta3(ic,angle)
144             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
145                pfaeta = pfaeta4(ic,angle)
146             else
147                pfaeta = cog(4,ic)
148             endif            
149    
150        else                      !X-view        else                      !X-view
151    
152           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
153              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
154           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
155              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
156           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
157              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
158           endif           else
159                pfaeta = cog(4,ic)
160             endif            
161                            
162        endif        endif
163                
 c      print*,'pfaeta ',pfaeta, angle  
   
164   100  return   100  return
165        end        end
166    
167  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
168        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
169  *--------------------------------------------------------------  *--------------------------------------------------------------
170  *     this function returns the average spatial resolution  *     this function returns the position (in strip units)
 *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  
171  *     it calls:  *     it calls:
172  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
173  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
174  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
175  *     according to the angle  *     according to the angle
176  *--------------------------------------------------------------  *--------------------------------------------------------------
177        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
178        include 'level1.f'        include 'level1.f'
179          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
180                
181        ris_eta = 0        pfaeta = 0
182    
183        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
184                
185           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
186           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              pfaeta = pfaeta2(ic,angle)+pfacorr(ic,angle)
187             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
188                pfaeta = pfaeta3(ic,angle)+pfacorr(ic,angle)
189             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
190                pfaeta = pfaeta4(ic,angle)+pfacorr(ic,angle)
191             else
192                pfaeta = cog(4,ic)
193             endif            
194    
195        else                      !X-view        else                      !X-view
196    
197           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
198              ris_eta = risx_eta2(angle)              pfaeta = pfaeta2(ic,angle)+pfacorr(ic,angle)
199           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
200              ris_eta = risx_eta3(angle)              pfaeta = pfaeta3(ic,angle)+pfacorr(ic,angle)
201           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
202              ris_eta = risx_eta4(angle)              pfaeta = pfaeta4(ic,angle)+pfacorr(ic,angle)
203           elseif(abs(angle).gt.21.)then           else
204              ris_eta = risx_eta4(21.)              pfaeta = cog(4,ic)
205           endif           endif            
206                            
207        endif        endif
208          
209     100  return
210          end
211    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
212    c      real function riseta(ic,angle)
213          real function riseta(iview,angle)
214    *--------------------------------------------------------------
215    *     this function returns the average spatial resolution
216    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
217    *     it calls:
218    *     - risxeta2(angle)
219    *     - risyeta2(angle)
220    *     - risxeta3(angle)
221    *     - risxeta4(angle)
222    *     according to the angle
223    *--------------------------------------------------------------
224          include 'commontracker.f'
225          include 'level1.f'
226          include 'calib.f'
227    
228  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        riseta = 0
229  c$$$     $     ,' -->',ris_eta  
230    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
231          if(mod(iview,2).eq.1)then !Y-view
232          
233    
234             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
235                riseta = risyeta2(angle)
236             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
237                riseta = risy_cog(angle) !ATTENZIONE!!
238             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
239                riseta = risy_cog(angle) !ATTENZIONE!!
240             else
241                riseta = risy_cog(angle)
242             endif            
243    
244          else                      !X-view
245    
246             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
247                riseta = risxeta2(angle)
248             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
249                riseta = risxeta3(angle)
250             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
251                riseta = risxeta4(angle)
252             else
253                riseta = risx_cog(angle)
254             endif            
255                
256          endif
257    
258    
259   100  return   100  return
# Line 104  c$$$     $     ,' -->',ris_eta Line 273  c$$$     $     ,' -->',ris_eta
273    
274        include 'commontracker.f'        include 'commontracker.f'
275        include 'level1.f'        include 'level1.f'
276  *      include 'calib.f'        include 'calib.f'
277        fbad_eta = 0        fbad_eta = 0
278    
279        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
280                
281           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
282                fbad_eta = fbad_cog(2,ic)
283             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
284                fbad_eta = fbad_cog(3,ic)
285             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
286                fbad_eta = fbad_cog(4,ic)
287             else
288                fbad_eta = fbad_cog(4,ic)
289             endif            
290    
291        else                      !X-view        else                      !X-view
292    
293           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
294              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
295           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
296              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
297           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
298              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
299           endif           else
300                fbad_eta = fbad_cog(4,ic)
301             endif            
302                            
303        endif        endif
304    
# Line 127  c$$$     $     ,' -->',ris_eta Line 306  c$$$     $     ,' -->',ris_eta
306        end        end
307    
308  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
309        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
310  *--------------------------------------------------------------  *--------------------------------------------------------------
311  *     this function returns  *     this function returns
# Line 150  c      real function pfaeta2(cog2,view,l Line 325  c      real function pfaeta2(cog2,view,l
325        real cog2,angle        real cog2,angle
326        integer iview,lad        integer iview,lad
327    
328  c      logical DEBUG        iview = VIEW(ic)            
329  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
330          cog2 = cog(2,ic)          
       iview = VIEW(ic)              !(1)  
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
331        pfaeta2=cog2        pfaeta2=cog2
332    
333  *     find angular bin  *     find angular bin
334  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
335        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
336           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
337              iangle=iang              iangle=iang
338              goto 98              goto 98
# Line 251  c$$$      endif Line 422  c$$$      endif
422        end        end
423    
424  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
425        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
426  *--------------------------------------------------------------  *--------------------------------------------------------------
427  *     this function returns  *     this function returns
# Line 274  c      real function pfaeta3(cog3,view,l Line 441  c      real function pfaeta3(cog3,view,l
441        real cog3,angle        real cog3,angle
442        integer iview,lad        integer iview,lad
443    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
444    
445          iview = VIEW(ic)            
446        iview = VIEW(ic)              !(1)        lad = nld(MAXS(ic),VIEW(ic))
447        lad = nld(MAXS(ic),VIEW(ic)) !(1)        cog3 = cog(3,ic)            
       cog3 = cog(3,ic)             !(1)  
448        pfaeta3=cog3        pfaeta3=cog3
449    
450  *     find angular bin  *     find angular bin
# Line 374  c$$$      endif Line 538  c$$$      endif
538        end        end
539    
540  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
541  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta4(cog4,view,lad,angle)  
       real function pfaeta4(ic,angle) !(1)  
542  *--------------------------------------------------------------  *--------------------------------------------------------------
543  *     this function returns  *     this function returns
544  *  *
# Line 397  c      real function pfaeta4(cog4,view,l Line 557  c      real function pfaeta4(cog4,view,l
557        real cog4,angle        real cog4,angle
558        integer iview,lad        integer iview,lad
559    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
560    
561        iview = VIEW(ic)             !(1)        iview = VIEW(ic)            
562        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
563        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
564        pfaeta4=cog4        pfaeta4=cog4
565    
566  *     find angular bin  *     find angular bin
# Line 498  c$$$      endif Line 656  c$$$      endif
656    
657    
658  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
659        real function cog0(ncog,ic)  c$$$      real function cog0(ncog,ic)
660  *-------------------------------------------------  c$$$*-------------------------------------------------
661  *     this function returns  c$$$*     this function returns
662  *  c$$$*
663  *     - the Center-Of-Gravity of the cluster IC  c$$$*     - the Center-Of-Gravity of the cluster IC
664  *     evaluated using NCOG strips,  c$$$*     evaluated using NCOG strips,
665  *     calculated relative to MAXS(IC)  c$$$*     calculated relative to MAXS(IC)
666  *      c$$$*    
667  *     - zero in case that not  enough strips  c$$$*     - zero in case that not  enough strips
668  *     have a positive signal  c$$$*     have a positive signal
669  *      c$$$*    
670  *     NOTE:  c$$$*     NOTE:
671  *     This is the old definition, used by Straulino.  c$$$*     This is the old definition, used by Straulino.
672  *     The new routine, according to Landi,  c$$$*     The new routine, according to Landi,
673  *     is COG(NCOG,IC)  c$$$*     is COG(NCOG,IC)
674  *-------------------------------------------------  c$$$*-------------------------------------------------
675    c$$$
676    c$$$
677        include 'commontracker.f'  c$$$      include 'commontracker.f'
678        include 'level1.f'  c$$$      include 'level1.f'
679          c$$$      
680  *     --> signal of the central strip  c$$$*     --> signal of the central strip
681        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
682    c$$$
683  *     signal of adjacent strips  c$$$*     signal of adjacent strips
684  *     --> left  c$$$*     --> left
685        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
686        if(  c$$$      if(
687       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
688       $     )  c$$$     $     )
689       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
690    c$$$
691        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
692        if(  c$$$      if(
693       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
694       $     )  c$$$     $     )
695       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
696    c$$$
697  *     --> right  c$$$*     --> right
698        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
699        if(  c$$$      if(
700       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
701       $     .or.  c$$$     $     .or.
702       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
703       $     )  c$$$     $     )
704       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
705    c$$$
706        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
707        if(  c$$$      if(
708       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
709       $     .or.  c$$$     $     .or.
710       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
711       $     )  c$$$     $     )
712       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
713          c$$$      
714  ************************************************************  c$$$************************************************************
715  *     COG computation  c$$$*     COG computation
716  ************************************************************  c$$$************************************************************
717    c$$$
718  c      print*,sl2,sl1,sc,sr1,sr2  c$$$c      print*,sl2,sl1,sc,sr1,sr2
719    c$$$
720        COG = 0.  c$$$      COG = 0.
721            c$$$        
722        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
723            c$$$        
724           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
725              COG = -sl1/(sl1+sc)          c$$$            COG = -sl1/(sl1+sc)        
726           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
727              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
728           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
729              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
730           else  c$$$         else
731              COG = 0.  c$$$            COG = 0.
732           endif  c$$$         endif
733            c$$$        
734        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
735    c$$$
736           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
737              COG = sr1/(sc+sr1)              c$$$            COG = sr1/(sc+sr1)            
738           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
739              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
740           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
741              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
742           else  c$$$         else
743              COG = 0.  c$$$            COG = 0.
744           endif  c$$$         endif
745    c$$$
746        endif  c$$$      endif
747    c$$$
748        COG0 = COG  c$$$      COG0 = COG
749    c$$$
750  c      print *,ncog,ic,cog,'/////////////'  c$$$c      print *,ncog,ic,cog,'/////////////'
751    c$$$
752        return  c$$$      return
753        end  c$$$      end
754    
755  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
756        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 614  c      print *,ncog,ic,cog,'//////////// Line 772  c      print *,ncog,ic,cog,'////////////
772        include 'calib.f'        include 'calib.f'
773        include 'level1.f'        include 'level1.f'
774                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
775    
776    
777        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 626  c      common/dbg/DEBUG Line 782  c      common/dbg/DEBUG
782  *     --> signal of the central strip  *     --> signal of the central strip
783           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
784  *     signal of adjacent strips  *     signal of adjacent strips
785           sl1 = 0                !left 1           sl1 = -9999.           !left 1
786           if(           if(
787       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
788       $        )       $        )
789       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
790                    
791           sl2 = 0                !left 2           sl2 = -9999.           !left 2
792           if(           if(
793       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
794       $        )       $        )
795       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
796                    
797           sr1 = 0                !right 1           sr1 = -9999.           !right 1
798           if(           if(
799       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
800       $        .or.       $        .or.
# Line 646  c      common/dbg/DEBUG Line 802  c      common/dbg/DEBUG
802       $        )       $        )
803       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
804                    
805           sr2 = 0                !right 2           sr2 = -9999.           !right 2
806           if(           if(
807       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
808       $        .or.       $        .or.
# Line 656  c      common/dbg/DEBUG Line 812  c      common/dbg/DEBUG
812                    
813           COG = 0.           COG = 0.
814                    
815    c         print*,'## ',sl2,sl1,sc,sr1,sr2
816    
817    c     ==============================================================
818           if(ncog.eq.1)then           if(ncog.eq.1)then
819              COG = 0.              COG = 0.
820                if(sr1.gt.sc)cog=1. !NEW
821                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
822    c     ==============================================================
823           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
824              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
825                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
826              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
827                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
828              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
829                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
830         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
831                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
832         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
833                endif
834    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
835    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
836    c     ==============================================================
837           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
838              COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
839    c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
840    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
841    c     ==============================================================
842           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
843              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
844                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sl1+sc+sr1).ne.0)
845              elseif(sl2.le.sr2)then       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
846                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)              elseif(sl2.lt.sr2)then
847                   if((sr2+sl1+sc+sr1).ne.0)
848         $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
849                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
850                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
851         $              .and.(sl2+sl1+sc+sr1).ne.0 )
852         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
853                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
854         $              .and.(sr2+sl1+sc+sr1).ne.0 )
855         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
856              endif              endif
857           else           else
858              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
859              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
860              COG = 0.              COG = 0.
861           endif           endif
862    
# Line 688  ' Line 870  '
870           iv=VIEW(ic)           iv=VIEW(ic)
871           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
872           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
873           istart = INDSTART(IC)           istart = INDSTART(IC)
874           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
875           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
876           COG = 0             COG = 0  
877             SGN = 0.
878           mu  = 0           mu  = 0
879           do i = istart,istop  c         print*,'-------'
880             do i = INDMAX(IC),istart,-1
881                ipos = i-INDMAX(ic)
882                cut  = incut*CLSIGMA(i)
883                if(CLSIGNAL(i).ge.cut)then              
884                   COG = COG + ipos*CLSIGNAL(i)
885                   SGN = SGN + CLSIGNAL(i)
886                   mu = mu + 1
887    c               print*,ipos,CLSIGNAL(i)
888                else
889                   goto 10
890                endif
891             enddo
892     10      continue
893             do i = INDMAX(IC)+1,istop
894              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
 cc            ivk  = nvk(MAXS(ic)+ipos)  
 cc            is   = nst(MAXS(ic)+ipos)  
 *            print*,'******************',istart,istop,ipos  
 *     $           ,MAXS(ic)+ipos,iv,ivk,is  
 cc            cut  = incut*SIGMA(iv,ivk,is)  
895              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
 c            if(SIGMA(iv,ivk,is).ne.CLSIGMA(i))  
 c     $           print*,'cog(0,ic) --> hai fatto qualche cazzata'  
896              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
897                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
898                   SGN = SGN + CLSIGNAL(i)
899                 mu = mu + 1                 mu = mu + 1
900  c               print*,ipos,CLSIGNAL(i),incut,cut  c               print*,ipos,CLSIGNAL(i)
901                else
902                   goto 20
903              endif              endif
904           enddo           enddo
905           if(DEDX(ic).le.0)then   20      continue
906              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)           if(SGN.le.0)then
907                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
908              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
909              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
910              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
911           else           else
912              COG=COG/DEDX(ic)              COG=COG/SGN
913           endif           endif
914  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'  c         print*,'-------'
 c     $        ,cog  
915                    
916        else        else
917                    
# Line 730  c     $        ,cog Line 922  c     $        ,cog
922    
923        endif        endif
924    
925  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
926    
927        return        return
928        end        end
929    
930  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
931    
932        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
933  *-------------------------------------------------------  *-------------------------------------------------------
934  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 751  c      print *,ncog,ic,cog,'//////////// Line 944  c      print *,ncog,ic,cog,'////////////
944        include 'calib.f'        include 'calib.f'
945    
946        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
947           f  = 4.           si = 8.4  !average good-strip noise
948           si = 8.4           f  = 4.   !average bad-strip noise: f*si
949           incut=incuty           incut=incuty
950        else                      !X-view        else                      !X-view
951           f  = 6.           si = 3.9  !average good-strip noise
952           si = 3.9           f  = 6.   !average bad-strip noise: f*si
953           incut=incutx           incut=incutx
954        endif        endif
955                
# Line 767  c      print *,ncog,ic,cog,'//////////// Line 960  c      print *,ncog,ic,cog,'////////////
960  *     --> signal of the central strip  *     --> signal of the central strip
961           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
962           fsc = 1           fsc = 1
963           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
964             fsc = clsigma(INDMAX(ic))/si
965  *     --> signal of adjacent strips  *     --> signal of adjacent strips
966           sl1 = 0                !left 1           sl1 = 0                !left 1
967           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 775  c      print *,ncog,ic,cog,'//////////// Line 969  c      print *,ncog,ic,cog,'////////////
969       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
970       $        )then       $        )then
971              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
972              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
973  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
974           endif           endif
975    
976           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 786  c            fsl1 = 0 Line 979  c            fsl1 = 0
979       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
980       $        )then       $        )then
981              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
982              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
983  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
984           endif           endif
985           sr1 = 0                !right 1           sr1 = 0                !right 1
986           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 798  c            fsl2 = 0 Line 990  c            fsl2 = 0
990       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
991       $        )then       $        )then
992              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
993              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
994  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
995           endif               endif    
996           sr2 = 0                !right 2           sr2 = 0                !right 2
997           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 810  c            fsr1 = 0 Line 1001  c            fsr1 = 0
1001       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1002       $        )then       $        )then
1003              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
1004              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
1005  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
1006           endif           endif
1007    
1008    
1009    
1010  ************************************************************  ************************************************************
1011  *     COG computation  *     COG2-3-4 computation
1012  ************************************************************  ************************************************************
1013    
1014  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
1015                    
1016           COG = 0.           vCOG = cog(ncog,ic)!0.
1017                    
1018           if(ncog.eq.2)then           if(ncog.eq.2)then
1019              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1020                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
1021                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1022                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1023              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
1024                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
1025                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1026                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1027              endif              endif
1028           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1029              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
1030              fbad_cog =              fbad_cog =
1031       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1032              fbad_cog =              fbad_cog =
1033       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1034           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1035              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1036                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1037                 fbad_cog =                 fbad_cog =
1038       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1039       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1040                 fbad_cog =                 fbad_cog =
1041       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1042       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1043              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1044                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1045                 fbad_cog =                 fbad_cog =
1046       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1047       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1048                 fbad_cog =                 fbad_cog =
1049       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1050       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1051              endif              endif
1052           else           else
1053              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1054              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1055              COG = 0.  c            COG = 0.
1056           endif           endif
1057                    
1058        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1059    *     =========================
1060    *     COG computation
1061    *     =========================
1062                    
1063           iv=VIEW(ic)           vCOG = cog(0,ic)
1064           istart=INDSTART(IC)  
1065           istop=TOTCLLENGTH           iv     = VIEW(ic)
1066           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1067           COG=0.           istop  = TOTCLLENGTH
1068           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1069           SDE=0.           SGN = 0.
1070           do i=istart,istop           SNU = 0.
1071              ipos=i-INDMAX(ic)           SDE = 0.
1072              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1073              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1074              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1075    c$$$            if(CLSIGNAL(i).gt.cut)then
1076    c$$$               COG = COG + ipos*CLSIGNAL(i)
1077    c$$$               SGN = SGN + CLSIGNAL(i)
1078    c$$$            else
1079    c$$$               goto 10
1080    c$$$            endif
1081    c$$$         enddo
1082    c$$$ 10      continue
1083    c$$$         do i=INDMAX(IC)+1,istop
1084    c$$$            ipos = i-INDMAX(ic)
1085    c$$$            cut  = incut*CLSIGMA(i)
1086    c$$$            if(CLSIGNAL(i).gt.cut)then
1087    c$$$               COG = COG + ipos*CLSIGNAL(i)
1088    c$$$               SGN = SGN + CLSIGNAL(i)
1089    c$$$            else
1090    c$$$               goto 20
1091    c$$$            endif
1092    c$$$         enddo
1093    c$$$ 20      continue
1094    c$$$         if(SGN.le.0)then
1095    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1096    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1097    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1098    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1099    c$$$         else
1100    c$$$            COG=COG/SGN
1101    c$$$         endif
1102    
1103             do i=INDMAX(IC),istart,-1
1104                ipos = i-INDMAX(ic)
1105                cut  = incut*CLSIGMA(i)
1106              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1107                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1108              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1109                   SDE = SDE + (ipos-vCOG)**2
1110                else
1111                   goto 10
1112                endif            
1113           enddo           enddo
1114           COG=COG/DEDX(ic)   10      continue
1115           do i=istart,istop           do i=INDMAX(IC)+1,istop
1116              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1117              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1118              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1119                 fs=1                 fs = clsigma(i)/si
1120                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1121                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1122                 SDE = SDE + (ipos-COG)**2              else
1123                   goto 20
1124              endif                          endif            
1125           enddo           enddo
1126           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1127             if(SDE.ne.0)then
1128                FBAD_COG=SNU/SDE
1129             else
1130                
1131             endif
1132    
1133        else        else
1134                                        
# Line 1028  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1259  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1259    
1260  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1261    
1262        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1263    
1264        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1265        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1100  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1331  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1331       +/       +/
1332    
1333        V(1)= abs(x)        V(1)= abs(x)
1334          if(V(1).gt.20.)V(1)=20.
1335    
1336        HQUADF = 0.        HQUADF = 0.
1337        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1114  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1346  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1346     20 CONTINUE     20 CONTINUE
1347        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1348                
1349        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1350    
1351        END        END
1352    
1353  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1354        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1355        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1356        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1357        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1190  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1422  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1422       +/       +/
1423    
1424        V(1) =  abs(x)        V(1) =  abs(x)
1425          if(V(1).gt.20.)V(1)=20.
1426    
1427        HQUADF = 0.        HQUADF = 0.
1428        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1429           HQDJ = 0.           HQDJ = 0.
# Line 1203  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1437  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1437     20 CONTINUE     20 CONTINUE
1438        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1439                
1440        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1441    
1442        END        END
1443  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1444        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1445        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1446        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1447        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1278  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1512  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1512       +/       +/
1513    
1514        V(1)=abs(x)        V(1)=abs(x)
1515          if(V(1).gt.20.)V(1)=20.
1516    
1517        HQUADF = 0.        HQUADF = 0.
1518        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1292  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1527  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1527     20 CONTINUE     20 CONTINUE
1528        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1529                
1530        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1531    
1532        END        END
1533  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1534        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1535        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1536        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1537        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1349  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1584  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1584       +/       +/
1585    
1586        v(1)= abs(x)        v(1)= abs(x)
1587          if(V(1).gt.20.)V(1)=20.
1588                
1589        HQUADF = 0.        HQUADF = 0.
1590        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1363  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1599  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1599     20 CONTINUE     20 CONTINUE
1600        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1601    
1602        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1603    
1604        END        END
1605  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1415  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1651  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1651       +/       +/
1652                
1653        V(1)=abs(x)        V(1)=abs(x)
1654          if(V(1).gt.20.)V(1)=20.
1655    
1656        HQUADF = 0.        HQUADF = 0.
1657        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1495  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1732  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1732       +/       +/
1733    
1734        V(1)=abs(x)        V(1)=abs(x)
1735          if(V(1).gt.20.)V(1)=20.
1736    
1737        HQUADF = 0.        HQUADF = 0.
1738        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1512  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1750  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1750        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1751    
1752        END        END
1753    
1754    
1755  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1756          real function pfacorr(ic,angle) !(1)
1757    *--------------------------------------------------------------
1758    *     this function returns the landi correction for this cluster
1759    *--------------------------------------------------------------
1760          include 'commontracker.f'
1761          include 'calib.f'
1762          include 'level1.f'
1763    
1764          real angle
1765          integer iview,lad
1766    
1767          iview = VIEW(ic)            
1768          lad = nld(MAXS(ic),VIEW(ic))
1769    
1770    *     find angular bin
1771    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1772          do iang=1,nangbin
1773             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1774                iangle=iang
1775                goto 98
1776             endif
1777          enddo
1778          if(DEBUG)
1779         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1780          if(angle.lt.angL(1))iang=1
1781          if(angle.gt.angR(nangbin))iang=nangbin
1782     98   continue                  !jump here if ok
1783    
1784          pfacorr = fcorr(iview,lad,iang)
1785    
1786          if(DEBUG)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1787    
1788    
1789     100  return
1790          end
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.23