/[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.26 by mocchiut, Thu Jan 16 15:29:56 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    
28          double precision NX, NY ! EM GCC4.7
29  c------------------------------------------------------------------------  c------------------------------------------------------------------------
30  c     variables used in the tracking procedure (mini and its subroutines)  c     variables used in the tracking procedure (mini and its subroutines)
31  c  c
# Line 44  c      DATA XGOOD,YGOOD/nplanes*1.,nplan Line 47  c      DATA XGOOD,YGOOD/nplanes*1.,nplan
47  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
48  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
49        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
50        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !"
51    
52  c$$$      DIMENSION DAL(5)                    !increment of vector alfa  c$$$      DIMENSION DAL(5)                    !increment of vector alfa
53        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 141  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 144  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
144        CHI2=0        CHI2=0
145    
146        if(TRKDEBUG) print*,'guess: ',al        if(TRKDEBUG) print*,'guess: ',al
147        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
148    
149  *      *    
150  *     -----------------------  *     -----------------------
# Line 251  c$$$            CHI2D(I)=CHI2D(I)*COST Line 254  c$$$            CHI2D(I)=CHI2D(I)*COST
254              ENDDO              ENDDO
255           ENDIF           ENDIF
256                
257           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)           if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5)
258    
259  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)  c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
260  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)  c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
# Line 418  c$$$     $        ISTEPMAX Line 421  c$$$     $        ISTEPMAX
421  *     ---------------------------------------------  *     ---------------------------------------------
422  *------------------------------------------------------------*  *------------------------------------------------------------*
423  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
424          IF(FACT.EQ.0)THEN
425             IFAIL=1
426             RETURN
427          ENDIF
428        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT        ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
429        ALTOL(1) = ALTOL(5)/DELETA1        ALTOL(1) = ALTOL(5)/DELETA1
430        ALTOL(2) = ALTOL(1)        ALTOL(2) = ALTOL(1)
# Line 554  c$$$               DAL(I)=0. Line 561  c$$$               DAL(I)=0.
561  *     ------------------------------------  *     ------------------------------------
562  *     Reduced chi^2  *     Reduced chi^2
563        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
   
