/[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.8 by pam-fi, Fri Feb 16 14:56:02 2007 UTC revision 1.19 by pam-fi, Tue Aug 21 15:21:22 2007 UTC
# Line 1  Line 1 
1    
2    
3          subroutine idtoc(ipfa,cpfa)
4          
5          integer ipfa
6          character*10 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    cc            print*,pfaeta2(ic,angle)
143             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
144                pfaeta = pfaeta3(ic,angle)
145             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
146                pfaeta = pfaeta4(ic,angle)
147             else
148                pfaeta = cog(4,ic)
149             endif            
150    
151        else                      !X-view        else                      !X-view
152    
153           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
154              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
155           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
156              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
157           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
158              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
159           endif           else
160                pfaeta = cog(4,ic)
161             endif            
162                            
163        endif        endif
164                
 c      print*,'pfaeta ',pfaeta, angle  
   
165   100  return   100  return
166        end        end
167    
168  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
169        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
170  *--------------------------------------------------------------  *--------------------------------------------------------------
171  *     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))  
172  *     it calls:  *     it calls:
173  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
174  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
175  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
176  *     according to the angle  *     according to the angle
177  *--------------------------------------------------------------  *--------------------------------------------------------------
178        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
179        include 'level1.f'        include 'level1.f'
180          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
181                
182        ris_eta = 0        pfaetal = 0
183    
184        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
185                
186           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
187           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
188    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
189             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
190                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
191             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
192                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
193             else
194                pfaetal = cog(4,ic)
195             endif            
196    
197        else                      !X-view        else                      !X-view
198    
199           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
200              ris_eta = risx_eta2(angle)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
201           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then  cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
202              ris_eta = risx_eta3(angle)           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
203           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then              pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
204              ris_eta = risx_eta4(angle)           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
205           elseif(abs(angle).gt.21.)then              pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
206              ris_eta = risx_eta4(21.)           else
207           endif              pfaetal = cog(4,ic)
208             endif            
209                            
210        endif        endif
211          
212     100  return
213          end
214    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
215    c      real function riseta(ic,angle)
216          real function riseta(iview,angle)
217    *--------------------------------------------------------------
218    *     this function returns the average spatial resolution
219    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
220    *     it calls:
221    *     - risxeta2(angle)
222    *     - risyeta2(angle)
223    *     - risxeta3(angle)
224    *     - risxeta4(angle)
225    *     according to the angle
226    *--------------------------------------------------------------
227          include 'commontracker.f'
228          include 'level1.f'
229          include 'calib.f'
230    
231  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        riseta = 0
232  c$$$     $     ,' -->',ris_eta  
233    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
234          if(mod(iview,2).eq.1)then !Y-view
235          
236    
237             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
238                riseta = risyeta2(angle)
239             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
240                riseta = risy_cog(angle) !ATTENZIONE!!
241             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
242                riseta = risy_cog(angle) !ATTENZIONE!!
243             else
244                riseta = risy_cog(angle)
245             endif            
246    
247          else                      !X-view
248    
249             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
250                riseta = risxeta2(angle)
251             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
252                riseta = risxeta3(angle)
253             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
254                riseta = risxeta4(angle)
255             else
256                riseta = risx_cog(angle)
257             endif            
258                
259          endif
260    
261    
262   100  return   100  return
# Line 104  c$$$     $     ,' -->',ris_eta Line 276  c$$$     $     ,' -->',ris_eta
276    
277        include 'commontracker.f'        include 'commontracker.f'
278        include 'level1.f'        include 'level1.f'
279  *      include 'calib.f'        include 'calib.f'
280        fbad_eta = 0        fbad_eta = 0
281    
282        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
283                
284           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
285                fbad_eta = fbad_cog(2,ic)
286             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
287                fbad_eta = fbad_cog(3,ic)
288             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
289                fbad_eta = fbad_cog(4,ic)
290             else
291                fbad_eta = fbad_cog(4,ic)
292             endif            
293    
294        else                      !X-view        else                      !X-view
295    
296           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
297              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
298           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
299              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
300           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
301              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
302           endif           else
303                fbad_eta = fbad_cog(4,ic)
304             endif            
305                            
306        endif        endif
307    
# Line 127  c$$$     $     ,' -->',ris_eta Line 309  c$$$     $     ,' -->',ris_eta
309        end        end
310    
311  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
312        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
313  *--------------------------------------------------------------  *--------------------------------------------------------------
314  *     this function returns  *     this function returns
# Line 150  c      real function pfaeta2(cog2,view,l Line 328  c      real function pfaeta2(cog2,view,l
328        real cog2,angle        real cog2,angle
329        integer iview,lad        integer iview,lad
330    
331  c      logical DEBUG        iview = VIEW(ic)            
332  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
333          cog2 = cog(2,ic)          
 c      print*,'## pfaeta2 ',ic,angle  
       iview = VIEW(ic)              !(1)  
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
334        pfaeta2=cog2        pfaeta2=cog2
335    
336  *     find angular bin  *     find angular bin
337  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
338        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
339           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
340              iangle=iang              iangle=iang
341              goto 98              goto 98
342           endif           endif
343        enddo        enddo
344        if(DEBUG)        if(DEBUG.EQ.1)
345       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
346        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
347        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
# Line 244  c$$$         pfaeta2=pfaeta2+1.   !temp Line 417  c$$$         pfaeta2=pfaeta2+1.   !temp
417  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
418  c$$$      endif  c$$$      endif
419    
420        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
421       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
422    
423    
# Line 252  c$$$      endif Line 425  c$$$      endif
425        end        end
426    
427  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
428        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
429  *--------------------------------------------------------------  *--------------------------------------------------------------
430  *     this function returns  *     this function returns
# Line 275  c      real function pfaeta3(cog3,view,l Line 444  c      real function pfaeta3(cog3,view,l
444        real cog3,angle        real cog3,angle
445        integer iview,lad        integer iview,lad
446    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
447    
448  c      print*,'## pfaeta3 ',ic,angle        iview = VIEW(ic)            
449          lad = nld(MAXS(ic),VIEW(ic))
450        iview = VIEW(ic)              !(1)        cog3 = cog(3,ic)            
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog3 = cog(3,ic)             !(1)  
451        pfaeta3=cog3        pfaeta3=cog3
452    
453  *     find angular bin  *     find angular bin
# Line 294  c         print*,'~~~~~~~~~~~~ ',iang,an Line 459  c         print*,'~~~~~~~~~~~~ ',iang,an
459              goto 98              goto 98
460           endif           endif
461        enddo        enddo
462        if(DEBUG)        if(DEBUG.EQ.1)
463       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
464        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
465        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
# Line 369  c$$$         pfaeta2=pfaeta2+1.   !temp Line 534  c$$$         pfaeta2=pfaeta2+1.   !temp
534  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
535  c$$$      endif  c$$$      endif
536    
537        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
538       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
539    
540   100  return   100  return
541        end        end
542    
543  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
544  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)  
545  *--------------------------------------------------------------  *--------------------------------------------------------------
546  *     this function returns  *     this function returns
547  *  *
# Line 399  c      real function pfaeta4(cog4,view,l Line 560  c      real function pfaeta4(cog4,view,l
560        real cog4,angle        real cog4,angle
561        integer iview,lad        integer iview,lad
562    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 c      print*,'## pfaeta4 ',ic,angle  
563    
564        iview = VIEW(ic)             !(1)        iview = VIEW(ic)            
565        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
566        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
567        pfaeta4=cog4        pfaeta4=cog4
568    
569  *     find angular bin  *     find angular bin
# Line 418  c         print*,'~~~~~~~~~~~~ ',iang,an Line 575  c         print*,'~~~~~~~~~~~~ ',iang,an
575              goto 98              goto 98
576           endif           endif
577        enddo        enddo
578        if(DEBUG)        if(DEBUG.EQ.1)
579       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
580        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
581        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
# Line 493  c$$$         pfaeta2=pfaeta2+1.   !temp Line 650  c$$$         pfaeta2=pfaeta2+1.   !temp
650  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
651  c$$$      endif  c$$$      endif
652    
653        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
654       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
655    
656   100  return   100  return
# Line 502  c$$$      endif Line 659  c$$$      endif
659    
660    
661  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
662        real function cog0(ncog,ic)  c$$$      real function cog0(ncog,ic)
663  *-------------------------------------------------  c$$$*-------------------------------------------------
664  *     this function returns  c$$$*     this function returns
665  *  c$$$*
666  *     - the Center-Of-Gravity of the cluster IC  c$$$*     - the Center-Of-Gravity of the cluster IC
667  *     evaluated using NCOG strips,  c$$$*     evaluated using NCOG strips,
668  *     calculated relative to MAXS(IC)  c$$$*     calculated relative to MAXS(IC)
669  *      c$$$*    
670  *     - zero in case that not  enough strips  c$$$*     - zero in case that not  enough strips
671  *     have a positive signal  c$$$*     have a positive signal
672  *      c$$$*    
673  *     NOTE:  c$$$*     NOTE:
674  *     This is the old definition, used by Straulino.  c$$$*     This is the old definition, used by Straulino.
675  *     The new routine, according to Landi,  c$$$*     The new routine, according to Landi,
676  *     is COG(NCOG,IC)  c$$$*     is COG(NCOG,IC)
677  *-------------------------------------------------  c$$$*-------------------------------------------------
678    c$$$
679    c$$$
680        include 'commontracker.f'  c$$$      include 'commontracker.f'
681        include 'level1.f'  c$$$      include 'level1.f'
682          c$$$      
683  *     --> signal of the central strip  c$$$*     --> signal of the central strip
684        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
685    c$$$
686  *     signal of adjacent strips  c$$$*     signal of adjacent strips
687  *     --> left  c$$$*     --> left
688        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
689        if(  c$$$      if(
690       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
691       $     )  c$$$     $     )
692       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
693    c$$$
694        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
695        if(  c$$$      if(
696       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
697       $     )  c$$$     $     )
698       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
699    c$$$
700  *     --> right  c$$$*     --> right
701        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
702        if(  c$$$      if(
703       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
704       $     .or.  c$$$     $     .or.
705       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
706       $     )  c$$$     $     )
707       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
708    c$$$
709        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
710        if(  c$$$      if(
711       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
712       $     .or.  c$$$     $     .or.
713       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
714       $     )  c$$$     $     )
715       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
716          c$$$      
717  ************************************************************  c$$$************************************************************
718  *     COG computation  c$$$*     COG computation
719  ************************************************************  c$$$************************************************************
720    c$$$
721  c      print*,sl2,sl1,sc,sr1,sr2  c$$$c      print*,sl2,sl1,sc,sr1,sr2
722    c$$$
723        COG = 0.  c$$$      COG = 0.
724            c$$$        
725        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
726            c$$$        
727           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
728              COG = -sl1/(sl1+sc)          c$$$            COG = -sl1/(sl1+sc)        
729           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
730              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
731           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
732              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
733           else  c$$$         else
734              COG = 0.  c$$$            COG = 0.
735           endif  c$$$         endif
736            c$$$        
737        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
738    c$$$
739           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
740              COG = sr1/(sc+sr1)              c$$$            COG = sr1/(sc+sr1)            
741           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
742              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
743           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
744              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
745           else  c$$$         else
746              COG = 0.  c$$$            COG = 0.
747           endif  c$$$         endif
748    c$$$
749        endif  c$$$      endif
750    c$$$
751        COG0 = COG  c$$$      COG0 = COG
752    c$$$
753  c      print *,ncog,ic,cog,'/////////////'  c$$$c      print *,ncog,ic,cog,'/////////////'
754    c$$$
755        return  c$$$      return
756        end  c$$$      end
757    
758  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
759        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 618  c      print *,ncog,ic,cog,'//////////// Line 775  c      print *,ncog,ic,cog,'////////////
775        include 'calib.f'        include 'calib.f'
776        include 'level1.f'        include 'level1.f'
777                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
778    
779    
780        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 630  c      common/dbg/DEBUG Line 785  c      common/dbg/DEBUG
785  *     --> signal of the central strip  *     --> signal of the central strip
786           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
787  *     signal of adjacent strips  *     signal of adjacent strips
788           sl1 = 0                !left 1           sl1 = -9999.           !left 1
789           if(           if(
790       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
791       $        )       $        )
792       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
793                    
794           sl2 = 0                !left 2           sl2 = -9999.           !left 2
795           if(           if(
796       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
797       $        )       $        )
798       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
799                    
800           sr1 = 0                !right 1           sr1 = -9999.           !right 1
801           if(           if(
802       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
803       $        .or.       $        .or.
# Line 650  c      common/dbg/DEBUG Line 805  c      common/dbg/DEBUG
805       $        )       $        )
806       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
807                    
808           sr2 = 0                !right 2           sr2 = -9999.           !right 2
809           if(           if(
810       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
811       $        .or.       $        .or.
# Line 662  c      common/dbg/DEBUG Line 817  c      common/dbg/DEBUG
817                    
818  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c         print*,'## ',sl2,sl1,sc,sr1,sr2
819    
820    c     ==============================================================
821           if(ncog.eq.1)then           if(ncog.eq.1)then
822              COG = 0.              COG = 0.
823                if(sr1.gt.sc)cog=1. !NEW
824                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
825    c     ==============================================================
826           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
827              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
828                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
829              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
830                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
831              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
832                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
833         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
834                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
835         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
836                endif
837    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
838    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
839    c     ==============================================================
840           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
841               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
842    c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
843    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
844    c     ==============================================================
845           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
846              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
847                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sl1+sc+sr1).ne.0)
848       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
849              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
850                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sl1+sc+sr1).ne.0)
851       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
852                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
853                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
854         $              .and.(sl2+sl1+sc+sr1).ne.0 )
855         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
856                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
857         $              .and.(sr2+sl1+sc+sr1).ne.0 )
858         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
859              endif              endif
860           else           else
861              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
# Line 696  ' Line 873  '
873           iv=VIEW(ic)           iv=VIEW(ic)
874           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
875           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
876           istart = INDSTART(IC)           istart = INDSTART(IC)
877           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
878           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
879           COG = 0             COG = 0  
880             SGN = 0.
881           mu  = 0           mu  = 0
882           do i = istart,istop  c         print*,'-------'
883             do i = INDMAX(IC),istart,-1
884                ipos = i-INDMAX(ic)
885                cut  = incut*CLSIGMA(i)
886                if(CLSIGNAL(i).ge.cut)then              
887                   COG = COG + ipos*CLSIGNAL(i)
888                   SGN = SGN + CLSIGNAL(i)
889                   mu = mu + 1
890    c               print*,ipos,CLSIGNAL(i)
891                else
892                   goto 10
893                endif
894             enddo
895     10      continue
896             do i = INDMAX(IC)+1,istop
897              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
898              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
899              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
900                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
901                   SGN = SGN + CLSIGNAL(i)
902                 mu = mu + 1                 mu = mu + 1
903    c               print*,ipos,CLSIGNAL(i)
904                else
905                   goto 20
906              endif              endif
907           enddo           enddo
908           if(SGNL(ic).le.0)then   20      continue
909              print*,'cog(0,ic) --> ic, dedx ',ic,SGNL(ic)           if(SGN.le.0)then
910                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
911              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
912              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
913              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
914           else           else
915              COG=COG/SGNL(ic)              COG=COG/SGN
916           endif           endif
917    c         print*,'-------'
918                    
919        else        else
920                    
# Line 734  c      print *,'## cog ',ncog,ic,cog,'// Line 931  c      print *,'## cog ',ncog,ic,cog,'//
931        end        end
932    
933  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
934    
935        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
936  *-------------------------------------------------------  *-------------------------------------------------------
937  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 749  c      print *,'## cog ',ncog,ic,cog,'// Line 947  c      print *,'## cog ',ncog,ic,cog,'//
947        include 'calib.f'        include 'calib.f'
948    
949        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
950           f  = 4.           si = 8.4  !average good-strip noise
951           si = 8.4           f  = 4.   !average bad-strip noise: f*si
952           incut=incuty           incut=incuty
953        else                      !X-view        else                      !X-view
954           f  = 6.           si = 3.9  !average good-strip noise
955           si = 3.9           f  = 6.   !average bad-strip noise: f*si
956           incut=incutx           incut=incutx
957        endif        endif
958                
# Line 765  c      print *,'## cog ',ncog,ic,cog,'// Line 963  c      print *,'## cog ',ncog,ic,cog,'//
963  *     --> signal of the central strip  *     --> signal of the central strip
964           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
965           fsc = 1           fsc = 1
966  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
967           if( CLBAD(INDMAX(ic)).eq.0 )fsc=f           fsc = clsigma(INDMAX(ic))/si
968  *     --> signal of adjacent strips  *     --> signal of adjacent strips
969           sl1 = 0                !left 1           sl1 = 0                !left 1
970           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 774  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 972  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
972       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
973       $        )then       $        )then
974              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
975  c            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
976              if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f              fsl1 = clsigma(INDMAX(ic)-1)/si
 c         else  
 c            fsl1 = 0  
