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 |
c------------------------------------------------------------------------ |
c------------------------------------------------------------------------ |
27 |
c variables used in the tracking procedure (mini and its subroutines) |
c variables used in the tracking procedure (mini and its subroutines) |
28 |
c |
c |
31 |
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 |
32 |
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 |
33 |
c------------------------------------------------------------------------ |
c------------------------------------------------------------------------ |
34 |
DATA ZINI/23.5/ !z coordinate of the reference plane |
DATA ZINI/23.5/ !!! ***PP*** to be changed !z coordinate of the reference plane |
35 |
|
|
36 |
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 |
37 |
|
|
38 |
DATA STEPAL/5*1.d-7/ !alpha vector step |
DATA STEPAL/5*1.d-7/ !alpha vector step |
39 |
DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization |
DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization |
41 |
* !the tracking procedure |
* !the tracking procedure |
42 |
DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess |
DATA STEPMAX/100./ !maximum number of steps in the trackin gprocess |
43 |
|
|
44 |
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 |
45 |
c DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !" |
c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !" |
46 |
DATA ALMAX/inf,inf,1.,inf,0.25e2/ !limits on alpha vector components |
DATA ALMAX/dinf,dinf,dinf,dinf,dinf/ !limits on alpha vector components |
47 |
DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !" |
DATA ALMIN/-dinf,-dinf,-dinf,-dinf,-dinf/ !" |
|
|
|
48 |
|
|
49 |
DIMENSION DAL(5) !increment of vector alfa |
DIMENSION DAL(5) !increment of vector alfa |
50 |
|
DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2 |
51 |
|
|
52 |
|
c elena-------- |
53 |
|
REAL*8 AVRESX,AVRESY |
54 |
|
c elena-------- |
55 |
|
|
56 |
INTEGER IFLAG |
INTEGER IFLAG |
57 |
c-------------------------------------------------------- |
c-------------------------------------------------------- |
58 |
c IFLAG =1 ---- chi2 derivatives computed by using |
c IFLAG =1 ---- chi2 derivatives computed by using |
63 |
c NB: the two metods gives equivalent results BUT |
c NB: the two metods gives equivalent results BUT |
64 |
c method 2 is faster!! |
c method 2 is faster!! |
65 |
c-------------------------------------------------------- |
c-------------------------------------------------------- |
66 |
DATA IFLAG/2/ |
DATA IFLAG/2/ |
67 |
|
|
68 |
|
c LOGICAL TRKDEBUG,TRKVERBOSE |
69 |
|
c COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
70 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
71 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
72 |
|
|
73 |
|
IF(IPRINT.EQ.1) THEN |
74 |
|
TRKVERBOSE = .TRUE. |
75 |
|
TRKDEBUG = .FALSE. |
76 |
|
ELSEIF(IPRINT.EQ.2)THEN |
77 |
|
TRKVERBOSE = .TRUE. |
78 |
|
TRKDEBUG = .TRUE. |
79 |
|
ELSE |
80 |
|
TRKVERBOSE = .FALSE. |
81 |
|
TRKDEBUG = .FALSE. |
82 |
|
ENDIF |
83 |
|
|
84 |
|
* ---------------------------------------------------------- |
85 |
|
* evaluate average spatial resolution |
86 |
|
* ---------------------------------------------------------- |
87 |
|
AVRESX = RESXAV |
88 |
|
AVRESY = RESYAV |
89 |
|
DO IP=1,6 |
90 |
|
IF( XGOOD(IP).EQ.1 )THEN |
91 |
|
NX=NX+1 |
92 |
|
AVRESX=AVRESX+RESX(IP) |
93 |
|
ENDIF |
94 |
|
IF(NX.NE.0)AVRESX=AVRESX/NX |
95 |
|
IF( YGOOD(IP).EQ.1 )THEN |
96 |
|
NY=NY+1 |
97 |
|
AVRESY=AVRESY+RESY(IP) |
98 |
|
ENDIF |
99 |
|
IF(NX.NE.0)AVRESY=AVRESY/NY |
100 |
|
ENDDO |
101 |
|
|
102 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
103 |
* define ALTOL(5) ---> tolerances on state vector |
* define ALTOL(5) ---> tolerances on state vector |
104 |
* |
* |
105 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
106 |
FACT=10. !scale factor to define |
* changed in order to evaluate energy-dependent |
107 |
!tolerance on alfa |
* tolerances on all 5 parameters |
108 |
ALTOL(1)=RESX(1)/FACT !al(1) = x |
FACT=100. !scale factor to define tolerance on alfa |
|
ALTOL(2)=RESY(1)/FACT !al(2) = y |
|
|
ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta) |
|
|
$ +RESY(1)**2)/44.51/FACT |
|
|
ALTOL(4)=ALTOL(3) !al(4)=phi |
|
109 |
c deflection error (see PDG) |
c deflection error (see PDG) |
110 |
DELETA1=0.01*RESX(1)/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.)) |
111 |
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) |
112 |
|
c$$$ ALTOL(1) = AVRESX/FACT !al(1) = x |
113 |
|
c$$$ ALTOL(2) = AVRESY/FACT !al(2) = y |
114 |
|
c$$$ ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta) |
115 |
|
c$$$ $ +AVRESY**2)/44.51/FACT |
116 |
|
c$$$ ALTOL(4) = ALTOL(3) !al(4)=phi |
117 |
|
c deflection error (see PDG) |
118 |
|
c$$$ DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.)) |
119 |
|
c$$$ DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36) |
120 |
* ---------------------------------------------------------- |
* ---------------------------------------------------------- |
121 |
* |
* |
122 |
ISTEP=0 !num. steps to minimize chi^2 |
ISTEP=0 !num. steps to minimize chi^2 |
123 |
JFAIL=0 !error flag |
JFAIL=0 !error flag |
124 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
|
125 |
IF(JFAIL.NE.0) THEN |
if(TRKDEBUG) print*,'guess: ',al |
126 |
IFAIL=1 |
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
127 |
if(DEBUG) |
|
|
$ PRINT *,'mini_2 ===> error on CHISQ computation !!!' |
|
|
RETURN |
|
|
ENDIF |
|
128 |
* |
* |
129 |
* ----------------------- |
* ----------------------- |
130 |
* START MINIMIZATION LOOP |
* START MINIMIZATION LOOP |
131 |
* ----------------------- |
* ----------------------- |
132 |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !! |
|
CHI2_P=CHI2 |
|
|
|
|
|
c print*,'@@@@@ ',istep,' - ',al |
|
133 |
|
|
134 |
cost=1e-7 |
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
135 |
|
IF(JFAIL.NE.0) THEN |
136 |
|
IFAIL=1 |
137 |
|
CHI2=-9999. |
138 |
|
if(TRKVERBOSE) |
139 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
140 |
|
RETURN |
141 |
|
ENDIF |
142 |
|
|
143 |
|
COST=1e-7 |
144 |
DO I=1,5 |
DO I=1,5 |
145 |
DO J=1,5 |
DO J=1,5 |
146 |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
CHI2DD(I,J)=CHI2DD(I,J)*COST |
147 |
ENDDO |
ENDDO |
148 |
CHI2D(I)=CHI2D(I)*COST |
CHI2D(I)=CHI2D(I)*COST |
149 |
ENDDO |
ENDDO |
150 |
|
|
151 |
|
IF(PFIXED.EQ.0.) THEN |
152 |
|
|
153 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
154 |
* track fitting with FREE deflection |
* track fitting with FREE deflection |
155 |
*------------------------------------------------------------* |
*------------------------------------------------------------* |
156 |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
157 |
IF(IFA.NE.0) THEN !not positive-defined |
IF(IFA.NE.0) THEN !not positive-defined |
158 |
if(DEBUG)then |
if(TRKVERBOSE)then |
159 |
PRINT *, |
PRINT *, |
160 |
$ 'MINI_HOUGH ==> '// |
$ '*** ERROR in mini ***'// |
161 |
$ '** ERROR ** on matrix inversion (not positive-defined)!!!' |
$ 'on matrix inversion (not pos-def)' |
162 |
$ ,DET |
$ ,DET |
163 |
endif |
endif |
164 |
IFAIL=1 |
IF(CHI2.EQ.0) CHI2=-9999. |
165 |
RETURN |
IF(CHI2.GT.0) CHI2=-CHI2 |
166 |
ENDIF |
IFAIL=1 |
167 |
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
RETURN |
168 |
|
ENDIF |
169 |
|
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
170 |
* ******************************************* |
* ******************************************* |
171 |
* find new value of AL-pha !* |
* find new value of AL-pha |
172 |
* !* |
* ******************************************* |
173 |
DO I=1,5 !* |
DO I=1,5 |
174 |
DAL(I)=0. !* |
DAL(I)=0. |
175 |
DO J=1,5 !* |
DO J=1,5 |
176 |
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) !* |
DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) |
177 |
COV(I,J)=CHI2DD(I,J) !* |
COV(I,J)=2.*COST*CHI2DD(I,J) |
178 |
ENDDO !* |
ENDDO |
179 |
ENDDO !* |
ENDDO |
180 |
DO I=1,5 !* |
DO I=1,5 |
181 |
AL(I)=AL(I)+DAL(I) !* |
AL(I)=AL(I)+DAL(I) |
182 |
ENDDO !* |
ENDDO |
183 |
|
*------------------------------------------------------------* |
184 |
|
* track fitting with FIXED deflection |
185 |
|
*------------------------------------------------------------* |
186 |
|
ELSE |
187 |
|
AL(5)=1./PFIXED |
188 |
|
DO I=1,4 |
189 |
|
CHI2D_R(I)=CHI2D(I) |
190 |
|
DO J=1,4 |
191 |
|
CHI2DD_R(I,J)=CHI2DD(I,J) |
192 |
|
ENDDO |
193 |
|
ENDDO |
194 |
|
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
195 |
|
IF(IFA.NE.0) THEN |
196 |
|
if(TRKVERBOSE)then |
197 |
|
PRINT *, |
198 |
|
$ '*** ERROR in mini ***'// |
199 |
|
$ 'on matrix inversion (not pos-def)' |
200 |
|
$ ,DET |
201 |
|
endif |
202 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
203 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
204 |
|
IFAIL=1 |
205 |
|
RETURN |
206 |
|
ENDIF |
207 |
|
CALL DSFINV(4,CHI2DD_R,4) |
208 |
* ******************************************* |
* ******************************************* |
209 |
|
* find new value of AL-pha |
210 |
|
* ******************************************* |
211 |
|
DO I=1,4 |
212 |
|
DAL(I)=0. |
213 |
|
DO J=1,4 |
214 |
|
DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) |
215 |
|
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
216 |
|
ENDDO |
217 |
|
ENDDO |
218 |
|
DAL(5)=0. |
219 |
|
DO I=1,4 |
220 |
|
AL(I)=AL(I)+DAL(I) |
221 |
|
ENDDO |
222 |
|
ENDIF |
223 |
|
|
224 |
|
if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5) |
225 |
|
|
226 |
|
*------------------------------------------------------------* |
227 |
|
* ---------------------------------------------------- * |
228 |
|
*------------------------------------------------------------* |
229 |
* check parameter bounds: |
* check parameter bounds: |
230 |
|
*------------------------------------------------------------* |
231 |
DO I=1,5 |
DO I=1,5 |
232 |
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 |
233 |
if(DEBUG)then |
if(TRKVERBOSE)then |
234 |
PRINT*,' **WARNING** ' |
PRINT*,' *** WARNING in mini *** ' |
235 |
PRINT*,'MINI_2 ==> AL(',I,') out of range' |
PRINT*,'MINI_2 ==> AL(',I,') out of range' |
236 |
PRINT*,' value: ',AL(I), |
PRINT*,' value: ',AL(I), |
237 |
$ ' limits: ',ALMIN(I),ALMAX(I) |
$ ' limits: ',ALMIN(I),ALMAX(I) |
238 |
print*,'istep ',istep |
print*,'istep ',istep |
239 |
endif |
endif |
240 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
241 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
242 |
IFAIL=1 |
IFAIL=1 |
243 |
RETURN |
RETURN |
244 |
ENDIF |
ENDIF |
245 |
ENDDO |
ENDDO |
246 |
* 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 |
|
247 |
* check number of steps: |
* check number of steps: |
248 |
IF(ISTEP.gt.ISTEPMAX) then |
*------------------------------------------------------------* |
249 |
IFAIL=1 |
IF(ISTEP.ge.ISTEPMAX) then |
250 |
if(DEBUG) |
c$$$ IFAIL=1 |
251 |
$ PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX |
c$$$ if(TRKVERBOSE) |
252 |
|
c$$$ $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=', |
253 |
|
c$$$ $ ISTEPMAX |
254 |
goto 11 |
goto 11 |
255 |
endif |
endif |
256 |
|
*------------------------------------------------------------* |
257 |
* --------------------------------------------- |
* --------------------------------------------- |
258 |
* evaluate deflection tolerance on the basis of |
* evaluate deflection tolerance on the basis of |
259 |
* estimated deflection |
* estimated deflection |
260 |
* --------------------------------------------- |
* --------------------------------------------- |
261 |
ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT |
*------------------------------------------------------------* |
262 |
|
c$$$ ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT |
263 |
|
ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT |
264 |
|
ALTOL(1) = ALTOL(5)/DELETA1 |
265 |
|
ALTOL(2) = ALTOL(1) |
266 |
|
ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51 |
267 |
|
ALTOL(4) = ALTOL(3) |
268 |
|
|
269 |
*---- check tolerances: |
*---- check tolerances: |
270 |
|
c$$$ DO I=1,5 |
271 |
|
c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step! |
272 |
|
c$$$ ENDDO |
273 |
|
c$$$ print*,'chi2 -- ',DCHI2 |
274 |
|
|
275 |
DO I=1,5 |
DO I=1,5 |
276 |
IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! |
IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step! |
277 |
ENDDO |
ENDDO |
278 |
|
|
279 |
|
* new estimate of chi^2: |
280 |
|
JFAIL=0 !error flag |
281 |
|
CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives |
282 |
|
IF(JFAIL.NE.0) THEN |
283 |
|
IFAIL=1 |
284 |
|
if(TRKVERBOSE)THEN |
285 |
|
CHI2=-9999. |
286 |
|
if(TRKVERBOSE) |
287 |
|
$ PRINT *,'*** ERROR in mini *** wrong CHISQ' |
288 |
|
ENDIF |
289 |
|
RETURN |
290 |
|
ENDIF |
291 |
|
COST=1e-7 |
292 |
|
DO I=1,5 |
293 |
|
DO J=1,5 |
294 |
|
CHI2DD(I,J)=CHI2DD(I,J)*COST |
295 |
|
ENDDO |
296 |
|
CHI2D(I)=CHI2D(I)*COST |
297 |
|
ENDDO |
298 |
|
IF(PFIXED.EQ.0.) THEN |
299 |
|
CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant |
300 |
|
IF(IFA.NE.0) THEN !not positive-defined |
301 |
|
if(TRKVERBOSE)then |
302 |
|
PRINT *, |
303 |
|
$ '*** ERROR in mini ***'// |
304 |
|
$ 'on matrix inversion (not pos-def)' |
305 |
|
$ ,DET |
306 |
|
endif |
307 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
308 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
309 |
|
IFAIL=1 |
310 |
|
RETURN |
311 |
|
ENDIF |
312 |
|
CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion |
313 |
|
DO I=1,5 |
314 |
|
DAL(I)=0. |
315 |
|
DO J=1,5 |
316 |
|
COV(I,J)=2.*COST*CHI2DD(I,J) |
317 |
|
ENDDO |
318 |
|
ENDDO |
319 |
|
ELSE |
320 |
|
DO I=1,4 |
321 |
|
CHI2D_R(I)=CHI2D(I) |
322 |
|
DO J=1,4 |
323 |
|
CHI2DD_R(I,J)=CHI2DD(I,J) |
324 |
|
ENDDO |
325 |
|
ENDDO |
326 |
|
CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA) |
327 |
|
IF(IFA.NE.0) THEN |
328 |
|
if(TRKVERBOSE)then |
329 |
|
PRINT *, |
330 |
|
$ '*** ERROR in mini ***'// |
331 |
|
$ 'on matrix inversion (not pos-def)' |
332 |
|
$ ,DET |
333 |
|
endif |
334 |
|
IF(CHI2.EQ.0) CHI2=-9999. |
335 |
|
IF(CHI2.GT.0) CHI2=-CHI2 |
336 |
|
IFAIL=1 |
337 |
|
RETURN |
338 |
|
ENDIF |
339 |
|
CALL DSFINV(4,CHI2DD_R,4) |
340 |
|
DO I=1,4 |
341 |
|
DAL(I)=0. |
342 |
|
DO J=1,4 |
343 |
|
COV(I,J)=2.*COST*CHI2DD_R(I,J) |
344 |
|
ENDDO |
345 |
|
ENDDO |
346 |
|
ENDIF |
347 |
|
***************************** |
348 |
|
|
349 |
* ------------------------------------ |
* ------------------------------------ |
350 |
* Number of Degree Of Freedom |
* Number of Degree Of Freedom |
354 |
$ +int(xgood(ip)) |
$ +int(xgood(ip)) |
355 |
$ +int(ygood(ip)) |
$ +int(ygood(ip)) |
356 |
enddo |
enddo |
357 |
ndof=ndof-5 |
if(pfixed.eq.0.) ndof=ndof-5 ! ***PP*** |
358 |
|
if(pfixed.ne.0.) ndof=ndof-4 ! ***PP*** |
359 |
|
if(ndof.le.0.) then |
360 |
|
ndof = 1 |
361 |
|
if(TRKVERBOSE) |
362 |
|
$ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' |
363 |
|
endif |
364 |
|
|
365 |
|
if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5) |
366 |
|
|
367 |
* ------------------------------------ |
* ------------------------------------ |
368 |
* Reduced chi^2 |
* Reduced chi^2 |
369 |
CHI2 = CHI2/dble(ndof) |
CHI2 = CHI2/dble(ndof) |
370 |
|
|
371 |
|
c print*,'mini2: chi2 ',chi2 |
|
11 CONTINUE |
|
372 |
|
|
373 |
101 CONTINUE |
11 CONTINUE |
374 |
|
|
375 |
c print*,'END MINI' |
NSTEP=ISTEP ! ***PP*** |
376 |
|
|
377 |
RETURN |
RETURN |
378 |
END |
END |
399 |
DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes) |
DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes) |
400 |
$ ,XV0(nplanes),YV0(nplanes) |
$ ,XV0(nplanes),YV0(nplanes) |
401 |
DIMENSION AL_P(5) |
DIMENSION AL_P(5) |
402 |
|
|
403 |
|
c LOGICAL TRKVERBOSE |
404 |
|
c COMMON/TRKD/TRKVERBOSE |
405 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
406 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
407 |
* |
* |
408 |
* chi^2 computation |
* chi^2 computation |
409 |
* |
* |
413 |
JFAIL=0 !error flag |
JFAIL=0 !error flag |
414 |
CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes |
CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes |
415 |
IF(JFAIL.NE.0) THEN |
IF(JFAIL.NE.0) THEN |
416 |
PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' |
IF(TRKVERBOSE) |
417 |
|
$ PRINT *,'CHISQ ==> error from trk routine POSXYZ !!' |
418 |
IFAIL=1 |
IFAIL=1 |
419 |
RETURN |
RETURN |
420 |
ENDIF |
ENDIF |
467 |
+ +((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 |
468 |
+ *( (1-XGOOD(I))*YGOOD(I) ) |
+ *( (1-XGOOD(I))*YGOOD(I) ) |
469 |
ENDDO |
ENDDO |
470 |
|
c print*,'CHISQ ',chi2 |
471 |
* ------------------------------------------------ |
* ------------------------------------------------ |
472 |
* |
* |
473 |
* calculation of derivatives (dX/dAL_fa and dY/dAL_fa) |
* calculation of derivatives (dX/dAL_fa and dY/dAL_fa) |
486 |
JFAIL=0 |
JFAIL=0 |
487 |
CALL POSXYZ(AL_P,JFAIL) |
CALL POSXYZ(AL_P,JFAIL) |
488 |
IF(JFAIL.NE.0) THEN |
IF(JFAIL.NE.0) THEN |
489 |
PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' |
IF(TRKVERBOSE) |
490 |
|
*23456789012345678901234567890123456789012345678901234567890123456789012 |
491 |
|
$ PRINT *,'CHISQ ==> error from trk routine POSXYZ' |
492 |
IFAIL=1 |
IFAIL=1 |
493 |
RETURN |
RETURN |
494 |
ENDIF |
ENDIF |
500 |
JFAIL=0 |
JFAIL=0 |
501 |
CALL POSXYZ(AL_P,JFAIL) |
CALL POSXYZ(AL_P,JFAIL) |
502 |
IF(JFAIL.NE.0) THEN |
IF(JFAIL.NE.0) THEN |
503 |
PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!' |
IF(TRKVERBOSE) |
504 |
|
$ PRINT *,'CHISQ ==> error from trk routine POSXYZ' |
505 |
IFAIL=1 |
IFAIL=1 |
506 |
RETURN |
RETURN |
507 |
ENDIF |
ENDIF |
532 |
|
|
533 |
COSTHE=DSQRT(1.-AL(3)**2) |
COSTHE=DSQRT(1.-AL(3)**2) |
534 |
IF(COSTHE.EQ.0.) THEN |
IF(COSTHE.EQ.0.) THEN |
535 |
PRINT *,'=== WARNING ===> COSTHE=0' |
IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0' |
536 |
STOP |
IFAIL=1 |
537 |
|
RETURN |
538 |
ENDIF |
ENDIF |
539 |
|
|
540 |
DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 |
DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3 |
618 |
|
|
619 |
include 'commontracker.f' !tracker general common |
include 'commontracker.f' !tracker general common |
620 |
include 'common_mini_2.f' !common for the tracking procedure |
include 'common_mini_2.f' !common for the tracking procedure |
621 |
|
|
622 |
|
c LOGICAL TRKVERBOSE |
623 |
|
c COMMON/TRKD/TRKVERBOSE |
624 |
|
LOGICAL TRKDEBUG,TRKVERBOSE |
625 |
|
COMMON/TRKD/TRKDEBUG,TRKVERBOSE |
626 |
c |
c |
627 |
DIMENSION AL_P(5) |
DIMENSION AL_P(5) |
628 |
* |
* |
642 |
VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
643 |
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)) |
644 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
645 |
|
|
646 |
|
c$$$ print*,'POSXY ',vout |
647 |
|
|
648 |
DO I=1,nplanes |
DO I=1,nplanes |
649 |
step=vout(3)-zv(i) |
step=vout(3)-zv(i) |
650 |
10 DO J=1,7 |
10 DO J=1,7 |
655 |
CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
656 |
IF(VOUT(3).GT.VECT(3)) THEN |
IF(VOUT(3).GT.VECT(3)) THEN |
657 |
IFAIL=1 |
IFAIL=1 |
658 |
if(WARNING) |
if(TRKVERBOSE) |
659 |
$ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' |
$ PRINT *,'posxy (grkuta): WARNING ===> backward track!!' |
660 |
if(WARNING)print*,'charge',charge |
c$$$ if(.TRUE.)print*,'charge',charge |
661 |
if(WARNING)print*,'vect',vect |
c$$$ if(.TRUE.)print*,'vect',vect |
662 |
if(WARNING)print*,'vout',vout |
c$$$ if(.TRUE.)print*,'vout',vout |
663 |
if(WARNING)print*,'step',step |
c$$$ if(.TRUE.)print*,'step',step |
664 |
|
if(TRKVERBOSE)print*,'charge',charge |
665 |
|
if(TRKVERBOSE)print*,'vect',vect |
666 |
|
if(TRKVERBOSE)print*,'vout',vout |
667 |
|
if(TRKVERBOSE)print*,'step',step |
668 |
RETURN |
RETURN |
669 |
ENDIF |
ENDIF |
670 |
Z=VOUT(3) |
Z=VOUT(3) |
733 |
|
|
734 |
return |
return |
735 |
end |
end |
736 |
|
|
737 |
|
|
738 |
|
*************************************************** |
739 |
|
* * |
740 |
|
* * |
741 |
|
* * |
742 |
|
* * |
743 |
|
* * |
744 |
|
* * |
745 |
|
************************************************** |
746 |
|
|
747 |
|
subroutine guess() |
748 |
|
|
749 |
|
c IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
750 |
|
|
751 |
|
include 'commontracker.f' !tracker general common |
752 |
|
include 'common_mini_2.f' !common for the tracking procedure |
753 |
|
|
754 |
|
REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES) |
755 |
|
REAL*4 CHI,XC,ZC,RADIUS |
756 |
|
* ---------------------------------------- |
757 |
|
* Y view |
758 |
|
* ---------------------------------------- |
759 |
|
* ---------------------------------------- |
760 |
|
* initial guess with a straigth line |
761 |
|
* ---------------------------------------- |
762 |
|
SZZ=0. |
763 |
|
SZY=0. |
764 |
|
SSY=0. |
765 |
|
SZ=0. |
766 |
|
S1=0. |
767 |
|
DO I=1,nplanes |
768 |
|
IF(YGOOD(I).EQ.1)THEN |
769 |
|
YY = YM(I) |
770 |
|
IF(XGOOD(I).EQ.0)THEN |
771 |
|
YY = (YM_A(I) + YM_B(I))/2 |
772 |
|
ENDIF |
773 |
|
SZZ=SZZ+ZM(I)*ZM(I) |
774 |
|
SZY=SZY+ZM(I)*YY |
775 |
|
SSY=SSY+YY |
776 |
|
SZ=SZ+ZM(I) |
777 |
|
S1=S1+1. |
778 |
|
ENDIF |
779 |
|
ENDDO |
780 |
|
DET=SZZ*S1-SZ*SZ |
781 |
|
AY=(SZY*S1-SZ*SSY)/DET |
782 |
|
BY=(SZZ*SSY-SZY*SZ)/DET |
783 |
|
Y0 = AY*ZINI+BY |
784 |
|
* ---------------------------------------- |
785 |
|
* X view |
786 |
|
* ---------------------------------------- |
787 |
|
* ---------------------------------------- |
788 |
|
* 1) initial guess with a circle |
789 |
|
* ---------------------------------------- |
790 |
|
NP=0 |
791 |
|
DO I=1,nplanes |
792 |
|
IF(XGOOD(I).EQ.1)THEN |
793 |
|
XX = XM(I) |
794 |
|
IF(YGOOD(I).EQ.0)THEN |
795 |
|
XX = (XM_A(I) + XM_B(I))/2 |
796 |
|
ENDIF |
797 |
|
NP=NP+1 |
798 |
|
XP(NP)=XX |
799 |
|
ZP(NP)=ZM(I) |
800 |
|
ENDIF |
801 |
|
ENDDO |
802 |
|
IFLAG=0 !no debug mode |
803 |
|
CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG) |
804 |
|
c print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG |
805 |
|
IF(IFLAG.NE.0)GOTO 10 !straigth fit |
806 |
|
if(CHI.gt.100)GOTO 10 !straigth fit |
807 |
|
ARG = RADIUS**2-(ZINI-ZC)**2 |
808 |
|
IF(ARG.LT.0)GOTO 10 !straigth fit |
809 |
|
DC = SQRT(ARG) |
810 |
|
IF(XC.GT.0)DC=-DC |
811 |
|
X0=XC+DC |
812 |
|
AX = -(ZINI-ZC)/DC |
813 |
|
DEF=100./(RADIUS*0.3*0.43) |
814 |
|
IF(XC.GT.0)DEF=-DEF |
815 |
|
|
816 |
|
IF(ABS(X0).GT.30)THEN |
817 |
|
PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS |
818 |
|
$ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF |
819 |
|
GOTO 10 !straigth fit |
820 |
|
ENDIF |
821 |
|
GOTO 20 !guess is ok |
822 |
|
|
823 |
|
* ---------------------------------------- |
824 |
|
* 2) initial guess with a straigth line |
825 |
|
* - if circle does not intersect reference plane |
826 |
|
* - if bad chi**2 |
827 |
|
* ---------------------------------------- |
828 |
|
10 CONTINUE |
829 |
|
SZZ=0. |
830 |
|
SZX=0. |
831 |
|
SSX=0. |
832 |
|
SZ=0. |
833 |
|
S1=0. |
834 |
|
DO I=1,nplanes |
835 |
|
IF(XGOOD(I).EQ.1)THEN |
836 |
|
XX = XM(I) |
837 |
|
IF(YGOOD(I).EQ.0)THEN |
838 |
|
XX = (XM_A(I) + XM_B(I))/2 |
839 |
|
ENDIF |
840 |
|
SZZ=SZZ+ZM(I)*ZM(I) |
841 |
|
SZX=SZX+ZM(I)*XX |
842 |
|
SSX=SSX+XX |
843 |
|
SZ=SZ+ZM(I) |
844 |
|
S1=S1+1. |
845 |
|
ENDIF |
846 |
|
ENDDO |
847 |
|
DET=SZZ*S1-SZ*SZ |
848 |
|
AX=(SZX*S1-SZ*SSX)/DET |
849 |
|
BX=(SZZ*SSX-SZX*SZ)/DET |
850 |
|
DEF = 0 |
851 |
|
X0 = AX*ZINI+BX |
852 |
|
|
853 |
|
20 CONTINUE |
854 |
|
* ---------------------------------------- |
855 |
|
* guess |
856 |
|
* ---------------------------------------- |
857 |
|
|
858 |
|
AL(1) = X0 |
859 |
|
AL(2) = Y0 |
860 |
|
tath = sqrt(AY**2+AX**2) |
861 |
|
AL(3) = tath/sqrt(1+tath**2) |
862 |
|
IF(AX.NE.0)THEN |
863 |
|
AL(4)= atan(AY/AX) |
864 |
|
ELSE |
865 |
|
AL(4) = acos(-1.)/2 |
866 |
|
IF(AY.LT.0)AL(4) = AL(4)+acos(-1.) |
867 |
|
ENDIF |
868 |
|
IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4) |
869 |
|
AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs |
870 |
|
AL(5) = DEF |
871 |
|
|
872 |
|
c print*,' guess: ',(al(i),i=1,5) |
873 |
|
|
874 |
|
end |