12 |
************************************************************************ |
************************************************************************ |
13 |
|
|
14 |
|
|
15 |
SUBROUTINE MINI_2(ISTEP,IFAIL) |
SUBROUTINE MINI2(ISTEP,IFAIL,IPRINT) |
16 |
|
|
17 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
18 |
|
|
22 |
c logical DEBUG |
c logical DEBUG |
23 |
c common/dbg/DEBUG |
c common/dbg/DEBUG |
24 |
|
|
25 |
parameter (inf=1.e8) !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 |
32 |
c ordering, but they maintain their Z coordinates. so plane number 1 is |
c ordering, but they maintain their Z coordinates. so plane number 1 is |
33 |
c the first one that a particle meets, and its Z coordinate is > 0 |
c the first one that a particle meets, and its Z coordinate is > 0 |
34 |
c------------------------------------------------------------------------ |
c------------------------------------------------------------------------ |
35 |
DATA ZINI/23.5/ !z coordinate of the reference plane |
DATA ZINI/23.5/ !!! ***PP*** to be changed !z coordinate of the reference plane |
36 |
|
|
37 |
DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking |
c DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking |
38 |
|
|
39 |
DATA STEPAL/5*1.d-7/ !alpha vector step |
DATA STEPAL/5*1.d-7/ !alpha vector step |
40 |
DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization |
DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization |
42 |
* !the tracking procedure |
* !the tracking procedure |
43 |
DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess |
DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess |
44 |
|
|
45 |
c DATA ALMAX/inf,inf,inf,inf,0.25e2/ !limits on alpha vector components |
c DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components |
46 |
c DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !" |
c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" |
47 |
DATA ALMAX/inf,inf,1.,inf,0.25e2/ !limits on alpha vector components |
DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components |
48 |
DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !" |
DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !" |
49 |
|
|
50 |
|
c$$$ DIMENSION DAL(5) !increment of vector alfa |
51 |
|
DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2 |
52 |
|
|
53 |
|
c elena-------- |
54 |
|
REAL*8 AVRESX,AVRESY |
55 |
|
c elena-------- |
56 |
|
|
|
|
|
|
DIMENSION DAL(5) !increment of vector alfa |
|
|
|
|
57 |
INTEGER IFLAG |
INTEGER IFLAG |
58 |
c-------------------------------------------------------- |
c-------------------------------------------------------- |
59 |
c IFLAG =1 ---- chi2 derivatives computed by using |
c IFLAG =1 ---- chi2 derivatives computed by using |
64 |
c NB: the two metods gives equivalent results BUT |
c NB: the two metods gives equivalent results BUT |
65 |
c method 2 is faster!! |
c method 2 is faster!! |
66 |
c-------------------------------------------------------- |
c-------------------------------------------------------- |
67 |
DATA IFLAG/2/ |
DATA IFLAG/2/ |
68 |
|
|
69 |
|
c LOGICAL TRKDEBUG,TRKVERBOSE |
70 |
|
c COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
71 |
|
LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT |
72 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
73 |
|
|
74 |
|
DIMENSION AL0(5) |
75 |
|
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) |
83 |
|
* |
84 |
|
STUDENT = .false. |
85 |
|
FIRSTSTEPS = .true. |
86 |
|
FIRSTSTUDENT = .true. |
87 |
|
IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true. |
88 |
|
|
89 |
|
IF(IPRINT.EQ.1) THEN |
90 |
|
TRKVERBOSE = .TRUE. |
91 |
|
TRKDEBUG = .FALSE. |
92 |
|
ELSEIF(IPRINT.EQ.2)THEN |
93 |
|
TRKVERBOSE = .TRUE. |
94 |
|
TRKDEBUG = .TRUE. |
95 |
|
ELSE |
96 |
|
TRKVERBOSE = .FALSE. |
97 |
|
TRKDEBUG = .FALSE. |
98 |
|
ENDIF |
99 |
|
|
100 |
|
* ---------------------------------------------------------- |
101 |
|
* evaluate average spatial resolution |
102 |
|
* ---------------------------------------------------------- |
103 |
|
AVRESX = RESXAV |
104 |
|
AVRESY = RESYAV |
105 |
|
NX = 0 !EM GCC4.7 |
106 |
|
NY = 0 !EM GCC4.7 |
107 |
|
DO IP=1,6 |
108 |
|
IF( XGOOD(IP).EQ.1 )THEN |
109 |
|
NX=NX+1!EM GCC4.7 |
110 |
|
AVRESX=AVRESX+RESX(IP) |
111 |
|
ENDIF |
112 |
|
IF( YGOOD(IP).EQ.1 )THEN |
113 |
|
NY=NY+1!EM GCC4.7 |
114 |
|
AVRESY=AVRESY+RESY(IP) |
115 |
|
ENDIF |
116 |
|
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 |
122 |
* |
* |
123 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
124 |
FACT=10. !scale factor to define |
* changed in order to evaluate energy-dependent |
125 |
!tolerance on alfa |
* tolerances on all 5 parameters |
126 |
ALTOL(1)=RESX(1)/FACT !al(1) = x |
cPP FACT=1.0e10 !scale factor to define tolerance on alfa |
127 |
ALTOL(2)=RESY(1)/FACT !al(2) = y |
c deflection error (see PDG) |
128 |
ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta) |
DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) |
129 |
$ +RESY(1)**2)/44.51/FACT |
DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) |
130 |
ALTOL(4)=ALTOL(3) !al(4)=phi |
c$$$ ALTOL(1) = AVRESX/FACT !al(1) = x |
131 |
|
c$$$ ALTOL(2) = AVRESY/FACT !al(2) = y |
132 |
|
c$$$ ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta) |
133 |
|
c$$$ $ +AVRESY**2)/44.51/FACT |
134 |
|
c$$$ ALTOL(4) = ALTOL(3) !al(4)=phi |
135 |
c deflection error (see PDG) |
c deflection error (see PDG) |
136 |
DELETA1=0.01*RESX(1)/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) |
c$$$ DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) |
137 |
DELETA2=0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) |
c$$$ DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) |
138 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
139 |
* |
* |
140 |
ISTEP=0 !num. steps to minimize chi^2 |
ISTEP=0 !num. steps to minimize chi^2 |
141 |
JFAIL=0 !error flag |
JFAIL=0 !error flag |
142 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
CHI2=0 |
143 |
IF(JFAIL.NE.0) THEN |
|
144 |
IFAIL=1 |
if(TRKDEBUG) print*,'guess: ',al |
145 |
if(DEBUG) |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5) |
146 |
$ PRINT *,'mini_2 ===> error on CHISQ computation !!!' |
|
|
RETURN |
|
|
ENDIF |
|
147 |
* |
* |
148 |
* ----------------------- |
* ----------------------- |
149 |
* START MINIMIZATION LOOP |
* START MINIMIZATION LOOP |
150 |
* ----------------------- |
* ----------------------- |
151 |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
|
CHI2_P=CHI2 |
|
|
|
|
|
c print*,'@@@@@ ',istep,' - ',al |
|
152 |
|
|
153 |
cost=1e-7 |
* ------------------------------- |
154 |
DO I=1,5 |
* **** Chi2+gaussian minimization |
155 |
DO J=1,5 |
* ------------------------------- |
156 |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
|
157 |
|
IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN |
158 |
|
|
159 |
|
IF(ISTEP.GE.3) FIRSTSTEPS = .false. |
160 |
|
|
161 |
|
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
162 |
|
IF(JFAIL.NE.0) THEN |
163 |
|
IFAIL=1 |
164 |
|
CHI2=-9999. |
165 |
|
if(TRKVERBOSE) |
166 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
167 |
|
RETURN |
168 |
|
ENDIF |
169 |
|
|
170 |
|
c COST=1e-5 |
171 |
|
COST=1. |
172 |
|
DO I=1,5 |
173 |
|
IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 |
174 |
ENDDO |
ENDDO |
175 |
CHI2D(I)=CHI2D(I)*COST |
DO I=1,5 |
176 |
ENDDO |
DO J=1,5 |
177 |
|
CHI2DD(I,J)=CHI2DD(I,J)*COST |
178 |
|
ENDDO |
179 |
|
c$$$ CHI2D(I)=CHI2D(I)*COST |
180 |
|
ENDDO |
181 |
|
|
182 |
|
IF(PFIXED.EQ.0.) THEN |
183 |
|
|
184 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
185 |
* track fitting with FREE deflection |
* track fitting with FREE deflection |
186 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
187 |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
188 |
IF(IFA.NE.0) THEN !not positive-defined |
IF(IFA.NE.0) THEN !not positive-defined |
189 |
if(DEBUG)then |
if(TRKVERBOSE)then |
190 |
PRINT *, |
PRINT *, |
191 |
$ 'MINI_HOUGH ==> '// |
$ '*** ERROR in mini ***'// |
192 |
$ '** ERROR ** on matrix inversion (not positive-defined)!!!' |
$ 'on matrix inversion (not pos-def)' |
193 |
$ ,DET |
$ ,DET |
194 |
endif |
endif |
195 |
IFAIL=1 |
IF(CHI2.EQ.0) CHI2=-9999. |
196 |
RETURN |
IF(CHI2.GT.0) CHI2=-CHI2 |
197 |
ENDIF |
IFAIL=1 |
198 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
RETURN |
199 |
|
ENDIF |
200 |
|
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
201 |
* ******************************************* |
* ******************************************* |
202 |
* find new value of AL-pha !* |
* find new value of AL-pha |
203 |
* !* |
* ******************************************* |
204 |
DO I=1,5 !* |
DO I=1,5 |
205 |
DAL(I)=0. !* |
DAL(I)=0. |
206 |
DO J=1,5 !* |
DO J=1,5 |
207 |
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) !* |
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) *COST |
208 |
COV(I,J)=CHI2DD(I,J) !* |
COV(I,J)=2.*COST*CHI2DD(I,J) |
209 |
ENDDO !* |
ENDDO |
210 |
ENDDO !* |
ENDDO |
211 |
DO I=1,5 !* |
DO I=1,5 |
212 |
AL(I)=AL(I)+DAL(I) !* |
AL(I)=AL(I)+DAL(I) |
213 |
ENDDO !* |
ENDDO |
214 |
|
*------------------------------------------------------------* |
215 |
|
* track fitting with FIXED deflection |
216 |
|
*------------------------------------------------------------* |
217 |
|
ELSE |
218 |
|
AL(5)=1./PFIXED |
219 |
|
DO I=1,4 |
220 |
|
CHI2D_R(I)=CHI2D(I) |
221 |
|
DO J=1,4 |
222 |
|
CHI2DD_R(I,J)=CHI2DD(I,J) |
223 |
|
ENDDO |
224 |
|
ENDDO |
225 |
|
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
226 |
|
IF(IFA.NE.0) THEN |
227 |
|
if(TRKVERBOSE)then |
228 |
|
PRINT *, |
229 |
|
$ '*** ERROR in mini ***'// |
230 |
|
$ 'on matrix inversion (not pos-def)' |
231 |
|
$ ,DET |
232 |
|
endif |
233 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
234 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
235 |
|
IFAIL=1 |
236 |
|
RETURN |
237 |
|
ENDIF |
238 |
|
CALL DSFINV(4,CHI2DD_R,4) |
239 |
* ******************************************* |
* ******************************************* |
240 |
|
* find new value of AL-pha |
241 |
|
* ******************************************* |
242 |
|
DO I=1,4 |
243 |
|
DAL(I)=0. |
244 |
|
DO J=1,4 |
245 |
|
DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) *COST |
246 |
|
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
247 |
|
ENDDO |
248 |
|
ENDDO |
249 |
|
DAL(5)=0. |
250 |
|
DO I=1,4 |
251 |
|
AL(I)=AL(I)+DAL(I) |
252 |
|
ENDDO |
253 |
|
ENDIF |
254 |
|
|
255 |
|
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,AL(5) |
256 |
|
|
257 |
|
c$$$ PRINT*,'DAL ',(DAL(K),K=1,5) |
258 |
|
c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5) |
259 |
|
|
260 |
|
|
261 |
|
ENDIF |
262 |
|
|
263 |
|
* ------------------------------- |
264 |
|
* **** Likelihood+Student minimization |
265 |
|
* ------------------------------- |
266 |
|
|
267 |
|
IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN |
268 |
|
|
269 |
|
IF(FIRSTSTUDENT) THEN |
270 |
|
FIRSTSTUDENT = .false. |
271 |
|
ISTEP = 1 |
272 |
|
ENDIF |
273 |
|
|
274 |
|
CALL CHISQSTT(1,JFAIL) |
275 |
|
DO I=1,5 |
276 |
|
DAL(I)=0. |
277 |
|
DO J=1,5 |
278 |
|
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) |
279 |
|
ENDDO |
280 |
|
ENDDO |
281 |
|
|
282 |
|
DO I=1,5 |
283 |
|
DO j=1,5 |
284 |
|
COV(I,J) = 2.*CHI2DD(I,J) |
285 |
|
ENDDO |
286 |
|
ENDDO |
287 |
|
|
288 |
|
CHI2TOLL = 1.E-3 |
289 |
|
ALPHA = 3.0 |
290 |
|
BETA = -0.4 |
291 |
|
E=1. |
292 |
|
EA=1. |
293 |
|
EB=1. |
294 |
|
EC=1. |
295 |
|
FA=1. |
296 |
|
FB=1. |
297 |
|
FC=1. |
298 |
|
SUCCESS_OLD = .FALSE. |
299 |
|
SUCCESS_NEW = .FALSE. |
300 |
|
|
301 |
|
CALL CHISQSTT(0,JFAIL) |
302 |
|
c$$$ PRINT*,CHI2 |
303 |
|
CHI2_NEW = CHI2 |
304 |
|
FC = CHI2 |
305 |
|
EC = 0. |
306 |
|
|
307 |
|
ICOUNT = 0 |
308 |
|
100 CONTINUE |
309 |
|
ICOUNT = ICOUNT+1 |
310 |
|
|
311 |
|
DO I=1,5 |
312 |
|
AL0(I)=AL(I) |
313 |
|
ENDDO |
314 |
|
DO I=1,5 |
315 |
|
AL(I)=AL(I)+E*DAL(I) |
316 |
|
ENDDO |
317 |
|
CALL CHISQSTT(0,JFAIL) |
318 |
|
CHI2_OLD = CHI2_NEW |
319 |
|
CHI2_NEW = CHI2 |
320 |
|
FA = FB |
321 |
|
FB = FC |
322 |
|
FC = CHI2 |
323 |
|
EA = EB |
324 |
|
EB = EC |
325 |
|
EC = E |
326 |
|
|
327 |
|
c$$$ PRINT*,E,CHI2_NEW |
328 |
|
|
329 |
|
IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success |
330 |
|
IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101 |
331 |
|
SUCCESS_OLD = SUCCESS_NEW |
332 |
|
SUCCESS_NEW = .TRUE. |
333 |
|
E = E*ALPHA |
334 |
|
ELSE ! failure |
335 |
|
SUCCESS_OLD = SUCCESS_NEW |
336 |
|
SUCCESS_NEW = .FALSE. |
337 |
|
CHI2_NEW = CHI2_OLD |
338 |
|
DO I=1,5 |
339 |
|
AL(I)=AL0(I) |
340 |
|
ENDDO |
341 |
|
IF(SUCCESS_OLD) THEN |
342 |
|
DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA) |
343 |
|
IF(DENOM.NE.0.) THEN |
344 |
|
E = EB - 0.5*( (EB-EA)**2*(FB-FC) |
345 |
|
$ - (EB-EC)**2*(FB-FA) ) / DENOM |
346 |
|
ELSE |
347 |
|
E = BETA*E |
348 |
|
ENDIF |
349 |
|
ELSE |
350 |
|
E = BETA*E |
351 |
|
ENDIF |
352 |
|
c$$$ E = BETA*E |
353 |
|
ENDIF |
354 |
|
IF(ICOUNT.GT.20) GOTO 101 |
355 |
|
GOTO 100 |
356 |
|
|
357 |
|
101 CONTINUE |
358 |
|
|
359 |
|
DO I=1,5 |
360 |
|
DAL(I)=E*DAL(I) |
361 |
|
ENDDO |
362 |
|
|
363 |
|
c$$$ print*,' ' |
364 |
|
c$$$ PRINT*,'DAL ',(DAL(K),K=1,5) |
365 |
|
c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5) |
366 |
|
c$$$ print*,'==== CHI2 ====' |
367 |
|
c$$$ print*,chi2 |
368 |
|
c$$$ print*,'==== CHI2d ====' |
369 |
|
c$$$ print*,(chi2d(i),i=1,5) |
370 |
|
c$$$ print*,'==== CHI2dd ====' |
371 |
|
c$$$ do j=1,5 |
372 |
|
c$$$ print*,(chi2dd(j,i),i=1,5) |
373 |
|
c$$$ enddo |
374 |
|
c$$$ print*,'================' |
375 |
|
c$$$ print*,' ' |
376 |
|
|
377 |
|
*========= FIN QUI ============= |
378 |
|
|
379 |
|
ENDIF |
380 |
|
|
381 |
|
|
382 |
|
|
383 |
|
|
384 |
|
|
385 |
|
*------------------------------------------------------------* |
386 |
|
* ---------------------------------------------------- * |
387 |
|
*------------------------------------------------------------* |
388 |
* check parameter bounds: |
* check parameter bounds: |
389 |
|
*------------------------------------------------------------* |
390 |
DO I=1,5 |
DO I=1,5 |
391 |
IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN |
IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN |
392 |
if(DEBUG)then |
if(TRKVERBOSE)then |
393 |
PRINT*,' **WARNING** ' |
PRINT*,' *** WARNING in mini *** ' |
394 |
PRINT*,'MINI_2 ==> AL(',I,') out of range' |
PRINT*,'MINI_2 ==> AL(',I,') out of range' |
395 |
PRINT*,' value: ',AL(I), |
PRINT*,' value: ',AL(I), |
396 |
$ ' limits: ',ALMIN(I),ALMAX(I) |
$ ' limits: ',ALMIN(I),ALMAX(I) |
397 |
print*,'istep ',istep |
print*,'istep ',istep |
398 |
endif |
endif |
399 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
400 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
401 |
IFAIL=1 |
IFAIL=1 |
402 |
RETURN |
RETURN |
403 |
ENDIF |
ENDIF |
404 |
ENDDO |
ENDDO |
405 |
* new estimate of chi^2: |
*------------------------------------------------------------* |
|
JFAIL=0 !error flag |
|
|
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
|
|
IF(JFAIL.NE.0) THEN |
|
|
IFAIL=1 |
|
|
if(DEBUG) |
|
|
$ PRINT *,'mini_2: ===> error on CHISQ computation !!!' |
|
|
RETURN |
|
|
ENDIF |
|
406 |
* check number of steps: |
* check number of steps: |
407 |
IF(ISTEP.gt.ISTEPMAX) then |
*------------------------------------------------------------* |
408 |
IFAIL=1 |
IF(ISTEP.ge.ISTEPMAX) then |
409 |
if(DEBUG) |
c$$$ IFAIL=1 |
410 |
$ PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX |
c$$$ if(TRKVERBOSE) |
411 |
|
c$$$ $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', |
412 |
|
c$$$ $ ISTEPMAX |
413 |
goto 11 |
goto 11 |
414 |
endif |
endif |
415 |
|
*------------------------------------------------------------* |
416 |
* --------------------------------------------- |
* --------------------------------------------- |
417 |
* evaluate deflection tolerance on the basis of |
* evaluate deflection tolerance on the basis of |
418 |
* estimated deflection |
* estimated deflection |
419 |
* --------------------------------------------- |
* --------------------------------------------- |
420 |
ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT |
*------------------------------------------------------------* |
421 |
|
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 |
427 |
|
ALTOL(1) = ALTOL(5)/DELETA1 |
428 |
|
ALTOL(2) = ALTOL(1) |
429 |
|
ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 |
430 |
|
ALTOL(4) = ALTOL(3) |
431 |
|
|
432 |
|
c$$$ print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step! |
433 |
|
|
434 |
*---- check tolerances: |
*---- check tolerances: |
435 |
|
c$$$ DO I=1,5 |
436 |
|
c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! |
437 |
|
c$$$ ENDDO |
438 |
|
c$$$ print*,'chi2 -- ',DCHI2 |
439 |
|
|
440 |
|
IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP*** |
441 |
DO I=1,5 |
DO I=1,5 |
442 |
IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! |
IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! |
443 |
ENDDO |
ENDDO |
444 |
|
|
445 |
|
***************************** |
446 |
|
* final estimate of chi^2 |
447 |
|
***************************** |
448 |
|
|
449 |
|
* ------------------------------- |
450 |
|
* **** Chi2+gaussian minimization |
451 |
|
* ------------------------------- |
452 |
|
|
453 |
|
IF(.NOT.STUDENT) THEN |
454 |
|
|
455 |
|
JFAIL=0 !error flag |
456 |
|
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
457 |
|
IF(JFAIL.NE.0) THEN |
458 |
|
IFAIL=1 |
459 |
|
if(TRKVERBOSE)THEN |
460 |
|
CHI2=-9999. |
461 |
|
if(TRKVERBOSE) |
462 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
463 |
|
ENDIF |
464 |
|
RETURN |
465 |
|
ENDIF |
466 |
|
c COST=1e-7 |
467 |
|
COST=1. |
468 |
|
DO I=1,5 |
469 |
|
IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 |
470 |
|
ENDDO |
471 |
|
DO I=1,5 |
472 |
|
DO J=1,5 |
473 |
|
CHI2DD(I,J)=CHI2DD(I,J)*COST |
474 |
|
ENDDO |
475 |
|
ENDDO |
476 |
|
IF(PFIXED.EQ.0.) THEN |
477 |
|
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
478 |
|
IF(IFA.NE.0) THEN !not positive-defined |
479 |
|
if(TRKVERBOSE)then |
480 |
|
PRINT *, |
481 |
|
$ '*** ERROR in mini ***'// |
482 |
|
$ 'on matrix inversion (not pos-def)' |
483 |
|
$ ,DET |
484 |
|
endif |
485 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
486 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
487 |
|
IFAIL=1 |
488 |
|
RETURN |
489 |
|
ENDIF |
490 |
|
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
491 |
|
DO I=1,5 |
492 |
|
c$$$ DAL(I)=0. |
493 |
|
DO J=1,5 |
494 |
|
COV(I,J)=2.*COST*CHI2DD(I,J) |
495 |
|
ENDDO |
496 |
|
ENDDO |
497 |
|
ELSE |
498 |
|
DO I=1,4 |
499 |
|
CHI2D_R(I)=CHI2D(I) |
500 |
|
DO J=1,4 |
501 |
|
CHI2DD_R(I,J)=CHI2DD(I,J) |
502 |
|
ENDDO |
503 |
|
ENDDO |
504 |
|
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
505 |
|
IF(IFA.NE.0) THEN |
506 |
|
if(TRKVERBOSE)then |
507 |
|
PRINT *, |
508 |
|
$ '*** ERROR in mini ***'// |
509 |
|
$ 'on matrix inversion (not pos-def)' |
510 |
|
$ ,DET |
511 |
|
endif |
512 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
513 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
514 |
|
IFAIL=1 |
515 |
|
RETURN |
516 |
|
ENDIF |
517 |
|
CALL DSFINV(4,CHI2DD_R,4) |
518 |
|
DO I=1,4 |
519 |
|
c$$$ DAL(I)=0. |
520 |
|
DO J=1,4 |
521 |
|
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
522 |
|
ENDDO |
523 |
|
ENDDO |
524 |
|
ENDIF |
525 |
|
|
526 |
|
ENDIF |
527 |
|
|
528 |
|
* ------------------------------- |
529 |
|
* **** Likelihood+student minimization |
530 |
|
* ------------------------------- |
531 |
|
|
532 |
|
IF(STUDENT) THEN |
533 |
|
CALL CHISQSTT(1,JFAIL) |
534 |
|
DO I=1,5 |
535 |
|
DO j=1,5 |
536 |
|
COV(I,J) = 2.*CHI2DD(I,J) |
537 |
|
ENDDO |
538 |
|
ENDDO |
539 |
|
ENDIF |
540 |
|
|
541 |
|
***************************** |
542 |
|
|
543 |
* ------------------------------------ |
* ------------------------------------ |
544 |
* Number of Degree Of Freedom |
* Number of Degree Of Freedom |
548 |
$ +int(xgood(ip)) |
$ +int(xgood(ip)) |
549 |
$ +int(ygood(ip)) |
$ +int(ygood(ip)) |
550 |
enddo |
enddo |
551 |
ndof=ndof-5 |
if(pfixed.eq.0.) ndof=ndof-5 ! ***PP*** |
552 |
|
if(pfixed.ne.0.) ndof=ndof-4 ! ***PP*** |
553 |
|
if(ndof.le.0.) then |
554 |
|
ndof = 1 |
555 |
|
if(TRKVERBOSE) |
556 |
|
$ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' |
557 |
|
endif |
558 |
|
|
559 |
* ------------------------------------ |
* ------------------------------------ |
560 |
* Reduced chi^2 |
* Reduced chi^2 |
561 |
CHI2 = CHI2/dble(ndof) |
CHI2 = CHI2/dble(ndof) |
562 |
|
c print*,'mini2: chi2 ',chi2 |
563 |
|
|
|
|
|
564 |
11 CONTINUE |
11 CONTINUE |
565 |
|
|
566 |
101 CONTINUE |
if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,AL(5) |
567 |
|
|
568 |
|
NSTEP=ISTEP ! ***PP*** |
569 |
|
|
570 |
c print*,'END MINI' |
c$$$ print*,'>>>>> NSTEP = ',NSTEP |
571 |
|
|
572 |
RETURN |
RETURN |
573 |
END |
END |
594 |
DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes) |
DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes) |
595 |
$ ,XV0(nplanes),YV0(nplanes) |
$ ,XV0(nplanes),YV0(nplanes) |
596 |
DIMENSION AL_P(5) |
DIMENSION AL_P(5) |
597 |
|
|
598 |
|
c LOGICAL TRKVERBOSE |
599 |
|
c COMMON/TRKD/TRKVERBOSE |
600 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
601 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
602 |
* |
* |
603 |
* chi^2 computation |
* chi^2 computation |
604 |
* |
* |
608 |
JFAIL=0 !error flag |
JFAIL=0 !error flag |
609 |
CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes |
CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes |
610 |
IF(JFAIL.NE.0) THEN |
IF(JFAIL.NE.0) THEN |
611 |
PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' |
IF(TRKVERBOSE) |
612 |
|
$ PRINT *,'CHISQ ==> error from trk routine POSXYZ !!' |
613 |
IFAIL=1 |
IFAIL=1 |
614 |
RETURN |
RETURN |
615 |
ENDIF |
ENDIF |
634 |
* measured position of the cluster. |
* measured position of the cluster. |
635 |
* --------------------------------------------------------- |
* --------------------------------------------------------- |
636 |
CHI2=0. |
CHI2=0. |
637 |
DO I=1,nplanes |
DO I=1,nplanes |
638 |
IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl |
IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl |
639 |
BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I)) |
BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I)) |
640 |
ALFA = XM_A(I) - BETA * YM_A(I) |
ALFA = XM_A(I) - BETA * YM_A(I) |
661 |
+ *( XGOOD(I)*(1-YGOOD(I)) ) |
+ *( XGOOD(I)*(1-YGOOD(I)) ) |
662 |
+ +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 |
+ +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 |
663 |
+ *( (1-XGOOD(I))*YGOOD(I) ) |
+ *( (1-XGOOD(I))*YGOOD(I) ) |
664 |
|
c$$$ print*,(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) ) |
665 |
|
c$$$ print*,(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) ) |
666 |
|
c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 |
667 |
|
c$$$ + *( XGOOD(I)*(1-YGOOD(I)) ) |
668 |
|
c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 |
669 |
|
c$$$ + *( (1-XGOOD(I))*YGOOD(I) ) |
670 |
|
c$$$ print*,XV(I),XM(I),XGOOD(I) |
671 |
|
c$$$ print*,YV(I),YM(I),YGOOD(I) |
672 |
ENDDO |
ENDDO |
673 |
|
c$$$ print*,'CHISQ ',chi2 |
674 |
* ------------------------------------------------ |
* ------------------------------------------------ |
675 |
* |
* |
676 |
* calculation of derivatives (dX/dAL_fa and dY/dAL_fa) |
* calculation of derivatives (dX/dAL_fa and dY/dAL_fa) |
689 |
JFAIL=0 |
JFAIL=0 |
690 |
CALL POSXYZ(AL_P,JFAIL) |
CALL POSXYZ(AL_P,JFAIL) |
691 |
IF(JFAIL.NE.0) THEN |
IF(JFAIL.NE.0) THEN |
692 |
PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' |
IF(TRKVERBOSE) |
693 |
|
*23456789012345678901234567890123456789012345678901234567890123456789012 |
694 |
|
$ PRINT *,'CHISQ ==> error from trk routine POSXYZ' |
695 |
IFAIL=1 |
IFAIL=1 |
696 |
RETURN |
RETURN |
697 |
ENDIF |
ENDIF |
703 |
JFAIL=0 |
JFAIL=0 |
704 |
CALL POSXYZ(AL_P,JFAIL) |
CALL POSXYZ(AL_P,JFAIL) |
705 |
IF(JFAIL.NE.0) THEN |
IF(JFAIL.NE.0) THEN |
706 |
PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' |
IF(TRKVERBOSE) |
707 |
|
$ PRINT *,'CHISQ ==> error from trk routine POSXYZ' |
708 |
IFAIL=1 |
IFAIL=1 |
709 |
RETURN |
RETURN |
710 |
ENDIF |
ENDIF |
735 |
|
|
736 |
COSTHE=DSQRT(1.-AL(3)**2) |
COSTHE=DSQRT(1.-AL(3)**2) |
737 |
IF(COSTHE.EQ.0.) THEN |
IF(COSTHE.EQ.0.) THEN |
738 |
PRINT *,'=== WARNING ===> COSTHE=0' |
IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' |
739 |
STOP |
IFAIL=1 |
740 |
|
RETURN |
741 |
ENDIF |
ENDIF |
742 |
|
|
743 |
DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 |
DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 |
790 |
RETURN |
RETURN |
791 |
END |
END |
792 |
|
|
793 |
|
****************************************************************************** |
794 |
|
* |
795 |
|
* routine to compute Likelihodd+Student and its derivatives |
796 |
|
* |
797 |
|
* (modified in respect to the previous one in order to include |
798 |
|
* single clusters. In this case the residual is evaluated by |
799 |
|
* calculating the distance between the track intersection and the |
800 |
|
* segment AB associated to the single cluster) |
801 |
|
* |
802 |
|
****************************************************************************** |
803 |
|
|
804 |
|
SUBROUTINE CHISQSTT(IFLAG,JFAIL) |
805 |
|
|
806 |
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
807 |
|
|
808 |
|
include 'commontracker.f' !tracker general common |
809 |
|
include 'common_mini_2.f' !common for the tracking procedure |
810 |
|
|
811 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
812 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
813 |
|
|
814 |
|
DIMENSION AL_P(5) |
815 |
|
DIMENSION VECTEMP(5) |
816 |
|
c$$$ DIMENSION U(5) ! BFGS |
817 |
|
|
818 |
|
DO I=1,5 |
819 |
|
AL_P(I)=AL(I) |
820 |
|
ENDDO |
821 |
|
JFAIL=0 !error flag |
822 |
|
CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes |
823 |
|
IF(JFAIL.NE.0) THEN |
824 |
|
IF(TRKVERBOSE) |
825 |
|
$ PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!' |
826 |
|
IFAIL=1 |
827 |
|
RETURN |
828 |
|
ENDIF |
829 |
|
|
830 |
|
DO I=1,nplanes |
831 |
|
DXDAL(I,1)=1. |
832 |
|
DYDAL(I,1)=0. |
833 |
|
DXDAL(I,2)=0. |
834 |
|
DYDAL(I,2)=1. |
835 |
|
COSTHE=DSQRT(1.-AL(3)**2) |
836 |
|
IF(COSTHE.EQ.0.) THEN |
837 |
|
IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' |
838 |
|
IFAIL=1 |
839 |
|
RETURN |
840 |
|
ENDIF |
841 |
|
DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 |
842 |
|
DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3 |
843 |
|
DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE |
844 |
|
DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE |
845 |
|
IF(AL(5).NE.0.) THEN |
846 |
|
DXDAL(I,5)= |
847 |
|
+ (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I)) |
848 |
|
+ *DCOS(AL(4))))/AL(5) |
849 |
|
DYDAL(I,5)= |
850 |
|
+ (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I)) |
851 |
|
+ *DSIN(AL(4))))/AL(5) |
852 |
|
ELSE |
853 |
|
DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 ) |
854 |
|
DYDAL(I,5)=0. |
855 |
|
ENDIF |
856 |
|
ENDDO |
857 |
|
|
858 |
|
IF(IFLAG.EQ.0) THEN ! function calulation |
859 |
|
CHI2=0. |
860 |
|
DO I=1,nplanes |
861 |
|
IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl |
862 |
|
BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I)) |
863 |
|
ALFA = XM_A(I) - BETA * YM_A(I) |
864 |
|
YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2) |
865 |
|
if(YM(I).lt.dmin1(YM_A(I),YM_B(I))) |
866 |
|
$ YM(I)=dmin1(YM_A(I),YM_B(I)) |
867 |
|
if(YM(I).gt.dmax1(YM_A(I),YM_B(I))) |
868 |
|
$ YM(I)=dmax1(YM_A(I),YM_B(I)) |
869 |
|
XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates |
870 |
|
ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl |
871 |
|
BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I)) |
872 |
|
ALFA = YM_A(I) - BETA * XM_A(I) |
873 |
|
XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2) |
874 |
|
if(XM(I).lt.dmin1(XM_A(I),XM_B(I))) |
875 |
|
$ XM(I)=dmin1(XM_A(I),XM_B(I)) |
876 |
|
if(XM(I).gt.dmax1(XM_A(I),XM_B(I))) |
877 |
|
$ XM(I)=dmax1(XM_A(I),XM_B(I)) |
878 |
|
YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates |
879 |
|
ENDIF |
880 |
|
TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/ |
881 |
|
$ (TAILX(I)*RESX(I)**2) ) |
882 |
|
TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/ |
883 |
|
$ (TAILY(I)*RESY(I)**2) ) |
884 |
|
CHI2=CHI2 |
885 |
|
$ +(TAILX(I)+1.0)*TERMX *( XGOOD(I) ) |
886 |
|
$ +(TAILY(I)+1.0)*TERMY *( YGOOD(I) ) |
887 |
|
ENDDO |
888 |
|
ENDIF |
889 |
|
|
890 |
|
IF(IFLAG.EQ.1) THEN ! derivative calulation |
891 |
|
DO I=1,5 |
892 |
|
CHI2DOLD(I)=CHI2D(I) |
893 |
|
ENDDO |
894 |
|
DO J=1,5 |
895 |
|
CHI2D(J)=0. |
896 |
|
DO I=1,nplanes |
897 |
|
CHI2D(J)=CHI2D(J) |
898 |
|
$ +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/ |
899 |
|
$ (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)* |
900 |
|
$ DXDAL(I,J) *XGOOD(I) |
901 |
|
$ +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/ |
902 |
|
$ (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)* |
903 |
|
$ DYDAL(I,J) *YGOOD(I) |
904 |
|
ENDDO |
905 |
|
ENDDO |
906 |
|
DO K=1,5 |
907 |
|
VECTEMP(K)=0. |
908 |
|
DO M=1,5 |
909 |
|
VECTEMP(K) = VECTEMP(K) + |
910 |
|
$ COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M)) |
911 |
|
ENDDO |
912 |
|
ENDDO |
913 |
|
DOWN1 = 0. |
914 |
|
DO K=1,5 |
915 |
|
DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K)) |
916 |
|
ENDDO |
917 |
|
IF(DOWN1.EQ.0.) THEN |
918 |
|
PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0' |
919 |
|
IFAIL=1 |
920 |
|
RETURN |
921 |
|
ENDIF |
922 |
|
DOWN2 = 0. |
923 |
|
DO K=1,5 |
924 |
|
DO M=1,5 |
925 |
|
DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K) |
926 |
|
ENDDO |
927 |
|
ENDDO |
928 |
|
IF(DOWN2.EQ.0.) THEN |
929 |
|
PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0' |
930 |
|
IFAIL=1 |
931 |
|
RETURN |
932 |
|
ENDIF |
933 |
|
c$$$ DO K=1,5 ! BFGS |
934 |
|
c$$$ U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2 |
935 |
|
c$$$ ENDDO |
936 |
|
DO I=1,5 |
937 |
|
DO J=1,5 |
938 |
|
CHI2DD(I,J) = COV(I,J)/2. |
939 |
|
$ +DAL(I)*DAL(J)/DOWN1 |
940 |
|
$ -VECTEMP(I)*VECTEMP(J)/DOWN2 |
941 |
|
c$$$ $ +DOWN2*U(I)*U(J) ! BFGS |
942 |
|
ENDDO |
943 |
|
ENDDO |
944 |
|
ENDIF |
945 |
|
|
946 |
|
RETURN |
947 |
|
END |
948 |
|
|
949 |
***************************************************************** |
***************************************************************** |
950 |
* |
* |
951 |
* Routine to compute the track intersection points |
* Routine to compute the track intersection points |
976 |
|
|
977 |
include 'commontracker.f' !tracker general common |
include 'commontracker.f' !tracker general common |
978 |
include 'common_mini_2.f' !common for the tracking procedure |
include 'common_mini_2.f' !common for the tracking procedure |
979 |
|
|
980 |
|
c LOGICAL TRKVERBOSE |
981 |
|
c COMMON/TRKD/TRKVERBOSE |
982 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
983 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
984 |
c |
c |
985 |
DIMENSION AL_P(5) |
DIMENSION AL_P(5) |
986 |
* |
* |
987 |
DO I=1,nplanes |
cpp DO I=1,nplanes |
988 |
ZV(I)=ZM(I) ! |
cpp ZV(I)=ZM(I) ! |
989 |
ENDDO |
cpp ENDDO |
990 |
* |
* |
991 |
* set parameters for GRKUTA |
* set parameters for GRKUTA |
992 |
* |
* |
1000 |
VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
1001 |
IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
1002 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
1003 |
|
|
1004 |
|
c$$$ print*,'POSXY (prima) ',vout |
1005 |
|
|
1006 |
DO I=1,nplanes |
DO I=1,nplanes |
1007 |
step=vout(3)-zv(i) |
c$$$ ipass = 0 ! TEST |
1008 |
|
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(WARNING) |
if(TRKVERBOSE) |
1030 |
$ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' |
$ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' |
1031 |
if(WARNING)print*,'charge',charge |
c$$$ if(.TRUE.)print*,'charge',charge |
1032 |
if(WARNING)print*,'vect',vect |
c$$$ if(.TRUE.)print*,'vect',vect |
1033 |
if(WARNING)print*,'vout',vout |
c$$$ if(.TRUE.)print*,'vout',vout |
1034 |
if(WARNING)print*,'step',step |
c$$$ if(.TRUE.)print*,'step',step |
1035 |
|
if(TRKVERBOSE)print*,'charge',charge |
1036 |
|
if(TRKVERBOSE)print*,'vect',vect |
1037 |
|
if(TRKVERBOSE)print*,'vout',vout |
1038 |
|
if(TRKVERBOSE)print*,'step',step |
1039 |
RETURN |
RETURN |
1040 |
ENDIF |
ENDIF |
1041 |
Z=VOUT(3) |
Z=VOUT(3) |
1049 |
GOTO 11 |
GOTO 11 |
1050 |
ENDIF |
ENDIF |
1051 |
|
|
1052 |
|
|
1053 |
* ----------------------------------------------- |
* ----------------------------------------------- |
1054 |
* evaluate track coordinates |
* evaluate track coordinates |
1055 |
100 XV(I)=VOUT(1) |
100 XV(I)=VOUT(1) |
1059 |
AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
1060 |
* ----------------------------------------------- |
* ----------------------------------------------- |
1061 |
|
|
1062 |
|
IF(TRACKMODE.EQ.1) THEN |
1063 |
|
* ----------------------------------------------- |
1064 |
|
* change of energy by bremsstrahlung for electrons |
1065 |
|
VOUT(7) = VOUT(7) * 0.997 !0.9968 |
1066 |
|
* ----------------------------------------------- |
1067 |
|
ENDIF |
1068 |
|
c$$$ PRINT *,'TRACKING -> END' ! TEST |
1069 |
|
|
1070 |
ENDDO |
ENDDO |
1071 |
|
|
1072 |
|
c$$$ print*,'POSXY (dopo) ',vout |
1073 |
|
|
1074 |
|
|
1075 |
RETURN |
RETURN |
1076 |
END |
END |
1077 |
|
|
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 |
1114 |
YGOOD(IP) = 0 |
YGOOD(IP) = 0 |
1115 |
|
DEDXTRK_X(IP) = 0 |
1116 |
|
DEDXTRK_Y(IP) = 0 |
1117 |
|
AXV(IP) = 0 |
1118 |
|
AYV(IP) = 0 |
1119 |
|
XV(IP) = -100 |
1120 |
|
YV(IP) = -100 |
1121 |
enddo |
enddo |
1122 |
|
|
1123 |
return |
return |
1124 |
end |
end |
1125 |
|
|
1126 |
|
|
1127 |
|
*************************************************** |
1128 |
|
* * |
1129 |
|
* * |
1130 |
|
* * |
1131 |
|
* * |
1132 |
|
* * |
1133 |
|
* * |
1134 |
|
************************************************** |
1135 |
|
|
1136 |
|
subroutine guess() |
1137 |
|
|
1138 |
|
c IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
1139 |
|
|
1140 |
|
include 'commontracker.f' !tracker general common |
1141 |
|
include 'common_mini_2.f' !common for the tracking procedure |
1142 |
|
|
1143 |
|
REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) |
1144 |
|
REAL*4 CHI,XC,ZC,RADIUS |
1145 |
|
* ---------------------------------------- |
1146 |
|
* Y view |
1147 |
|
* ---------------------------------------- |
1148 |
|
* ---------------------------------------- |
1149 |
|
* initial guess with a straigth line |
1150 |
|
* ---------------------------------------- |
1151 |
|
SZZ=0. |
1152 |
|
SZY=0. |
1153 |
|
SSY=0. |
1154 |
|
SZ=0. |
1155 |
|
S1=0. |
1156 |
|
DO I=1,nplanes |
1157 |
|
IF(YGOOD(I).EQ.1)THEN |
1158 |
|
YY = REAL(YM(I))!EM GCC4.7 |
1159 |
|
IF(XGOOD(I).EQ.0)THEN |
1160 |
|
YY = REAL((YM_A(I) + YM_B(I))/2.)!EM GCC4.7 |
1161 |
|
ENDIF |
1162 |
|
SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7 |
1163 |
|
SZY=SZY+REAL(ZM(I)*YY)!EM GCC4.7 |
1164 |
|
SSY=SSY+YY |
1165 |
|
SZ=SZ+REAL(ZM(I))!EM GCC4.7 |
1166 |
|
S1=S1+1. |
1167 |
|
ENDIF |
1168 |
|
ENDDO |
1169 |
|
DET=SZZ*S1-SZ*SZ |
1170 |
|
AY=(SZY*S1-SZ*SSY)/DET |
1171 |
|
BY=(SZZ*SSY-SZY*SZ)/DET |
1172 |
|
Y0 = REAL(AY*ZINI+BY)!EM GCC4.7 |
1173 |
|
* ---------------------------------------- |
1174 |
|
* X view |
1175 |
|
* ---------------------------------------- |
1176 |
|
* ---------------------------------------- |
1177 |
|
* 1) initial guess with a circle |
1178 |
|
* ---------------------------------------- |
1179 |
|
NP=0 |
1180 |
|
DO I=1,nplanes |
1181 |
|
IF(XGOOD(I).EQ.1)THEN |
1182 |
|
XX = REAL(XM(I))!EM GCC4.7 |
1183 |
|
IF(YGOOD(I).EQ.0)THEN |
1184 |
|
XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7 |
1185 |
|
ENDIF |
1186 |
|
NP=NP+1 |
1187 |
|
XP(NP)=XX |
1188 |
|
ZP(NP)=REAL(ZM(I))!EM GCC4.7 |
1189 |
|
ENDIF |
1190 |
|
ENDDO |
1191 |
|
IFLAG=0 !no debug mode |
1192 |
|
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
1193 |
|
|
1194 |
|
c$$$ print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
1195 |
|
c$$$ print*,' XP ',(xp(i),i=1,np) |
1196 |
|
c$$$ print*,' ZP ',(zp(i),i=1,np) |
1197 |
|
c$$$ print*,' AP ',(ap(i),i=1,np) |
1198 |
|
c$$$ print*,' XP ',(rp(i),i=1,np) |
1199 |
|
|
1200 |
|
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
1201 |
|
c if(CHI.gt.100)GOTO 10 !straigth fit |
1202 |
|
ARG = REAL(RADIUS**2-(ZINI-ZC)**2)!EM GCC4.7 |
1203 |
|
IF(ARG.LT.0)GOTO 10 !straigth fit |
1204 |
|
DC = SQRT(ARG) |
1205 |
|
IF(XC.GT.0)DC=-DC |
1206 |
|
X0=XC+DC |
1207 |
|
AX = REAL(-(ZINI-ZC)/DC)!EM GCC4.7 |
1208 |
|
DEF=100./(RADIUS*0.3*0.43) |
1209 |
|
IF(XC.GT.0)DEF=-DEF |
1210 |
|
|
1211 |
|
|
1212 |
|
|
1213 |
|
IF(ABS(X0).GT.30)THEN |
1214 |
|
c$$$ PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
1215 |
|
c$$$ $ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
1216 |
|
GOTO 10 !straigth fit |
1217 |
|
ENDIF |
1218 |
|
GOTO 20 !guess is ok |
1219 |
|
|
1220 |
|
* ---------------------------------------- |
1221 |
|
* 2) initial guess with a straigth line |
1222 |
|
* - if circle does not intersect reference plane |
1223 |
|
* - if bad chi**2 |
1224 |
|
* ---------------------------------------- |
1225 |
|
10 CONTINUE |
1226 |
|
SZZ=0. |
1227 |
|
SZX=0. |
1228 |
|
SSX=0. |
1229 |
|
SZ=0. |
1230 |
|
S1=0. |
1231 |
|
DO I=1,nplanes |
1232 |
|
IF(XGOOD(I).EQ.1)THEN |
1233 |
|
XX = REAL(XM(I))!EM GCC4.7 |
1234 |
|
IF(YGOOD(I).EQ.0)THEN |
1235 |
|
XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7 |
1236 |
|
ENDIF |
1237 |
|
SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7 |
1238 |
|
SZX=SZX+REAL(ZM(I)*XX)!EM GCC4.7 |
1239 |
|
SSX=SSX+XX |
1240 |
|
SZ=SZ+REAL(ZM(I))!EM GCC4.7 |
1241 |
|
S1=S1+1. |
1242 |
|
ENDIF |
1243 |
|
ENDDO |
1244 |
|
DET=SZZ*S1-SZ*SZ |
1245 |
|
AX=(SZX*S1-SZ*SSX)/DET |
1246 |
|
BX=(SZZ*SSX-SZX*SZ)/DET |
1247 |
|
DEF = 0 |
1248 |
|
X0 = REAL(AX*ZINI+BX)!EM GCC4.7 |
1249 |
|
|
1250 |
|
20 CONTINUE |
1251 |
|
* ---------------------------------------- |
1252 |
|
* guess |
1253 |
|
* ---------------------------------------- |
1254 |
|
|
1255 |
|
AL(1) = X0 |
1256 |
|
AL(2) = Y0 |
1257 |
|
tath = sqrt(AY**2+AX**2) |
1258 |
|
AL(3) = tath/sqrt(1+tath**2) |
1259 |
|
|
1260 |
|
AL(4)=0. |
1261 |
|
IF( AX.NE.0.OR.AY.NE.0. ) THEN |
1262 |
|
AL(4) = ASIN(AY/SQRT(AX**2+AY**2)) |
1263 |
|
IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4) |
1264 |
|
IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4) |
1265 |
|
ENDIF |
1266 |
|
IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0) |
1267 |
|
IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0) |
1268 |
|
|
1269 |
|
AL(5) = DEF |
1270 |
|
|
1271 |
|
c print*,' guess: ',(al(i),i=1,5) |
1272 |
|
|
1273 |
|
end |