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

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

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

revision 1.22 by bongi, Thu Nov 20 15:06:27 2008 UTC revision 1.23 by pam-fi, Fri Dec 5 08:26:47 2008 UTC
# Line 630  c$$$      ENDDO Line 630  c$$$      ENDDO
630  *     measured position of the cluster.  *     measured position of the cluster.
631  *     ---------------------------------------------------------  *     ---------------------------------------------------------
632        CHI2=0.        CHI2=0.
633          DO I=1,nplanes        
634        DO I=1,nplanes                                                     IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
635           IF( XGOOD(I).NE.YGOOD(I) ) THEN ! singlet              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
636              IF(XGOOD(I).EQ.1) THEN              ALFA = XM_A(I) - BETA * YM_A(I)
637                 Z =              YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
638       $         ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) -              if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
639       $           (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/       $           YM(I)=dmin1(YM_A(I),YM_B(I))
640       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
641       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )       $           YM(I)=dmax1(YM_A(I),YM_B(I))
642                 ZM(I) = Z              XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
643                 ZV(I) = Z           ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
644                 XV(I) = XV_A(I)+(XV_B(I)-XV_A(I))*              BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
645       $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))              ALFA = YM_A(I) - BETA * XM_A(I)
646                 Y =              XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
647       $         ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
648       $           (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(YM_A(I)-YM_B(I)) )/       $           XM(I)=dmin1(XM_A(I),XM_B(I))
649       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
650       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )       $           XM(I)=dmax1(XM_A(I),XM_B(I))
651                 YM(I) = Y              YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
652                 YV(I) = Y           ENDIF
653                 XM(I) = XM_A(I)+(XM_B(I)-XM_A(I))*           CHI2=CHI2
654       $                 (Y-YM_A(I))/(YM_B(I)-YM_A(I))       +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
655         +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
656                 CHI2=CHI2+(XV(I)-XM(I))**2/RESX(I)**2       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
657         +                                       *( XGOOD(I)*(1-YGOOD(I)) )
658              ENDIF       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
659              IF(YGOOD(I).EQ.1) THEN       +                                       *( (1-XGOOD(I))*YGOOD(I) )
660                 Z =  c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
661       $         ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) -  c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
662       $           (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/  c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
663       $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -  c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
664       $           (ZV_A(I)-ZV_B(I))*(XM_A(I)-XM_B(I)) )  c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
665                 ZM(I) = Z  c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
666                 ZV(I) = Z  c$$$         print*,XV(I),XM(I),XGOOD(I)
667                 YV(I) = YV_A(I)+(YV_B(I)-YV_A(I))*  c$$$         print*,YV(I),YM(I),YGOOD(I)
      $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))  
                X =  
      $         ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(XV_A(I)-XV_B(I)) -  
      $           (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(XM_A(I)-XM_B(I)) )/  
      $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -  
      $           (ZV_A(I)-ZV_B(I))*(XM_A(I)-XM_B(I)) )  
                XM(I) = X  
                XV(I) = X  
                YM(I) = YM_A(I)+(YM_B(I)-YM_A(I))*  
      $                 (X-XM_A(I))/(XM_B(I)-XM_A(I))  
   
                CHI2=CHI2+(YV(I)-YM(I))**2/RESY(I)**2  
   
             ENDIF  
          ELSEIF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.1)THEN !Y-cl  
             CHI2=CHI2  
      +           +(XV(I)-XM(I))**2/RESX(i)**2  
      +           +(YV(I)-YM(I))**2/RESY(i)**2  
          ENDIF  
       ENDDO  
       DO I=1,nplanes  
          XV0(I)=XV(I)  
          YV0(I)=YV(I)  
668        ENDDO        ENDDO
   
 c$$$      DO I=1,nplanes          
 c$$$         IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl  
 c$$$            BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))  
 c$$$            ALFA = XM_A(I) - BETA * YM_A(I)  
 c$$$            YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)  
 c$$$            if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))  
 c$$$     $           YM(I)=dmin1(YM_A(I),YM_B(I))  
 c$$$            if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))  
 c$$$     $           YM(I)=dmax1(YM_A(I),YM_B(I))  
 c$$$            XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates  
 c$$$         ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl  
 c$$$            BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))  
 c$$$            ALFA = YM_A(I) - BETA * XM_A(I)  
 c$$$            XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)  
 c$$$            if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))  
 c$$$     $           XM(I)=dmin1(XM_A(I),XM_B(I))  
 c$$$            if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))  
 c$$$     $           XM(I)=dmax1(XM_A(I),XM_B(I))  
 c$$$            YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates  
 c$$$         ENDIF  
 c$$$         CHI2=CHI2  
 c$$$     +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )  
 c$$$     +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )  
 c$$$     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2  
 c$$$     +                                       *( XGOOD(I)*(1-YGOOD(I)) )  
 c$$$     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2  
 c$$$     +                                       *( (1-XGOOD(I))*YGOOD(I) )  
 c$$$      ENDDO  
   
