/[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.22 by pam-fi, Tue Sep 4 09:47:49 2007 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             tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001
74    
75          elseif(mod(iview,2).eq.1)then
76    c     =================================================
77    c     Y view
78    c     =================================================        
79    c     here bbb is the x component of the m.filed
80             angy = ang
81             bx   = bbb
82             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
83    
84          endif      
85          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
86    
87          return
88          end
89    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
90          real function fieldcorr(iview,bbb)
91    
92          include 'commontracker.f'
93    
94          fieldcorr = 0.
95    
96          if(mod(iview,2).eq.0)then
97    
98    c     =================================================
99    c     X view
100    c     =================================================
101    c     here bbb is the y component of the m.field
102             by   = bbb
103             if(iview.eq.12) by = -1. * bbb
104             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
105    
106          elseif(mod(iview,2).eq.1)then
107    c     =================================================
108    c     Y view
109    c     =================================================        
110    c     here bbb is the x component of the m.filed
111             bx   = bbb
112             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
113    
114          endif      
115          
116          return
117          end
118    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
119    
120          subroutine applypfa(PFAtt,ic,ang,corr,res)
121    *---------------------------------------------------------------
122    *     this subroutine calculate the coordinate of cluster ic (in
123    *     strip units), relative to the strip with the maximum signal,
124    *     and its spatial resolution (in cm), applying PFAtt.
125    *     ang is the effective angle, relative to the sensor
126    *---------------------------------------------------------------
127    
128          character*4 PFAtt
129          include 'commontracker.f'
130          include 'level1.f'
131    
132          corr = 0
133          res  = 0
134    
135          if(ic.le.0)return
136    
137          iview   = VIEW(ic)
138    
139          if(mod(iview,2).eq.0)then
140    c     =================================================
141    c     X view
142    c     =================================================
143    
144             res = RESXAV
145    
146             if(PFAtt.eq.'COG1')then
147    
148                corr   = 0
149                res = 1e-4*pitchX/sqrt(12.)!!res
150    
151             elseif(PFAtt.eq.'COG2')then
152    
153                corr    = cog(2,ic)            
154                res = risx_cog(abs(ang))!TEMPORANEO              
155                res = res*fbad_cog(2,ic)
156    
157             elseif(PFAtt.eq.'COG3')then
158    
159                corr    = cog(3,ic)            
160                res = risx_cog(abs(ang))!TEMPORANEO                      
161                res = res*fbad_cog(3,ic)
162    
163             elseif(PFAtt.eq.'COG4')then
164    
165                corr    = cog(4,ic)            
166                res = risx_cog(abs(ang))!TEMPORANEO                      
167                res = res*fbad_cog(4,ic)
168    
169             elseif(PFAtt.eq.'ETA2')then
170    
171                corr  = pfaeta2(ic,ang)          
172                res = risxeta2(abs(ang))
173                res = res*fbad_cog(2,ic)
174    
175             elseif(PFAtt.eq.'ETA3')then                        
176    
177                corr  = pfaeta3(ic,ang)          
178                res = risxeta3(abs(ang))                      
179                res = res*fbad_cog(3,ic)              
180    
181             elseif(PFAtt.eq.'ETA4')then                        
182    
183                corr  = pfaeta4(ic,ang)            
184                res = risxeta4(abs(ang))                      
185                res = res*fbad_cog(4,ic)              
186    
187             elseif(PFAtt.eq.'ETA')then  
188    
189                corr  = pfaeta(ic,ang)            
190    c            res = riseta(ic,ang)                    
191                res = riseta(iview,ang)                    
192                res = res*fbad_eta(ic,ang)            
193    
194             elseif(PFAtt.eq.'ETAL')then  
195    
196                corr  = pfaetal(ic,ang)            
197                res = riseta(iview,ang)                    
198                res = res*fbad_eta(ic,ang)            
199    
200             elseif(PFAtt.eq.'COG')then          
201    
202                corr  = cog(0,ic)            
203                res = risx_cog(abs(ang))                    
204                res = res*fbad_cog(0,ic)
205    
206             else
207                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
208             endif
209    
210    
211    *     ======================================
212    *     temporary patch for saturated clusters
213    *     ======================================
214             if( nsatstrips(ic).gt.0 )then
215                corr  = cog(4,ic)            
216                res = pitchX*1e-4/sqrt(12.)
217    cc            cc=cog(4,ic)
218    c$$$            print*,ic,' *** ',cc
219    c$$$            print*,ic,' *** ',res
220             endif
221    
222    
223          elseif(mod(iview,2).eq.1)then
224    c     =================================================
225    c     Y view
226    c     =================================================
227    
228             res = RESYAV
229    
230             if(PFAtt.eq.'COG1')then
231    
232                corr  = 0  
233                res = 1e-4*pitchY/sqrt(12.)!res  
234    
235             elseif(PFAtt.eq.'COG2')then
236    
237                corr    = cog(2,ic)
238                res = risy_cog(abs(ang))!TEMPORANEO
239                res = res*fbad_cog(2,ic)
240    
241             elseif(PFAtt.eq.'COG3')then
242    
243                corr    = cog(3,ic)
244                res = risy_cog(abs(ang))!TEMPORANEO
245                res = res*fbad_cog(3,ic)
246    
247             elseif(PFAtt.eq.'COG4')then
248    
249                corr    = cog(4,ic)
250                res = risy_cog(abs(ang))!TEMPORANEO
251                res = res*fbad_cog(4,ic)
252    
253             elseif(PFAtt.eq.'ETA2')then
254    
255                corr    = pfaeta2(ic,ang)          
256                res = risyeta2(abs(ang))              
257                res = res*fbad_cog(2,ic)
258    
259             elseif(PFAtt.eq.'ETA3')then                      
260    
261                corr    = pfaeta3(ic,ang)
262                res = res*fbad_cog(3,ic)  
263    
264             elseif(PFAtt.eq.'ETA4')then  
265    
266                corr    = pfaeta4(ic,ang)
267                res = res*fbad_cog(4,ic)
268    
269             elseif(PFAtt.eq.'ETA')then
270    
271                corr    = pfaeta(ic,ang)
272    c            res = riseta(ic,ang)  
273                res = riseta(iview,ang)  
274                res = res*fbad_eta(ic,ang)
275    
276             elseif(PFAtt.eq.'ETAL')then
277    
278                corr    = pfaetal(ic,ang)
279                res = riseta(iview,ang)  
280                res = res*fbad_eta(ic,ang)
281    
282             elseif(PFAtt.eq.'COG')then
283    
284                corr    = cog(0,ic)            
285                res = risy_cog(abs(ang))
286                res = res*fbad_cog(0,ic)
287    
288             else
289                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
290             endif
291    
292    
293    *     ======================================
294    *     temporary patch for saturated clusters
295    *     ======================================
296             if( nsatstrips(ic).gt.0 )then
297                corr    = cog(4,ic)            
298                res = pitchY*1e-4/sqrt(12.)
299    cc            cc=cog(4,ic)
300    c$$$            print*,ic,' *** ',cc
301    c$$$            print*,ic,' *** ',res
302             endif
303            
304          endif
305          end
306    
307    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
308        integer function npfastrips(ic,angle)        integer function npfastrips(ic,angle)
309  *--------------------------------------------------------------  *--------------------------------------------------------------
310  *     thid function returns the number of strips used  *     thid function returns the number of strips used
# Line 369  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 622  c      if(mod(int(VIEW(ic)),2).eq.1)then
622        enddo        enddo
623        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
624       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
625        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
626        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
627   98   continue                  !jump here if ok   98   continue                  !jump here if ok
628    
629    
# Line 487  c         print*,'~~~~~~~~~~~~ ',iang,an Line 740  c         print*,'~~~~~~~~~~~~ ',iang,an
740        enddo        enddo
741        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
742       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
743        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
744        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
745   98   continue                  !jump here if ok   98   continue                  !jump here if ok
746    
747    
# Line 603  c         print*,'~~~~~~~~~~~~ ',iang,an Line 856  c         print*,'~~~~~~~~~~~~ ',iang,an
856        enddo        enddo
857        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
858       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
859        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
860        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
861   98   continue                  !jump here if ok   98   continue                  !jump here if ok
862    
863    
# Line 749  c         print*,'## ',sl2,sl1,sc,sr1,sr Line 1002  c         print*,'## ',sl2,sl1,sc,sr1,sr
1002  c     ==============================================================  c     ==============================================================
1003           if(ncog.eq.1)then           if(ncog.eq.1)then
1004              COG = 0.              COG = 0.
1005              if(sr1.gt.sc)cog=1. !NEW              if(sr1.gt.sc)cog=1.
1006              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1007  c     ==============================================================  c     ==============================================================
1008           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
1009                COG = 0.
1010              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1011                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1012              elseif(sl1.lt.sr1)then              elseif(sl1.lt.sr1)then
1013                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                                         if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
1014              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1015                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1016       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1017                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1018       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1019              endif              endif
1020  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1021  c     $           ,' : ',sl2,sl1,sc,sr1,sr2  c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1022  c     ==============================================================  c     ==============================================================
1023           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1024               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
1025  c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)              sss = sc
1026                if( sl1.ne.-9999. )COG = COG-sl1
1027                if( sl1.ne.-9999. )sss = sss+sl1
1028                if( sr1.ne.-9999. )COG = COG+sr1
1029                if( sr1.ne.-9999. )sss = sss+sr1
1030                if(sss.ne.0)COG=COG/sss
1031    
1032    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1033    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1034  c     $            ,' : ',sl2,sl1,sc,sr1,sr2  c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1035  c     ==============================================================  c     ==============================================================
1036           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1037    
1038                COG = 0
1039                sss = sc
1040                if( sl1.ne.-9999. )COG = COG-sl1
1041                if( sl1.ne.-9999. )sss = sss+sl1
1042                if( sr1.ne.-9999. )COG = COG+sr1
1043                if( sr1.ne.-9999. )sss = sss+sr1
1044              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1045                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1046       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1047              elseif(sl2.lt.sr2)then              elseif(sl2.lt.sr2)then
1048                 if((sr2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1049       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1050              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1051                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1052       $              .and.(sl2+sl1+sc+sr1).ne.0 )       $              .and.(sl2+sss).ne.0 )
1053       $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW       $              cog = (cog-2*sl2)/(sl2+sss)
1054                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1055       $              .and.(sr2+sl1+sc+sr1).ne.0 )       $              .and.(sr2+sss).ne.0 )
1056       $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW                     $              cog = (2*sr2+cog)/(sr2+sss)              
1057              endif              endif
1058    c     ==============================================================
1059             elseif(ncog.eq.5)then
1060                COG = 0
1061                sss = sc
1062                if( sl1.ne.-9999. )COG = COG-sl1
1063                if( sl1.ne.-9999. )sss = sss+sl1
1064                if( sr1.ne.-9999. )COG = COG+sr1
1065                if( sr1.ne.-9999. )sss = sss+sr1
1066                if( sl2.ne.-9999. )COG = COG-2*sl2
1067                if( sl2.ne.-9999. )sss = sss+sl2
1068                if( sr2.ne.-9999. )COG = COG+2*sr2
1069                if( sr2.ne.-9999. )sss = sss+sr2
1070                if(sss.ne.0)COG=COG/sss
1071           else           else
1072              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1073       $           ,' not implemented'       $           ,' not implemented'
# Line 1679  c$$$ Line 1961  c$$$
1961        enddo        enddo
1962        if(DEBUG.eq.1)        if(DEBUG.eq.1)
1963       $     print*,'pfacorr *** warning *** angle out of range: ',angle       $     print*,'pfacorr *** warning *** angle out of range: ',angle
1964        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
1965        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
1966   98   continue                  !jump here if ok   98   continue                  !jump here if ok
1967    
1968        pfacorr = fcorr(iview,lad,iang)        pfacorr = fcorr(iview,lad,iang)
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23