977           endif           endif
978    
979           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 786  c            fsl1 = 0 Line 982  c            fsl1 = 0
982       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
983       $        )then       $        )then
984              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
985  c            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
986              if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f              fsl2 = clsigma(INDMAX(ic)-2)/si
 c         else  
 c            fsl2 = 0  
987           endif           endif
988           sr1 = 0                !right 1           sr1 = 0                !right 1
989           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 799  c            fsl2 = 0 Line 993  c            fsl2 = 0
993       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
994       $        )then       $        )then
995              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
996  c            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
997              if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f              fsr1 = clsigma(INDMAX(ic)+1)/si
 c         else  
 c            fsr1 = 0  
998           endif               endif    
999           sr2 = 0                !right 2           sr2 = 0                !right 2
1000           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 812  c            fsr1 = 0 Line 1004  c            fsr1 = 0
1004       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1005       $        )then       $        )then
1006              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
1007  c            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
1008              if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f              fsr2 = clsigma(INDMAX(ic)+2)/si
 c         else  
 c            fsr2 = 0  
1009           endif           endif
1010    
1011    
1012    
1013  ************************************************************  ************************************************************
1014  *     COG computation  *     COG2-3-4 computation
1015  ************************************************************  ************************************************************
1016    
1017  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
1018                    
1019           COG = 0.           vCOG = cog(ncog,ic)!0.
1020                    
1021           if(ncog.eq.2)then           if(ncog.eq.2)then
1022              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1023                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
1024                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1025                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1026              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
1027                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
1028                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1029                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1030              endif              endif
1031           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1032              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
1033              fbad_cog =              fbad_cog =
1034       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1035              fbad_cog =              fbad_cog =
1036       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1037           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1038              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1039                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1040                 fbad_cog =                 fbad_cog =
1041       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1042       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1043                 fbad_cog =                 fbad_cog =
1044       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1045       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1046              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1047                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1048                 fbad_cog =                 fbad_cog =
1049       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1050       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1051                 fbad_cog =                 fbad_cog =
1052       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1053       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1054              endif              endif
1055           else           else
1056              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1057              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1058              COG = 0.  c            COG = 0.
1059           endif           endif
1060                    
1061        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1062    *     =========================
1063    *     COG computation
1064    *     =========================
1065                    
1066           iv=VIEW(ic)           vCOG = cog(0,ic)
1067           istart=INDSTART(IC)  
1068           istop=TOTCLLENGTH           iv     = VIEW(ic)
1069           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1070           COG=0.           istop  = TOTCLLENGTH
1071           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1072           SDE=0.           SGN = 0.
1073           do i=istart,istop           SNU = 0.
1074              ipos=i-INDMAX(ic)           SDE = 0.
1075              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1076              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1077              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1078    c$$$            if(CLSIGNAL(i).gt.cut)then
1079    c$$$               COG = COG + ipos*CLSIGNAL(i)
1080    c$$$               SGN = SGN + CLSIGNAL(i)
1081    c$$$            else
1082    c$$$               goto 10
1083    c$$$            endif
1084    c$$$         enddo
1085    c$$$ 10      continue
1086    c$$$         do i=INDMAX(IC)+1,istop
1087    c$$$            ipos = i-INDMAX(ic)
1088    c$$$            cut  = incut*CLSIGMA(i)
1089    c$$$            if(CLSIGNAL(i).gt.cut)then
1090    c$$$               COG = COG + ipos*CLSIGNAL(i)
1091    c$$$               SGN = SGN + CLSIGNAL(i)
1092    c$$$            else
1093    c$$$               goto 20
1094    c$$$            endif
1095    c$$$         enddo
1096    c$$$ 20      continue
1097    c$$$         if(SGN.le.0)then
1098    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1099    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1100    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1101    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1102    c$$$         else
1103    c$$$            COG=COG/SGN
1104    c$$$         endif
1105    
1106             do i=INDMAX(IC),istart,-1
1107                ipos = i-INDMAX(ic)
1108                cut  = incut*CLSIGMA(i)
1109              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1110                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1111              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1112                   SDE = SDE + (ipos-vCOG)**2
1113                else
1114                   goto 10
1115                endif            
1116           enddo           enddo
1117           COG=COG/SGNL(ic)   10      continue
1118           do i=istart,istop           do i=INDMAX(IC)+1,istop
1119              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1120              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1121              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1122                 fs=1                 fs = clsigma(i)/si
1123                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1124                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1125                 SDE = SDE + (ipos-COG)**2              else
1126                   goto 20
1127              endif                          endif            
1128           enddo           enddo
1129           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1130             if(SDE.ne.0)then
1131                FBAD_COG=SNU/SDE
1132             else
1133                
1134             endif
1135    
1136        else        else
1137                                        
# Line 1031  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1262  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1262    
1263  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1264    
1265        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1266    
1267        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1268        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1103  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1334  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1334       +/       +/
1335    
1336        V(1)= abs(x)        V(1)= abs(x)
1337          if(V(1).gt.20.)V(1)=20.
1338    
1339        HQUADF = 0.        HQUADF = 0.
1340        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1117  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1349  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1349     20 CONTINUE     20 CONTINUE
1350        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1351                
1352        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1353    
1354        END        END
1355    
1356  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1357        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1358        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1359        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1360        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1193  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1425  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1425       +/       +/
1426    
1427        V(1) =  abs(x)        V(1) =  abs(x)
1428          if(V(1).gt.20.)V(1)=20.
1429    
1430        HQUADF = 0.        HQUADF = 0.
1431        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1432           HQDJ = 0.           HQDJ = 0.
# Line 1206  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1440  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1440     20 CONTINUE     20 CONTINUE
1441        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1442                
1443        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1444    
1445        END        END
1446  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1447        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1448        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1449        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1450        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1281  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1515  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1515       +/       +/
1516    
1517        V(1)=abs(x)        V(1)=abs(x)
1518          if(V(1).gt.20.)V(1)=20.
1519    
1520        HQUADF = 0.        HQUADF = 0.
1521        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1295  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1530  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1530     20 CONTINUE     20 CONTINUE
1531        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1532                
1533        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1534    
1535        END        END
1536  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1537        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1538        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1539        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1540        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1352  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1587  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1587       +/       +/
1588    
1589        v(1)= abs(x)        v(1)= abs(x)
1590          if(V(1).gt.20.)V(1)=20.
1591                
1592        HQUADF = 0.        HQUADF = 0.
1593        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1366  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1602  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1602     20 CONTINUE     20 CONTINUE
1603        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1604    
1605        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1606    
1607        END        END
1608  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1418  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1654  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1654       +/       +/
1655                
1656        V(1)=abs(x)        V(1)=abs(x)
1657          if(V(1).gt.20.)V(1)=20.
1658    
1659        HQUADF = 0.        HQUADF = 0.
1660        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1498  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1735  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1735       +/       +/
1736    
1737        V(1)=abs(x)        V(1)=abs(x)
1738          if(V(1).gt.20.)V(1)=20.
1739    
1740        HQUADF = 0.        HQUADF = 0.
1741        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1515  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1753  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1753        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1754    
1755        END        END
1756    
1757    
1758  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1759          real function pfacorr(ic,angle)
1760    *--------------------------------------------------------------
1761    *     this function returns the landi correction for this cluster
1762    *--------------------------------------------------------------
1763          include 'commontracker.f'
1764          include 'calib.f'
1765          include 'level1.f'
1766    
1767          real angle
1768          integer iview,lad
1769    
1770          iview = VIEW(ic)            
1771          lad = nld(MAXS(ic),VIEW(ic))
1772    
1773    *     find angular bin
1774    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1775          do iang=1,nangbin
1776             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1777                iangle=iang
1778                goto 98
1779             endif
1780          enddo
1781          if(DEBUG.eq.1)
1782         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1783          if(angle.lt.angL(1))iang=1
1784          if(angle.gt.angR(nangbin))iang=nangbin
1785     98   continue                  !jump here if ok
1786    
1787          pfacorr = fcorr(iview,lad,iang)
1788    
1789          if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1790    
1791    
1792     100  return
1793          end
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.23