/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/functionspfa.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/functionspfa.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by pam-fi, Tue May 30 16:30:37 2006 UTC revision 1.14 by pam-fi, Tue Aug 7 13:56:29 2007 UTC
# Line 1  Line 1 
1    
2    
3          subroutine idtoc(ipfa,cpfa)
4          
5          integer ipfa
6          character*4 cpfa
7    
8          CPFA='COG4'
9          if(ipfa.eq.0)CPFA='ETA'
10          if(ipfa.eq.2)CPFA='ETA2'
11          if(ipfa.eq.3)CPFA='ETA3'
12          if(ipfa.eq.4)CPFA='ETA4'
13          if(ipfa.eq.10)CPFA='COG'
14          if(ipfa.eq.11)CPFA='COG1'
15          if(ipfa.eq.12)CPFA='COG2'
16          if(ipfa.eq.13)CPFA='COG3'
17          if(ipfa.eq.14)CPFA='COG4'
18          
19          end
20    
21  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
22  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
23  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms
# Line 6  Line 26 
26  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
27    
28    
29          integer function npfastrips(ic,PFA,angle)
30    *--------------------------------------------------------------
31    *     thid function returns the number of strips used
32    *     to evaluate the position of a cluster, according to the p.f.a.
33    *--------------------------------------------------------------
34          include 'commontracker.f'
35          include 'level1.f'
36          include 'calib.f'
37    
38          character*4 usedPFA,PFA
39    
40    
41          usedPFA=PFA
42    
43          npfastrips=0
44    
45          if(usedPFA.eq.'COG1')npfastrips=1
46          if(usedPFA.eq.'COG2')npfastrips=2
47          if(usedPFA.eq.'COG3')npfastrips=3
48          if(usedPFA.eq.'COG4')npfastrips=4
49          if(usedPFA.eq.'ETA2')npfastrips=2
50          if(usedPFA.eq.'ETA3')npfastrips=3
51          if(usedPFA.eq.'ETA4')npfastrips=4
52    *     ----------------------------------------------------------------
53          if(usedPFA.eq.'ETA')then
54    c         print*,VIEW(ic),angle
55             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
56                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
57                   npfastrips=2
58                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
59                   npfastrips=3
60                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
61                   npfastrips=4
62                else
63                   npfastrips=4
64    c               usedPFA='COG'
65                endif                        
66             else                   !X-view
67                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
68                   npfastrips=2
69                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
70                   npfastrips=3
71                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
72                   npfastrips=4
73                else
74                   npfastrips=4
75    c               usedPFA='COG'
76                endif                        
77             endif
78          endif
79    *     ----------------------------------------------------------------
80          if(usedPFA.eq.'COG')then
81    
82             iv=VIEW(ic)
83             if(mod(iv,2).eq.1)incut=incuty
84             if(mod(iv,2).eq.0)incut=incutx
85             istart = INDSTART(IC)
86             istop  = TOTCLLENGTH
87             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
88             mu  = 0
89             do i = INDMAX(IC),istart,-1
90                ipos = i-INDMAX(ic)
91                cut  = incut*CLSIGMA(i)
92                if(CLSIGNAL(i).ge.cut)then
93                   mu = mu + 1
94                   print*,i,mu
95                else
96                   goto 10
97                endif
98             enddo
99     10      continue
100             do i = INDMAX(IC)+1,istop
101                ipos = i-INDMAX(ic)
102                cut  = incut*CLSIGMA(i)
103                if(CLSIGNAL(i).ge.cut)then
104                   mu = mu + 1
105                   print*,i,mu
106                else
107                   goto 20
108                endif
109             enddo
110     20      continue
111             npfastrips=mu
112    
113          endif
114    *     ----------------------------------------------------------------
115    
116    c      print*,pfastrips
117    
118          return
119          end
120    
121  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
122        real function pfa_eta(ic,angle)        real function pfaeta(ic,angle)
123  *--------------------------------------------------------------  *--------------------------------------------------------------
124  *     this function returns the position (in strip units)  *     this function returns the position (in strip units)
125  *     it calls:  *     it calls:
126  *     - pfa_eta2(ic,angle)  *     - pfaeta2(ic,angle)
127  *     - pfa_eta3(ic,angle)  *     - pfaeta3(ic,angle)
128  *     - pfa_eta4(ic,angle)  *     - pfaeta4(ic,angle)
129  *     according to the angle  *     according to the angle
130  *--------------------------------------------------------------  *--------------------------------------------------------------
131        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
132        include 'level1.f'        include 'level1.f'
133          include 'calib.f'
134                
135        pfa_eta = 0        pfaeta = 0
136    
137        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
138                
139           pfa_eta = pfa_eta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
140                  pfaeta = pfaeta2(ic,angle)
141             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
142                pfaeta = pfaeta3(ic,angle)
143             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
144                pfaeta = pfaeta4(ic,angle)
145             else
146                pfaeta = cog(4,ic)
147             endif            
148    
149        else                      !X-view        else                      !X-view
150    
151           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
152              pfa_eta = pfa_eta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
153           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
154              pfa_eta = pfa_eta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
155           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
156              pfa_eta = pfa_eta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
157           endif           else
158                pfaeta = cog(4,ic)
159             endif            
160                            
161        endif        endif
162          
   
