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

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

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

revision 1.2 by pam-fi, Tue May 30 16:30:37 2006 UTC revision 1.7 by pam-fi, Thu Jan 11 10:20:58 2007 UTC
# Line 7  Line 7 
7    
8    
9  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
10        real function pfa_eta(ic,angle)        real function pfaeta(ic,angle)
11  *--------------------------------------------------------------  *--------------------------------------------------------------
12  *     this function returns the position (in strip units)  *     this function returns the position (in strip units)
13  *     it calls:  *     it calls:
14  *     - pfa_eta2(ic,angle)  *     - pfaeta2(ic,angle)
15  *     - pfa_eta3(ic,angle)  *     - pfaeta3(ic,angle)
16  *     - pfa_eta4(ic,angle)  *     - pfaeta4(ic,angle)
17  *     according to the angle  *     according to the angle
18  *--------------------------------------------------------------  *--------------------------------------------------------------
19        include 'commontracker.f'        include 'commontracker.f'
20  c      include 'calib.f'  c      include 'calib.f'
21        include 'level1.f'        include 'level1.f'
22                
23        pfa_eta = 0        pfaeta = 0
24    
25        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
26                
27           pfa_eta = pfa_eta2(ic,angle)           pfaeta = pfaeta2(ic,angle)
28        
29        else                      !X-view        else                      !X-view
30    
31           if(abs(angle).le.10.)then           if(abs(angle).le.10.)then
32              pfa_eta = pfa_eta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
33           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then
34              pfa_eta = pfa_eta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
35           elseif(abs(angle).gt.15.)then           elseif(abs(angle).gt.15.)then
36              pfa_eta = pfa_eta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
37           endif           endif
38                            
39        endif        endif
40          
41    c      print*,'pfaeta ',pfaeta, angle
42    
43   100  return   100  return
44        end        end
# Line 46  c      include 'calib.f' Line 47  c      include 'calib.f'
47        real function ris_eta(ic,angle)        real function ris_eta(ic,angle)
48  *--------------------------------------------------------------  *--------------------------------------------------------------
49  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
50  *     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
51  *     it calls:  *     it calls:
52  *     - risx_eta2(angle)  *     - risx_eta2(angle)
53  *     - risy_eta2(angle)  *     - risy_eta2(angle)
# Line 129  c$$$     $     ,' -->',ris_eta Line 130  c$$$     $     ,' -->',ris_eta
130  c*****************************************************  c*****************************************************
131  cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
132  c*****************************************************  c*****************************************************
133  c      real function pfa_eta2(cog2,view,lad,angle)  c      real function pfaeta2(cog2,view,lad,angle)
134        real function pfa_eta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
135  *--------------------------------------------------------------  *--------------------------------------------------------------
136  *     this function returns  *     this function returns
137  *  *
# Line 152  c      real function pfa_eta2(cog2,view, Line 153  c      real function pfa_eta2(cog2,view,
153  c      logical DEBUG  c      logical DEBUG
154  c      common/dbg/DEBUG  c      common/dbg/DEBUG
155    
156    c      print*,'## pfaeta2 ',ic,angle
157        iview = VIEW(ic)              !(1)        iview = VIEW(ic)              !(1)
158        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic)) !(1)
159        cog2 = cog(2,ic)             !(1)        cog2 = cog(2,ic)             !(1)
160        pfa_eta2=cog2        pfaeta2=cog2
161    
162  *     find angular bin  *     find angular bin
163  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 167  c         print*,'~~~~~~~~~~~~ ',iang,an Line 169  c         print*,'~~~~~~~~~~~~ ',iang,an
169           endif           endif
170        enddo        enddo
171        if(DEBUG)        if(DEBUG)
172       $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
173        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
174        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
175   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 177  c$$$*     find extremes of interpolation Line 179  c$$$*     find extremes of interpolation
179  c$$$      iflag=0  c$$$      iflag=0
180  c$$$*     --------------------------------  c$$$*     --------------------------------
181  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
182  c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2  c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2
183  c$$$*         goto 100  c$$$*         goto 100
184  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
185  c$$$*     non salto piu`, ma scalo di 1 o -1  c$$$*     non salto piu`, ma scalo di 1 o -1
# Line 230  c            print*,'-----',x1,x2,y1,y2 Line 232  c            print*,'-----',x1,x2,y1,y2
232        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
233        BB=y1-AA*x1        BB=y1-AA*x1
234    
235        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
236        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
237    
238  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
239  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
240  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
241  c$$$      endif  c$$$      endif
242  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
243  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
244  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
245  c$$$      endif  c$$$      endif
246    
247        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'
248       $     ,cog2-iadd,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
249    
250    
251   100  return   100  return
# Line 253  c$$$      endif Line 255  c$$$      endif
255  c*****************************************************  c*****************************************************
256  cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
257  c*****************************************************  c*****************************************************
258  c      real function pfa_eta3(cog3,view,lad,angle)  c      real function pfaeta3(cog3,view,lad,angle)
259        real function pfa_eta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
260  *--------------------------------------------------------------  *--------------------------------------------------------------
261  *     this function returns  *     this function returns
262  *  *
# Line 276  c      real function pfa_eta3(cog3,view, Line 278  c      real function pfa_eta3(cog3,view,
278  c      logical DEBUG  c      logical DEBUG
279  c      common/dbg/DEBUG  c      common/dbg/DEBUG
280    
281    c      print*,'## pfaeta3 ',ic,angle
282    
283        iview = VIEW(ic)              !(1)        iview = VIEW(ic)              !(1)
284        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic)) !(1)
285        cog3 = cog(3,ic)             !(1)        cog3 = cog(3,ic)             !(1)
286        pfa_eta3=cog3        pfaeta3=cog3
287    
288  *     find angular bin  *     find angular bin
289  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 292  c         print*,'~~~~~~~~~~~~ ',iang,an Line 295  c         print*,'~~~~~~~~~~~~ ',iang,an
295           endif           endif
296        enddo        enddo
297        if(DEBUG)        if(DEBUG)
298       $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
299        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
300        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
301   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 354  c            print*,'-----',x1,x2,y1,y2 Line 357  c            print*,'-----',x1,x2,y1,y2
357        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
358        BB=y1-AA*x1        BB=y1-AA*x1
359    
360        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
361        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
362    
363  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
364  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
365  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
366  c$$$      endif  c$$$      endif
367  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
368  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
369  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
370  c$$$      endif  c$$$      endif
371    
372        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'
373       $     ,cog3-iadd,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
374    
375   100  return   100  return
376        end        end
# Line 376  c$$$      endif Line 379  c$$$      endif
379  c*****************************************************  c*****************************************************
380  cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  cccccc 02/02/2006 modified by Elena Vannuccini --> (1)
381  c*****************************************************  c*****************************************************
382  c      real function pfa_eta4(cog4,view,lad,angle)  c      real function pfaeta4(cog4,view,lad,angle)
383        real function pfa_eta4(ic,angle) !(1)        real function pfaeta4(ic,angle) !(1)
384  *--------------------------------------------------------------  *--------------------------------------------------------------
385  *     this function returns  *     this function returns
386  *  *
# Line 399  c      real function pfa_eta4(cog4,view, Line 402  c      real function pfa_eta4(cog4,view,
402  c      logical DEBUG  c      logical DEBUG
403  c      common/dbg/DEBUG  c      common/dbg/DEBUG
404    
405        iview = VIEW(ic)              !(1)  c      print*,'## pfaeta4 ',ic,angle
406    
407          iview = VIEW(ic)             !(1)
408        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic)) !(1)
409        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)               !(1)
410        pfa_eta4=cog4        pfaeta4=cog4
411    
412  *     find angular bin  *     find angular bin
413  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 414  c         print*,'~~~~~~~~~~~~ ',iang,an Line 419  c         print*,'~~~~~~~~~~~~ ',iang,an
419           endif           endif
420        enddo        enddo
421        if(DEBUG)        if(DEBUG)
422       $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
423        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
424        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
425   98   continue                  !jump here if ok   98   continue                  !jump here if ok
# Line 476  c            print*,'-----',x1,x2,y1,y2 Line 481  c            print*,'-----',x1,x2,y1,y2
481        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
482        BB=y1-AA*x1        BB=y1-AA*x1
483    
484        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
485        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
486    
487  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
488  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
489  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
490  c$$$      endif  c$$$      endif
491  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
492  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
493  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
494  c$$$      endif  c$$$      endif
495    
496        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'
497       $     ,cog4-iadd,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
498    
499   100  return   100  return
500        end        end
# Line 655  c      common/dbg/DEBUG Line 660  c      common/dbg/DEBUG
660                    
661           COG = 0.           COG = 0.
662                    
663           if(ncog.eq.2)then  c         print*,'## ',sl2,sl1,sc,sr1,sr2
664    
665             if(ncog.eq.1)then
666                COG = 0.
667             elseif(ncog.eq.2)then
668              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
669                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
670              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
671                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)            
672              endif              endif
673           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
674              COG = (sr1-sl1)/(sl1+sc+sr1)               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)
675           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
676              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
677                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sl1+sc+sr1).ne.0)
678         $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
679              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
680                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)                  if((sl2+sl1+sc+sr1).ne.0)
681         $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
682              endif              endif
683           else           else
684              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
685              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
686              COG = 0.              COG = 0.
687           endif           endif
688    
# Line 686  ' Line 697  '
697           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
698           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
699    
700           istart=INDSTART(IC)           istart = INDSTART(IC)
701           istop=TOTCLLENGTH           istop  = TOTCLLENGTH
702           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
703           COG=0             COG = 0  
704           mu=0           mu  = 0
705           do i=istart,istop           do i = istart,istop
706              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
707              ivk=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
 *            print*,'******************',istart,istop,ipos  
 *     $           ,MAXS(ic)+ipos,iv,ivk,is  
             cut=incut*SIGMA(iv,ivk,is)  
708              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
709                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
710                 mu = mu + 1                 mu = mu + 1
 c               print*,ipos,CLSIGNAL(i),incut,cut  
711              endif              endif
712           enddo           enddo
713           COG=COG/DEDX(ic)           if(DEDX(ic).le.0)then
714  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)
715  c     $        ,cog              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
716                print*,(CLSIGNAL(i),i=istart,istop)
717                print*,'cog(0,ic) --> NOT EVALUATED '
718             else
719                COG=COG/DEDX(ic)
720             endif
721                    
722        else        else
723                    
# Line 717  c     $        ,cog Line 728  c     $        ,cog
728    
729        endif        endif
730    
731  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
732    
733        return        return
734        end        end
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23