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 |
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 |
73 |
|
|
74 |
DIMENSION AL0(5) |
DIMENSION AL0(5) |
75 |
LOGICAL SUCCESS_NEW,SUCCESS_OLD |
LOGICAL SUCCESS_NEW,SUCCESS_OLD |
76 |
|
|
77 |
|
c$$$ PRINT*,'==========' ! TEST |
78 |
|
c$$$ PRINT*,'START MINI' ! TEST |
79 |
|
c$$$ PRINT*,'==========' ! TEST |
80 |
|
|
81 |
* |
* |
82 |
* define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student) |
* define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student) |
83 |
* |
* |
102 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
103 |
AVRESX = RESXAV |
AVRESX = RESXAV |
104 |
AVRESY = RESYAV |
AVRESY = RESYAV |
105 |
|
NX = 0 !EM GCC4.7 |
106 |
|
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 |
NX=NX+1!EM GCC4.7 |
110 |
AVRESX=AVRESX+RESX(IP) |
AVRESX=AVRESX+RESX(IP) |
111 |
ENDIF |
ENDIF |
|
IF(NX.NE.0)AVRESX=AVRESX/NX |
|
112 |
IF( YGOOD(IP).EQ.1 )THEN |
IF( YGOOD(IP).EQ.1 )THEN |
113 |
NY=NY+1 |
NY=NY+1!EM GCC4.7 |
114 |
AVRESY=AVRESY+RESY(IP) |
AVRESY=AVRESY+RESY(IP) |
115 |
ENDIF |
ENDIF |
|
IF(NX.NE.0)AVRESY=AVRESY/NY |
|
116 |
ENDDO |
ENDDO |
117 |
|
IF(NX.NE.0.0)AVRESX=AVRESX/NX |
118 |
|
IF(NY.NE.0.0)AVRESY=AVRESY/NY |
119 |
|
|
120 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
121 |
* define ALTOL(5) ---> tolerances on state vector |
* define ALTOL(5) ---> tolerances on state vector |
141 |
JFAIL=0 !error flag |
JFAIL=0 !error flag |
142 |
CHI2=0 |
CHI2=0 |
143 |
|
|
144 |
if(TRKDEBUG) print*,'guess: ',al |
if(TRKDEBUG) print*,'mini : guess ',al |
145 |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
if(TRKDEBUG) print*,'mini : step ',istep,chi2,AL(5) |
146 |
|
|
147 |
* |
* |
148 |
* ----------------------- |
* ----------------------- |
252 |
ENDDO |
ENDDO |
253 |
ENDIF |
ENDIF |
254 |
|
|
255 |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
if(TRKDEBUG) print*,'mini : 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) |
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) |
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*,'mini : -ok- ',istep,chi2,AL(5) |
567 |
|
|
568 |
NSTEP=ISTEP ! ***PP*** |
NSTEP=ISTEP ! ***PP*** |
569 |
|
|
1004 |
c$$$ print*,'POSXY (prima) ',vout |
c$$$ print*,'POSXY (prima) ',vout |
1005 |
|
|
1006 |
DO I=1,nplanes |
DO I=1,nplanes |
1007 |
cpp step=vout(3)-zv(i) |
c$$$ ipass = 0 ! TEST |
1008 |
step=vout(3)-zm(i) |
c$$$ PRINT *,'TRACKING -> START PLANE: ',I ! TEST |
1009 |
|
cPPP step=vout(3)-zm(i) |
1010 |
|
cPP step=(zm(i)-vout(3))/VOUT(6) |
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 |
1017 |
|
IFAIL=1 |
1018 |
|
if(TRKVERBOSE) |
1019 |
|
$ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' |
1020 |
|
RETURN |
1021 |
|
ENDIF |
1022 |
|
step=(zm(i)-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 |
1026 |
|
c$$$ PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST |
1027 |
IF(VOUT(3).GT.VECT(3)) THEN |
IF(VOUT(3).GT.VECT(3)) THEN |
1028 |
IFAIL=1 |
IFAIL=1 |
1029 |
if(TRKVERBOSE) |
if(TRKVERBOSE) |
1065 |
VOUT(7) = VOUT(7) * 0.997 !0.9968 |
VOUT(7) = VOUT(7) * 0.997 !0.9968 |
1066 |
* ----------------------------------------------- |
* ----------------------------------------------- |
1067 |
ENDIF |
ENDIF |
1068 |
|
c$$$ PRINT *,'TRACKING -> END' ! TEST |
1069 |
|
|
1070 |
ENDDO |
ENDDO |
1071 |
|
|
1104 |
YM(IP) = -100. !0. |
YM(IP) = -100. !0. |
1105 |
XM_A(IP) = -100. !0. |
XM_A(IP) = -100. !0. |
1106 |
YM_A(IP) = -100. !0. |
YM_A(IP) = -100. !0. |
1107 |
c ZM_A(IP) = 0 |
ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position |
1108 |
XM_B(IP) = -100. !0. |
XM_B(IP) = -100. !0. |
1109 |
YM_B(IP) = -100. !0. |
YM_B(IP) = -100. !0. |
1110 |
c ZM_B(IP) = 0 |
ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position |
1111 |
RESX(IP) = 1000. !3.d-4 |
RESX(IP) = 1000. !3.d-4 |
1112 |
RESY(IP) = 1000. !12.d-4 |
RESY(IP) = 1000. !12.d-4 |
1113 |
XGOOD(IP) = 0 |
XGOOD(IP) = 0 |
1142 |
|
|
1143 |
REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) |
REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) |
1144 |
REAL*4 CHI,XC,ZC,RADIUS |
REAL*4 CHI,XC,ZC,RADIUS |
1145 |
|
|
1146 |
|
c$$$ DO I=1,nplanes |
1147 |
|
c$$$ print *,i,' - ',XGOOD(I),YGOOD(I) |
1148 |
|
c$$$ print *,i,' - ',xm(i),ym(i),zm(i) |
1149 |
|
c$$$ print *,i,' A ',xm_a(i),ym_a(i),zm_a(i) |
1150 |
|
c$$$ print *,i,' B ',xm_b(i),ym_b(i),zm_b(i) |
1151 |
|
c$$$ enddo |
1152 |
|
|
1153 |
|
|
1154 |
* ---------------------------------------- |
* ---------------------------------------- |
1155 |
* Y view |
* Y view |
1156 |
* ---------------------------------------- |
* ---------------------------------------- |
1164 |
S1=0. |
S1=0. |
1165 |
DO I=1,nplanes |
DO I=1,nplanes |
1166 |
IF(YGOOD(I).EQ.1)THEN |
IF(YGOOD(I).EQ.1)THEN |
1167 |
YY = YM(I) |
YY = REAL(YM(I))!EM GCC4.7 |
1168 |
IF(XGOOD(I).EQ.0)THEN |
IF(XGOOD(I).EQ.0)THEN |
1169 |
YY = (YM_A(I) + YM_B(I))/2 |
YY = REAL((YM_A(I) + YM_B(I))/2.)!EM GCC4.7 |
1170 |
ENDIF |
ENDIF |
1171 |
SZZ=SZZ+ZM(I)*ZM(I) |
SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7 |
1172 |
SZY=SZY+ZM(I)*YY |
SZY=SZY+REAL(ZM(I)*YY)!EM GCC4.7 |
1173 |
SSY=SSY+YY |
SSY=SSY+YY |
1174 |
SZ=SZ+ZM(I) |
SZ=SZ+REAL(ZM(I))!EM GCC4.7 |
1175 |
S1=S1+1. |
S1=S1+1. |
1176 |
ENDIF |
ENDIF |
1177 |
ENDDO |
ENDDO |
1178 |
DET=SZZ*S1-SZ*SZ |
DET=SZZ*S1-SZ*SZ |
1179 |
AY=(SZY*S1-SZ*SSY)/DET |
AY=(SZY*S1-SZ*SSY)/DET |
1180 |
BY=(SZZ*SSY-SZY*SZ)/DET |
BY=(SZZ*SSY-SZY*SZ)/DET |
1181 |
Y0 = AY*ZINI+BY |
Y0 = REAL(AY*ZINI+BY)!EM GCC4.7 |
1182 |
* ---------------------------------------- |
* ---------------------------------------- |
1183 |
* X view |
* X view |
1184 |
* ---------------------------------------- |
* ---------------------------------------- |
1188 |
NP=0 |
NP=0 |
1189 |
DO I=1,nplanes |
DO I=1,nplanes |
1190 |
IF(XGOOD(I).EQ.1)THEN |
IF(XGOOD(I).EQ.1)THEN |
1191 |
XX = XM(I) |
XX = REAL(XM(I))!EM GCC4.7 |
1192 |
IF(YGOOD(I).EQ.0)THEN |
IF(YGOOD(I).EQ.0)THEN |
1193 |
XX = (XM_A(I) + XM_B(I))/2 |
XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7 |
1194 |
ENDIF |
ENDIF |
1195 |
NP=NP+1 |
NP=NP+1 |
1196 |
XP(NP)=XX |
XP(NP)=XX |
1197 |
ZP(NP)=ZM(I) |
ZP(NP)=REAL(ZM(I))!EM GCC4.7 |
1198 |
ENDIF |
ENDIF |
1199 |
ENDDO |
ENDDO |
1200 |
IFLAG=0 !no debug mode |
IFLAG=0 !no debug mode |
1208 |
|
|
1209 |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
1210 |
c if(CHI.gt.100)GOTO 10 !straigth fit |
c if(CHI.gt.100)GOTO 10 !straigth fit |
1211 |
ARG = RADIUS**2-(ZINI-ZC)**2 |
ARG = REAL(RADIUS**2-(ZINI-ZC)**2)!EM GCC4.7 |
1212 |
IF(ARG.LT.0)GOTO 10 !straigth fit |
IF(ARG.LT.0)GOTO 10 !straigth fit |
1213 |
DC = SQRT(ARG) |
DC = SQRT(ARG) |
1214 |
IF(XC.GT.0)DC=-DC |
IF(XC.GT.0)DC=-DC |
1215 |
X0=XC+DC |
X0=XC+DC |
1216 |
AX = -(ZINI-ZC)/DC |
AX = REAL(-(ZINI-ZC)/DC)!EM GCC4.7 |
1217 |
DEF=100./(RADIUS*0.3*0.43) |
DEF=100./(RADIUS*0.3*0.43) |
1218 |
IF(XC.GT.0)DEF=-DEF |
IF(XC.GT.0)DEF=-DEF |
1219 |
|
|
1239 |
S1=0. |
S1=0. |
1240 |
DO I=1,nplanes |
DO I=1,nplanes |
1241 |
IF(XGOOD(I).EQ.1)THEN |
IF(XGOOD(I).EQ.1)THEN |
1242 |
XX = XM(I) |
XX = REAL(XM(I))!EM GCC4.7 |
1243 |
IF(YGOOD(I).EQ.0)THEN |
IF(YGOOD(I).EQ.0)THEN |
1244 |
XX = (XM_A(I) + XM_B(I))/2 |
XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7 |
1245 |
ENDIF |
ENDIF |
1246 |
SZZ=SZZ+ZM(I)*ZM(I) |
SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7 |
1247 |
SZX=SZX+ZM(I)*XX |
SZX=SZX+REAL(ZM(I)*XX)!EM GCC4.7 |
1248 |
SSX=SSX+XX |
SSX=SSX+XX |
1249 |
SZ=SZ+ZM(I) |
SZ=SZ+REAL(ZM(I))!EM GCC4.7 |
1250 |
S1=S1+1. |
S1=S1+1. |
1251 |
ENDIF |
ENDIF |
1252 |
ENDDO |
ENDDO |
1254 |
AX=(SZX*S1-SZ*SSX)/DET |
AX=(SZX*S1-SZ*SSX)/DET |
1255 |
BX=(SZZ*SSX-SZX*SZ)/DET |
BX=(SZZ*SSX-SZX*SZ)/DET |
1256 |
DEF = 0 |
DEF = 0 |
1257 |
X0 = AX*ZINI+BX |
X0 = REAL(AX*ZINI+BX)!EM GCC4.7 |
1258 |
|
|
1259 |
20 CONTINUE |
20 CONTINUE |
1260 |
* ---------------------------------------- |
* ---------------------------------------- |
1265 |
AL(2) = Y0 |
AL(2) = Y0 |
1266 |
tath = sqrt(AY**2+AX**2) |
tath = sqrt(AY**2+AX**2) |
1267 |
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 |
|
1268 |
|
|
1269 |
AL(4)=0. |
AL(4)=0. |
1270 |
IF( AX.NE.0.OR.AY.NE.0. ) THEN |
IF( AX.NE.0.OR.AY.NE.0. ) THEN |