--- DarthVader/TrackerLevel2/src/F77/mini.f 2008/04/11 13:44:38 1.21 +++ DarthVader/TrackerLevel2/src/F77/mini.f 2008/11/20 15:06:27 1.22 @@ -630,42 +630,95 @@ * measured position of the cluster. * --------------------------------------------------------- CHI2=0. - 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) + + 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) 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 * ------------------------------------------------ * @@ -966,6 +1019,7 @@ * ***************************************************************** +cPPP --- new --- (with singlets in 3D) SUBROUTINE POSXYZ(AL_P,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) @@ -979,6 +1033,7 @@ COMMON/TRKD/TRKDEBUG,TRKVERBOSE c DIMENSION AL_P(5) + LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER * cpp DO I=1,nplanes cpp ZV(I)=ZM(I) ! @@ -1000,22 +1055,40 @@ c$$$ print*,'POSXY (prima) ',vout DO I=1,nplanes -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) + 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!!!' 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 - step=(zm(i)-vect(3))/VOUT(6) +cPP step=(zm(i)-vect(3))/VOUT(6) + step=(ZNEXT-vect(3))/VOUT(6) 11 continue CALL GRKUTA(CHARGE,STEP,VECT,VOUT) c$$$ ipass = ipass + 1 ! TEST @@ -1032,29 +1105,85 @@ 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.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)) + 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)) 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 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.) -* ----------------------------------------------- + 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 + +* ----------------------------------------------- + IF(TRACKMODE.EQ.1) THEN * ----------------------------------------------- * change of energy by bremsstrahlung for electrons