163   100  return   100  return
164        end        end
165    
166  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
167        real function ris_eta(ic,angle)  c      real function riseta(ic,angle)
168          real function riseta(iview,angle)
169  *--------------------------------------------------------------  *--------------------------------------------------------------
170  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
171  *     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
172  *     it calls:  *     it calls:
173  *     - risx_eta2(angle)  *     - risxeta2(angle)
174  *     - risy_eta2(angle)  *     - risyeta2(angle)
175  *     - risx_eta3(angle)  *     - risxeta3(angle)
176  *     - risx_eta4(angle)  *     - risxeta4(angle)
177  *     according to the angle  *     according to the angle
178  *--------------------------------------------------------------  *--------------------------------------------------------------
179        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
180        include 'level1.f'        include 'level1.f'
181          include 'calib.f'
182    
183  c$$$      logical DEBUG        riseta = 0
 c$$$      common/dbg/DEBUG  
         
       ris_eta = 0  
184    
185        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
186          if(mod(iview,2).eq.1)then !Y-view
187                
188           ris_eta = risy_eta2(angle)  
189           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
190                riseta = risyeta2(angle)
191             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
192                riseta = risy_cog(angle) !ATTENZIONE!!
193             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
194                riseta = risy_cog(angle) !ATTENZIONE!!
195             else
196                riseta = risy_cog(angle)
197             endif            
198    
199        else                      !X-view        else                      !X-view
200    
201           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
202              ris_eta = risx_eta2(angle)              riseta = risxeta2(angle)
203           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
204              ris_eta = risx_eta3(angle)              riseta = risxeta3(angle)
205           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
206              ris_eta = risx_eta4(angle)              riseta = risxeta4(angle)
207           elseif(abs(angle).gt.21.)then           else
208              ris_eta = risx_eta4(21.)              riseta = risx_cog(angle)
209           endif           endif            
210                            
211        endif        endif
212    
213  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'  cc      print*,'---- ',riseta,iview,angle
 c$$$     $     ,' -->',ris_eta  
   
214    
215   100  return   100  return
216        end        end
# Line 103  c$$$     $     ,' -->',ris_eta Line 229  c$$$     $     ,' -->',ris_eta
229    
230        include 'commontracker.f'        include 'commontracker.f'
231        include 'level1.f'        include 'level1.f'
232  *      include 'calib.f'        include 'calib.f'
233        fbad_eta = 0        fbad_eta = 0
234    
235        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
236                
237           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
238                fbad_eta = fbad_cog(2,ic)
239             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
240                fbad_eta = fbad_cog(3,ic)
241             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
242                fbad_eta = fbad_cog(4,ic)
243             else
244                fbad_eta = fbad_cog(4,ic)
245             endif            
246    
247        else                      !X-view        else                      !X-view
248    
249           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
250              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
251           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
252              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
253           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
254              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
255           endif           else
256                fbad_eta = fbad_cog(4,ic)
257             endif            
258                            
259        endif        endif
260    
# Line 126  c$$$     $     ,' -->',ris_eta Line 262  c$$$     $     ,' -->',ris_eta
262        end        end
263    
264  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
265  c*****************************************************        real function pfaeta2(ic,angle) !(1)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta2(cog2,view,lad,angle)  
       real function pfa_eta2(ic,angle) !(1)  
