/[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.21 by pam-fi, Fri Apr 11 13:44:38 2008 UTC revision 1.22 by bongi, Thu Nov 20 15:06:27 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           IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl        DO I=1,nplanes                                          
635              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))           IF( XGOOD(I).NE.YGOOD(I) ) THEN ! singlet
636              ALFA = XM_A(I) - BETA * YM_A(I)              IF(XGOOD(I).EQ.1) THEN
637              YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)                 Z =
638              if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))       $         ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) -
639       $           YM(I)=dmin1(YM_A(I),YM_B(I))       $           (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/
640              if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -
641       $           YM(I)=dmax1(YM_A(I),YM_B(I))       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )
642              XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates                 ZM(I) = Z
643           ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl                 ZV(I) = Z
644              BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))                 XV(I) = XV_A(I)+(XV_B(I)-XV_A(I))*
645              ALFA = YM_A(I) - BETA * XM_A(I)       $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))
646              XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)                 Y =
647              if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))       $         ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(YV_A(I)-YV_B(I)) -
648       $           XM(I)=dmin1(XM_A(I),XM_B(I))       $           (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(YM_A(I)-YM_B(I)) )/
649              if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -
650       $           XM(I)=dmax1(XM_A(I),XM_B(I))       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )
651              YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates                 YM(I) = Y
652           ENDIF                 YV(I) = Y
653           CHI2=CHI2                 XM(I) = XM_A(I)+(XM_B(I)-XM_A(I))*
654       +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )       $                 (Y-YM_A(I))/(YM_B(I)-YM_A(I))
655       +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )  
656       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2                 CHI2=CHI2+(XV(I)-XM(I))**2/RESX(I)**2
657       +                                       *( XGOOD(I)*(1-YGOOD(I)) )  
658       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2              ENDIF
659       +                                       *( (1-XGOOD(I))*YGOOD(I) )              IF(YGOOD(I).EQ.1) THEN
660  c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )                 Z =
661  c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )       $         ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) -
662  c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2       $           (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/
663  c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )       $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -
664  c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2       $           (ZV_A(I)-ZV_B(I))*(XM_A(I)-XM_B(I)) )
665  c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )                 ZM(I) = Z
666  c$$$         print*,XV(I),XM(I),XGOOD(I)                 ZV(I) = Z
667  c$$$         print*,YV(I),YM(I),YGOOD(I)                 YV(I) = YV_A(I)+(YV_B(I)-YV_A(I))*
668         $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))
669                   X =
670         $         ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(XV_A(I)-XV_B(I)) -
671         $           (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(XM_A(I)-XM_B(I)) )/
672         $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -
673         $           (ZV_A(I)-ZV_B(I))*(XM_A(I)-XM_B(I)) )
674                   XM(I) = X
675                   XV(I) = X
676                   YM(I) = YM_A(I)+(YM_B(I)-YM_A(I))*
677         $                 (X-XM_A(I))/(XM_B(I)-XM_A(I))
678    
679                   CHI2=CHI2+(YV(I)-YM(I))**2/RESY(I)**2
680    
681                ENDIF
682             ELSEIF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.1)THEN !Y-cl
683                CHI2=CHI2
684         +           +(XV(I)-XM(I))**2/RESX(i)**2
685         +           +(YV(I)-YM(I))**2/RESY(i)**2
686             ENDIF
687          ENDDO
688          DO I=1,nplanes
689             XV0(I)=XV(I)
690             YV0(I)=YV(I)
691        ENDDO        ENDDO
692    
693    c$$$      DO I=1,nplanes        
694    c$$$         IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
695    c$$$            BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
696    c$$$            ALFA = XM_A(I) - BETA * YM_A(I)
697    c$$$            YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
698    c$$$            if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
699    c$$$     $           YM(I)=dmin1(YM_A(I),YM_B(I))
700    c$$$            if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
701    c$$$     $           YM(I)=dmax1(YM_A(I),YM_B(I))
702    c$$$            XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
703    c$$$         ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
704    c$$$            BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
705    c$$$            ALFA = YM_A(I) - BETA * XM_A(I)
706    c$$$            XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
707    c$$$            if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
708    c$$$     $           XM(I)=dmin1(XM_A(I),XM_B(I))
709    c$$$            if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
710    c$$$     $           XM(I)=dmax1(XM_A(I),XM_B(I))
711    c$$$            YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
712    c$$$         ENDIF
713    c$$$         CHI2=CHI2
714    c$$$     +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
715    c$$$     +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
716    c$$$     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
717    c$$$     +                                       *( XGOOD(I)*(1-YGOOD(I)) )
718    c$$$     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
719    c$$$     +                                       *( (1-XGOOD(I))*YGOOD(I) )
720    c$$$      ENDDO
721    
722  c$$$      print*,'CHISQ ',chi2  c$$$      print*,'CHISQ ',chi2
723  *     ------------------------------------------------  *     ------------------------------------------------
724  *      *    
# Line 966  c$$$     $              +DOWN2*U(I)*U(J) Line 1019  c$$$     $              +DOWN2*U(I)*U(J)
1019  *      *    
1020  *****************************************************************  *****************************************************************
1021    
1022    cPPP --- new --- (with singlets in 3D)
1023        SUBROUTINE POSXYZ(AL_P,IFAIL)        SUBROUTINE POSXYZ(AL_P,IFAIL)
1024    
1025        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# Line 979  c      COMMON/TRKD/TRKVERBOSE Line 1033  c      COMMON/TRKD/TRKVERBOSE
1033        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
1034  c        c      
1035        DIMENSION AL_P(5)        DIMENSION AL_P(5)
1036          LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER
1037  *      *    
1038  cpp      DO I=1,nplanes  cpp      DO I=1,nplanes
1039  cpp         ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
# Line 1000  cpp      ENDDO       Line 1055  cpp      ENDDO      
1055  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1056    
1057        DO I=1,nplanes        DO I=1,nplanes
1058  c$$$         ipass = 0 ! TEST           IF(XGOOD(I).EQ.YGOOD(I)) SINGLET = .false.
1059  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST           IF(XGOOD(I).NE.YGOOD(I)) SINGLET = .true.
1060  cPPP         step=vout(3)-zm(i)           ZNEXT = ZM(I)
1061  cPP         step=(zm(i)-vout(3))/VOUT(6)           IF(SINGLET) THEN
1062                IF(ZM_A(I).GT.ZM_B(I)+TOLL) THEN
1063                   ZNEXT = ZM_A(I)
1064                   ZNEXT2 = ZM_B(I)
1065                   SINGLET_FIRST = .true.
1066                   ZDEGENER = .false.
1067                ELSEIF(ZM_B(I).GT.ZM_A(I)+TOLL) THEN
1068                   ZNEXT = ZM_B(I)
1069                   ZNEXT2 = ZM_A(I)
1070                   SINGLET_FIRST = .true.
1071                   ZDEGENER = .false.
1072                ELSE
1073                   ZNEXT = 0.5*(ZM_A(I)+ZM_B(I))
1074                   SINGLET_FIRST = .false.
1075                   ZDEGENER = .true.
1076                ENDIF
1077             ENDIF
1078    c$$$         IF(SINGLET) PRINT*,'SINGLET!!!'
1079   10      DO J=1,7   10      DO J=1,7
1080              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1081              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1082           ENDDO           ENDDO
1083  cPPP         step=vect(3)-zm(i)  
1084           IF(VOUT(6).GE.0.) THEN           IF(VOUT(6).GE.0.) THEN
1085              IFAIL=1              IFAIL=1
1086              if(TRKVERBOSE)                    if(TRKVERBOSE)      
1087       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1088              RETURN              RETURN
1089           ENDIF           ENDIF
1090           step=(zm(i)-vect(3))/VOUT(6)  cPP         step=(zm(i)-vect(3))/VOUT(6)
1091             step=(ZNEXT-vect(3))/VOUT(6)
1092   11      continue   11      continue
1093           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1094  c$$$         ipass = ipass + 1 ! TEST  c$$$         ipass = ipass + 1 ! TEST
# Line 1032  c$$$            if(.TRUE.)print*,'step', Line 1105  c$$$            if(.TRUE.)print*,'step',
1105              if(TRKVERBOSE)print*,'vect',vect              if(TRKVERBOSE)print*,'vect',vect
1106              if(TRKVERBOSE)print*,'vout',vout              if(TRKVERBOSE)print*,'vout',vout
1107              if(TRKVERBOSE)print*,'step',step              if(TRKVERBOSE)print*,'step',step
1108                if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
1109              RETURN              RETURN
1110           ENDIF           ENDIF
1111           Z=VOUT(3)           Z=VOUT(3)
1112           IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100           IF(Z.LE.ZNEXT+TOLL.AND.Z.GE.ZNEXT-TOLL) GOTO 100
1113           IF(Z.GT.ZM(I)+TOLL) GOTO 10                 IF(Z.GT.ZNEXT+TOLL) GOTO 10
1114           IF(Z.LE.ZM(I)-TOLL) THEN           IF(Z.LE.ZNEXT-TOLL) THEN
1115              STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZNEXT-VECT(3))/(Z-VECT(3))
1116              DO J=1,7              DO J=1,7
1117                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
1118              ENDDO              ENDDO
1119              GOTO 11              GOTO 11
1120           ENDIF           ENDIF
1121    c$$$         IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100
1122    c$$$         IF(Z.GT.ZM(I)+TOLL) GOTO 10      
1123    c$$$         IF(Z.LE.ZM(I)-TOLL) THEN
1124    c$$$            STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
1125    c$$$            DO J=1,7
1126    c$$$               VECT(J)=VECTINI(J)
1127    c$$$            ENDDO
1128    c$$$            GOTO 11
1129    c$$$         ENDIF
1130    
1131    
1132  *     -----------------------------------------------  *     -----------------------------------------------
1133  *        evaluate track coordinates  *        evaluate track coordinates
  100     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.)  
 *     -----------------------------------------------  