669  c$$$      print*,'CHISQ ',chi2  c$$$      print*,'CHISQ ',chi2
670  *     ------------------------------------------------  *     ------------------------------------------------
671  *      *    
# Line 1019  c$$$     $              +DOWN2*U(I)*U(J) Line 966  c$$$     $              +DOWN2*U(I)*U(J)
966  *      *    
967  *****************************************************************  *****************************************************************
968    
 cPPP --- new --- (with singlets in 3D)  
969        SUBROUTINE POSXYZ(AL_P,IFAIL)        SUBROUTINE POSXYZ(AL_P,IFAIL)
970    
971        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# Line 1033  c      COMMON/TRKD/TRKVERBOSE Line 979  c      COMMON/TRKD/TRKVERBOSE
979        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
980  c        c      
981        DIMENSION AL_P(5)        DIMENSION AL_P(5)
       LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER  
982  *      *    
983  cpp      DO I=1,nplanes  cpp      DO I=1,nplanes
984  cpp         ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
# Line 1055  cpp      ENDDO       Line 1000  cpp      ENDDO      
1000  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1001    
1002        DO I=1,nplanes        DO I=1,nplanes
1003           IF(XGOOD(I).EQ.YGOOD(I)) SINGLET = .false.  c$$$         ipass = 0 ! TEST
1004           IF(XGOOD(I).NE.YGOOD(I)) SINGLET = .true.  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1005           ZNEXT = ZM(I)  cPPP         step=vout(3)-zm(i)
1006           IF(SINGLET) THEN  cPP         step=(zm(i)-vout(3))/VOUT(6)
             IF(ZM_A(I).GT.ZM_B(I)+TOLL) THEN  
                ZNEXT = ZM_A(I)  
                ZNEXT2 = ZM_B(I)  
                SINGLET_FIRST = .true.  
                ZDEGENER = .false.  
             ELSEIF(ZM_B(I).GT.ZM_A(I)+TOLL) THEN  
                ZNEXT = ZM_B(I)  
                ZNEXT2 = ZM_A(I)  
                SINGLET_FIRST = .true.  
                ZDEGENER = .false.  
             ELSE  
                ZNEXT = 0.5*(ZM_A(I)+ZM_B(I))  
                SINGLET_FIRST = .false.  
                ZDEGENER = .true.  
             ENDIF  
          ENDIF  
 c$$$         IF(SINGLET) PRINT*,'SINGLET!!!'  
1007   10      DO J=1,7   10      DO J=1,7
1008              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1009              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1010           ENDDO           ENDDO
1011    cPPP         step=vect(3)-zm(i)
1012           IF(VOUT(6).GE.0.) THEN           IF(VOUT(6).GE.0.) THEN
1013              IFAIL=1              IFAIL=1
1014              if(TRKVERBOSE)                    if(TRKVERBOSE)      
1015       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1016              RETURN              RETURN
1017           ENDIF           ENDIF
1018  cPP         step=(zm(i)-vect(3))/VOUT(6)           step=(zm(i)-vect(3))/VOUT(6)
          step=(ZNEXT-vect(3))/VOUT(6)  
1019   11      continue   11      continue
1020           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1021  c$$$         ipass = ipass + 1 ! TEST  c$$$         ipass = ipass + 1 ! TEST
# Line 1105  c$$$            if(.TRUE.)print*,'step', Line 1032  c$$$            if(.TRUE.)print*,'step',
1032              if(TRKVERBOSE)print*,'vect',vect              if(TRKVERBOSE)print*,'vect',vect
1033              if(TRKVERBOSE)print*,'vout',vout              if(TRKVERBOSE)print*,'vout',vout
1034              if(TRKVERBOSE)print*,'step',step              if(TRKVERBOSE)print*,'step',step
             if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1  