266  *--------------------------------------------------------------  *--------------------------------------------------------------
267  *     this function returns  *     this function returns
268  *  *
# Line 149  c      real function pfa_eta2(cog2,view, Line 281  c      real function pfa_eta2(cog2,view,
281        real cog2,angle        real cog2,angle
282        integer iview,lad        integer iview,lad
283    
284  c      logical DEBUG        iview = VIEW(ic)            
285  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
286          cog2 = cog(2,ic)          
287        iview = VIEW(ic)              !(1)        pfaeta2=cog2
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfa_eta2=cog2  
288    
289  *     find angular bin  *     find angular bin
290  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
291        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
292           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
293              iangle=iang              iangle=iang
294              goto 98              goto 98
295           endif           endif
296        enddo        enddo
297        if(DEBUG)        if(DEBUG)
298       $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
299        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
300        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
301   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 177  c$$$*     find extremes of interpolation Line 305  c$$$*     find extremes of interpolation
305  c$$$      iflag=0  c$$$      iflag=0
306  c$$$*     --------------------------------  c$$$*     --------------------------------
307  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
308  c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2  c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2
309  c$$$*         goto 100  c$$$*         goto 100
310  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
311  c$$$*     non salto piu`, ma scalo di 1 o -1  c$$$*     non salto piu`, ma scalo di 1 o -1
# Line 230  c            print*,'-----',x1,x2,y1,y2 Line 358  c            print*,'-----',x1,x2,y1,y2
358        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
359        BB=y1-AA*x1        BB=y1-AA*x1
360    
361        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
362        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
363    
364  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
365  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
366  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
367  c$$$      endif  c$$$      endif
368  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
369  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
370  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
371  c$$$      endif  c$$$      endif
372    
373        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'
374       $     ,cog2-iadd,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
375    
376    
377   100  return   100  return
378        end        end
379    
380  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
381  c*****************************************************        real function pfaeta3(ic,angle) !(1)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta3(cog3,view,lad,angle)  
       real function pfa_eta3(ic,angle) !(1)  
382  *--------------------------------------------------------------  *--------------------------------------------------------------
383  *     this function returns  *     this function returns
384  *  *
# Line 273  c      real function pfa_eta3(cog3,view, Line 397  c      real function pfa_eta3(cog3,view,
397        real cog3,angle        real cog3,angle
398        integer iview,lad        integer iview,lad
399    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
400    
401        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
402        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
403        cog3 = cog(3,ic)             !(1)        cog3 = cog(3,ic)            
404        pfa_eta3=cog3        pfaeta3=cog3
405    
406  *     find angular bin  *     find angular bin
407  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 292  c         print*,'~~~~~~~~~~~~ ',iang,an Line 413  c         print*,'~~~~~~~~~~~~ ',iang,an
413           endif           endif
414        enddo        enddo
415        if(DEBUG)        if(DEBUG)
416       $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
417        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
418        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
419   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 354  c            print*,'-----',x1,x2,y1,y2 Line 475  c            print*,'-----',x1,x2,y1,y2
475        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
476        BB=y1-AA*x1        BB=y1-AA*x1
477    
478        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
479        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
480    
481  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
482  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
483  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
484  c$$$      endif  c$$$      endif
485  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
486  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
487  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
488  c$$$      endif  c$$$      endif
489    
490        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'
491       $     ,cog3-iadd,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
492    
493   100  return   100  return
494        end        end
495    
496  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
497  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta4(cog4,view,lad,angle)  
       real function pfa_eta4(ic,angle) !(1)  