1134    
1135     100     IF(SINGLET.AND.(.NOT.ZDEGENER).AND.SINGLET_FIRST) THEN
1136                IF(ZM_A(I).GT.ZM_B(I)) THEN
1137                   XV_A(I) = VOUT(1)
1138                   YV_A(I) = VOUT(2)
1139                   ZV_A(I) = VOUT(3)
1140                ELSE
1141                   XV_B(I) = VOUT(1)
1142                   YV_B(I) = VOUT(2)
1143                   ZV_B(I) = VOUT(3)
1144                ENDIF
1145                AXVUP = DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1146                AYVUP = DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1147                ZNEXT = ZNEXT2
1148                SINGLET_FIRST = .false.
1149                GOTO 10
1150             ENDIF
1151    
1152             IF(SINGLET.AND.(.NOT.ZDEGENER).AND.(.NOT.SINGLET_FIRST)) THEN
1153                IF(ZM_A(I).LT.ZM_B(I)) THEN
1154                   XV_A(I) = VOUT(1)
1155                   YV_A(I) = VOUT(2)
1156                   ZV_A(I) = VOUT(3)
1157                ELSE
1158                   XV_B(I) = VOUT(1)
1159                   YV_B(I) = VOUT(2)
1160                   ZV_B(I) = VOUT(3)
1161                ENDIF
1162                AXV(I)=0.5*(DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)+AXVUP)
1163                AYV(I)=0.5*(DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)+AYVUP)
1164             ENDIF
1165    
1166             IF(SINGLET.AND.ZDEGENER) THEN
1167                XV_A(I) = VOUT(1)
1168                YV_A(I) = VOUT(2)
1169                ZV_A(I) = VOUT(3)+0.1
1170                XV_B(I) = VOUT(1)
1171                YV_B(I) = VOUT(2)
1172                ZV_B(I) = VOUT(3)-0.1
1173                AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1174                AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1175             ENDIF
1176    
1177             IF(.NOT.SINGLET) THEN
1178                XV(I)=VOUT(1)
1179                YV(I)=VOUT(2)
1180                ZV(I)=VOUT(3)
1181                AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1182                AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1183             ENDIF
1184    
1185    *     -----------------------------------------------
1186            
1187           IF(TRACKMODE.EQ.1) THEN           IF(TRACKMODE.EQ.1) THEN
1188  *     -----------------------------------------------  *     -----------------------------------------------
1189  *        change of energy by bremsstrahlung for electrons  *        change of energy by bremsstrahlung for electrons
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23