/[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.6 by pam-fi, Fri Dec 1 10:43:18 2006 UTC revision 1.9 by pam-fi, Fri Apr 27 10:39:58 2007 UTC
# Line 1  Line 1 
1    
2  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
3  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
4  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms
# Line 6  Line 7 
7  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
8    
9    
10          integer function npfastrips(ic,PFA,angle)
11    *--------------------------------------------------------------
12    *     thid function returns the number of strips used
13    *     to evaluate the position of a cluster, according to the p.f.a.
14    *--------------------------------------------------------------
15          include 'commontracker.f'
16          include 'level1.f'
17          include 'calib.f'
18    
19          character*4 usedPFA,PFA
20    
21    
22          usedPFA=PFA
23    
24          npfastrips=0
25    
26          if(usedPFA.eq.'COG1')npfastrips=1
27          if(usedPFA.eq.'COG2')npfastrips=2
28          if(usedPFA.eq.'COG3')npfastrips=3
29          if(usedPFA.eq.'COG4')npfastrips=4
30          if(usedPFA.eq.'ETA2')npfastrips=2
31          if(usedPFA.eq.'ETA3')npfastrips=3
32          if(usedPFA.eq.'ETA4')npfastrips=4
33    *     ----------------------------------------------------------------
34          if(usedPFA.eq.'ETA')then
35    c         print*,VIEW(ic),angle
36             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
37                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
38                   npfastrips=2
39                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
40                   npfastrips=3
41                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
42                   npfastrips=4
43                else
44                   npfastrips=4
45    c               usedPFA='COG'
46                endif                        
47             else                   !X-view
48                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
49                   npfastrips=2
50                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
51                   npfastrips=3
52                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
53                   npfastrips=4
54                else
55                   npfastrips=4
56    c               usedPFA='COG'
57                endif                        
58             endif
59          endif
60    *     ----------------------------------------------------------------
61          if(usedPFA.eq.'COG')then
62    
63             iv=VIEW(ic)
64             if(mod(iv,2).eq.1)incut=incuty
65             if(mod(iv,2).eq.0)incut=incutx
66             istart = INDSTART(IC)
67             istop  = TOTCLLENGTH
68             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
69             mu  = 0
70             do i = INDMAX(IC),istart,-1
71                ipos = i-INDMAX(ic)
72                cut  = incut*CLSIGMA(i)
73                if(CLSIGNAL(i).ge.cut)then
74                   mu = mu + 1
75                   print*,i,mu
76                else
77                   goto 10
78                endif
79             enddo
80     10      continue
81             do i = INDMAX(IC)+1,istop
82                ipos = i-INDMAX(ic)
83                cut  = incut*CLSIGMA(i)
84                if(CLSIGNAL(i).ge.cut)then
85                   mu = mu + 1
86                   print*,i,mu
87                else
88                   goto 20
89                endif
90             enddo
91     20      continue
92             npfastrips=mu
93    
94          endif
95    *     ----------------------------------------------------------------
96    
97    c      print*,pfastrips
98    
99          return
100          end
101    
102  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
103        real function pfaeta(ic,angle)        real function pfaeta(ic,angle)
104  *--------------------------------------------------------------  *--------------------------------------------------------------
# Line 17  Line 110 
110  *     according to the angle  *     according to the angle
111  *--------------------------------------------------------------  *--------------------------------------------------------------
112        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
113        include 'level1.f'        include 'level1.f'
114          include 'calib.f'
115                
116        pfaeta = 0        pfaeta = 0
117    
118        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
119                
120           pfaeta = pfaeta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
121                  pfaeta = pfaeta2(ic,angle)
122             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
123                pfaeta = pfaeta3(ic,angle)
124             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
125                pfaeta = pfaeta4(ic,angle)
126             else
127                pfaeta = cog(4,ic)
128             endif            
129    
130        else                      !X-view        else                      !X-view
131    
132           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
133              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
134           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
135              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
136           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
137              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
138           endif           else
139                pfaeta = cog(4,ic)
140             endif            
141                            
142        endif        endif
143                
 c      print*,'pfaeta ',pfaeta, angle  
   
144   100  return   100  return
145        end        end
146    
# Line 56  c      print*,'pfaeta ',pfaeta, angle Line 157  c      print*,'pfaeta ',pfaeta, angle
157  *     according to the angle  *     according to the angle
158  *--------------------------------------------------------------  *--------------------------------------------------------------
159        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
160        include 'level1.f'        include 'level1.f'
161          include 'calib.f'
162    
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
         
163        ris_eta = 0        ris_eta = 0
164    
165        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
166                
167           ris_eta = risy_eta2(angle)  
168           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
169                ris_eta = risy_eta2(angle)
170             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
171                ris_eta = risy_cog(angle) !ATTENZIONE!!
172             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
173                ris_eta = risy_cog(angle) !ATTENZIONE!!
174             else
175                ris_eta = risy_cog(angle)
176             endif            
177    
178        else                      !X-view        else                      !X-view
179    
180           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
181              ris_eta = risx_eta2(angle)              ris_eta = risx_eta2(angle)
182           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
183              ris_eta = risx_eta3(angle)              ris_eta = risx_eta3(angle)
184           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
185              ris_eta = risx_eta4(angle)              ris_eta = risx_eta4(angle)
186           elseif(abs(angle).gt.21.)then           else
187              ris_eta = risx_eta4(21.)              ris_eta = risx_cog(angle)
188           endif           endif            
189                            
190        endif        endif
191    
 c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'  
 c$$$     $     ,' -->',ris_eta  
   
   
192   100  return   100  return
193        end        end
194    
# Line 104  c$$$     $     ,' -->',ris_eta Line 206  c$$$     $     ,' -->',ris_eta
206    
207        include 'commontracker.f'        include 'commontracker.f'
208        include 'level1.f'        include 'level1.f'
209  *      include 'calib.f'        include 'calib.f'
210        fbad_eta = 0        fbad_eta = 0
211    
212        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
213                
214           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
215                fbad_eta = fbad_cog(2,ic)
216             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
217                fbad_eta = fbad_cog(3,ic)
218             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
219                fbad_eta = fbad_cog(4,ic)
220             else
221                fbad_eta = fbad_cog(4,ic)
222             endif            
223    
224        else                      !X-view        else                      !X-view
225    
226           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
227              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
228           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
229              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
230           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
231              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
232           endif           else
233                fbad_eta = fbad_cog(4,ic)
234             endif            
235                            
236        endif        endif
237    
# Line 127  c$$$     $     ,' -->',ris_eta Line 239  c$$$     $     ,' -->',ris_eta
239        end        end
240    
241  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
242        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
243  *--------------------------------------------------------------  *--------------------------------------------------------------
244  *     this function returns  *     this function returns
# Line 150  c      real function pfaeta2(cog2,view,l Line 258  c      real function pfaeta2(cog2,view,l
258        real cog2,angle        real cog2,angle
259        integer iview,lad        integer iview,lad
260    
261  c      logical DEBUG        iview = VIEW(ic)            
262  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
263          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)  
264        pfaeta2=cog2        pfaeta2=cog2
265    
266  *     find angular bin  *     find angular bin
267  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
268        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
269           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
270              iangle=iang              iangle=iang
271              goto 98              goto 98
# Line 252  c$$$      endif Line 355  c$$$      endif
355        end        end
356    
357  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
358        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
359  *--------------------------------------------------------------  *--------------------------------------------------------------
360  *     this function returns  *     this function returns
# Line 275  c      real function pfaeta3(cog3,view,l Line 374  c      real function pfaeta3(cog3,view,l
374        real cog3,angle        real cog3,angle
375        integer iview,lad        integer iview,lad
376    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
377    
378  c      print*,'## pfaeta3 ',ic,angle        iview = VIEW(ic)            
379          lad = nld(MAXS(ic),VIEW(ic))
380        iview = VIEW(ic)              !(1)        cog3 = cog(3,ic)            
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog3 = cog(3,ic)             !(1)  
381        pfaeta3=cog3        pfaeta3=cog3
382    
383  *     find angular bin  *     find angular bin
# Line 376  c$$$      endif Line 471  c$$$      endif
471        end        end
472    
473  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
474  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)  
475  *--------------------------------------------------------------  *--------------------------------------------------------------
476  *     this function returns  *     this function returns
477  *  *
# Line 399  c      real function pfaeta4(cog4,view,l Line 490  c      real function pfaeta4(cog4,view,l
490        real cog4,angle        real cog4,angle
491        integer iview,lad        integer iview,lad
492    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 c      print*,'## pfaeta4 ',ic,angle  
493    
494        iview = VIEW(ic)             !(1)        iview = VIEW(ic)            
495        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
496        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
497        pfaeta4=cog4        pfaeta4=cog4
498    
499  *     find angular bin  *     find angular bin
# Line 618  c      print *,ncog,ic,cog,'//////////// Line 705  c      print *,ncog,ic,cog,'////////////
705        include 'calib.f'        include 'calib.f'
706        include 'level1.f'        include 'level1.f'
707                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
708    
709    
710        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 682  c         print*,'## ',sl2,sl1,sc,sr1,sr Line 767  c         print*,'## ',sl2,sl1,sc,sr1,sr
767              endif              endif
768           else           else
769              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
770              print*,'                 (NCOG must be 0-4)'       $           ,' not implemented'
771              COG = 0.              COG = 0.
772           endif           endif
773    
# Line 696  ' Line 781  '
781           iv=VIEW(ic)           iv=VIEW(ic)
782           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
783           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
784           istart = INDSTART(IC)           istart = INDSTART(IC)
785           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
786           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
787           COG = 0             COG = 0  
788             SGN = 0.
789           mu  = 0           mu  = 0
790           do i = istart,istop  c         print*,'-------'
791             do i = INDMAX(IC),istart,-1
792                ipos = i-INDMAX(ic)
793                cut  = incut*CLSIGMA(i)
794                if(CLSIGNAL(i).ge.cut)then              
795                   COG = COG + ipos*CLSIGNAL(i)
796                   SGN = SGN + CLSIGNAL(i)
797                   mu = mu + 1
798                   print*,ipos,CLSIGNAL(i)
799                else
800                   goto 10
801                endif
802             enddo
803     10      continue
804             do i = INDMAX(IC)+1,istop
805              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
 cc            ivk  = nvk(MAXS(ic)+ipos)  
 cc            is   = nst(MAXS(ic)+ipos)  
 *            print*,'******************',istart,istop,ipos  
 *     $           ,MAXS(ic)+ipos,iv,ivk,is  
 cc            cut  = incut*SIGMA(iv,ivk,is)  