498  *--------------------------------------------------------------  *--------------------------------------------------------------
499  *     this function returns  *     this function returns
500  *  *
# Line 396  c      real function pfa_eta4(cog4,view, Line 513  c      real function pfa_eta4(cog4,view,
513        real cog4,angle        real cog4,angle
514        integer iview,lad        integer iview,lad
515    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
516    
517        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
518        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
519        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
520        pfa_eta4=cog4        pfaeta4=cog4
521    
522  *     find angular bin  *     find angular bin
523  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 414  c         print*,'~~~~~~~~~~~~ ',iang,an Line 529  c         print*,'~~~~~~~~~~~~ ',iang,an
529           endif           endif
530        enddo        enddo
531        if(DEBUG)        if(DEBUG)
532       $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
533        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
534        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
535   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 476  c            print*,'-----',x1,x2,y1,y2 Line 591  c            print*,'-----',x1,x2,y1,y2
591        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
592        BB=y1-AA*x1        BB=y1-AA*x1
593    
594        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
595        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
596    
597  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
598  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
599  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
600  c$$$      endif  c$$$      endif
601  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
602  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
603  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
604  c$$$      endif  c$$$      endif
605    
606        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'
607       $     ,cog4-iadd,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
608    
609   100  return   100  return
610        end        end
# Line 497  c$$$      endif Line 612  c$$$      endif
612    
613    
614  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
615        real function cog0(ncog,ic)  c$$$      real function cog0(ncog,ic)
616  *-------------------------------------------------  c$$$*-------------------------------------------------
617  *     this function returns  c$$$*     this function returns
618  *  c$$$*
619  *     - the Center-Of-Gravity of the cluster IC  c$$$*     - the Center-Of-Gravity of the cluster IC
620  *     evaluated using NCOG strips,  c$$$*     evaluated using NCOG strips,
621  *     calculated relative to MAXS(IC)  c$$$*     calculated relative to MAXS(IC)
622  *      c$$$*    
623  *     - zero in case that not  enough strips  c$$$*     - zero in case that not  enough strips
624  *     have a positive signal  c$$$*     have a positive signal
625  *      c$$$*    
626  *     NOTE:  c$$$*     NOTE:
627  *     This is the old definition, used by Straulino.  c$$$*     This is the old definition, used by Straulino.
628  *     The new routine, according to Landi,  c$$$*     The new routine, according to Landi,
629  *     is COG(NCOG,IC)  c$$$*     is COG(NCOG,IC)
630  *-------------------------------------------------  c$$$*-------------------------------------------------
631    c$$$
632    c$$$
633        include 'commontracker.f'  c$$$      include 'commontracker.f'
634        include 'level1.f'  c$$$      include 'level1.f'
635          c$$$      
636  *     --> signal of the central strip  c$$$*     --> signal of the central strip
637        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
638    c$$$
639  *     signal of adjacent strips  c$$$*     signal of adjacent strips
640  *     --> left  c$$$*     --> left
641        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
642        if(  c$$$      if(
643       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
644       $     )  c$$$     $     )
645       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
646    c$$$
647        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
648        if(  c$$$      if(
649       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
650       $     )  c$$$     $     )
651       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
652    c$$$
653  *     --> right  c$$$*     --> right
654        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
655        if(  c$$$      if(
656       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
657       $     .or.  c$$$     $     .or.
658       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
659       $     )  c$$$     $     )
660       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
661    c$$$
662        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
663        if(  c$$$      if(
664       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
665       $     .or.  c$$$     $     .or.
666       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
667       $     )  c$$$     $     )
668       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
669          c$$$      
670  ************************************************************  c$$$************************************************************
671  *     COG computation  c$$$*     COG computation
672  ************************************************************  c$$$************************************************************
673    c$$$
674  c      print*,sl2,sl1,sc,sr1,sr2  c$$$c      print*,sl2,sl1,sc,sr1,sr2
675    c$$$
676        COG = 0.  c$$$      COG = 0.
677            c$$$        
678        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
679            c$$$        
680           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
681              COG = -sl1/(sl1+sc)          c$$$            COG = -sl1/(sl1+sc)        
682           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
683              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
684           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
685              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
686           else  c$$$         else
687              COG = 0.  c$$$            COG = 0.
688           endif  c$$$         endif
689            c$$$        
690        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
691    c$$$
692           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
693              COG = sr1/(sc+sr1)              c$$$            COG = sr1/(sc+sr1)            
694           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
695              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
696           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
697              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
698           else  c$$$         else
699              COG = 0.  c$$$            COG = 0.
700           endif  c$$$         endif
701    c$$$
702        endif  c$$$      endif
703    c$$$
704        COG0 = COG  c$$$      COG0 = COG
705    c$$$
706  c      print *,ncog,ic,cog,'/////////////'  c$$$c      print *,ncog,ic,cog,'/////////////'
707    c$$$
708        return  c$$$      return
709        end  c$$$      end
710    
711  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
712        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 613  c      print *,ncog,ic,cog,'//////////// Line 728  c      print *,ncog,ic,cog,'////////////
728        include 'calib.f'        include 'calib.f'
729        include 'level1.f'        include 'level1.f'
730                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
731    
732    
733        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 625  c      common/dbg/DEBUG Line 738  c      common/dbg/DEBUG
738  *     --> signal of the central strip  *     --> signal of the central strip
739           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
740  *     signal of adjacent strips  *     signal of adjacent strips
741           sl1 = 0                !left 1           sl1 = -9999.           !left 1
742           if(           if(
743       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
744       $        )       $        )
745       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
746                    
747           sl2 = 0                !left 2           sl2 = -9999.           !left 2
748           if(           if(
749       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
750       $        )       $        )
751       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
752                    
753           sr1 = 0                !right 1           sr1 = -9999.           !right 1
754           if(           if(
755       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
756       $        .or.       $        .or.
# Line 645  c      common/dbg/DEBUG Line 758  c      common/dbg/DEBUG
758       $        )       $        )
759       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
760                    
761           sr2 = 0                !right 2           sr2 = -9999.           !right 2
762           if(           if(
763       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
764       $        .or.       $        .or.
# Line 655  c      common/dbg/DEBUG Line 768  c      common/dbg/DEBUG
768                    
769           COG = 0.           COG = 0.
770                    
771           if(ncog.eq.2)then  c         print*,'## ',sl2,sl1,sc,sr1,sr2
772    
773    c     ==============================================================
774             if(ncog.eq.1)then
775                COG = 0.
776                if(sr1.gt.sc)cog=1. !NEW
777                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
778    c     ==============================================================
779             elseif(ncog.eq.2)then
780              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
781                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
782              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
783                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
784              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
785                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
786         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
787                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
788         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
789                endif
790    c$$$            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
791    c$$$     $           ,VIEW(ic),LADDER(ic)
792    c$$$     $           ,' : ',sl2,sl1,sc,sr1,sr2
793    c     ==============================================================
794           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
795              COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
796    c$$$             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
797    c$$$     $           ,VIEW(ic),LADDER(ic)
798    c$$$     $            ,' : ',sl2,sl1,sc,sr1,sr2
799    c     ==============================================================
800           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
801              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
802                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sl1+sc+sr1).ne.0)
803              elseif(sl2.le.sr2)then       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
804                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)              elseif(sl2.lt.sr2)then
805                   if((sr2+sl1+sc+sr1).ne.0)
806         $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
807                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
808                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
809         $              .and.(sl2+sl1+sc+sr1).ne.0 )
810         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
811                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
812         $              .and.(sr2+sl1+sc+sr1).ne.0 )
813         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
814              endif              endif
815    c$$$            if(cog==0)print*,'Strange cluster (4) - @maxs ',MAXS(ic)
816    c$$$     $           ,VIEW(ic),LADDER(ic)
817    c$$$     $           ,' : ',sl2,sl1,sc,sr1,sr2
818    c     ==============================================================
819           else           else
820              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
821              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
822              COG = 0.              COG = 0.
823           endif           endif
824    
# Line 685  ' Line 832  '
832           iv=VIEW(ic)           iv=VIEW(ic)
833           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
834           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
835             istart = INDSTART(IC)
836           istart=INDSTART(IC)           istop  = TOTCLLENGTH
          istop=TOTCLLENGTH  
