/[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.12 by pam-fi, Wed May 23 15:31:19 2007 UTC
# Line 1  Line 1 
1    
2    
3          subroutine idtoc(ipfa,cpfa)
4          
5          integer ipfa
6          character*4 cpfa
7    
8          CPFA='COG4'
9          if(ipfa.eq.0)CPFA='ETA'
10          if(ipfa.eq.2)CPFA='ETA2'
11          if(ipfa.eq.3)CPFA='ETA3'
12          if(ipfa.eq.4)CPFA='ETA4'
13          if(ipfa.eq.10)CPFA='COG'
14          if(ipfa.eq.11)CPFA='COG1'
15          if(ipfa.eq.12)CPFA='COG2'
16          if(ipfa.eq.13)CPFA='COG3'
17          if(ipfa.eq.14)CPFA='COG4'
18          
19          end
20    
21  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
22  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
23  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms
# Line 6  Line 26 
26  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
27    
28    
29          integer function npfastrips(ic,PFA,angle)
30    *--------------------------------------------------------------
31    *     thid function returns the number of strips used
32    *     to evaluate the position of a cluster, according to the p.f.a.
33    *--------------------------------------------------------------
34          include 'commontracker.f'
35          include 'level1.f'
36          include 'calib.f'
37    
38          character*4 usedPFA,PFA
39    
40    
41          usedPFA=PFA
42    
43          npfastrips=0
44    
45          if(usedPFA.eq.'COG1')npfastrips=1
46          if(usedPFA.eq.'COG2')npfastrips=2
47          if(usedPFA.eq.'COG3')npfastrips=3
48          if(usedPFA.eq.'COG4')npfastrips=4
49          if(usedPFA.eq.'ETA2')npfastrips=2
50          if(usedPFA.eq.'ETA3')npfastrips=3
51          if(usedPFA.eq.'ETA4')npfastrips=4
52    *     ----------------------------------------------------------------
53          if(usedPFA.eq.'ETA')then
54    c         print*,VIEW(ic),angle
55             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
56                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
57                   npfastrips=2
58                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
59                   npfastrips=3
60                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
61                   npfastrips=4
62                else
63                   npfastrips=4
64    c               usedPFA='COG'
65                endif                        
66             else                   !X-view
67                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
68                   npfastrips=2
69                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
70                   npfastrips=3
71                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
72                   npfastrips=4
73                else
74                   npfastrips=4
75    c               usedPFA='COG'
76                endif                        
77             endif
78          endif
79    *     ----------------------------------------------------------------
80          if(usedPFA.eq.'COG')then
81    
82             iv=VIEW(ic)
83             if(mod(iv,2).eq.1)incut=incuty
84             if(mod(iv,2).eq.0)incut=incutx
85             istart = INDSTART(IC)
86             istop  = TOTCLLENGTH
87             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
88             mu  = 0
89             do i = INDMAX(IC),istart,-1
90                ipos = i-INDMAX(ic)
91                cut  = incut*CLSIGMA(i)
92                if(CLSIGNAL(i).ge.cut)then
93                   mu = mu + 1
94                   print*,i,mu
95                else
96                   goto 10
97                endif
98             enddo
99     10      continue
100             do i = INDMAX(IC)+1,istop
101                ipos = i-INDMAX(ic)
102                cut  = incut*CLSIGMA(i)
103                if(CLSIGNAL(i).ge.cut)then
104                   mu = mu + 1
105                   print*,i,mu
106                else
107                   goto 20
108                endif
109             enddo
110     20      continue
111             npfastrips=mu
112    
113          endif
114    *     ----------------------------------------------------------------
115    
116    c      print*,pfastrips
117    
118          return
119          end
120    
121  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
122        real function pfaeta(ic,angle)        real function pfaeta(ic,angle)
123  *--------------------------------------------------------------  *--------------------------------------------------------------
# Line 17  Line 129 
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        pfaeta = 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           pfaeta = pfaeta2(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              pfaeta = pfaeta2(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              pfaeta = pfaeta3(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              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
157           endif           else
158                pfaeta = cog(4,ic)
159             endif            
160                            
161        endif        endif
162                
 c      print*,'pfaeta ',pfaeta, angle  
   
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 pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
# Line 56  c      print*,'pfaeta ',pfaeta, angle Line 177  c      print*,'pfaeta ',pfaeta, 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 = risy_eta2(angle)
191             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
192                riseta = risy_cog(angle) !ATTENZIONE!!
193             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
194                riseta = risy_cog(angle) !ATTENZIONE!!
195             else
196                riseta = risy_cog(angle)
197             endif            
198    
199        else                      !X-view        else                      !X-view
200    
201           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
202              ris_eta = risx_eta2(angle)              riseta = risx_eta2(angle)
203           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
204              ris_eta = risx_eta3(angle)              riseta = risx_eta3(angle)
205           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
206              ris_eta = risx_eta4(angle)              riseta = risx_eta4(angle)
207           elseif(abs(angle).gt.21.)then           else
208              ris_eta = risx_eta4(21.)              riseta = risx_cog(angle)
209           endif           endif            
210                            
211        endif        endif
212    
213  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        print*,'---- ',riseta,iview,angle
 c$$$     $     ,' -->',ris_eta  
   
214    
215   100  return   100  return
216        end        end
# Line 104  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 127  c$$$     $     ,' -->',ris_eta Line 262  c$$$     $     ,' -->',ris_eta
262        end        end
263    
264  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
265        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
266  *--------------------------------------------------------------  *--------------------------------------------------------------
267  *     this function returns  *     this function returns
# Line 150  c      real function pfaeta2(cog2,view,l Line 281  c      real function pfaeta2(cog2,view,l
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)          
 c      print*,'## pfaeta2 ',ic,angle  
       iview = VIEW(ic)              !(1)  
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
287        pfaeta2=cog2        pfaeta2=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
# Line 252  c$$$      endif Line 378  c$$$      endif
378        end        end
379    
380  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
381        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
382  *--------------------------------------------------------------  *--------------------------------------------------------------
383  *     this function returns  *     this function returns
# Line 275  c      real function pfaeta3(cog3,view,l Line 397  c      real function pfaeta3(cog3,view,l
397        real cog3,angle        real cog3,angle
398        integer iview,lad        integer iview,lad
399    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
400    
401  c      print*,'## pfaeta3 ',ic,angle        iview = VIEW(ic)            
402          lad = nld(MAXS(ic),VIEW(ic))
403        iview = VIEW(ic)              !(1)        cog3 = cog(3,ic)            
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog3 = cog(3,ic)             !(1)  
404        pfaeta3=cog3        pfaeta3=cog3
405    
406  *     find angular bin  *     find angular bin
# Line 376  c$$$      endif Line 494  c$$$      endif
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 pfaeta4(cog4,view,lad,angle)  
       real function pfaeta4(ic,angle) !(1)  
498  *--------------------------------------------------------------  *--------------------------------------------------------------
499  *     this function returns  *     this function returns
500  *  *
# Line 399  c      real function pfaeta4(cog4,view,l Line 513  c      real function pfaeta4(cog4,view,l
513        real cog4,angle        real cog4,angle
514        integer iview,lad        integer iview,lad
515    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 c      print*,'## pfaeta4 ',ic,angle  
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        pfaeta4=cog4        pfaeta4=cog4
521    
522  *     find angular bin  *     find angular bin
# Line 618  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 677  c         print*,'## ',sl2,sl1,sc,sr1,sr Line 785  c         print*,'## ',sl2,sl1,sc,sr1,sr
785                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sl1+sc+sr1).ne.0)
786       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
787              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
788                  if((sl2+sl1+sc+sr1).ne.0)                  if((sr2+sl1+sc+sr1).ne.0)
789       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
790              endif              endif
791           else           else
792              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
# Line 696  ' Line 804  '
804           iv=VIEW(ic)           iv=VIEW(ic)
805           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
806           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
807           istart = INDSTART(IC)           istart = INDSTART(IC)
808           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
809           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
810           COG = 0             COG = 0  
811             SGN = 0.
812           mu  = 0           mu  = 0
813           do i = istart,istop  c         print*,'-------'
814             do i = INDMAX(IC),istart,-1
815                ipos = i-INDMAX(ic)
816                cut  = incut*CLSIGMA(i)
817                if(CLSIGNAL(i).ge.cut)then              
818                   COG = COG + ipos*CLSIGNAL(i)
819                   SGN = SGN + CLSIGNAL(i)
820                   mu = mu + 1
821    c               print*,ipos,CLSIGNAL(i)
822                else
823                   goto 10
824                endif
825             enddo
826     10      continue
827             do i = INDMAX(IC)+1,istop
828              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
829              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
830              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
831                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
832                   SGN = SGN + CLSIGNAL(i)
833                 mu = mu + 1                 mu = mu + 1
834    c               print*,ipos,CLSIGNAL(i)
835                else
836                   goto 20
837              endif              endif
838           enddo           enddo
839           if(SGNL(ic).le.0)then   20      continue
840              print*,'cog(0,ic) --> ic, dedx ',ic,SGNL(ic)           if(SGN.le.0)then
841                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
842              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
843              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
844              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
845           else           else
846              COG=COG/SGNL(ic)              COG=COG/SGN
847           endif           endif
848    c         print*,'-------'
849                    
850        else        else
851                    
# Line 734  c      print *,'## cog ',ncog,ic,cog,'// Line 862  c      print *,'## cog ',ncog,ic,cog,'//
862        end        end
863    
864  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
865    
866        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
867  *-------------------------------------------------------  *-------------------------------------------------------
868  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 749  c      print *,'## cog ',ncog,ic,cog,'// Line 878  c      print *,'## cog ',ncog,ic,cog,'//
878        include 'calib.f'        include 'calib.f'
879    
880        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
881           f  = 4.           si = 8.4  !average good-strip noise
882           si = 8.4           f  = 4.   !average bad-strip noise: f*si
883           incut=incuty           incut=incuty
884        else                      !X-view        else                      !X-view
885           f  = 6.           si = 3.9  !average good-strip noise
886           si = 3.9           f  = 6.   !average bad-strip noise: f*si
887           incut=incutx           incut=incutx
888        endif        endif
889                
# Line 765  c      print *,'## cog ',ncog,ic,cog,'// Line 894  c      print *,'## cog ',ncog,ic,cog,'//
894  *     --> signal of the central strip  *     --> signal of the central strip
895           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
896           fsc = 1           fsc = 1
897  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
898           if( CLBAD(INDMAX(ic)).eq.0 )fsc=f           fsc = clsigma(INDMAX(ic))/si
899  *     --> signal of adjacent strips  *     --> signal of adjacent strips
900           sl1 = 0                !left 1           sl1 = 0                !left 1
901           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 774  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 903  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
903       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
904       $        )then       $        )then
905              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
906  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
907              if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f              fsl1 = clsigma(INDMAX(ic)-1)/si
 c         else  
 c            fsl1 = 0  