564  c      print*,'mini2: chi2 ',chi2  c      print*,'mini2: chi2 ',chi2
565    
566   11   CONTINUE         11   CONTINUE      
567    
568        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,AL(5)
569    
570        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
571    
# Line 630  c$$$      ENDDO Line 636  c$$$      ENDDO
636  *     measured position of the cluster.  *     measured position of the cluster.
637  *     ---------------------------------------------------------  *     ---------------------------------------------------------
638        CHI2=0.        CHI2=0.
639          DO I=1,nplanes        
640        DO I=1,nplanes                                                     IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
641           IF( XGOOD(I).NE.YGOOD(I) ) THEN ! singlet              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
642              IF(XGOOD(I).EQ.1) THEN              ALFA = XM_A(I) - BETA * YM_A(I)
643                 Z =              YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
644       $         ( (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)))
645       $           (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))
646       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
647       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )       $           YM(I)=dmax1(YM_A(I),YM_B(I))
648                 ZM(I) = Z              XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
649                 ZV(I) = Z           ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
650                 XV(I) = XV_A(I)+(XV_B(I)-XV_A(I))*              BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
651       $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))              ALFA = YM_A(I) - BETA * XM_A(I)
652                 Y =              XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
653       $         ( (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)))
654       $           (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))
655       $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -              if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
656       $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )       $           XM(I)=dmax1(XM_A(I),XM_B(I))
657                 YM(I) = Y              YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
658                 YV(I) = Y           ENDIF
659                 XM(I) = XM_A(I)+(XM_B(I)-XM_A(I))*           CHI2=CHI2
660       $                 (Y-YM_A(I))/(YM_B(I)-YM_A(I))       +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
661         +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
662                 CHI2=CHI2+(XV(I)-XM(I))**2/RESX(I)**2       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
663         +                                       *( XGOOD(I)*(1-YGOOD(I)) )
664              ENDIF       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
665              IF(YGOOD(I).EQ.1) THEN       +                                       *( (1-XGOOD(I))*YGOOD(I) )
666                 Z =  c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
667       $         ( (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) )
668       $           (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
669       $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -  c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
670       $           (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
671                 ZM(I) = Z  c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
672                 ZV(I) = Z  c$$$         print*,XV(I),XM(I),XGOOD(I)
673                 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  
674        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  
   
675  c$$$      print*,'CHISQ ',chi2  c$$$      print*,'CHISQ ',chi2
676  *     ------------------------------------------------  *     ------------------------------------------------
677  *      *    
# Line 1019  c$$$     $              +DOWN2*U(I)*U(J) Line 972  c$$$     $              +DOWN2*U(I)*U(J)
972  *      *    
973  *****************************************************************  *****************************************************************
974    
 cPPP --- new --- (with singlets in 3D)  
975        SUBROUTINE POSXYZ(AL_P,IFAIL)        SUBROUTINE POSXYZ(AL_P,IFAIL)
976    
977        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# Line 1033  c      COMMON/TRKD/TRKVERBOSE Line 985  c      COMMON/TRKD/TRKVERBOSE
985        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
986  c        c      
987        DIMENSION AL_P(5)        DIMENSION AL_P(5)
       LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER  
988  *      *    
989  cpp      DO I=1,nplanes  cpp      DO I=1,nplanes
990  cpp         ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
# Line 1055  cpp      ENDDO       Line 1006  cpp      ENDDO      
1006  c$$$      print*,'POSXY (prima) ',vout  c$$$      print*,'POSXY (prima) ',vout
1007    
1008        DO I=1,nplanes        DO I=1,nplanes
1009           IF(XGOOD(I).EQ.YGOOD(I)) SINGLET = .false.  c$$$         ipass = 0 ! TEST
1010           IF(XGOOD(I).NE.YGOOD(I)) SINGLET = .true.  c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1011           ZNEXT = ZM(I)  cPPP         step=vout(3)-zm(i)
1012           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!!!'  
1013   10      DO J=1,7   10      DO J=1,7
1014              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1015              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1016           ENDDO           ENDDO
1017    cPPP         step=vect(3)-zm(i)
1018           IF(VOUT(6).GE.0.) THEN           IF(VOUT(6).GE.0.) THEN
1019              IFAIL=1              IFAIL=1
1020              if(TRKVERBOSE)                    if(TRKVERBOSE)      
1021       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1022              RETURN              RETURN
1023           ENDIF           ENDIF
1024  cPP         step=(zm(i)-vect(3))/VOUT(6)           step=(zm(i)-vect(3))/VOUT(6)
          step=(ZNEXT-vect(3))/VOUT(6)  
1025   11      continue   11      continue
1026           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1027  c$$$         ipass = ipass + 1 ! TEST  c$$$         ipass = ipass + 1 ! TEST
# Line 1105  c$$$            if(.TRUE.)print*,'step', Line 1038  c$$$            if(.TRUE.)print*,'step',
1038              if(TRKVERBOSE)print*,'vect',vect              if(TRKVERBOSE)print*,'vect',vect
1039              if(TRKVERBOSE)print*,'vout',vout              if(TRKVERBOSE)print*,'vout',vout
1040              if(TRKVERBOSE)print*,'step',step              if(TRKVERBOSE)print*,'step',step
             if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1  
1041              RETURN              RETURN
1042           ENDIF           ENDIF
1043           Z=VOUT(3)           Z=VOUT(3)
1044           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
1045           IF(Z.GT.ZNEXT+TOLL) GOTO 10           IF(Z.GT.ZM(I)+TOLL) GOTO 10      
1046           IF(Z.LE.ZNEXT-TOLL) THEN           IF(Z.LE.ZM(I)-TOLL) THEN
1047              STEP=STEP*(ZNEXT-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
1048              DO J=1,7              DO J=1,7
1049                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
1050              ENDDO              ENDDO
1051              GOTO 11              GOTO 11
1052           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  
1053    
1054    
1055  *     -----------------------------------------------  *     -----------------------------------------------
1056  *        evaluate track coordinates  *        evaluate track coordinates
1057     100     XV(I)=VOUT(1)
1058   100     IF(SINGLET.AND.(.NOT.ZDEGENER).AND.SINGLET_FIRST) THEN           YV(I)=VOUT(2)
1059              IF(ZM_A(I).GT.ZM_B(I)) THEN           ZV(I)=VOUT(3)
1060                 XV_A(I) = VOUT(1)           AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1061                 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  
   
1062  *     -----------------------------------------------  *     -----------------------------------------------
1063            
1064           IF(TRACKMODE.EQ.1) THEN           IF(TRACKMODE.EQ.1) THEN
1065  *     -----------------------------------------------  *     -----------------------------------------------
1066  *        change of energy by bremsstrahlung for electrons  *        change of energy by bremsstrahlung for electrons
# Line 1260  c$$$      print*,'POSXY (dopo) ',vout Line 1137  c$$$      print*,'POSXY (dopo) ',vout
1137    
1138        subroutine guess()        subroutine guess()
1139    
1140  c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z) ! EM GCC4.7
1141    
1142        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
1143        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
1144                
1145        REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)        REAL*8 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) ! EM GCC4.7
1146        REAL*4 CHI,XC,ZC,RADIUS        REAL*4 CHI,XC,ZC,RADIUS
1147  *     ----------------------------------------  *     ----------------------------------------
1148  *     Y view  *     Y view
# Line 1381  c$$$     $     ,' - CHI ',CHI,' - X0,AX, Line 1258  c$$$     $     ,' - CHI ',CHI,' - X0,AX,
1258        AL(2) = Y0        AL(2) = Y0
1259        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1260        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  
1261    
1262        AL(4)=0.        AL(4)=0.
1263        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.26
  ViewVC Help
Powered by ViewVC 1.1.23