837           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
838           COG=0             COG = 0  
839           mu=0           SGN = 0.
840           do i=istart,istop           mu  = 0
841              ipos=i-INDMAX(ic)  c         print*,'-------'
842              ivk=nvk(MAXS(ic)+ipos)           do i = INDMAX(IC),istart,-1
843              is=nst(MAXS(ic)+ipos)              ipos = i-INDMAX(ic)
844  *            print*,'******************',istart,istop,ipos              cut  = incut*CLSIGMA(i)
845  *     $           ,MAXS(ic)+ipos,iv,ivk,is              if(CLSIGNAL(i).ge.cut)then              
846              cut=incut*SIGMA(iv,ivk,is)                 COG = COG + ipos*CLSIGNAL(i)
847                   SGN = SGN + CLSIGNAL(i)
848                   mu = mu + 1
849    c               print*,ipos,CLSIGNAL(i)
850                else
851                   goto 10
852                endif
853             enddo
854     10      continue
855             do i = INDMAX(IC)+1,istop
856                ipos = i-INDMAX(ic)
857                cut  = incut*CLSIGMA(i)
858              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
859                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
860                   SGN = SGN + CLSIGNAL(i)
861                 mu = mu + 1                 mu = mu + 1
862  c               print*,ipos,CLSIGNAL(i),incut,cut  c               print*,ipos,CLSIGNAL(i)
863                else
864                   goto 20
865              endif              endif
866           enddo           enddo
867           COG=COG/DEDX(ic)   20      continue
868  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'           if(SGN.le.0)then
869  c     $        ,cog              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
870                print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
871                print*,(CLSIGNAL(i),i=istart,istop)
872    c            print*,'cog(0,ic) --> NOT EVALUATED '
873             else
874                COG=COG/SGN
875             endif
876    c         print*,'-------'
877                    
878        else        else
879                    
# Line 717  c     $        ,cog Line 884  c     $        ,cog
884    
885        endif        endif
886    
887  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
888    
889        return        return
890        end        end
891    
892  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
893    
894        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
895  *-------------------------------------------------------  *-------------------------------------------------------
896  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 906  c      print *,ncog,ic,cog,'////////////
906        include 'calib.f'        include 'calib.f'
907    
908        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
909           f  = 4.           si = 8.4  !average good-strip noise
910           si = 8.4           f  = 4.   !average bad-strip noise: f*si
911           incut=incuty           incut=incuty
912        else                      !X-view        else                      !X-view
913           f  = 6.           si = 3.9  !average good-strip noise
914           si = 3.9           f  = 6.   !average bad-strip noise: f*si
915           incut=incutx           incut=incutx
916        endif        endif
917                
# Line 754  c      print *,ncog,ic,cog,'//////////// Line 922  c      print *,ncog,ic,cog,'////////////
922  *     --> signal of the central strip  *     --> signal of the central strip
923           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
924           fsc = 1           fsc = 1
925           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
926             fsc = clsigma(INDMAX(ic))/si
927  *     --> signal of adjacent strips  *     --> signal of adjacent strips
928           sl1 = 0                !left 1           sl1 = 0                !left 1
929           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 762  c      print *,ncog,ic,cog,'//////////// Line 931  c      print *,ncog,ic,cog,'////////////
931       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
932       $        )then       $        )then
933              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
934              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
935  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
936           endif           endif
937    
938           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 773  c            fsl1 = 0 Line 941  c            fsl1 = 0
941       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
942       $        )then       $        )then
943              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
944              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
945  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
946           endif           endif
947           sr1 = 0                !right 1           sr1 = 0                !right 1
948           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 785  c            fsl2 = 0 Line 952  c            fsl2 = 0
952       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
953       $        )then       $        )then
954              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
955              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
956  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
957           endif               endif    
958           sr2 = 0                !right 2           sr2 = 0                !right 2
959           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 797  c            fsr1 = 0 Line 963  c            fsr1 = 0
963       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
964       $        )then       $        )then
965              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
966              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
967  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
968           endif           endif
969    
970    
971    
972  ************************************************************  ************************************************************
973  *     COG computation  *     COG2-3-4 computation
974  ************************************************************  ************************************************************
975    
976  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
977                    
978           COG = 0.           vCOG = cog(ncog,ic)!0.
979                    
980           if(ncog.eq.2)then           if(ncog.eq.2)then
981              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
982                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
983                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
984                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
985              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
986                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
987                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
988                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
989              endif              endif
990           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
991              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
992              fbad_cog =              fbad_cog =
993       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
994              fbad_cog =              fbad_cog =
995       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
996           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
997              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
998                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
999                 fbad_cog =                 fbad_cog =
1000       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1001       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1002                 fbad_cog =                 fbad_cog =
1003       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1004       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1005              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1006                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1007                 fbad_cog =                 fbad_cog =
1008       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1009       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1010                 fbad_cog =                 fbad_cog =
1011       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1012       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1013              endif              endif
1014           else           else
1015              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1016              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1017              COG = 0.  c            COG = 0.
1018           endif           endif
1019                    
1020        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1021    *     =========================
1022    *     COG computation
1023    *     =========================
1024                    
1025           iv=VIEW(ic)           vCOG = cog(0,ic)
1026           istart=INDSTART(IC)  
1027           istop=TOTCLLENGTH           iv     = VIEW(ic)
1028           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1029           COG=0.           istop  = TOTCLLENGTH
1030           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1031           SDE=0.           SGN = 0.
1032           do i=istart,istop           SNU = 0.
1033              ipos=i-INDMAX(ic)           SDE = 0.
1034              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1035              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1036              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1037    c$$$            if(CLSIGNAL(i).gt.cut)then
1038    c$$$               COG = COG + ipos*CLSIGNAL(i)
1039    c$$$               SGN = SGN + CLSIGNAL(i)
1040    c$$$            else
1041    c$$$               goto 10
1042    c$$$            endif
1043    c$$$         enddo
1044    c$$$ 10      continue
1045    c$$$         do i=INDMAX(IC)+1,istop
1046    c$$$            ipos = i-INDMAX(ic)
1047    c$$$            cut  = incut*CLSIGMA(i)
1048    c$$$            if(CLSIGNAL(i).gt.cut)then
1049    c$$$               COG = COG + ipos*CLSIGNAL(i)
1050    c$$$               SGN = SGN + CLSIGNAL(i)
1051    c$$$            else
1052    c$$$               goto 20
1053    c$$$            endif
1054    c$$$         enddo
1055    c$$$ 20      continue
1056    c$$$         if(SGN.le.0)then
1057    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1058    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1059    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1060    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1061    c$$$         else
1062    c$$$            COG=COG/SGN
1063    c$$$         endif
1064    
1065             do i=INDMAX(IC),istart,-1
1066                ipos = i-INDMAX(ic)
1067                cut  = incut*CLSIGMA(i)
1068              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1069                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1070              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1071                   SDE = SDE + (ipos-vCOG)**2
1072                else
1073                   goto 10
1074                endif            
1075           enddo           enddo
1076           COG=COG/DEDX(ic)   10      continue
1077           do i=istart,istop           do i=INDMAX(IC)+1,istop
1078              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1079              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1080              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1081                 fs=1                 fs = clsigma(i)/si
1082                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1083                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1084                 SDE = SDE + (ipos-COG)**2              else
1085                   goto 20
1086              endif                          endif            
1087           enddo           enddo
1088           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1089             if(SDE.ne.0)then
1090                FBAD_COG=SNU/SDE
1091             else
1092                
1093             endif
1094    
1095        else        else
1096                                        
# Line 1015  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1221  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1221    
1222  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1223    
1224        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1225    
1226        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1227        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1087  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1293  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1293       +/       +/
1294    
1295        V(1)= abs(x)        V(1)= abs(x)
1296          if(V(1).gt.20.)V(1)=20.
1297    
1298        HQUADF = 0.        HQUADF = 0.
1299        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1101  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1308  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1308     20 CONTINUE     20 CONTINUE
1309        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1310                
1311        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1312    
1313        END        END
1314    
1315  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1316        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1317        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1318        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1319        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1177  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1384  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1384       +/       +/
1385    
1386        V(1) =  abs(x)        V(1) =  abs(x)
1387          if(V(1).gt.20.)V(1)=20.
1388    
1389        HQUADF = 0.        HQUADF = 0.
1390        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1391           HQDJ = 0.           HQDJ = 0.
# Line 1190  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1399  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1399     20 CONTINUE     20 CONTINUE
1400        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1401                
1402        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1403    
1404        END        END
1405  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1406        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1407        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1408        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1409        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1265  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1474  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1474       +/       +/
1475    
1476        V(1)=abs(x)        V(1)=abs(x)
1477          if(V(1).gt.20.)V(1)=20.
1478    
1479        HQUADF = 0.        HQUADF = 0.
1480        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1279  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1489  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1489     20 CONTINUE     20 CONTINUE
1490        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1491                
1492        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1493    
1494        END        END
1495  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1496        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1497        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1498        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1499        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1336  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1546  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1546       +/       +/
1547    
1548        v(1)= abs(x)        v(1)= abs(x)
1549          if(V(1).gt.20.)V(1)=20.
1550                
1551        HQUADF = 0.        HQUADF = 0.
1552        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1350  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1561  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1561     20 CONTINUE     20 CONTINUE
1562        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1563    
1564        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1565    
1566        END        END
1567  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1402  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1613  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1613       +/       +/
1614                
1615        V(1)=abs(x)        V(1)=abs(x)
1616          if(V(1).gt.20.)V(1)=20.
1617    
1618        HQUADF = 0.        HQUADF = 0.
1619        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1482  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1694  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1694       +/       +/
1695    
1696        V(1)=abs(x)        V(1)=abs(x)
1697          if(V(1).gt.20.)V(1)=20.
1698    
1699        HQUADF = 0.        HQUADF = 0.
1700        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1499  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1712  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1712        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1713    
1714        END        END
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.23