1035              RETURN              RETURN
1036           ENDIF           ENDIF
1037           Z=VOUT(3)           Z=VOUT(3)
1038           IF(Z.LE.ZNEXT+TOLL.AND.Z.GE.ZNEXT-TOLL) GOTO 100           IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100
1039           IF(Z.GT.ZNEXT+TOLL) GOTO 10           IF(Z.GT.ZM(I)+TOLL) GOTO 10      
1040           IF(Z.LE.ZNEXT-TOLL) THEN           IF(Z.LE.ZM(I)-TOLL) THEN
1041              STEP=STEP*(ZNEXT-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
1042              DO J=1,7              DO J=1,7
1043                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
1044              ENDDO              ENDDO
1045              GOTO 11              GOTO 11
1046           ENDIF           ENDIF
 c$$$         IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100  
 c$$$         IF(Z.GT.ZM(I)+TOLL) GOTO 10        
 c$$$         IF(Z.LE.ZM(I)-TOLL) THEN  
 c$$$            STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))  
 c$$$            DO J=1,7  
 c$$$               VECT(J)=VECTINI(J)  
 c$$$            ENDDO  
 c$$$            GOTO 11  
 c$$$         ENDIF  
1047    
1048    
1049  *     -----------------------------------------------  *     -----------------------------------------------
1050  *        evaluate track coordinates  *        evaluate track coordinates
1051     100     XV(I)=VOUT(1)
1052   100     IF(SINGLET.AND.(.NOT.ZDEGENER).AND.SINGLET_FIRST) THEN           YV(I)=VOUT(2)
1053              IF(ZM_A(I).GT.ZM_B(I)) THEN           ZV(I)=VOUT(3)
1054                 XV_A(I) = VOUT(1)           AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1055                 YV_A(I) = VOUT(2)           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
                ZV_A(I) = VOUT(3)  
             ELSE  
                XV_B(I) = VOUT(1)  
                YV_B(I) = VOUT(2)  
                ZV_B(I) = VOUT(3)  
             ENDIF  
             AXVUP = DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)  
             AYVUP = DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)  
             ZNEXT = ZNEXT2  
             SINGLET_FIRST = .false.  
             GOTO 10  
          ENDIF  
   
          IF(SINGLET.AND.(.NOT.ZDEGENER).AND.(.NOT.SINGLET_FIRST)) THEN  
             IF(ZM_A(I).LT.ZM_B(I)) THEN  
                XV_A(I) = VOUT(1)  
                YV_A(I) = VOUT(2)  
                ZV_A(I) = VOUT(3)  
             ELSE  
                XV_B(I) = VOUT(1)  
                YV_B(I) = VOUT(2)  
                ZV_B(I) = VOUT(3)  
             ENDIF  
             AXV(I)=0.5*(DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)+AXVUP)  
             AYV(I)=0.5*(DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)+AYVUP)  
          ENDIF  
   
          IF(SINGLET.AND.ZDEGENER) THEN  
             XV_A(I) = VOUT(1)  
             YV_A(I) = VOUT(2)  
             ZV_A(I) = VOUT(3)+0.1  
             XV_B(I) = VOUT(1)  
             YV_B(I) = VOUT(2)  
             ZV_B(I) = VOUT(3)-0.1  
             AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)  
             AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)  
          ENDIF  
   
          IF(.NOT.SINGLET) THEN  
             XV(I)=VOUT(1)  
             YV(I)=VOUT(2)  
             ZV(I)=VOUT(3)  
             AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)  
             AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)  
          ENDIF  
   
1056  *     -----------------------------------------------  *     -----------------------------------------------
1057            
1058           IF(TRACKMODE.EQ.1) THEN           IF(TRACKMODE.EQ.1) THEN
1059  *     -----------------------------------------------  *     -----------------------------------------------
1060  *        change of energy by bremsstrahlung for electrons  *        change of energy by bremsstrahlung for electrons
# Line 1381  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1252  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1252        AL(2) = Y0        AL(2) = Y0
1253        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1254        AL(3) = tath/sqrt(1+tath**2)        AL(3) = tath/sqrt(1+tath**2)
 c$$$      IF(AX.NE.0)THEN  
 c$$$         AL(4)= atan(AY/AX)  
 c$$$      ELSE  
 c$$$         AL(4) = acos(-1.)/2  
 c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  
 c$$$  ENDIF  
 c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)  
 c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.  
   
 c$$$      AL(4) =  0.  
 c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN  
 c$$$         AL(4)= atan(AY/AX)  
 c$$$      ELSEIF(AY.EQ.0)THEN  
 c$$$         AL(4) = 0.  
 c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)          
 c$$$      ELSEIF(AX.EQ.0)THEN  
 c$$$         AL(4) = acos(-1.)/2  
 c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  
 c$$$      ENDIF  
 c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)  
 c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.  
   
 c$$$      AL(4)=0.  
 c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN  
 c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))  
 c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)  
 c$$$      ENDIF  
1255    
1256        AL(4)=0.        AL(4)=0.
1257        IF( AX.NE.0.OR.AY.NE.0. ) THEN        IF( AX.NE.0.OR.AY.NE.0. ) THEN
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23