/[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.27 by mocchiut, Fri Jan 17 12:56:51 2014 UTC
# Line 23  c      logical DEBUG Line 23  c      logical DEBUG
23  c      common/dbg/DEBUG  c      common/dbg/DEBUG
24    
25        parameter (dinf=1.d15)      !just a huge number...        parameter (dinf=1.d15)      !just a huge number...
26          parameter (dinfneg=-dinf)   ! just a huge negative number...
27  c------------------------------------------------------------------------  c------------------------------------------------------------------------
28  c     variables used in the tracking procedure (mini and its subroutines)  c     variables used in the tracking procedure (mini and its subroutines)
29  c  c
# Line 44  c      DATA XGOOD,YGOOD/nplanes*1.,nplan Line 45  c      DATA XGOOD,YGOOD/nplanes*1.,nplan
45  c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components  c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
46  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
47        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
48        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !"
49    
50  c$$$      DIMENSION DAL(5)                    !increment of vector alfa  c$$$      DIMENSION DAL(5)                    !increment of vector alfa
51        DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2        DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2
# Line 101  c$$$      PRINT*,'==========' ! TEST Line 102  c$$$      PRINT*,'==========' ! TEST
102  *     ----------------------------------------------------------  *     ----------------------------------------------------------
103        AVRESX = RESXAV        AVRESX = RESXAV
104        AVRESY = RESYAV        AVRESY = RESYAV
105        NX = 0.0        NX = 0 !EM GCC4.7
106        NY = 0.0        NY = 0 !EM GCC4.7
107        DO IP=1,6        DO IP=1,6
108           IF( XGOOD(IP).EQ.1 )THEN           IF( XGOOD(IP).EQ.1 )THEN
109              NX=NX+1.0              NX=NX+1!EM GCC4.7
110              AVRESX=AVRESX+RESX(IP)              AVRESX=AVRESX+RESX(IP)
111           ENDIF           ENDIF
112           IF( YGOOD(IP).EQ.1 )THEN           IF( YGOOD(IP).EQ.1 )THEN
113              NY=NY+1.0              NY=NY+1!EM GCC4.7
114              AVRESY=AVRESY+RESY(IP)              AVRESY=AVRESY+RESY(IP)
115           ENDIF           ENDIF
116        ENDDO        ENDDO
# Line 141  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 142  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
142        CHI2=0        CHI2=0
143    
144        if(TRKDEBUG) print*,'guess: ',al        if(TRKDEBUG) print*,'guess: ',al
145        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
146    
147  *      *    
148  *     -----------------------  *     -----------------------
# Line 251  c$$$            CHI2D(I)=CHI2D(I)*COST Line 252  c$$$            CHI2D(I)=CHI2D(I)*COST
252              ENDDO              ENDDO
253           ENDIF           ENDIF
254                
255           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
256    
257  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
258  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
# Line 418  c$$$     $        ISTEPMAX Line 419  c$$$     $        ISTEPMAX
419  *     ---------------------------------------------  *     ---------------------------------------------
420  *------------------------------------------------------------*  *------------------------------------------------------------*
421  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
422          IF(FACT.EQ.0)THEN
423             IFAIL=1
424             RETURN
425          ENDIF
426        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
427        ALTOL(1) = ALTOL(5)/DELETA1        ALTOL(1) = ALTOL(5)/DELETA1
428        ALTOL(2) = ALTOL(1)        ALTOL(2) = ALTOL(1)
# Line 554  c$$$               DAL(I)=0. Line 559  c$$$               DAL(I)=0.
559  *     ------------------------------------  *     ------------------------------------
560  *     Reduced chi^2  *     Reduced chi^2
561        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
   
562  c      print*,'mini2: chi2 ',chi2  c      print*,'mini2: chi2 ',chi2
563    
564   11   CONTINUE         11   CONTINUE      
565    
566        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,AL(5)
567    
568        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
569    
# Line 630  c$$$      ENDDO Line 634  c$$$      ENDDO
634  *     measured position of the cluster.  *     measured position of the cluster.
635  *     ---------------------------------------------------------  *     ---------------------------------------------------------
636        CHI2=0.        CHI2=0.
637          DO I=1,nplanes        
638        DO I=1,nplanes                                                     IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
639           IF( XGOOD(I).NE.YGOOD(I) ) THEN ! singlet              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
640              IF(XGOOD(I).EQ.1) THEN              ALFA = XM_A(I) - BETA * YM_A(I)
641                 Z =              YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
642       $         ( (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)))
643       $           (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))
644       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
645       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )       $           YM(I)=dmax1(YM_A(I),YM_B(I))
646                 ZM(I) = Z              XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
647                 ZV(I) = Z           ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
648                 XV(I) = XV_A(I)+(XV_B(I)-XV_A(I))*              BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
649       $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))              ALFA = YM_A(I) - BETA * XM_A(I)
650                 Y =              XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
651       $         ( (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)))
652       $           (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))
653       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
654       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )       $           XM(I)=dmax1(XM_A(I),XM_B(I))
655                 YM(I) = Y              YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
656                 YV(I) = Y           ENDIF
657                 XM(I) = XM_A(I)+(XM_B(I)-XM_A(I))*           CHI2=CHI2
658       $                 (Y-YM_A(I))/(YM_B(I)-YM_A(I))       +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
659         +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
660                 CHI2=CHI2+(XV(I)-XM(I))**2/RESX(I)**2       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
661         +                                       *( XGOOD(I)*(1-YGOOD(I)) )
662              ENDIF       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
663              IF(YGOOD(I).EQ.1) THEN       +                                       *( (1-XGOOD(I))*YGOOD(I) )
664                 Z =  c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
665       $         ( (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) )
666       $           (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
667       $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -  c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
668       $           (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
669                 ZM(I) = Z  c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
670                 ZV(I) = Z  c$$$         print*,XV(I),XM(I),XGOOD(I)
671                 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  
672        ENDDO        ENDDO
       DO I=1,nplanes  
          XV0(I)=XV(I)  
          YV0(I)=YV(I)  
       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  
   
673  c$$$      print*,'CHISQ ',chi2  c$$$      print*,'CHISQ ',chi2
674  *     ------------------------------------------------  *     ------------------------------------------------
675  *      *    
# Line 1019  c$$$     $              +DOWN2*U(I)*U(J) Line 970  c$$$     $              +DOWN2*U(I)*U(J)
970  *      *    
971  *****************************************************************  *****************************************************************
972    
 cPPP --- new --- (with singlets in 3D)  
973        SUBROUTINE POSXYZ(AL_P,IFAIL)        SUBROUTINE POSXYZ(AL_P,IFAIL)
974    
975        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# Line 1033  c      COMMON/TRKD/TRKVERBOSE Line 983  c      COMMON/TRKD/TRKVERBOSE
983        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
984  c        c      
985        DIMENSION AL_P(5)        DIMENSION AL_P(5)
       LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER  
986  *      *    
987  cpp      DO I=1,nplanes  cpp      DO I=1,nplanes
988  cpp         ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
# Line 1055  cpp      ENDDO       Line 1004  cpp      ENDDO      
1004  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1005    
1006        DO I=1,nplanes        DO I=1,nplanes
1007           IF(XGOOD(I).EQ.YGOOD(I)) SINGLET = .false.  c$$$         ipass = 0 ! TEST
1008           IF(XGOOD(I).NE.YGOOD(I)) SINGLET = .true.  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1009           ZNEXT = ZM(I)  cPPP         step=vout(3)-zm(i)
1010           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!!!'  
1011   10      DO J=1,7   10      DO J=1,7
1012              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1013              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1014           ENDDO           ENDDO
1015    cPPP         step=vect(3)-zm(i)
1016           IF(VOUT(6).GE.0.) THEN           IF(VOUT(6).GE.0.) THEN
1017              IFAIL=1              IFAIL=1
1018              if(TRKVERBOSE)                    if(TRKVERBOSE)      
1019       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1020              RETURN              RETURN
1021           ENDIF           ENDIF
1022  cPP         step=(zm(i)-vect(3))/VOUT(6)           step=(zm(i)-vect(3))/VOUT(6)
          step=(ZNEXT-vect(3))/VOUT(6)  
1023   11      continue   11      continue
1024           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1025  c$$$         ipass = ipass + 1 ! TEST  c$$$         ipass = ipass + 1 ! TEST
# Line 1105  c$$$            if(.TRUE.)print*,'step', Line 1036  c$$$            if(.TRUE.)print*,'step',
1036              if(TRKVERBOSE)print*,'vect',vect              if(TRKVERBOSE)print*,'vect',vect
1037              if(TRKVERBOSE)print*,'vout',vout              if(TRKVERBOSE)print*,'vout',vout
1038              if(TRKVERBOSE)print*,'step',step              if(TRKVERBOSE)print*,'step',step
             if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1  
1039              RETURN              RETURN
1040           ENDIF           ENDIF
1041           Z=VOUT(3)           Z=VOUT(3)
1042           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
1043           IF(Z.GT.ZNEXT+TOLL) GOTO 10           IF(Z.GT.ZM(I)+TOLL) GOTO 10      
1044           IF(Z.LE.ZNEXT-TOLL) THEN           IF(Z.LE.ZM(I)-TOLL) THEN
1045              STEP=STEP*(ZNEXT-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
1046              DO J=1,7              DO J=1,7
1047                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
1048              ENDDO              ENDDO
1049              GOTO 11              GOTO 11
1050           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  
1051    
1052    
1053  *     -----------------------------------------------  *     -----------------------------------------------
1054  *        evaluate track coordinates  *        evaluate track coordinates
1055     100     XV(I)=VOUT(1)
1056   100     IF(SINGLET.AND.(.NOT.ZDEGENER).AND.SINGLET_FIRST) THEN           YV(I)=VOUT(2)
1057              IF(ZM_A(I).GT.ZM_B(I)) THEN           ZV(I)=VOUT(3)
1058                 XV_A(I) = VOUT(1)           AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1059                 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  
   
1060  *     -----------------------------------------------  *     -----------------------------------------------
1061            
1062           IF(TRACKMODE.EQ.1) THEN           IF(TRACKMODE.EQ.1) THEN
1063  *     -----------------------------------------------  *     -----------------------------------------------
1064  *        change of energy by bremsstrahlung for electrons  *        change of energy by bremsstrahlung for electrons
# Line 1280  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1155  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1155        S1=0.        S1=0.
1156        DO I=1,nplanes        DO I=1,nplanes
1157           IF(YGOOD(I).EQ.1)THEN           IF(YGOOD(I).EQ.1)THEN
1158              YY = YM(I)              YY = REAL(YM(I))!EM GCC4.7
1159              IF(XGOOD(I).EQ.0)THEN              IF(XGOOD(I).EQ.0)THEN
1160                 YY = (YM_A(I) + YM_B(I))/2                 YY = REAL((YM_A(I) + YM_B(I))/2.)!EM GCC4.7
1161              ENDIF              ENDIF
1162              SZZ=SZZ+ZM(I)*ZM(I)              SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7
1163              SZY=SZY+ZM(I)*YY              SZY=SZY+REAL(ZM(I)*YY)!EM GCC4.7
1164              SSY=SSY+YY              SSY=SSY+YY
1165              SZ=SZ+ZM(I)              SZ=SZ+REAL(ZM(I))!EM GCC4.7
1166              S1=S1+1.              S1=S1+1.
1167           ENDIF           ENDIF
1168        ENDDO        ENDDO
1169        DET=SZZ*S1-SZ*SZ        DET=SZZ*S1-SZ*SZ
1170        AY=(SZY*S1-SZ*SSY)/DET        AY=(SZY*S1-SZ*SSY)/DET
1171        BY=(SZZ*SSY-SZY*SZ)/DET        BY=(SZZ*SSY-SZY*SZ)/DET
1172        Y0 = AY*ZINI+BY        Y0 = REAL(AY*ZINI+BY)!EM GCC4.7
1173  *     ----------------------------------------  *     ----------------------------------------
1174  *     X view  *     X view
1175  *     ----------------------------------------  *     ----------------------------------------
# Line 1304  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1179  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1179        NP=0        NP=0
1180        DO I=1,nplanes        DO I=1,nplanes
1181           IF(XGOOD(I).EQ.1)THEN           IF(XGOOD(I).EQ.1)THEN
1182              XX = XM(I)              XX = REAL(XM(I))!EM GCC4.7
1183              IF(YGOOD(I).EQ.0)THEN              IF(YGOOD(I).EQ.0)THEN
1184                 XX = (XM_A(I) + XM_B(I))/2                 XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7
1185              ENDIF              ENDIF
1186              NP=NP+1              NP=NP+1
1187              XP(NP)=XX              XP(NP)=XX
1188              ZP(NP)=ZM(I)              ZP(NP)=REAL(ZM(I))!EM GCC4.7
1189           ENDIF           ENDIF
1190        ENDDO        ENDDO
1191        IFLAG=0                   !no debug mode        IFLAG=0                   !no debug mode
# Line 1324  c$$$      print*,' XP ',(rp(i),i=1,np) Line 1199  c$$$      print*,' XP ',(rp(i),i=1,np)
1199    
1200        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
1201  c      if(CHI.gt.100)GOTO 10 !straigth fit  c      if(CHI.gt.100)GOTO 10 !straigth fit
1202        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = REAL(RADIUS**2-(ZINI-ZC)**2)!EM GCC4.7
1203        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
1204        DC = SQRT(ARG)              DC = SQRT(ARG)      
1205        IF(XC.GT.0)DC=-DC        IF(XC.GT.0)DC=-DC
1206        X0=XC+DC        X0=XC+DC
1207        AX = -(ZINI-ZC)/DC        AX = REAL(-(ZINI-ZC)/DC)!EM GCC4.7
1208        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
1209        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
1210                
# Line 1355  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1230  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1230        S1=0.        S1=0.
1231        DO I=1,nplanes        DO I=1,nplanes
1232           IF(XGOOD(I).EQ.1)THEN           IF(XGOOD(I).EQ.1)THEN
1233              XX = XM(I)              XX = REAL(XM(I))!EM GCC4.7
1234              IF(YGOOD(I).EQ.0)THEN              IF(YGOOD(I).EQ.0)THEN
1235                 XX = (XM_A(I) + XM_B(I))/2                 XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7
1236              ENDIF              ENDIF
1237              SZZ=SZZ+ZM(I)*ZM(I)              SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7
1238              SZX=SZX+ZM(I)*XX              SZX=SZX+REAL(ZM(I)*XX)!EM GCC4.7
1239              SSX=SSX+XX              SSX=SSX+XX
1240              SZ=SZ+ZM(I)              SZ=SZ+REAL(ZM(I))!EM GCC4.7
1241              S1=S1+1.              S1=S1+1.
1242           ENDIF           ENDIF
1243        ENDDO        ENDDO
# Line 1370  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1245  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1245        AX=(SZX*S1-SZ*SSX)/DET        AX=(SZX*S1-SZ*SSX)/DET
1246        BX=(SZZ*SSX-SZX*SZ)/DET        BX=(SZZ*SSX-SZX*SZ)/DET
1247        DEF = 0        DEF = 0
1248        X0  = AX*ZINI+BX        X0  = REAL(AX*ZINI+BX)!EM GCC4.7
1249    
1250   20   CONTINUE   20   CONTINUE
1251  *     ----------------------------------------  *     ----------------------------------------
# Line 1381  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1256  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1256        AL(2) = Y0        AL(2) = Y0
1257        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1258        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  
1259    
1260        AL(4)=0.        AL(4)=0.
1261        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.27
  ViewVC Help
Powered by ViewVC 1.1.23