806              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
 c            if(SIGMA(iv,ivk,is).ne.CLSIGMA(i))  
 c     $           print*,'cog(0,ic) --> hai fatto qualche cazzata'  
807              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
808                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
809                   SGN = SGN + CLSIGNAL(i)
810                 mu = mu + 1                 mu = mu + 1
811  c               print*,ipos,CLSIGNAL(i),incut,cut                 print*,ipos,CLSIGNAL(i)
812                else
813                   goto 20
814              endif              endif
815           enddo           enddo
816           if(DEDX(ic).le.0)then   20      continue
817              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)           if(SGN.le.0)then
818    c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN
819              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
820              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
821              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
822           else           else
823              COG=COG/DEDX(ic)              COG=COG/SGN
824           endif           endif
825  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'  c         print*,'-------'
 c     $        ,cog  
826                    
827        else        else
828                    
# Line 744  c      print *,'## cog ',ncog,ic,cog,'// Line 839  c      print *,'## cog ',ncog,ic,cog,'//
839        end        end
840    
841  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
842    
843        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
844  *-------------------------------------------------------  *-------------------------------------------------------
845  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 759  c      print *,'## cog ',ncog,ic,cog,'// Line 855  c      print *,'## cog ',ncog,ic,cog,'//
855        include 'calib.f'        include 'calib.f'
856    
857        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
858           f  = 4.           si = 8.4  !average good-strip noise
859           si = 8.4           f  = 4.   !average bad-strip noise: f*si
860           incut=incuty           incut=incuty
861        else                      !X-view        else                      !X-view
862           f  = 6.           si = 3.9  !average good-strip noise
863           si = 3.9           f  = 6.   !average bad-strip noise: f*si
864           incut=incutx           incut=incutx
865        endif        endif
866                
# Line 775  c      print *,'## cog ',ncog,ic,cog,'// Line 871  c      print *,'## cog ',ncog,ic,cog,'//
871  *     --> signal of the central strip  *     --> signal of the central strip
872           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
873           fsc = 1           fsc = 1
874           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
875             fsc = clsigma(INDMAX(ic))/si
876  *     --> signal of adjacent strips  *     --> signal of adjacent strips
877           sl1 = 0                !left 1           sl1 = 0                !left 1
878           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 783  c      print *,'## cog ',ncog,ic,cog,'// Line 880  c      print *,'## cog ',ncog,ic,cog,'//
880       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
881       $        )then       $        )then
882              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
883              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
884  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
885           endif           endif
886    
887           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 794  c            fsl1 = 0 Line 890  c            fsl1 = 0
890       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
891       $        )then       $        )then
892              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
893              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
894  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
895           endif           endif
896           sr1 = 0                !right 1           sr1 = 0                !right 1
897           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 806  c            fsl2 = 0 Line 901  c            fsl2 = 0
901       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
902       $        )then       $        )then
903              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
904              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
905  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
906           endif               endif    
907           sr2 = 0                !right 2           sr2 = 0                !right 2
908           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 818  c            fsr1 = 0 Line 912  c            fsr1 = 0
912       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
913       $        )then       $        )then
914              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
915              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
916  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
917           endif           endif
918    
919    
920    
921  ************************************************************  ************************************************************
922  *     COG computation  *     COG2-3-4 computation
923  ************************************************************  ************************************************************
924    
925  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
926                    
927           COG = 0.           vCOG = cog(ncog,ic)!0.
928                    
929           if(ncog.eq.2)then           if(ncog.eq.2)then
930              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
931                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
932                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
933                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
934              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
935                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
936                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
937                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
938              endif              endif
939           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
940              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
941              fbad_cog =              fbad_cog =
942       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
943              fbad_cog =              fbad_cog =
944       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
945           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
946              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
947                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
948                 fbad_cog =                 fbad_cog =
949       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
950       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
951                 fbad_cog =                 fbad_cog =
952       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
953       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
954              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
955                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
956                 fbad_cog =                 fbad_cog =
957       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
958       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
959                 fbad_cog =                 fbad_cog =
960       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
961       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
962              endif              endif
963           else           else
964              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
965              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
966              COG = 0.  c            COG = 0.
967           endif           endif
968                    
969        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
970    *     =========================
971    *     COG computation
972    *     =========================
973                    
974           iv=VIEW(ic)           vCOG = cog(0,ic)
975           istart=INDSTART(IC)  
976           istop=TOTCLLENGTH           iv     = VIEW(ic)
977           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
978           COG=0.           istop  = TOTCLLENGTH
979           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
980           SDE=0.           SGN = 0.
981           do i=istart,istop           SNU = 0.
982              ipos=i-INDMAX(ic)           SDE = 0.
983              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
984              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
985              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
986    c$$$            if(CLSIGNAL(i).gt.cut)then
987    c$$$               COG = COG + ipos*CLSIGNAL(i)
988    c$$$               SGN = SGN + CLSIGNAL(i)
989    c$$$            else
990    c$$$               goto 10
991    c$$$            endif
992    c$$$         enddo
993    c$$$ 10      continue
994    c$$$         do i=INDMAX(IC)+1,istop
995    c$$$            ipos = i-INDMAX(ic)
996    c$$$            cut  = incut*CLSIGMA(i)
997    c$$$            if(CLSIGNAL(i).gt.cut)then
998    c$$$               COG = COG + ipos*CLSIGNAL(i)
999    c$$$               SGN = SGN + CLSIGNAL(i)
1000    c$$$            else
1001    c$$$               goto 20
1002    c$$$            endif
1003    c$$$         enddo
1004    c$$$ 20      continue
1005    c$$$         if(SGN.le.0)then
1006    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1007    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1008    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1009    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1010    c$$$         else
1011    c$$$            COG=COG/SGN
1012    c$$$         endif
1013    
1014             do i=INDMAX(IC),istart,-1
1015                ipos = i-INDMAX(ic)
1016                cut  = incut*CLSIGMA(i)
1017              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1018                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1019              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1020                   SDE = SDE + (ipos-vCOG)**2
1021                else
1022                   goto 10
1023                endif            
1024           enddo           enddo
1025           COG=COG/DEDX(ic)   10      continue
1026           do i=istart,istop           do i=INDMAX(IC)+1,istop
1027              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1028              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1029              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1030                 fs=1                 fs = clsigma(i)/si
1031                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1032                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1033                 SDE = SDE + (ipos-COG)**2              else
1034                   goto 20
1035              endif                          endif            
1036           enddo           enddo
1037           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1038             if(SDE.ne.0)then
1039                FBAD_COG=SNU/SDE
1040             else
1041                
1042             endif
1043    
1044        else        else
1045                                        
# Line 1108  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1242  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1242       +/       +/
1243    
1244        V(1)= abs(x)        V(1)= abs(x)
1245          if(V(1).gt.20.)V(1)=20.
1246    
1247        HQUADF = 0.        HQUADF = 0.
1248        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1198  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1333  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1333       +/       +/
1334    
1335        V(1) =  abs(x)        V(1) =  abs(x)
1336          if(V(1).gt.20.)V(1)=20.
1337    
1338        HQUADF = 0.        HQUADF = 0.
1339        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1340           HQDJ = 0.           HQDJ = 0.
# Line 1286  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1423  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1423       +/       +/
1424    
1425        V(1)=abs(x)        V(1)=abs(x)
1426          if(V(1).gt.20.)V(1)=20.
1427    
1428        HQUADF = 0.        HQUADF = 0.
1429        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1357  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1495  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1495       +/       +/
1496    
1497        v(1)= abs(x)        v(1)= abs(x)
1498          if(V(1).gt.20.)V(1)=20.
1499                
1500        HQUADF = 0.        HQUADF = 0.
1501        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1423  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1562  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1562       +/       +/
1563                
1564        V(1)=abs(x)        V(1)=abs(x)
1565          if(V(1).gt.20.)V(1)=20.
1566    
1567        HQUADF = 0.        HQUADF = 0.
1568        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1503  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1643  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1643       +/       +/
1644    
1645        V(1)=abs(x)        V(1)=abs(x)
1646          if(V(1).gt.20.)V(1)=20.
1647    
1648        HQUADF = 0.        HQUADF = 0.
1649        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1520  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1661  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1661        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1662    
1663        END        END
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.23