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 |
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 |
52 |
|
|
53 |
c elena-------- |
c elena-------- |
68 |
|
|
69 |
c LOGICAL TRKDEBUG,TRKVERBOSE |
c LOGICAL TRKDEBUG,TRKVERBOSE |
70 |
c COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
c COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
71 |
LOGICAL TRKDEBUG,TRKVERBOSE |
LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT |
72 |
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
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 |
IF(IPRINT.EQ.1) THEN |
90 |
TRKVERBOSE = .TRUE. |
TRKVERBOSE = .TRUE. |
91 |
TRKDEBUG = .FALSE. |
TRKDEBUG = .FALSE. |
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 |
123 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
124 |
* changed in order to evaluate energy-dependent |
* changed in order to evaluate energy-dependent |
125 |
* tolerances on all 5 parameters |
* tolerances on all 5 parameters |
126 |
FACT=100. !scale factor to define tolerance on alfa |
cPP FACT=1.0e10 !scale factor to define tolerance on alfa |
127 |
c deflection error (see PDG) |
c deflection error (see PDG) |
128 |
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.)) |
129 |
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) |
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 |
* ----------------------- |
* ----------------------- |
150 |
* ----------------------- |
* ----------------------- |
151 |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
152 |
|
|
153 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
* ------------------------------- |
154 |
IF(JFAIL.NE.0) THEN |
* **** Chi2+gaussian minimization |
155 |
IFAIL=1 |
* ------------------------------- |
156 |
CHI2=-9999. |
|
157 |
if(TRKVERBOSE) |
IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN |
158 |
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
|
159 |
RETURN |
IF(ISTEP.GE.3) FIRSTSTEPS = .false. |
160 |
ENDIF |
|
161 |
|
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
162 |
COST=1e-5 |
IF(JFAIL.NE.0) THEN |
163 |
DO I=1,5 |
IFAIL=1 |
164 |
DO J=1,5 |
CHI2=-9999. |
165 |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
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 |
175 |
|
DO I=1,5 |
176 |
|
DO J=1,5 |
177 |
|
CHI2DD(I,J)=CHI2DD(I,J)*COST |
178 |
|
ENDDO |
179 |
|
c$$$ CHI2D(I)=CHI2D(I)*COST |
180 |
ENDDO |
ENDDO |
|
CHI2D(I)=CHI2D(I)*COST |
|
|
ENDDO |
|
181 |
|
|
182 |
IF(PFIXED.EQ.0.) THEN |
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(TRKVERBOSE)then |
if(TRKVERBOSE)then |
190 |
PRINT *, |
PRINT *, |
191 |
$ '*** ERROR in mini ***'// |
$ '*** ERROR in mini ***'// |
192 |
$ 'on matrix inversion (not pos-def)' |
$ 'on matrix inversion (not pos-def)' |
193 |
$ ,DET |
$ ,DET |
194 |
endif |
endif |
195 |
IF(CHI2.EQ.0) CHI2=-9999. |
IF(CHI2.EQ.0) CHI2=-9999. |
196 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IF(CHI2.GT.0) CHI2=-CHI2 |
197 |
IFAIL=1 |
IFAIL=1 |
198 |
RETURN |
RETURN |
199 |
ENDIF |
ENDIF |
200 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
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)=2.*COST*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 |
* track fitting with FIXED deflection |
216 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
217 |
ELSE |
ELSE |
218 |
AL(5)=1./PFIXED |
AL(5)=1./PFIXED |
219 |
DO I=1,4 |
DO I=1,4 |
220 |
CHI2D_R(I)=CHI2D(I) |
CHI2D_R(I)=CHI2D(I) |
221 |
DO J=1,4 |
DO J=1,4 |
222 |
CHI2DD_R(I,J)=CHI2DD(I,J) |
CHI2DD_R(I,J)=CHI2DD(I,J) |
223 |
|
ENDDO |
224 |
ENDDO |
ENDDO |
225 |
ENDDO |
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
226 |
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
IF(IFA.NE.0) THEN |
227 |
IF(IFA.NE.0) THEN |
if(TRKVERBOSE)then |
228 |
if(TRKVERBOSE)then |
PRINT *, |
229 |
PRINT *, |
$ '*** ERROR in mini ***'// |
230 |
$ '*** ERROR in mini ***'// |
$ 'on matrix inversion (not pos-def)' |
231 |
$ 'on matrix inversion (not pos-def)' |
$ ,DET |
232 |
$ ,DET |
endif |
233 |
endif |
IF(CHI2.EQ.0) CHI2=-9999. |
234 |
IF(CHI2.EQ.0) CHI2=-9999. |
IF(CHI2.GT.0) CHI2=-CHI2 |
235 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IFAIL=1 |
236 |
IFAIL=1 |
RETURN |
237 |
RETURN |
ENDIF |
238 |
ENDIF |
CALL DSFINV(4,CHI2DD_R,4) |
|
CALL DSFINV(4,CHI2DD_R,4) |
|
239 |
* ******************************************* |
* ******************************************* |
240 |
* find new value of AL-pha |
* find new value of AL-pha |
241 |
* ******************************************* |
* ******************************************* |
242 |
DO I=1,4 |
DO I=1,4 |
243 |
DAL(I)=0. |
DAL(I)=0. |
244 |
DO J=1,4 |
DO J=1,4 |
245 |
DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) |
DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) *COST |
246 |
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
247 |
|
ENDDO |
248 |
ENDDO |
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*,'mini : 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 |
ENDDO |
314 |
DAL(5)=0. |
DO I=1,5 |
315 |
DO I=1,4 |
AL(I)=AL(I)+E*DAL(I) |
|
AL(I)=AL(I)+DAL(I) |
|
316 |
ENDDO |
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 |
ENDIF |
380 |
|
|
381 |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
|
382 |
|
|
383 |
|
|
384 |
|
|
385 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
386 |
* ---------------------------------------------------- * |
* ---------------------------------------------------- * |
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) |
429 |
ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 |
ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 |
430 |
ALTOL(4) = ALTOL(3) |
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 |
c$$$ DO I=1,5 |
436 |
c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! |
c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! |
437 |
c$$$ ENDDO |
c$$$ ENDDO |
438 |
c$$$ print*,'chi2 -- ',DCHI2 |
c$$$ print*,'chi2 -- ',DCHI2 |
439 |
|
|
440 |
IF(ISTEP.LT.3) GOTO 10 ! ***PP*** |
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 |
* new estimate of chi^2: |
***************************** |
446 |
JFAIL=0 !error flag |
* final estimate of chi^2 |
447 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
***************************** |
448 |
IF(JFAIL.NE.0) THEN |
|
449 |
IFAIL=1 |
* ------------------------------- |
450 |
if(TRKVERBOSE)THEN |
* **** Chi2+gaussian minimization |
451 |
CHI2=-9999. |
* ------------------------------- |
452 |
if(TRKVERBOSE) |
|
453 |
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
IF(.NOT.STUDENT) THEN |
454 |
ENDIF |
|
455 |
RETURN |
JFAIL=0 !error flag |
456 |
ENDIF |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
457 |
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 |
|
458 |
IFAIL=1 |
IFAIL=1 |
459 |
RETURN |
if(TRKVERBOSE)THEN |
460 |
|
CHI2=-9999. |
461 |
|
if(TRKVERBOSE) |
462 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
463 |
|
ENDIF |
464 |
|
RETURN |
465 |
ENDIF |
ENDIF |
466 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
c COST=1e-7 |
467 |
DO I=1,5 |
COST=1. |
468 |
DAL(I)=0. |
DO I=1,5 |
469 |
DO J=1,5 |
IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2 |
470 |
COV(I,J)=2.*COST*CHI2DD(I,J) |
ENDDO |
471 |
ENDDO |
DO I=1,5 |
472 |
ENDDO |
DO J=1,5 |
473 |
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) |
|
474 |
ENDDO |
ENDDO |
475 |
ENDDO |
ENDDO |
476 |
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
IF(PFIXED.EQ.0.) THEN |
477 |
IF(IFA.NE.0) THEN |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
478 |
if(TRKVERBOSE)then |
IF(IFA.NE.0) THEN !not positive-defined |
479 |
PRINT *, |
if(TRKVERBOSE)then |
480 |
$ '*** ERROR in mini ***'// |
PRINT *, |
481 |
$ 'on matrix inversion (not pos-def)' |
$ '*** ERROR in mini ***'// |
482 |
$ ,DET |
$ 'on matrix inversion (not pos-def)' |
483 |
endif |
$ ,DET |
484 |
IF(CHI2.EQ.0) CHI2=-9999. |
endif |
485 |
IF(CHI2.GT.0) CHI2=-CHI2 |
IF(CHI2.EQ.0) CHI2=-9999. |
486 |
IFAIL=1 |
IF(CHI2.GT.0) CHI2=-CHI2 |
487 |
RETURN |
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 |
ENDIF |
525 |
CALL DSFINV(4,CHI2DD_R,4) |
|
526 |
DO I=1,4 |
ENDIF |
527 |
DAL(I)=0. |
|
528 |
DO J=1,4 |
* ------------------------------- |
529 |
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
* **** 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 |
ENDDO |
538 |
ENDDO |
ENDDO |
539 |
ENDIF |
ENDIF |
540 |
|
|
541 |
***************************** |
***************************** |
542 |
|
|
543 |
* ------------------------------------ |
* ------------------------------------ |
556 |
$ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' |
$ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' |
557 |
endif |
endif |
558 |
|
|
|
if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5) |
|
|
|
|
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*,'mini : -ok- ',istep,chi2,AL(5) |
567 |
|
|
568 |
NSTEP=ISTEP ! ***PP*** |
NSTEP=ISTEP ! ***PP*** |
569 |
|
|
570 |
|
c$$$ print*,'>>>>> NSTEP = ',NSTEP |
571 |
|
|
572 |
RETURN |
RETURN |
573 |
END |
END |
574 |
|
|
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) |
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 |
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 |
* |
* |
1004 |
c$$$ print*,'POSXY (prima) ',vout |
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(TRKVERBOSE) |
if(TRKVERBOSE) |
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 |
c$$$ print*,'POSXY (dopo) ',vout |
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 |
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 |
1201 |
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
1202 |
c print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
|
1203 |
|
c$$$ print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
1204 |
|
c$$$ print*,' XP ',(xp(i),i=1,np) |
1205 |
|
c$$$ print*,' ZP ',(zp(i),i=1,np) |
1206 |
|
c$$$ print*,' AP ',(ap(i),i=1,np) |
1207 |
|
c$$$ print*,' XP ',(rp(i),i=1,np) |
1208 |
|
|
1209 |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
1210 |
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 |
|
|
1220 |
|
|
1221 |
|
|
1222 |
IF(ABS(X0).GT.30)THEN |
IF(ABS(X0).GT.30)THEN |
1223 |
c$$$ PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
c$$$ PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
1224 |
c$$$ $ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
c$$$ $ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
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 |