--- DarthVader/TrackerLevel2/src/F77/mini.f 2008/11/20 15:06:27 1.22 +++ DarthVader/TrackerLevel2/src/F77/mini.f 2008/12/05 08:26:47 1.23 @@ -630,95 +630,42 @@ * measured position of the cluster. * --------------------------------------------------------- CHI2=0. - - DO I=1,nplanes - IF( XGOOD(I).NE.YGOOD(I) ) THEN ! singlet - IF(XGOOD(I).EQ.1) THEN - Z = - $ ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) - - $ (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/ - $ ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) - - $ (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) ) - ZM(I) = Z - ZV(I) = Z - XV(I) = XV_A(I)+(XV_B(I)-XV_A(I))* - $ (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I)) - Y = - $ ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(YV_A(I)-YV_B(I)) - - $ (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(YM_A(I)-YM_B(I)) )/ - $ ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) - - $ (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) ) - YM(I) = Y - YV(I) = Y - XM(I) = XM_A(I)+(XM_B(I)-XM_A(I))* - $ (Y-YM_A(I))/(YM_B(I)-YM_A(I)) - - CHI2=CHI2+(XV(I)-XM(I))**2/RESX(I)**2 - - ENDIF - IF(YGOOD(I).EQ.1) THEN - Z = - $ ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) - - $ (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_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)) ) - ZM(I) = Z - ZV(I) = Z - YV(I) = YV_A(I)+(YV_B(I)-YV_A(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) + DO I=1,nplanes + IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl + BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I)) + ALFA = XM_A(I) - BETA * YM_A(I) + YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2) + if(YM(I).lt.dmin1(YM_A(I),YM_B(I))) + $ YM(I)=dmin1(YM_A(I),YM_B(I)) + if(YM(I).gt.dmax1(YM_A(I),YM_B(I))) + $ YM(I)=dmax1(YM_A(I),YM_B(I)) + XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates + ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl + BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I)) + ALFA = YM_A(I) - BETA * XM_A(I) + XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2) + if(XM(I).lt.dmin1(XM_A(I),XM_B(I))) + $ XM(I)=dmin1(XM_A(I),XM_B(I)) + if(XM(I).gt.dmax1(XM_A(I),XM_B(I))) + $ XM(I)=dmax1(XM_A(I),XM_B(I)) + YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates + ENDIF + CHI2=CHI2 + + +(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) ) + + +(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) ) + + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 + + *( XGOOD(I)*(1-YGOOD(I)) ) + + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 + + *( (1-XGOOD(I))*YGOOD(I) ) +c$$$ print*,(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) ) +c$$$ print*,(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) ) +c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 +c$$$ + *( XGOOD(I)*(1-YGOOD(I)) ) +c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 +c$$$ + *( (1-XGOOD(I))*YGOOD(I) ) +c$$$ print*,XV(I),XM(I),XGOOD(I) +c$$$ print*,YV(I),YM(I),YGOOD(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 - c$$$ print*,'CHISQ ',chi2 * ------------------------------------------------ * @@ -1019,7 +966,6 @@ * ***************************************************************** -cPPP --- new --- (with singlets in 3D) SUBROUTINE POSXYZ(AL_P,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) @@ -1033,7 +979,6 @@ COMMON/TRKD/TRKDEBUG,TRKVERBOSE c DIMENSION AL_P(5) - LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER * cpp DO I=1,nplanes cpp ZV(I)=ZM(I) ! @@ -1055,40 +1000,22 @@ c$$$ print*,'POSXY (prima) ',vout DO I=1,nplanes - IF(XGOOD(I).EQ.YGOOD(I)) SINGLET = .false. - IF(XGOOD(I).NE.YGOOD(I)) SINGLET = .true. - ZNEXT = ZM(I) - IF(SINGLET) THEN - 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!!!' +c$$$ ipass = 0 ! TEST +c$$$ PRINT *,'TRACKING -> START PLANE: ',I ! TEST +cPPP step=vout(3)-zm(i) +cPP step=(zm(i)-vout(3))/VOUT(6) 10 DO J=1,7 VECT(J)=VOUT(J) VECTINI(J)=VOUT(J) ENDDO - +cPPP step=vect(3)-zm(i) IF(VOUT(6).GE.0.) THEN IFAIL=1 if(TRKVERBOSE) $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' RETURN ENDIF -cPP step=(zm(i)-vect(3))/VOUT(6) - step=(ZNEXT-vect(3))/VOUT(6) + step=(zm(i)-vect(3))/VOUT(6) 11 continue CALL GRKUTA(CHARGE,STEP,VECT,VOUT) c$$$ ipass = ipass + 1 ! TEST @@ -1105,85 +1032,29 @@ if(TRKVERBOSE)print*,'vect',vect if(TRKVERBOSE)print*,'vout',vout if(TRKVERBOSE)print*,'step',step - if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 RETURN ENDIF Z=VOUT(3) - IF(Z.LE.ZNEXT+TOLL.AND.Z.GE.ZNEXT-TOLL) GOTO 100 - IF(Z.GT.ZNEXT+TOLL) GOTO 10 - IF(Z.LE.ZNEXT-TOLL) THEN - STEP=STEP*(ZNEXT-VECT(3))/(Z-VECT(3)) + IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100 + IF(Z.GT.ZM(I)+TOLL) GOTO 10 + IF(Z.LE.ZM(I)-TOLL) THEN + STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3)) DO J=1,7 VECT(J)=VECTINI(J) ENDDO GOTO 11 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 * ----------------------------------------------- * evaluate track coordinates - - 100 IF(SINGLET.AND.(.NOT.ZDEGENER).AND.SINGLET_FIRST) THEN - IF(ZM_A(I).GT.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 - 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 - + 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.) * ----------------------------------------------- - + IF(TRACKMODE.EQ.1) THEN * ----------------------------------------------- * change of energy by bremsstrahlung for electrons @@ -1381,33 +1252,6 @@ AL(2) = Y0 tath = sqrt(AY**2+AX**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 AL(4)=0. IF( AX.NE.0.OR.AY.NE.0. ) THEN