/[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.20 by pam-fi, Mon Aug 27 12:57:15 2007 UTC revision 1.23 by pam-fi, Wed Mar 5 17:00:20 2008 UTC
# Line 4  Line 4 
4  *            *          
5  *     subroutine idtoc(ipfa,cpfa)  *     subroutine idtoc(ipfa,cpfa)
6  *  *
7    *     subroutine applypfa(PFAtt,ic,ang,corr,res)
8    *
9  *     integer function npfastrips(ic,angle)  *     integer function npfastrips(ic,angle)
10  *  *
11  *     real function pfaeta(ic,angle)  *     real function pfaeta(ic,angle)
# Line 26  Line 28 
28  *  *
29  *     real function pfacorr(ic,angle)  *     real function pfacorr(ic,angle)
30  *  *
31    *     real function effectiveangle(ang,iview,bbb)
32    *     real function fieldcorr(iview,bbb)
33    *
34  *     NB - The angle is the "effective angle", which is relative  *     NB - The angle is the "effective angle", which is relative
35  *          to the sensor and it takes into account the magnetic field  *          to the sensor and it takes into account the magnetic field
36  *  *
# Line 49  Line 54 
54        if(ipfa.eq.14)CPFA='COG4'        if(ipfa.eq.14)CPFA='COG4'
55                
56        end        end
57    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
58          real function effectiveangle(ang,iview,bbb)
59    
60          include 'commontracker.f'
61    
62          effectiveangle = 0.
63    
64          if(mod(iview,2).eq.0)then
65    c     =================================================
66    c     X view
67    c     =================================================
68    c     here bbb is the y component of the m.field
69             angx = ang
70             by   = bbb
71             if(iview.eq.12) angx = -1. * ang
72             if(iview.eq.12) by   = -1. * bbb
73    cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
74             tgtemp  = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001
75    
76          elseif(mod(iview,2).eq.1)then
77    c     =================================================
78    c     Y view
79    c     =================================================        
80    c     here bbb is the x component of the m.filed
81             angy = ang
82             bx   = bbb
83             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
84    
85          endif      
86          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
87    
88          return
89          end
90    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
91          real function fieldcorr(iview,bbb)
92    
93          include 'commontracker.f'
94    
95          fieldcorr = 0.
96    
97          if(mod(iview,2).eq.0)then
98    
99    c     =================================================
100    c     X view
101    c     =================================================
102    c     here bbb is the y component of the m.field
103             by   = bbb
104             if(iview.eq.12) by = -1. * bbb
105             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
106    
107          elseif(mod(iview,2).eq.1)then
108    c     =================================================
109    c     Y view
110    c     =================================================        
111    c     here bbb is the x component of the m.filed
112             bx   = bbb
113             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
114    
115          endif      
116          
117          return
118          end
119    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
120    
121          subroutine applypfa(PFAtt,ic,ang,corr,res)
122    *---------------------------------------------------------------
123    *     this subroutine calculate the coordinate of cluster ic (in
124    *     strip units), relative to the strip with the maximum signal,
125    *     and its spatial resolution (in cm), applying PFAtt.
126    *     ang is the effective angle, relative to the sensor
127    *---------------------------------------------------------------
128    
129          character*4 PFAtt
130          include 'commontracker.f'
131          include 'level1.f'
132    
133          corr = 0
134          res  = 0
135    
136          if(ic.le.0)return
137    
138          iview   = VIEW(ic)
139    
140          if(mod(iview,2).eq.0)then
141    c     =================================================
142    c     X view
143    c     =================================================
144    
145             res = RESXAV
146    
147             if(PFAtt.eq.'COG1')then
148    
149                corr   = 0
150                res = 1e-4*pitchX/sqrt(12.)!!res
151    
152             elseif(PFAtt.eq.'COG2')then
153    
154                corr    = cog(2,ic)            
155                res = risx_cog(abs(ang))!TEMPORANEO              
156                res = res*fbad_cog(2,ic)
157    
158             elseif(PFAtt.eq.'COG3')then
159    
160                corr    = cog(3,ic)            
161                res = risx_cog(abs(ang))!TEMPORANEO                      
162                res = res*fbad_cog(3,ic)
163    
164             elseif(PFAtt.eq.'COG4')then
165    
166                corr    = cog(4,ic)            
167                res = risx_cog(abs(ang))!TEMPORANEO                      
168                res = res*fbad_cog(4,ic)
169    
170             elseif(PFAtt.eq.'ETA2')then
171    
172                corr  = pfaeta2(ic,ang)          
173                res = risxeta2(abs(ang))
174                res = res*fbad_cog(2,ic)
175    
176             elseif(PFAtt.eq.'ETA3')then                        
177    
178                corr  = pfaeta3(ic,ang)          
179                res = risxeta3(abs(ang))                      
180                res = res*fbad_cog(3,ic)              
181    
182             elseif(PFAtt.eq.'ETA4')then                        
183    
184                corr  = pfaeta4(ic,ang)            
185                res = risxeta4(abs(ang))                      
186                res = res*fbad_cog(4,ic)              
187    
188             elseif(PFAtt.eq.'ETA')then  
189    
190                corr  = pfaeta(ic,ang)            
191    c            res = riseta(ic,ang)                    
192                res = riseta(iview,ang)                    
193                res = res*fbad_eta(ic,ang)            
194    
195             elseif(PFAtt.eq.'ETAL')then  
196    
197                corr  = pfaetal(ic,ang)            
198                res = riseta(iview,ang)                    
199                res = res*fbad_eta(ic,ang)            
200    
201             elseif(PFAtt.eq.'COG')then          
202    
203                corr  = cog(0,ic)            
204                res = risx_cog(abs(ang))                    
205                res = res*fbad_cog(0,ic)
206    
207             else
208                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
209             endif
210    
211    
212    *     ======================================
213    *     temporary patch for saturated clusters
214    *     ======================================
215             if( nsatstrips(ic).gt.0 )then
216                corr  = cog(4,ic)            
217                res = pitchX*1e-4/sqrt(12.)
218    cc            cc=cog(4,ic)
219    c$$$            print*,ic,' *** ',cc
220    c$$$            print*,ic,' *** ',res
221             endif
222    
223    
224          elseif(mod(iview,2).eq.1)then
225    c     =================================================
226    c     Y view
227    c     =================================================
228    
229             res = RESYAV
230    
231             if(PFAtt.eq.'COG1')then
232    
233                corr  = 0  
234                res = 1e-4*pitchY/sqrt(12.)!res  
235    
236             elseif(PFAtt.eq.'COG2')then
237    
238                corr    = cog(2,ic)
239                res = risy_cog(abs(ang))!TEMPORANEO
240                res = res*fbad_cog(2,ic)
241    
242             elseif(PFAtt.eq.'COG3')then
243    
244                corr    = cog(3,ic)
245                res = risy_cog(abs(ang))!TEMPORANEO
246                res = res*fbad_cog(3,ic)
247    
248             elseif(PFAtt.eq.'COG4')then
249    
250                corr    = cog(4,ic)
251                res = risy_cog(abs(ang))!TEMPORANEO
252                res = res*fbad_cog(4,ic)
253    
254             elseif(PFAtt.eq.'ETA2')then
255    
256                corr    = pfaeta2(ic,ang)          
257                res = risyeta2(abs(ang))              
258                res = res*fbad_cog(2,ic)
259    
260             elseif(PFAtt.eq.'ETA3')then                      
261    
262                corr    = pfaeta3(ic,ang)
263                res = res*fbad_cog(3,ic)  
264    
265             elseif(PFAtt.eq.'ETA4')then  
266    
267                corr    = pfaeta4(ic,ang)
268                res = res*fbad_cog(4,ic)
269    
270             elseif(PFAtt.eq.'ETA')then
271    
272                corr    = pfaeta(ic,ang)
273    c            res = riseta(ic,ang)  
274                res = riseta(iview,ang)  
275                res = res*fbad_eta(ic,ang)
276    
277             elseif(PFAtt.eq.'ETAL')then
278    
279                corr    = pfaetal(ic,ang)
280                res = riseta(iview,ang)  
281                res = res*fbad_eta(ic,ang)
282    
283             elseif(PFAtt.eq.'COG')then
284    
285                corr    = cog(0,ic)            
286                res = risy_cog(abs(ang))
287                res = res*fbad_cog(0,ic)
288    
289             else
290                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
291             endif
292    
293    
294    *     ======================================
295    *     temporary patch for saturated clusters
296    *     ======================================
297             if( nsatstrips(ic).gt.0 )then
298                corr    = cog(4,ic)            
299                res = pitchY*1e-4/sqrt(12.)
300    cc            cc=cog(4,ic)
301    c$$$            print*,ic,' *** ',cc
302    c$$$            print*,ic,' *** ',res
303             endif
304            
305          endif
306          end
307    
308    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
309        integer function npfastrips(ic,angle)        integer function npfastrips(ic,angle)
310  *--------------------------------------------------------------  *--------------------------------------------------------------
311  *     thid function returns the number of strips used  *     thid function returns the number of strips used
# Line 359  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 613  c      if(mod(int(VIEW(ic)),2).eq.1)then
613        cog2 = cog(2,ic)                  cog2 = cog(2,ic)          
614        pfaeta2=cog2        pfaeta2=cog2
615    
616    *     ----------------
617  *     find angular bin  *     find angular bin
618    *     ----------------
619  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
620        do iang=1,nangbin        do iang=1,nangbin
621           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
# Line 369  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 625  c      if(mod(int(VIEW(ic)),2).eq.1)then
625        enddo        enddo
626        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
627       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
628        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
629        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
630   98   continue                  !jump here if ok   98   continue                  !jump here if ok
631    
632    *     -------------
633    *     within +/-0.5
634    *     -------------
635    
636  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  
 c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2  
 c$$$*         goto 100  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
637        iadd=0        iadd=0
638   10   continue   10   continue
639        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
640           cog2 = cog2 + 1           cog2 = cog2 + 1
641           iadd = iadd + 1           iadd = iadd + 1
642             if(iadd>iaddmax)goto 111
643           goto 10           goto 10
644        endif        endif
645   20   continue   20   continue
646        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
647           cog2 = cog2 - 1           cog2 = cog2 - 1
648           iadd = iadd - 1           iadd = iadd - 1
649             if(iadd<-1*iaddmax)goto 111
650           goto 20           goto 20
651        endif        endif
652          goto 1111
653     111  continue
654          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
655          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
656          cog2=0
657     1111 continue
658    
659  *     --------------------------------  *     --------------------------------
660  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 473  c$$$      endif Line 719  c$$$      endif
719    
720        iview = VIEW(ic)                    iview = VIEW(ic)            
721        lad = nld(MAXS(ic),VIEW(ic))        lad = nld(MAXS(ic),VIEW(ic))
722        cog3 = cog(3,ic)                    cog3 = cog(3,ic)
723          cc = cog3
724          cog3 = cc
725        pfaeta3=cog3        pfaeta3=cog3
726    
727    *     ----------------
728  *     find angular bin  *     find angular bin
729    *     ----------------
730  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
731        do iang=1,nangbin        do iang=1,nangbin
732  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 487  c         print*,'~~~~~~~~~~~~ ',iang,an Line 737  c         print*,'~~~~~~~~~~~~ ',iang,an
737        enddo        enddo
738        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
739       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
740        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
741        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
742   98   continue                  !jump here if ok   98   continue                  !jump here if ok
743    
744    *     -------------
745    *     within +/-0.5
746    *     -------------
747    
748  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
749        iadd=0        iadd=0
750   10   continue   10   continue
751        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
752           cog3 = cog3 + 1           cog3   = cog3 + 1.
753           iadd = iadd + 1           iadd = iadd + 1
754             if(iadd>iaddmax) goto 111
755           goto 10           goto 10
756        endif        endif
757   20   continue   20   continue
758        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
759           cog3 = cog3 - 1           cog3 = cog3 - 1.
760           iadd = iadd - 1           iadd = iadd - 1
761             if(iadd<-1*iaddmax) goto 111
762           goto 20           goto 20
763        endif        endif
764          goto 1111
765     111  continue
766          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
767          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
768          cog3=0      
769     1111 continue
770    
771  *     --------------------------------  *     --------------------------------
772  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 551  c            print*,'-----',x1,x2,y1,y2 Line 792  c            print*,'-----',x1,x2,y1,y2
792        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
793        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
794    
 c$$$      if(iflag.eq.1)then  
 c$$$         pfaeta2=pfaeta2-1.   !temp  
 c$$$         cog2=cog2-1.           !temp  
 c$$$      endif  
 c$$$      if(iflag.eq.-1)then  
 c$$$         pfaeta2=pfaeta2+1.   !temp  
 c$$$         cog2=cog2+1.           !temp  
 c$$$      endif  
795    
796        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
797       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
# Line 592  c$$$      endif Line 825  c$$$      endif
825        cog4=cog(4,ic)        cog4=cog(4,ic)
826        pfaeta4=cog4        pfaeta4=cog4
827    
828    *     ----------------
829  *     find angular bin  *     find angular bin
830    *     ----------------
831  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
832        do iang=1,nangbin        do iang=1,nangbin
833  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 603  c         print*,'~~~~~~~~~~~~ ',iang,an Line 838  c         print*,'~~~~~~~~~~~~ ',iang,an
838        enddo        enddo
839        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
840       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
841        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
842        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
843   98   continue                  !jump here if ok   98   continue                  !jump here if ok
844    
845    *     -------------
846    *     within +/-0.5
847    *     -------------
848    
849  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
850        iadd=0        iadd=0
851   10   continue   10   continue
852        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
853           cog4 = cog4 + 1           cog4 = cog4 + 1
854           iadd = iadd + 1           iadd = iadd + 1
855             if(iadd>iaddmax)goto 111
856           goto 10           goto 10
857        endif        endif
858   20   continue   20   continue
859        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
860           cog4 = cog4 - 1           cog4 = cog4 - 1
861           iadd = iadd - 1           iadd = iadd - 1
862             if(iadd<-1*iaddmax)goto 111
863           goto 20           goto 20
864        endif        endif
865          goto 1111
866     111  continue
867          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
868          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
869          cog4=0
870     1111 continue
871    
872  *     --------------------------------  *     --------------------------------
873  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 744  c$$$      endif Line 970  c$$$      endif
970                    
971           COG = 0.           COG = 0.
972                    
973  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
974    
975  c     ==============================================================  c     ==============================================================
976           if(ncog.eq.1)then           if(ncog.eq.1)then
977              COG = 0.              COG = 0.
978              if(sr1.gt.sc)cog=1. !NEW              if(sr1.gt.sc)cog=1.
979              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
980  c     ==============================================================  c     ==============================================================
981           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
982                COG = 0.
983              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
984                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
985              elseif(sl1.lt.sr1)then              elseif(sl1.lt.sr1)then
986                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                                         if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
987              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
988                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
989       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
990                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
991       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
992              endif              endif
993  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
994  c     $           ,' : ',sl2,sl1,sc,sr1,sr2  c     $           ,' : ',sl2,sl1,sc,sr1,sr2
995  c     ==============================================================  c     ==============================================================
996           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
997               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
998  c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)              sss = sc
999                if( sl1.ne.-9999. )COG = COG-sl1
1000                if( sl1.ne.-9999. )sss = sss+sl1
1001                if( sr1.ne.-9999. )COG = COG+sr1
1002                if( sr1.ne.-9999. )sss = sss+sr1
1003                if(sss.ne.0)COG=COG/sss
1004    
1005    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1006    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1007  c     $            ,' : ',sl2,sl1,sc,sr1,sr2  c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1008  c     ==============================================================  c     ==============================================================
1009           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1010    
1011                COG = 0
1012                sss = sc
1013                if( sl1.ne.-9999. )COG = COG-sl1
1014                if( sl1.ne.-9999. )sss = sss+sl1
1015                if( sr1.ne.-9999. )COG = COG+sr1
1016                if( sr1.ne.-9999. )sss = sss+sr1
1017              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1018                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1019       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1020              elseif(sl2.lt.sr2)then              elseif(sl2.lt.sr2)then
1021                 if((sr2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1022       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1023              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1024                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1025       $              .and.(sl2+sl1+sc+sr1).ne.0 )       $              .and.(sl2+sss).ne.0 )
1026       $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW       $              cog = (cog-2*sl2)/(sl2+sss)
1027                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1028       $              .and.(sr2+sl1+sc+sr1).ne.0 )       $              .and.(sr2+sss).ne.0 )
1029       $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW                     $              cog = (2*sr2+cog)/(sr2+sss)              
1030              endif              endif
1031    c     ==============================================================
1032             elseif(ncog.eq.5)then
1033                COG = 0
1034                sss = sc
1035                if( sl1.ne.-9999. )COG = COG-sl1
1036                if( sl1.ne.-9999. )sss = sss+sl1
1037                if( sr1.ne.-9999. )COG = COG+sr1
1038                if( sr1.ne.-9999. )sss = sss+sr1
1039                if( sl2.ne.-9999. )COG = COG-2*sl2
1040                if( sl2.ne.-9999. )sss = sss+sl2
1041                if( sr2.ne.-9999. )COG = COG+2*sr2
1042                if( sr2.ne.-9999. )sss = sss+sr2
1043                if(sss.ne.0)COG=COG/sss
1044           else           else
1045              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1046       $           ,' not implemented'       $           ,' not implemented'
# Line 856  c         print*,'-------' Line 1111  c         print*,'-------'
1111    
1112  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1113    
1114          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1115             if(DEBUG.eq.1)
1116         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1117             if(DEBUG.eq.1)
1118         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1119          endif
1120    
1121    
1122        return        return
1123        end        end
1124    
# Line 1679  c$$$ Line 1942  c$$$
1942        enddo        enddo
1943        if(DEBUG.eq.1)        if(DEBUG.eq.1)
1944       $     print*,'pfacorr *** warning *** angle out of range: ',angle       $     print*,'pfacorr *** warning *** angle out of range: ',angle
1945        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
1946        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
1947   98   continue                  !jump here if ok   98   continue                  !jump here if ok
1948    
1949        pfacorr = fcorr(iview,lad,iang)        pfacorr = fcorr(iview,lad,iang)
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23