43 |
|
|
44 |
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 |
45 |
c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" |
c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" |
46 |
DATA ALMAX/dinf,dinf,dinf,dinf,dinf/ !limits on alpha vector components |
DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components |
47 |
DATA ALMIN/-dinf,-dinf,-dinf,-dinf,-dinf/ !" |
DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" |
48 |
|
|
49 |
DIMENSION DAL(5) !increment of vector alfa |
c$$$ DIMENSION DAL(5) !increment of vector alfa |
50 |
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 |
51 |
|
|
52 |
c elena-------- |
c elena-------- |
67 |
|
|
68 |
c LOGICAL TRKDEBUG,TRKVERBOSE |
c LOGICAL TRKDEBUG,TRKVERBOSE |
69 |
c COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
c COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
70 |
LOGICAL TRKDEBUG,TRKVERBOSE |
LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT |
71 |
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
72 |
|
|
73 |
|
DIMENSION AL0(5) |
74 |
|
LOGICAL SUCCESS_NEW,SUCCESS_OLD |
75 |
|
|
76 |
|
c$$$ PRINT*,'==========' ! TEST |
77 |
|
c$$$ PRINT*,'START MINI' ! TEST |
78 |
|
c$$$ PRINT*,'==========' ! TEST |
79 |
|
|
80 |
|
* |
81 |
|
* define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student) |
82 |
|
* |
83 |
|
STUDENT = .false. |
84 |
|
FIRSTSTEPS = .true. |
85 |
|
FIRSTSTUDENT = .true. |
86 |
|
IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true. |
87 |
|
|
88 |
IF(IPRINT.EQ.1) THEN |
IF(IPRINT.EQ.1) THEN |
89 |
TRKVERBOSE = .TRUE. |
TRKVERBOSE = .TRUE. |
90 |
TRKDEBUG = .FALSE. |
TRKDEBUG = .FALSE. |
101 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
102 |
AVRESX = RESXAV |
AVRESX = RESXAV |
103 |
AVRESY = RESYAV |
AVRESY = RESYAV |
104 |
|
NX = 0.0 |
105 |
|
NY = 0.0 |
106 |
DO IP=1,6 |
DO IP=1,6 |
107 |
IF( XGOOD(IP).EQ.1 )THEN |
IF( XGOOD(IP).EQ.1 )THEN |
108 |
NX=NX+1 |
NX=NX+1.0 |
109 |
AVRESX=AVRESX+RESX(IP) |
AVRESX=AVRESX+RESX(IP) |
110 |
ENDIF |
ENDIF |
|
IF(NX.NE.0)AVRESX=AVRESX/NX |
|
111 |
IF( YGOOD(IP).EQ.1 )THEN |
IF( YGOOD(IP).EQ.1 )THEN |
112 |
NY=NY+1 |
NY=NY+1.0 |
113 |
AVRESY=AVRESY+RESY(IP) |
AVRESY=AVRESY+RESY(IP) |
114 |
ENDIF |
ENDIF |
|
IF(NX.NE.0)AVRESY=AVRESY/NY |
|
115 |
ENDDO |
ENDDO |
116 |
|
IF(NX.NE.0.0)AVRESX=AVRESX/NX |
117 |
|
IF(NY.NE.0.0)AVRESY=AVRESY/NY |
118 |
|
|
119 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
120 |
* define ALTOL(5) ---> tolerances on state vector |
* define ALTOL(5) ---> tolerances on state vector |
122 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
123 |
* changed in order to evaluate energy-dependent |
* changed in order to evaluate energy-dependent |
124 |
* tolerances on all 5 parameters |
* tolerances on all 5 parameters |
125 |
FACT=100. !scale factor to define tolerance on alfa |
cPP FACT=1.0e10 !scale factor to define tolerance on alfa |
126 |
c deflection error (see PDG) |
c deflection error (see PDG) |
127 |
DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) |
DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) |
128 |
DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) |
DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) |
138 |
* |
* |
139 |
ISTEP=0 !num. steps to minimize chi^2 |
ISTEP=0 !num. steps to minimize chi^2 |
140 |
JFAIL=0 !error flag |
JFAIL=0 !error flag |
141 |
|
CHI2=0 |
142 |
|
|
143 |
if(TRKDEBUG) print*,'guess: ',al |
if(TRKDEBUG) print*,'guess: ',al |
144 |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
149 |
* ----------------------- |
* ----------------------- |
150 |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
151 |
|
|
152 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
* ------------------------------- |
153 |
IF(JFAIL.NE.0) THEN |
* **** Chi2+gaussian minimization |
154 |
IFAIL=1 |
* ------------------------------- |
155 |
CHI2=-9999. |
|
156 |
if(TRKVERBOSE) |
IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN |
157 |
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
|
158 |
RETURN |
IF(ISTEP.GE.3) FIRSTSTEPS = .false. |
159 |
ENDIF |
|
160 |
|
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
161 |
COST=1e-7 |
IF(JFAIL.NE.0) THEN |
162 |
DO I=1,5 |
IFAIL=1 |
163 |
DO J=1,5 |
CHI2=-9999. |
164 |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
if(TRKVERBOSE) |
165 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
166 |
|
RETURN |
167 |
|
ENDIF |
168 |
|
|
169 |
|
c COST=1e-5 |
170 |
|
COST=1. |
171 |
|
DO I=1,5 |
172 |
|
IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 |
173 |
|
ENDDO |
174 |
|
DO I=1,5 |
175 |
|
DO J=1,5 |
176 |
|
CHI2DD(I,J)=CHI2DD(I,J)*COST |
177 |
|
ENDDO |
178 |
|
c$$$ CHI2D(I)=CHI2D(I)*COST |
179 |
ENDDO |
ENDDO |
|
CHI2D(I)=CHI2D(I)*COST |
|
|
ENDDO |
|
180 |
|
|
181 |
IF(PFIXED.EQ.0.) THEN |
IF(PFIXED.EQ.0.) THEN |
182 |
|
|
183 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
184 |
* track fitting with FREE deflection |
* track fitting with FREE deflection |
185 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
186 |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
187 |
IF(IFA.NE.0) THEN !not positive-defined |
IF(IFA.NE.0) THEN !not positive-defined |
188 |
if(TRKVERBOSE)then |
if(TRKVERBOSE)then |
189 |
PRINT *, |
PRINT *, |
190 |
$ '*** ERROR in mini ***'// |
$ '*** ERROR in mini ***'// |
191 |
$ 'on matrix inversion (not pos-def)' |
$ 'on matrix inversion (not pos-def)' |
192 |
$ ,DET |
$ ,DET |
193 |
endif |
endif |
194 |
IF(CHI2.EQ.0) CHI2=-9999. |
IF(CHI2.EQ.0) CHI2=-9999. |
195 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IF(CHI2.GT.0) CHI2=-CHI2 |
196 |
IFAIL=1 |
IFAIL=1 |
197 |
RETURN |
RETURN |
198 |
ENDIF |
ENDIF |
199 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
200 |
* ******************************************* |
* ******************************************* |
201 |
* find new value of AL-pha |
* find new value of AL-pha |
202 |
* ******************************************* |
* ******************************************* |
203 |
DO I=1,5 |
DO I=1,5 |
204 |
DAL(I)=0. |
DAL(I)=0. |
205 |
DO J=1,5 |
DO J=1,5 |
206 |
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) |
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) *COST |
207 |
COV(I,J)=2.*COST*CHI2DD(I,J) |
COV(I,J)=2.*COST*CHI2DD(I,J) |
208 |
ENDDO |
ENDDO |
209 |
ENDDO |
ENDDO |
210 |
DO I=1,5 |
DO I=1,5 |
211 |
AL(I)=AL(I)+DAL(I) |
AL(I)=AL(I)+DAL(I) |
212 |
ENDDO |
ENDDO |
213 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
214 |
* track fitting with FIXED deflection |
* track fitting with FIXED deflection |
215 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
216 |
ELSE |
ELSE |
217 |
AL(5)=1./PFIXED |
AL(5)=1./PFIXED |
218 |
DO I=1,4 |
DO I=1,4 |
219 |
CHI2D_R(I)=CHI2D(I) |
CHI2D_R(I)=CHI2D(I) |
220 |
DO J=1,4 |
DO J=1,4 |
221 |
CHI2DD_R(I,J)=CHI2DD(I,J) |
CHI2DD_R(I,J)=CHI2DD(I,J) |
222 |
|
ENDDO |
223 |
ENDDO |
ENDDO |
224 |
ENDDO |
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
225 |
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
IF(IFA.NE.0) THEN |
226 |
IF(IFA.NE.0) THEN |
if(TRKVERBOSE)then |
227 |
if(TRKVERBOSE)then |
PRINT *, |
228 |
PRINT *, |
$ '*** ERROR in mini ***'// |
229 |
$ '*** ERROR in mini ***'// |
$ 'on matrix inversion (not pos-def)' |
230 |
$ 'on matrix inversion (not pos-def)' |
$ ,DET |
231 |
$ ,DET |
endif |
232 |
endif |
IF(CHI2.EQ.0) CHI2=-9999. |
233 |
IF(CHI2.EQ.0) CHI2=-9999. |
IF(CHI2.GT.0) CHI2=-CHI2 |
234 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IFAIL=1 |
235 |
IFAIL=1 |
RETURN |
236 |
RETURN |
ENDIF |
237 |
ENDIF |
CALL DSFINV(4,CHI2DD_R,4) |
|
CALL DSFINV(4,CHI2DD_R,4) |
|
238 |
* ******************************************* |
* ******************************************* |
239 |
* find new value of AL-pha |
* find new value of AL-pha |
240 |
* ******************************************* |
* ******************************************* |
241 |
DO I=1,4 |
DO I=1,4 |
242 |
DAL(I)=0. |
DAL(I)=0. |
243 |
DO J=1,4 |
DO J=1,4 |
244 |
DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) |
DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) *COST |
245 |
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
246 |
|
ENDDO |
247 |
ENDDO |
ENDDO |
248 |
|
DAL(5)=0. |
249 |
|
DO I=1,4 |
250 |
|
AL(I)=AL(I)+DAL(I) |
251 |
|
ENDDO |
252 |
|
ENDIF |
253 |
|
|
254 |
|
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
255 |
|
|
256 |
|
c$$$ PRINT*,'DAL ',(DAL(K),K=1,5) |
257 |
|
c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5) |
258 |
|
|
259 |
|
|
260 |
|
ENDIF |
261 |
|
|
262 |
|
* ------------------------------- |
263 |
|
* **** Likelihood+Student minimization |
264 |
|
* ------------------------------- |
265 |
|
|
266 |
|
IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN |
267 |
|
|
268 |
|
IF(FIRSTSTUDENT) THEN |
269 |
|
FIRSTSTUDENT = .false. |
270 |
|
ISTEP = 1 |
271 |
|
ENDIF |
272 |
|
|
273 |
|
CALL CHISQSTT(1,JFAIL) |
274 |
|
DO I=1,5 |
275 |
|
DAL(I)=0. |
276 |
|
DO J=1,5 |
277 |
|
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) |
278 |
|
ENDDO |
279 |
|
ENDDO |
280 |
|
|
281 |
|
DO I=1,5 |
282 |
|
DO j=1,5 |
283 |
|
COV(I,J) = 2.*CHI2DD(I,J) |
284 |
|
ENDDO |
285 |
|
ENDDO |
286 |
|
|
287 |
|
CHI2TOLL = 1.E-3 |
288 |
|
ALPHA = 3.0 |
289 |
|
BETA = -0.4 |
290 |
|
E=1. |
291 |
|
EA=1. |
292 |
|
EB=1. |
293 |
|
EC=1. |
294 |
|
FA=1. |
295 |
|
FB=1. |
296 |
|
FC=1. |
297 |
|
SUCCESS_OLD = .FALSE. |
298 |
|
SUCCESS_NEW = .FALSE. |
299 |
|
|
300 |
|
CALL CHISQSTT(0,JFAIL) |
301 |
|
c$$$ PRINT*,CHI2 |
302 |
|
CHI2_NEW = CHI2 |
303 |
|
FC = CHI2 |
304 |
|
EC = 0. |
305 |
|
|
306 |
|
ICOUNT = 0 |
307 |
|
100 CONTINUE |
308 |
|
ICOUNT = ICOUNT+1 |
309 |
|
|
310 |
|
DO I=1,5 |
311 |
|
AL0(I)=AL(I) |
312 |
ENDDO |
ENDDO |
313 |
DAL(5)=0. |
DO I=1,5 |
314 |
DO I=1,4 |
AL(I)=AL(I)+E*DAL(I) |
|
AL(I)=AL(I)+DAL(I) |
|
315 |
ENDDO |
ENDDO |
316 |
|
CALL CHISQSTT(0,JFAIL) |
317 |
|
CHI2_OLD = CHI2_NEW |
318 |
|
CHI2_NEW = CHI2 |
319 |
|
FA = FB |
320 |
|
FB = FC |
321 |
|
FC = CHI2 |
322 |
|
EA = EB |
323 |
|
EB = EC |
324 |
|
EC = E |
325 |
|
|
326 |
|
c$$$ PRINT*,E,CHI2_NEW |
327 |
|
|
328 |
|
IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success |
329 |
|
IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101 |
330 |
|
SUCCESS_OLD = SUCCESS_NEW |
331 |
|
SUCCESS_NEW = .TRUE. |
332 |
|
E = E*ALPHA |
333 |
|
ELSE ! failure |
334 |
|
SUCCESS_OLD = SUCCESS_NEW |
335 |
|
SUCCESS_NEW = .FALSE. |
336 |
|
CHI2_NEW = CHI2_OLD |
337 |
|
DO I=1,5 |
338 |
|
AL(I)=AL0(I) |
339 |
|
ENDDO |
340 |
|
IF(SUCCESS_OLD) THEN |
341 |
|
DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA) |
342 |
|
IF(DENOM.NE.0.) THEN |
343 |
|
E = EB - 0.5*( (EB-EA)**2*(FB-FC) |
344 |
|
$ - (EB-EC)**2*(FB-FA) ) / DENOM |
345 |
|
ELSE |
346 |
|
E = BETA*E |
347 |
|
ENDIF |
348 |
|
ELSE |
349 |
|
E = BETA*E |
350 |
|
ENDIF |
351 |
|
c$$$ E = BETA*E |
352 |
|
ENDIF |
353 |
|
IF(ICOUNT.GT.20) GOTO 101 |
354 |
|
GOTO 100 |
355 |
|
|
356 |
|
101 CONTINUE |
357 |
|
|
358 |
|
DO I=1,5 |
359 |
|
DAL(I)=E*DAL(I) |
360 |
|
ENDDO |
361 |
|
|
362 |
|
c$$$ print*,' ' |
363 |
|
c$$$ PRINT*,'DAL ',(DAL(K),K=1,5) |
364 |
|
c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5) |
365 |
|
c$$$ print*,'==== CHI2 ====' |
366 |
|
c$$$ print*,chi2 |
367 |
|
c$$$ print*,'==== CHI2d ====' |
368 |
|
c$$$ print*,(chi2d(i),i=1,5) |
369 |
|
c$$$ print*,'==== CHI2dd ====' |
370 |
|
c$$$ do j=1,5 |
371 |
|
c$$$ print*,(chi2dd(j,i),i=1,5) |
372 |
|
c$$$ enddo |
373 |
|
c$$$ print*,'================' |
374 |
|
c$$$ print*,' ' |
375 |
|
|
376 |
|
*========= FIN QUI ============= |
377 |
|
|
378 |
ENDIF |
ENDIF |
379 |
|
|
380 |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
|
381 |
|
|
382 |
|
|
383 |
|
|
384 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
385 |
* ---------------------------------------------------- * |
* ---------------------------------------------------- * |
424 |
ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 |
ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 |
425 |
ALTOL(4) = ALTOL(3) |
ALTOL(4) = ALTOL(3) |
426 |
|
|
427 |
|
c$$$ print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step! |
428 |
|
|
429 |
*---- check tolerances: |
*---- check tolerances: |
430 |
c$$$ DO I=1,5 |
c$$$ DO I=1,5 |
431 |
c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! |
c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! |
432 |
c$$$ ENDDO |
c$$$ ENDDO |
433 |
c$$$ print*,'chi2 -- ',DCHI2 |
c$$$ print*,'chi2 -- ',DCHI2 |
434 |
|
|
435 |
|
IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP*** |
436 |
DO I=1,5 |
DO I=1,5 |
437 |
IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! |
IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! |
438 |
ENDDO |
ENDDO |
439 |
|
|
440 |
* new estimate of chi^2: |
***************************** |
441 |
JFAIL=0 !error flag |
* final estimate of chi^2 |
442 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
***************************** |
443 |
IF(JFAIL.NE.0) THEN |
|
444 |
IFAIL=1 |
* ------------------------------- |
445 |
if(TRKVERBOSE)THEN |
* **** Chi2+gaussian minimization |
446 |
CHI2=-9999. |
* ------------------------------- |
447 |
if(TRKVERBOSE) |
|
448 |
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
IF(.NOT.STUDENT) THEN |
449 |
ENDIF |
|
450 |
RETURN |
JFAIL=0 !error flag |
451 |
ENDIF |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
452 |
COST=1e-7 |
IF(JFAIL.NE.0) THEN |
|
DO I=1,5 |
|
|
DO J=1,5 |
|
|
CHI2DD(I,J)=CHI2DD(I,J)*COST |
|
|
ENDDO |
|
|
CHI2D(I)=CHI2D(I)*COST |
|
|
ENDDO |
|
|
IF(PFIXED.EQ.0.) THEN |
|
|
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
|
|
IF(IFA.NE.0) THEN !not positive-defined |
|
|
if(TRKVERBOSE)then |
|
|
PRINT *, |
|
|
$ '*** ERROR in mini ***'// |
|
|
$ 'on matrix inversion (not pos-def)' |
|
|
$ ,DET |
|
|
endif |
|
|
IF(CHI2.EQ.0) CHI2=-9999. |
|
|
IF(CHI2.GT.0) CHI2=-CHI2 |
|
453 |
IFAIL=1 |
IFAIL=1 |
454 |
RETURN |
if(TRKVERBOSE)THEN |
455 |
|
CHI2=-9999. |
456 |
|
if(TRKVERBOSE) |
457 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
458 |
|
ENDIF |
459 |
|
RETURN |
460 |
ENDIF |
ENDIF |
461 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
c COST=1e-7 |
462 |
DO I=1,5 |
COST=1. |
463 |
DAL(I)=0. |
DO I=1,5 |
464 |
DO J=1,5 |
IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 |
465 |
COV(I,J)=2.*COST*CHI2DD(I,J) |
ENDDO |
466 |
ENDDO |
DO I=1,5 |
467 |
ENDDO |
DO J=1,5 |
468 |
ELSE |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
|
DO I=1,4 |
|
|
CHI2D_R(I)=CHI2D(I) |
|
|
DO J=1,4 |
|
|
CHI2DD_R(I,J)=CHI2DD(I,J) |
|
469 |
ENDDO |
ENDDO |
470 |
ENDDO |
ENDDO |
471 |
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
IF(PFIXED.EQ.0.) THEN |
472 |
IF(IFA.NE.0) THEN |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
473 |
if(TRKVERBOSE)then |
IF(IFA.NE.0) THEN !not positive-defined |
474 |
PRINT *, |
if(TRKVERBOSE)then |
475 |
$ '*** ERROR in mini ***'// |
PRINT *, |
476 |
$ 'on matrix inversion (not pos-def)' |
$ '*** ERROR in mini ***'// |
477 |
$ ,DET |
$ 'on matrix inversion (not pos-def)' |
478 |
endif |
$ ,DET |
479 |
IF(CHI2.EQ.0) CHI2=-9999. |
endif |
480 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IF(CHI2.EQ.0) CHI2=-9999. |
481 |
IFAIL=1 |
IF(CHI2.GT.0) CHI2=-CHI2 |
482 |
RETURN |
IFAIL=1 |
483 |
|
RETURN |
484 |
|
ENDIF |
485 |
|
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
486 |
|
DO I=1,5 |
487 |
|
c$$$ DAL(I)=0. |
488 |
|
DO J=1,5 |
489 |
|
COV(I,J)=2.*COST*CHI2DD(I,J) |
490 |
|
ENDDO |
491 |
|
ENDDO |
492 |
|
ELSE |
493 |
|
DO I=1,4 |
494 |
|
CHI2D_R(I)=CHI2D(I) |
495 |
|
DO J=1,4 |
496 |
|
CHI2DD_R(I,J)=CHI2DD(I,J) |
497 |
|
ENDDO |
498 |
|
ENDDO |
499 |
|
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
500 |
|
IF(IFA.NE.0) THEN |
501 |
|
if(TRKVERBOSE)then |
502 |
|
PRINT *, |
503 |
|
$ '*** ERROR in mini ***'// |
504 |
|
$ 'on matrix inversion (not pos-def)' |
505 |
|
$ ,DET |
506 |
|
endif |
507 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
508 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
509 |
|
IFAIL=1 |
510 |
|
RETURN |
511 |
|
ENDIF |
512 |
|
CALL DSFINV(4,CHI2DD_R,4) |
513 |
|
DO I=1,4 |
514 |
|
c$$$ DAL(I)=0. |
515 |
|
DO J=1,4 |
516 |
|
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
517 |
|
ENDDO |
518 |
|
ENDDO |
519 |
ENDIF |
ENDIF |
520 |
CALL DSFINV(4,CHI2DD_R,4) |
|
521 |
DO I=1,4 |
ENDIF |
522 |
DAL(I)=0. |
|
523 |
DO J=1,4 |
* ------------------------------- |
524 |
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
* **** Likelihood+student minimization |
525 |
|
* ------------------------------- |
526 |
|
|
527 |
|
IF(STUDENT) THEN |
528 |
|
CALL CHISQSTT(1,JFAIL) |
529 |
|
DO I=1,5 |
530 |
|
DO j=1,5 |
531 |
|
COV(I,J) = 2.*CHI2DD(I,J) |
532 |
ENDDO |
ENDDO |
533 |
ENDDO |
ENDDO |
534 |
ENDIF |
ENDIF |
535 |
|
|
536 |
***************************** |
***************************** |
537 |
|
|
538 |
* ------------------------------------ |
* ------------------------------------ |
551 |
$ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' |
$ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' |
552 |
endif |
endif |
553 |
|
|
|
if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5) |
|
|
|
|
554 |
* ------------------------------------ |
* ------------------------------------ |
555 |
* Reduced chi^2 |
* Reduced chi^2 |
556 |
CHI2 = CHI2/dble(ndof) |
CHI2 = CHI2/dble(ndof) |
559 |
|
|
560 |
11 CONTINUE |
11 CONTINUE |
561 |
|
|
562 |
|
if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5) |
563 |
|
|
564 |
NSTEP=ISTEP ! ***PP*** |
NSTEP=ISTEP ! ***PP*** |
565 |
|
|
566 |
|
c$$$ print*,'>>>>> NSTEP = ',NSTEP |
567 |
|
|
568 |
RETURN |
RETURN |
569 |
END |
END |
570 |
|
|
630 |
* measured position of the cluster. |
* measured position of the cluster. |
631 |
* --------------------------------------------------------- |
* --------------------------------------------------------- |
632 |
CHI2=0. |
CHI2=0. |
633 |
DO I=1,nplanes |
DO I=1,nplanes |
634 |
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 |
635 |
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)) |
636 |
ALFA = XM_A(I) - BETA * YM_A(I) |
ALFA = XM_A(I) - BETA * YM_A(I) |
657 |
+ *( XGOOD(I)*(1-YGOOD(I)) ) |
+ *( XGOOD(I)*(1-YGOOD(I)) ) |
658 |
+ +((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 |
659 |
+ *( (1-XGOOD(I))*YGOOD(I) ) |
+ *( (1-XGOOD(I))*YGOOD(I) ) |
660 |
|
c$$$ print*,(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) ) |
661 |
|
c$$$ print*,(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) ) |
662 |
|
c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 |
663 |
|
c$$$ + *( XGOOD(I)*(1-YGOOD(I)) ) |
664 |
|
c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 |
665 |
|
c$$$ + *( (1-XGOOD(I))*YGOOD(I) ) |
666 |
|
c$$$ print*,XV(I),XM(I),XGOOD(I) |
667 |
|
c$$$ print*,YV(I),YM(I),YGOOD(I) |
668 |
ENDDO |
ENDDO |
669 |
c print*,'CHISQ ',chi2 |
c$$$ print*,'CHISQ ',chi2 |
670 |
* ------------------------------------------------ |
* ------------------------------------------------ |
671 |
* |
* |
672 |
* calculation of derivatives (dX/dAL_fa and dY/dAL_fa) |
* calculation of derivatives (dX/dAL_fa and dY/dAL_fa) |
786 |
RETURN |
RETURN |
787 |
END |
END |
788 |
|
|
789 |
|
****************************************************************************** |
790 |
|
* |
791 |
|
* routine to compute Likelihodd+Student and its derivatives |
792 |
|
* |
793 |
|
* (modified in respect to the previous one in order to include |
794 |
|
* single clusters. In this case the residual is evaluated by |
795 |
|
* calculating the distance between the track intersection and the |
796 |
|
* segment AB associated to the single cluster) |
797 |
|
* |
798 |
|
****************************************************************************** |
799 |
|
|
800 |
|
SUBROUTINE CHISQSTT(IFLAG,JFAIL) |
801 |
|
|
802 |
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
803 |
|
|
804 |
|
include 'commontracker.f' !tracker general common |
805 |
|
include 'common_mini_2.f' !common for the tracking procedure |
806 |
|
|
807 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
808 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
809 |
|
|
810 |
|
DIMENSION AL_P(5) |
811 |
|
DIMENSION VECTEMP(5) |
812 |
|
c$$$ DIMENSION U(5) ! BFGS |
813 |
|
|
814 |
|
DO I=1,5 |
815 |
|
AL_P(I)=AL(I) |
816 |
|
ENDDO |
817 |
|
JFAIL=0 !error flag |
818 |
|
CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes |
819 |
|
IF(JFAIL.NE.0) THEN |
820 |
|
IF(TRKVERBOSE) |
821 |
|
$ PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!' |
822 |
|
IFAIL=1 |
823 |
|
RETURN |
824 |
|
ENDIF |
825 |
|
|
826 |
|
DO I=1,nplanes |
827 |
|
DXDAL(I,1)=1. |
828 |
|
DYDAL(I,1)=0. |
829 |
|
DXDAL(I,2)=0. |
830 |
|
DYDAL(I,2)=1. |
831 |
|
COSTHE=DSQRT(1.-AL(3)**2) |
832 |
|
IF(COSTHE.EQ.0.) THEN |
833 |
|
IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' |
834 |
|
IFAIL=1 |
835 |
|
RETURN |
836 |
|
ENDIF |
837 |
|
DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 |
838 |
|
DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3 |
839 |
|
DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE |
840 |
|
DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE |
841 |
|
IF(AL(5).NE.0.) THEN |
842 |
|
DXDAL(I,5)= |
843 |
|
+ (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I)) |
844 |
|
+ *DCOS(AL(4))))/AL(5) |
845 |
|
DYDAL(I,5)= |
846 |
|
+ (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I)) |
847 |
|
+ *DSIN(AL(4))))/AL(5) |
848 |
|
ELSE |
849 |
|
DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 ) |
850 |
|
DYDAL(I,5)=0. |
851 |
|
ENDIF |
852 |
|
ENDDO |
853 |
|
|
854 |
|
IF(IFLAG.EQ.0) THEN ! function calulation |
855 |
|
CHI2=0. |
856 |
|
DO I=1,nplanes |
857 |
|
IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl |
858 |
|
BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I)) |
859 |
|
ALFA = XM_A(I) - BETA * YM_A(I) |
860 |
|
YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2) |
861 |
|
if(YM(I).lt.dmin1(YM_A(I),YM_B(I))) |
862 |
|
$ YM(I)=dmin1(YM_A(I),YM_B(I)) |
863 |
|
if(YM(I).gt.dmax1(YM_A(I),YM_B(I))) |
864 |
|
$ YM(I)=dmax1(YM_A(I),YM_B(I)) |
865 |
|
XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates |
866 |
|
ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl |
867 |
|
BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I)) |
868 |
|
ALFA = YM_A(I) - BETA * XM_A(I) |
869 |
|
XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2) |
870 |
|
if(XM(I).lt.dmin1(XM_A(I),XM_B(I))) |
871 |
|
$ XM(I)=dmin1(XM_A(I),XM_B(I)) |
872 |
|
if(XM(I).gt.dmax1(XM_A(I),XM_B(I))) |
873 |
|
$ XM(I)=dmax1(XM_A(I),XM_B(I)) |
874 |
|
YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates |
875 |
|
ENDIF |
876 |
|
TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/ |
877 |
|
$ (TAILX(I)*RESX(I)**2) ) |
878 |
|
TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/ |
879 |
|
$ (TAILY(I)*RESY(I)**2) ) |
880 |
|
CHI2=CHI2 |
881 |
|
$ +(TAILX(I)+1.0)*TERMX *( XGOOD(I) ) |
882 |
|
$ +(TAILY(I)+1.0)*TERMY *( YGOOD(I) ) |
883 |
|
ENDDO |
884 |
|
ENDIF |
885 |
|
|
886 |
|
IF(IFLAG.EQ.1) THEN ! derivative calulation |
887 |
|
DO I=1,5 |
888 |
|
CHI2DOLD(I)=CHI2D(I) |
889 |
|
ENDDO |
890 |
|
DO J=1,5 |
891 |
|
CHI2D(J)=0. |
892 |
|
DO I=1,nplanes |
893 |
|
CHI2D(J)=CHI2D(J) |
894 |
|
$ +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/ |
895 |
|
$ (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)* |
896 |
|
$ DXDAL(I,J) *XGOOD(I) |
897 |
|
$ +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/ |
898 |
|
$ (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)* |
899 |
|
$ DYDAL(I,J) *YGOOD(I) |
900 |
|
ENDDO |
901 |
|
ENDDO |
902 |
|
DO K=1,5 |
903 |
|
VECTEMP(K)=0. |
904 |
|
DO M=1,5 |
905 |
|
VECTEMP(K) = VECTEMP(K) + |
906 |
|
$ COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M)) |
907 |
|
ENDDO |
908 |
|
ENDDO |
909 |
|
DOWN1 = 0. |
910 |
|
DO K=1,5 |
911 |
|
DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K)) |
912 |
|
ENDDO |
913 |
|
IF(DOWN1.EQ.0.) THEN |
914 |
|
PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0' |
915 |
|
IFAIL=1 |
916 |
|
RETURN |
917 |
|
ENDIF |
918 |
|
DOWN2 = 0. |
919 |
|
DO K=1,5 |
920 |
|
DO M=1,5 |
921 |
|
DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K) |
922 |
|
ENDDO |
923 |
|
ENDDO |
924 |
|
IF(DOWN2.EQ.0.) THEN |
925 |
|
PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0' |
926 |
|
IFAIL=1 |
927 |
|
RETURN |
928 |
|
ENDIF |
929 |
|
c$$$ DO K=1,5 ! BFGS |
930 |
|
c$$$ U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2 |
931 |
|
c$$$ ENDDO |
932 |
|
DO I=1,5 |
933 |
|
DO J=1,5 |
934 |
|
CHI2DD(I,J) = COV(I,J)/2. |
935 |
|
$ +DAL(I)*DAL(J)/DOWN1 |
936 |
|
$ -VECTEMP(I)*VECTEMP(J)/DOWN2 |
937 |
|
c$$$ $ +DOWN2*U(I)*U(J) ! BFGS |
938 |
|
ENDDO |
939 |
|
ENDDO |
940 |
|
ENDIF |
941 |
|
|
942 |
|
RETURN |
943 |
|
END |
944 |
|
|
945 |
***************************************************************** |
***************************************************************** |
946 |
* |
* |
947 |
* Routine to compute the track intersection points |
* Routine to compute the track intersection points |
980 |
c |
c |
981 |
DIMENSION AL_P(5) |
DIMENSION AL_P(5) |
982 |
* |
* |
983 |
DO I=1,nplanes |
cpp DO I=1,nplanes |
984 |
ZV(I)=ZM(I) ! |
cpp ZV(I)=ZM(I) ! |
985 |
ENDDO |
cpp ENDDO |
986 |
* |
* |
987 |
* set parameters for GRKUTA |
* set parameters for GRKUTA |
988 |
* |
* |
997 |
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)) |
998 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
999 |
|
|
1000 |
c$$$ print*,'POSXY ',vout |
c$$$ print*,'POSXY (prima) ',vout |
1001 |
|
|
1002 |
DO I=1,nplanes |
DO I=1,nplanes |
1003 |
step=vout(3)-zv(i) |
c$$$ ipass = 0 ! TEST |
1004 |
|
c$$$ PRINT *,'TRACKING -> START PLANE: ',I ! TEST |
1005 |
|
cPPP step=vout(3)-zm(i) |
1006 |
|
cPP step=(zm(i)-vout(3))/VOUT(6) |
1007 |
10 DO J=1,7 |
10 DO J=1,7 |
1008 |
VECT(J)=VOUT(J) |
VECT(J)=VOUT(J) |
1009 |
VECTINI(J)=VOUT(J) |
VECTINI(J)=VOUT(J) |
1010 |
ENDDO |
ENDDO |
1011 |
|
cPPP step=vect(3)-zm(i) |
1012 |
|
IF(VOUT(6).GE.0.) THEN |
1013 |
|
IFAIL=1 |
1014 |
|
if(TRKVERBOSE) |
1015 |
|
$ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' |
1016 |
|
RETURN |
1017 |
|
ENDIF |
1018 |
|
step=(zm(i)-vect(3))/VOUT(6) |
1019 |
11 continue |
11 continue |
1020 |
CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
1021 |
|
c$$$ ipass = ipass + 1 ! TEST |
1022 |
|
c$$$ PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST |
1023 |
IF(VOUT(3).GT.VECT(3)) THEN |
IF(VOUT(3).GT.VECT(3)) THEN |
1024 |
IFAIL=1 |
IFAIL=1 |
1025 |
if(TRKVERBOSE) |
if(TRKVERBOSE) |
1045 |
GOTO 11 |
GOTO 11 |
1046 |
ENDIF |
ENDIF |
1047 |
|
|
1048 |
|
|
1049 |
* ----------------------------------------------- |
* ----------------------------------------------- |
1050 |
* evaluate track coordinates |
* evaluate track coordinates |
1051 |
100 XV(I)=VOUT(1) |
100 XV(I)=VOUT(1) |
1055 |
AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
1056 |
* ----------------------------------------------- |
* ----------------------------------------------- |
1057 |
|
|
1058 |
|
IF(TRACKMODE.EQ.1) THEN |
1059 |
|
* ----------------------------------------------- |
1060 |
|
* change of energy by bremsstrahlung for electrons |
1061 |
|
VOUT(7) = VOUT(7) * 0.997 !0.9968 |
1062 |
|
* ----------------------------------------------- |
1063 |
|
ENDIF |
1064 |
|
c$$$ PRINT *,'TRACKING -> END' ! TEST |
1065 |
|
|
1066 |
ENDDO |
ENDDO |
1067 |
|
|
1068 |
|
c$$$ print*,'POSXY (dopo) ',vout |
1069 |
|
|
1070 |
|
|
1071 |
RETURN |
RETURN |
1072 |
END |
END |
1073 |
|
|
1108 |
RESY(IP) = 1000. !12.d-4 |
RESY(IP) = 1000. !12.d-4 |
1109 |
XGOOD(IP) = 0 |
XGOOD(IP) = 0 |
1110 |
YGOOD(IP) = 0 |
YGOOD(IP) = 0 |
1111 |
|
DEDXTRK_X(IP) = 0 |
1112 |
|
DEDXTRK_Y(IP) = 0 |
1113 |
|
AXV(IP) = 0 |
1114 |
|
AYV(IP) = 0 |
1115 |
|
XV(IP) = -100 |
1116 |
|
YV(IP) = -100 |
1117 |
enddo |
enddo |
1118 |
|
|
1119 |
return |
return |
1186 |
ENDDO |
ENDDO |
1187 |
IFLAG=0 !no debug mode |
IFLAG=0 !no debug mode |
1188 |
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
1189 |
c print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
|
1190 |
|
c$$$ print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
1191 |
|
c$$$ print*,' XP ',(xp(i),i=1,np) |
1192 |
|
c$$$ print*,' ZP ',(zp(i),i=1,np) |
1193 |
|
c$$$ print*,' AP ',(ap(i),i=1,np) |
1194 |
|
c$$$ print*,' XP ',(rp(i),i=1,np) |
1195 |
|
|
1196 |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
1197 |
if(CHI.gt.100)GOTO 10 !straigth fit |
c if(CHI.gt.100)GOTO 10 !straigth fit |
1198 |
ARG = RADIUS**2-(ZINI-ZC)**2 |
ARG = RADIUS**2-(ZINI-ZC)**2 |
1199 |
IF(ARG.LT.0)GOTO 10 !straigth fit |
IF(ARG.LT.0)GOTO 10 !straigth fit |
1200 |
DC = SQRT(ARG) |
DC = SQRT(ARG) |
1204 |
DEF=100./(RADIUS*0.3*0.43) |
DEF=100./(RADIUS*0.3*0.43) |
1205 |
IF(XC.GT.0)DEF=-DEF |
IF(XC.GT.0)DEF=-DEF |
1206 |
|
|
1207 |
|
|
1208 |
|
|
1209 |
IF(ABS(X0).GT.30)THEN |
IF(ABS(X0).GT.30)THEN |
1210 |
PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
c$$$ PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
1211 |
$ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
c$$$ $ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
1212 |
GOTO 10 !straigth fit |
GOTO 10 !straigth fit |
1213 |
ENDIF |
ENDIF |
1214 |
GOTO 20 !guess is ok |
GOTO 20 !guess is ok |
1252 |
AL(2) = Y0 |
AL(2) = Y0 |
1253 |
tath = sqrt(AY**2+AX**2) |
tath = sqrt(AY**2+AX**2) |
1254 |
AL(3) = tath/sqrt(1+tath**2) |
AL(3) = tath/sqrt(1+tath**2) |
1255 |
IF(AX.NE.0)THEN |
c$$$ IF(AX.NE.0)THEN |
1256 |
AL(4)= atan(AY/AX) |
c$$$ AL(4)= atan(AY/AX) |
1257 |
ELSE |
c$$$ ELSE |
1258 |
AL(4) = acos(-1.)/2 |
c$$$ AL(4) = acos(-1.)/2 |
1259 |
IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) |
c$$$ IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) |
1260 |
|
c$$$ ENDIF |
1261 |
|
c$$$ IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) |
1262 |
|
c$$$ AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys. |
1263 |
|
|
1264 |
|
c$$$ AL(4) = 0. |
1265 |
|
c$$$ IF(AX.NE.0.AND.AY.NE.0)THEN |
1266 |
|
c$$$ AL(4)= atan(AY/AX) |
1267 |
|
c$$$ ELSEIF(AY.EQ.0)THEN |
1268 |
|
c$$$ AL(4) = 0. |
1269 |
|
c$$$ IF(AX.LT.0)AL(4) = AL(4)+acos(-1.) |
1270 |
|
c$$$ ELSEIF(AX.EQ.0)THEN |
1271 |
|
c$$$ AL(4) = acos(-1.)/2 |
1272 |
|
c$$$ IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) |
1273 |
|
c$$$ ENDIF |
1274 |
|
c$$$ IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) |
1275 |
|
c$$$ AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys. |
1276 |
|
|
1277 |
|
c$$$ AL(4)=0. |
1278 |
|
c$$$ IF( AX.NE.0.OR.AY.NE.0. ) THEN |
1279 |
|
c$$$ AL(4) = ASIN(AY/SQRT(AX**2+AY**2)) |
1280 |
|
c$$$ IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4) |
1281 |
|
c$$$ ENDIF |
1282 |
|
|
1283 |
|
AL(4)=0. |
1284 |
|
IF( AX.NE.0.OR.AY.NE.0. ) THEN |
1285 |
|
AL(4) = ASIN(AY/SQRT(AX**2+AY**2)) |
1286 |
|
IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4) |
1287 |
|
IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4) |
1288 |
ENDIF |
ENDIF |
1289 |
IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) |
IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0) |
1290 |
AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs |
IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0) |
1291 |
|
|
1292 |
AL(5) = DEF |
AL(5) = DEF |
1293 |
|
|
1294 |
c print*,' guess: ',(al(i),i=1,5) |
c print*,' guess: ',(al(i),i=1,5) |