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