908           endif           endif
909    
910           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 786  c            fsl1 = 0 Line 913  c            fsl1 = 0
913       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
914       $        )then       $        )then
915              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
916  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
917              if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f              fsl2 = clsigma(INDMAX(ic)-2)/si
 c         else  
 c            fsl2 = 0  
918           endif           endif
919           sr1 = 0                !right 1           sr1 = 0                !right 1
920           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 799  c            fsl2 = 0 Line 924  c            fsl2 = 0
924       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
925       $        )then       $        )then
926              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
927  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
928              if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f              fsr1 = clsigma(INDMAX(ic)+1)/si
 c         else  
 c            fsr1 = 0  
929           endif               endif    
930           sr2 = 0                !right 2           sr2 = 0                !right 2
931           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 812  c            fsr1 = 0 Line 935  c            fsr1 = 0
935       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
936       $        )then       $        )then
937              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
938  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
939              if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f              fsr2 = clsigma(INDMAX(ic)+2)/si
 c         else  
 c            fsr2 = 0  
940           endif           endif
941    
942    
943    
944  ************************************************************  ************************************************************
945  *     COG computation  *     COG2-3-4 computation
946  ************************************************************  ************************************************************
947    
948  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
949                    
950           COG = 0.           vCOG = cog(ncog,ic)!0.
951                    
952           if(ncog.eq.2)then           if(ncog.eq.2)then
953              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
954                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
955                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
956                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
957              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
958                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
959                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
960                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
961              endif              endif
962           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
963              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
964              fbad_cog =              fbad_cog =
965       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
966              fbad_cog =              fbad_cog =
967       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
968           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
969              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
970                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
971                 fbad_cog =                 fbad_cog =
972       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
973       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
974                 fbad_cog =                 fbad_cog =
975       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
976       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
977              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
978                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
979                 fbad_cog =                 fbad_cog =
980       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
981       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
982                 fbad_cog =                 fbad_cog =
983       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
984       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
985              endif              endif
986           else           else
987              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
988              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
989              COG = 0.  c            COG = 0.
990           endif           endif
991                    
992        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
993    *     =========================
994    *     COG computation
995    *     =========================
996                    
997           iv=VIEW(ic)           vCOG = cog(0,ic)
998           istart=INDSTART(IC)  
999           istop=TOTCLLENGTH           iv     = VIEW(ic)
1000           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1001           COG=0.           istop  = TOTCLLENGTH
1002           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1003           SDE=0.           SGN = 0.
1004           do i=istart,istop           SNU = 0.
1005              ipos=i-INDMAX(ic)           SDE = 0.
1006              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1007              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1008              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1009    c$$$            if(CLSIGNAL(i).gt.cut)then
1010    c$$$               COG = COG + ipos*CLSIGNAL(i)
1011    c$$$               SGN = SGN + CLSIGNAL(i)
1012    c$$$            else
1013    c$$$               goto 10
1014    c$$$            endif
1015    c$$$         enddo
1016    c$$$ 10      continue
1017    c$$$         do i=INDMAX(IC)+1,istop
1018    c$$$            ipos = i-INDMAX(ic)
1019    c$$$            cut  = incut*CLSIGMA(i)
1020    c$$$            if(CLSIGNAL(i).gt.cut)then
1021    c$$$               COG = COG + ipos*CLSIGNAL(i)
1022    c$$$               SGN = SGN + CLSIGNAL(i)
1023    c$$$            else
1024    c$$$               goto 20
1025    c$$$            endif
1026    c$$$         enddo
1027    c$$$ 20      continue
1028    c$$$         if(SGN.le.0)then
1029    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1030    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1031    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1032    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1033    c$$$         else
1034    c$$$            COG=COG/SGN
1035    c$$$         endif
1036    
1037             do i=INDMAX(IC),istart,-1
1038                ipos = i-INDMAX(ic)
1039                cut  = incut*CLSIGMA(i)
1040              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1041                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1042              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1043                   SDE = SDE + (ipos-vCOG)**2
1044                else
1045                   goto 10
1046                endif            
1047           enddo           enddo
1048           COG=COG/SGNL(ic)   10      continue
1049           do i=istart,istop           do i=INDMAX(IC)+1,istop
1050              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1051              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1052              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1053                 fs=1                 fs = clsigma(i)/si
1054                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1055                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1056                 SDE = SDE + (ipos-COG)**2              else
1057                   goto 20
1058              endif                          endif            
1059           enddo           enddo
1060           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1061             if(SDE.ne.0)then
1062                FBAD_COG=SNU/SDE
1063             else
1064                
1065             endif
1066    
1067        else        else
1068                                        
# Line 1103  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1265  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1265       +/       +/
1266    
1267        V(1)= abs(x)        V(1)= abs(x)
1268          if(V(1).gt.20.)V(1)=20.
1269    
1270        HQUADF = 0.        HQUADF = 0.
1271        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1193  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1356  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1356       +/       +/
1357    
1358        V(1) =  abs(x)        V(1) =  abs(x)
1359          if(V(1).gt.20.)V(1)=20.
1360    
1361        HQUADF = 0.        HQUADF = 0.
1362        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1363           HQDJ = 0.           HQDJ = 0.
# Line 1281  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1446  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1446       +/       +/
1447    
1448        V(1)=abs(x)        V(1)=abs(x)
1449          if(V(1).gt.20.)V(1)=20.
1450    
1451        HQUADF = 0.        HQUADF = 0.
1452        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1352  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1518  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1518       +/       +/
1519    
1520        v(1)= abs(x)        v(1)= abs(x)
1521          if(V(1).gt.20.)V(1)=20.
1522                
1523        HQUADF = 0.        HQUADF = 0.
1524        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1418  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1585  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1585       +/       +/
1586                
1587        V(1)=abs(x)        V(1)=abs(x)
1588          if(V(1).gt.20.)V(1)=20.
1589    
1590        HQUADF = 0.        HQUADF = 0.
1591        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1498  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1666  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1666       +/       +/
1667    
1668        V(1)=abs(x)        V(1)=abs(x)
1669          if(V(1).gt.20.)V(1)=20.
1670    
1671        HQUADF = 0.        HQUADF = 0.
1672        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1515  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1684  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1684        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1685    
1686        END        END
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23