/[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.11 by pam-fi, Tue May 15 16:22:18 2007 UTC revision 1.15 by pam-fi, Fri Aug 17 13:28:02 2007 UTC
# Line 170  c      real function riseta(ic,angle) Line 170  c      real function riseta(ic,angle)
170  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
171  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
172  *     it calls:  *     it calls:
173  *     - risx_eta2(angle)  *     - risxeta2(angle)
174  *     - risy_eta2(angle)  *     - risyeta2(angle)
175  *     - risx_eta3(angle)  *     - risxeta3(angle)
176  *     - risx_eta4(angle)  *     - risxeta4(angle)
177  *     according to the angle  *     according to the angle
178  *--------------------------------------------------------------  *--------------------------------------------------------------
179        include 'commontracker.f'        include 'commontracker.f'
# Line 187  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 187  c      if(mod(int(VIEW(ic)),2).eq.1)then
187                
188    
189           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
190              riseta = risy_eta2(angle)              riseta = risyeta2(angle)
191           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
192              riseta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
193           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
# Line 199  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 199  c      if(mod(int(VIEW(ic)),2).eq.1)then
199        else                      !X-view        else                      !X-view
200    
201           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
202              riseta = risx_eta2(angle)              riseta = risxeta2(angle)
203           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
204              riseta = risx_eta3(angle)              riseta = risxeta3(angle)
205           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
206              riseta = risx_eta4(angle)              riseta = risxeta4(angle)
207           else           else
208              riseta = risx_cog(angle)              riseta = risx_cog(angle)
209           endif                       endif            
210                            
211        endif        endif
212    
       print*,'---- ',riseta,iview,angle  
213    
214   100  return   100  return
215        end        end
# Line 612  c$$$      endif Line 611  c$$$      endif
611    
612    
613  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
614        real function cog0(ncog,ic)  c$$$      real function cog0(ncog,ic)
615  *-------------------------------------------------  c$$$*-------------------------------------------------
616  *     this function returns  c$$$*     this function returns
617  *  c$$$*
618  *     - the Center-Of-Gravity of the cluster IC  c$$$*     - the Center-Of-Gravity of the cluster IC
619  *     evaluated using NCOG strips,  c$$$*     evaluated using NCOG strips,
620  *     calculated relative to MAXS(IC)  c$$$*     calculated relative to MAXS(IC)
621  *      c$$$*    
622  *     - zero in case that not  enough strips  c$$$*     - zero in case that not  enough strips
623  *     have a positive signal  c$$$*     have a positive signal
624  *      c$$$*    
625  *     NOTE:  c$$$*     NOTE:
626  *     This is the old definition, used by Straulino.  c$$$*     This is the old definition, used by Straulino.
627  *     The new routine, according to Landi,  c$$$*     The new routine, according to Landi,
628  *     is COG(NCOG,IC)  c$$$*     is COG(NCOG,IC)
629  *-------------------------------------------------  c$$$*-------------------------------------------------
630    c$$$
631    c$$$
632        include 'commontracker.f'  c$$$      include 'commontracker.f'
633        include 'level1.f'  c$$$      include 'level1.f'
634          c$$$      
635  *     --> signal of the central strip  c$$$*     --> signal of the central strip
636        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
637    c$$$
638  *     signal of adjacent strips  c$$$*     signal of adjacent strips
639  *     --> left  c$$$*     --> left
640        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
641        if(  c$$$      if(
642       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
643       $     )  c$$$     $     )
644       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
645    c$$$
646        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
647        if(  c$$$      if(
648       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
649       $     )  c$$$     $     )
650       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
651    c$$$
652  *     --> right  c$$$*     --> right
653        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
654        if(  c$$$      if(
655       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
656       $     .or.  c$$$     $     .or.
657       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
658       $     )  c$$$     $     )
659       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
660    c$$$
661        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
662        if(  c$$$      if(
663       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
664       $     .or.  c$$$     $     .or.
665       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
666       $     )  c$$$     $     )
667       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
668          c$$$      
669  ************************************************************  c$$$************************************************************
670  *     COG computation  c$$$*     COG computation
671  ************************************************************  c$$$************************************************************
672    c$$$
673  c      print*,sl2,sl1,sc,sr1,sr2  c$$$c      print*,sl2,sl1,sc,sr1,sr2
674    c$$$
675        COG = 0.  c$$$      COG = 0.
676            c$$$        
677        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
678            c$$$        
679           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
680              COG = -sl1/(sl1+sc)          c$$$            COG = -sl1/(sl1+sc)        
681           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
682              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
683           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
684              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
685           else  c$$$         else
686              COG = 0.  c$$$            COG = 0.
687           endif  c$$$         endif
688            c$$$        
689        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
690    c$$$
691           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
692              COG = sr1/(sc+sr1)              c$$$            COG = sr1/(sc+sr1)            
693           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
694              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
695           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
696              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
697           else  c$$$         else
698              COG = 0.  c$$$            COG = 0.
699           endif  c$$$         endif
700    c$$$
701        endif  c$$$      endif
702    c$$$
703        COG0 = COG  c$$$      COG0 = COG
704    c$$$
705  c      print *,ncog,ic,cog,'/////////////'  c$$$c      print *,ncog,ic,cog,'/////////////'
706    c$$$
707        return  c$$$      return
708        end  c$$$      end
709    
710  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
711        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 737  c      print *,ncog,ic,cog,'////////////
737  *     --> signal of the central strip  *     --> signal of the central strip
738           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
739  *     signal of adjacent strips  *     signal of adjacent strips
740           sl1 = 0                !left 1           sl1 = -9999.           !left 1
741           if(           if(
742       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
743       $        )       $        )
744       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
745                    
746           sl2 = 0                !left 2           sl2 = -9999.           !left 2
747           if(           if(
748       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
749       $        )       $        )
750       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
751                    
752           sr1 = 0                !right 1           sr1 = -9999.           !right 1
753           if(           if(
754       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
755       $        .or.       $        .or.
# Line 758  c      print *,ncog,ic,cog,'//////////// Line 757  c      print *,ncog,ic,cog,'////////////
757       $        )       $        )
758       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
759                    
760           sr2 = 0                !right 2           sr2 = -9999.           !right 2
761           if(           if(
762       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
763       $        .or.       $        .or.
# Line 770  c      print *,ncog,ic,cog,'//////////// Line 769  c      print *,ncog,ic,cog,'////////////
769                    
770  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c         print*,'## ',sl2,sl1,sc,sr1,sr2
771    
772    c     ==============================================================
773           if(ncog.eq.1)then           if(ncog.eq.1)then
774              COG = 0.              COG = 0.
775                if(sr1.gt.sc)cog=1. !NEW
776                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
777    c     ==============================================================
778           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
779              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
780                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
781              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
782                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
783              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
784                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
785         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
786                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
787         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
788                endif
789    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
790    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
791    c     ==============================================================
792           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
793               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
794    c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
795    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
796    c     ==============================================================
797           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
798              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
799                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sl1+sc+sr1).ne.0)
800       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
801              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
802                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sl1+sc+sr1).ne.0)
803       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
804                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
805                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
806         $              .and.(sl2+sl1+sc+sr1).ne.0 )
807         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
808                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
809         $              .and.(sr2+sl1+sc+sr1).ne.0 )
810         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
811              endif              endif
812           else           else
813              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
# Line 818  c         print*,'-------' Line 839  c         print*,'-------'
839                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
840                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
841                 mu = mu + 1                 mu = mu + 1
842                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
843              else              else
844                 goto 10                 goto 10
845              endif              endif
# Line 831  c         print*,'-------' Line 852  c         print*,'-------'
852                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
853                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
854                 mu = mu + 1                 mu = mu + 1
855                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
856              else              else
857                 goto 20                 goto 20
858              endif              endif
859           enddo           enddo
860   20      continue   20      continue
861           if(SGN.le.0)then           if(SGN.le.0)then
862  c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
863              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
864              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
865  c            print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
# Line 1193  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1214  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1214    
1215  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1216    
1217        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1218    
1219        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1220        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1280  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1301  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1301     20 CONTINUE     20 CONTINUE
1302        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1303                
1304        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1305    
1306        END        END
1307    
1308  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1309        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1310        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1311        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1312        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1371  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1392  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1392     20 CONTINUE     20 CONTINUE
1393        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1394                
1395        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1396    
1397        END        END
1398  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1399        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1400        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1401        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1402        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1461  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1482  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1482     20 CONTINUE     20 CONTINUE
1483        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1484                
1485        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1486    
1487        END        END
1488  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1489        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1490        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1491        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1492        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1533  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1554  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1554     20 CONTINUE     20 CONTINUE
1555        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1556    
1557        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1558    
1559        END        END
1560  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.23