/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/mini.f
ViewVC logotype

Contents of /DarthVader/TrackerLevel2/src/F77/mini.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Fri May 19 13:15:56 2006 UTC (18 years, 6 months ago) by mocchiut
Branch: DarthVader
CVS Tags: v0r01, start
Changes since 1.1: +0 -0 lines
Imported sources

1 ************************************************************************
2 *
3 * subroutine to evaluate the vector alfa (AL)
4 * which minimizes CHI^2
5 *
6 * - modified from mini.f in order to call differente chi^2 routine.
7 * The new one includes also single clusters: in this case
8 * the residual is defined as the distance between the track and the
9 * segment AB associated to the single cluster.
10 *
11 *
12 ************************************************************************
13
14
15 SUBROUTINE MINI_2(ISTEP,IFAIL)
16
17 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
18
19 include 'commontracker.f' !tracker general common
20 include 'common_mini_2.f' !common for the tracking procedure
21
22 logical DEBUG
23 common/dbg/DEBUG
24
25 parameter (inf=1.e8) !just a huge number...
26 c------------------------------------------------------------------------
27 c variables used in the tracking procedure (mini and its subroutines)
28 c
29 c N.B.: in mini & C. (and in the following block of variables too)
30 c the plane ordering is reversed in respect of normal
31 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
33 c------------------------------------------------------------------------
34 DATA ZINI/23.5/ !z coordinate of the reference plane
35
36 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
39 DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization
40 DATA TOLL/1.d-8/ !tolerance in reaching the next plane during
41 * !the tracking procedure
42 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
45 c DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !"
46 DATA ALMAX/inf,inf,1.,inf,0.25e2/ !limits on alpha vector components
47 DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !"
48
49
50 DIMENSION DAL(5) !increment of vector alfa
51
52 INTEGER IFLAG
53 c--------------------------------------------------------
54 c IFLAG =1 ---- chi2 derivatives computed by using
55 c incremental ratios and posxyz.f
56 c IFLAG =2 ---- the approximation of Golden is used
57 c (see chisq.f)
58 c
59 c NB: the two metods gives equivalent results BUT
60 c method 2 is faster!!
61 c--------------------------------------------------------
62 DATA IFLAG/2/
63
64 * ----------------------------------------------------------
65 * define ALTOL(5) ---> tolerances on state vector
66 *
67 * ----------------------------------------------------------
68 FACT=10. !scale factor to define
69 !tolerance on alfa
70 ALTOL(1)=RESX(1)/FACT !al(1) = x
71 ALTOL(2)=RESY(1)/FACT !al(2) = y
72 ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta)
73 $ +RESY(1)**2)/44.51/FACT
74 ALTOL(4)=ALTOL(3) !al(4)=phi
75 c deflection error (see PDG)
76 DELETA1=0.01*RESX(1)/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
77 DELETA2=0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
78 * ----------------------------------------------------------
79 *
80 ISTEP=0 !num. steps to minimize chi^2
81 JFAIL=0 !error flag
82 CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
83 IF(JFAIL.NE.0) THEN
84 IFAIL=1
85 if(DEBUG)
86 $ PRINT *,'mini_2 ===> error on CHISQ computation !!!'
87 RETURN
88 ENDIF
89 *
90 * -----------------------
91 * START MINIMIZATION LOOP
92 * -----------------------
93 10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !!
94 CHI2_P=CHI2
95
96 c print*,'@@@@@ ',istep,' - ',al
97
98 cost=1e-7
99 DO I=1,5
100 DO J=1,5
101 CHI2DD(I,J)=CHI2DD(I,J)*COST
102 ENDDO
103 CHI2D(I)=CHI2D(I)*COST
104 ENDDO
105 *------------------------------------------------------------*
106 * track fitting with FREE deflection
107 *------------------------------------------------------------*
108 CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
109 IF(IFA.NE.0) THEN !not positive-defined
110 if(DEBUG)then
111 PRINT *,
112 $ 'MINI_HOUGH ==> '//
113 $ '** ERROR ** on matrix inversion (not positive-defined)!!!'
114 $ ,DET
115 endif
116 IFAIL=1
117 RETURN
118 ENDIF
119 CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion
120 * *******************************************
121 * find new value of AL-pha !*
122 * !*
123 DO I=1,5 !*
124 DAL(I)=0. !*
125 DO J=1,5 !*
126 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) !*
127 COV(I,J)=CHI2DD(I,J) !*
128 ENDDO !*
129 ENDDO !*
130 DO I=1,5 !*
131 AL(I)=AL(I)+DAL(I) !*
132 ENDDO !*
133 * *******************************************
134 * check parameter bounds:
135 DO I=1,5
136 IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN
137 if(DEBUG)then
138 PRINT*,' **WARNING** '
139 PRINT*,'MINI_2 ==> AL(',I,') out of range'
140 PRINT*,' value: ',AL(I),
141 $ ' limits: ',ALMIN(I),ALMAX(I)
142 print*,'istep ',istep
143 endif
144 IFAIL=1
145 RETURN
146 ENDIF
147 ENDDO
148 * new estimate of chi^2:
149 JFAIL=0 !error flag
150 CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
151 IF(JFAIL.NE.0) THEN
152 IFAIL=1
153 if(DEBUG)
154 $ PRINT *,'mini_2: ===> error on CHISQ computation !!!'
155 RETURN
156 ENDIF
157 * check number of steps:
158 IF(ISTEP.gt.ISTEPMAX) then
159 IFAIL=1
160 if(DEBUG)
161 $ PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX
162 goto 11
163 endif
164 * ---------------------------------------------
165 * evaluate deflection tolerance on the basis of
166 * estimated deflection
167 * ---------------------------------------------
168 ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
169 *---- check tolerances:
170 DO I=1,5
171 IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
172 ENDDO
173
174
175 * ------------------------------------
176 * Number of Degree Of Freedom
177 ndof=0
178 do ip=1,nplanes
179 ndof=ndof
180 $ +int(xgood(ip))
181 $ +int(ygood(ip))
182 enddo
183 ndof=ndof-5
184 * ------------------------------------
185 * Reduced chi^2
186 CHI2 = CHI2/dble(ndof)
187
188
189 11 CONTINUE
190
191 101 CONTINUE
192
193 c print*,'END MINI'
194
195 RETURN
196 END
197
198 ******************************************************************************
199 *
200 * routine to compute chi^2 and its derivatives
201 *
202 *
203 * (modified in respect to the previous one in order to include
204 * single clusters. In this case the residual is evaluated by
205 * calculating the distance between the track intersection and the
206 * segment AB associated to the single cluster)
207 *
208 ******************************************************************************
209
210 SUBROUTINE CHISQ(IFLAG,IFAIL)
211
212 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
213
214 include 'commontracker.f' !tracker general common
215 include 'common_mini_2.f' !common for the tracking procedure
216
217 DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)
218 $ ,XV0(nplanes),YV0(nplanes)
219 DIMENSION AL_P(5)
220 *
221 * chi^2 computation
222 *
223 DO I=1,5
224 AL_P(I)=AL(I)
225 ENDDO
226 JFAIL=0 !error flag
227 CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes
228 IF(JFAIL.NE.0) THEN
229 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'
230 IFAIL=1
231 RETURN
232 ENDIF
233 DO I=1,nplanes
234 XV0(I)=XV(I)
235 YV0(I)=YV(I)
236 ENDDO
237 * ------------------------------------------------
238 c$$$ CHI2=0.
239 c$$$ DO I=1,nplanes
240 c$$$ CHI2=CHI2
241 c$$$ + +(XV(I)-XM(I))**2/RESX(i)**2 *XGOOD(I)*YGOOD(I)
242 c$$$ + +(YV(I)-YM(I))**2/RESY(i)**2 *YGOOD(I)*XGOOD(I)
243 c$$$ ENDDO
244 * ---------------------------------------------------------
245 * For planes with only a X or Y-cl included, instead of
246 * a X-Y couple, the residual for chi^2 calculation is
247 * evaluated by finding the point x-y, along the segment AB,
248 * closest to the track.
249 * The X or Y coordinate, respectivelly for X and Y-cl, is
250 * then assigned to XM or YM, which is then considered the
251 * measured position of the cluster.
252 * ---------------------------------------------------------
253 CHI2=0.
254 DO I=1,nplanes
255 IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
256 BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
257 ALFA = XM_A(I) - BETA * YM_A(I)
258 YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
259 if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
260 $ YM(I)=dmin1(YM_A(I),YM_B(I))
261 if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
262 $ YM(I)=dmax1(YM_A(I),YM_B(I))
263 XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
264 ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
265 BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
266 ALFA = YM_A(I) - BETA * XM_A(I)
267 XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
268 if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
269 $ XM(I)=dmin1(XM_A(I),XM_B(I))
270 if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
271 $ XM(I)=dmax1(XM_A(I),XM_B(I))
272 YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
273 ENDIF
274 CHI2=CHI2
275 + +(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) )
276 + +(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) )
277 + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
278 + *( XGOOD(I)*(1-YGOOD(I)) )
279 + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
280 + *( (1-XGOOD(I))*YGOOD(I) )
281 ENDDO
282 * ------------------------------------------------
283 *
284 * calculation of derivatives (dX/dAL_fa and dY/dAL_fa)
285 *
286 * //////////////////////////////////////////////////
287 * METHOD 1 -- incremental ratios
288 * //////////////////////////////////////////////////
289
290 IF(IFLAG.EQ.1) THEN
291
292 DO J=1,5
293 DO JJ=1,5
294 AL_P(JJ)=AL(JJ)
295 ENDDO
296 AL_P(J)=AL_P(J)+STEPAL(J)/2.
297 JFAIL=0
298 CALL POSXYZ(AL_P,JFAIL)
299 IF(JFAIL.NE.0) THEN
300 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'
301 IFAIL=1
302 RETURN
303 ENDIF
304 DO I=1,nplanes
305 XV2(I)=XV(I)
306 YV2(I)=YV(I)
307 ENDDO
308 AL_P(J)=AL_P(J)-STEPAL(J)
309 JFAIL=0
310 CALL POSXYZ(AL_P,JFAIL)
311 IF(JFAIL.NE.0) THEN
312 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'
313 IFAIL=1
314 RETURN
315 ENDIF
316 DO I=1,nplanes
317 XV1(I)=XV(I)
318 YV1(I)=YV(I)
319 ENDDO
320 DO I=1,nplanes
321 DXDAL(I,J)=(XV2(I)-XV1(I))/STEPAL(J)
322 DYDAL(I,J)=(YV2(I)-YV1(I))/STEPAL(J)
323 ENDDO
324 ENDDO
325
326 ENDIF
327
328 * //////////////////////////////////////////////////
329 * METHOD 2 -- Bob Golden
330 * //////////////////////////////////////////////////
331
332 IF(IFLAG.EQ.2) THEN
333
334 DO I=1,nplanes
335 DXDAL(I,1)=1.
336 DYDAL(I,1)=0.
337
338 DXDAL(I,2)=0.
339 DYDAL(I,2)=1.
340
341 COSTHE=DSQRT(1.-AL(3)**2)
342 IF(COSTHE.EQ.0.) THEN
343 PRINT *,'=== WARNING ===> COSTHE=0'
344 STOP
345 ENDIF
346
347 DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
348 DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
349
350 DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
351 DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
352
353 IF(AL(5).NE.0.) THEN
354 DXDAL(I,5)=
355 + (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
356 + *DCOS(AL(4))))/AL(5)
357 DYDAL(I,5)=
358 + (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
359 + *DSIN(AL(4))))/AL(5)
360 ELSE
361 DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
362 DYDAL(I,5)=0.
363 ENDIF
364
365 ENDDO
366 ENDIF
367 *
368 * x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
369 * >>> CHI2D evaluation
370 *
371 DO J=1,5
372 CHI2D(J)=0.
373 DO I=1,nplanes
374 CHI2D(J)=CHI2D(J)
375 + +2.*(XV0(I)-XM(I))/RESX(i)**2*DXDAL(I,J) *XGOOD(I)
376 + +2.*(YV0(I)-YM(I))/RESY(i)**2*DYDAL(I,J) *YGOOD(I)
377 ENDDO
378 ENDDO
379 *
380 * >>> CHI2DD evaluation
381 *
382 DO I=1,5
383 DO J=1,5
384 CHI2DD(I,J)=0.
385 DO K=1,nplanes
386 CHI2DD(I,J)=CHI2DD(I,J)
387 + +2.*DXDAL(K,I)*DXDAL(K,J)/RESX(k)**2 *XGOOD(K)
388 + +2.*DYDAL(K,I)*DYDAL(K,J)/RESY(k)**2 *YGOOD(K)
389 ENDDO
390 ENDDO
391 ENDDO
392 * x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
393
394 RETURN
395 END
396
397
398 *****************************************************************
399 *
400 * Routine to compute the track intersection points
401 * on the tracking-system planes, given the track parameters
402 *
403 * The routine is based on GRKUTA, which computes the
404 * trajectory of a charged particle in a magnetic field
405 * by solving the equatins of motion with Runge-Kuta method.
406 *
407 * Variables that have to be assigned when the subroutine
408 * is called are:
409 *
410 * ZM(1,NPLANES) ----> z coordinates of the planes
411 * AL_P(1,5) ----> track-parameter vector
412 *
413 * -----------------------------------------------------------
414 * NB !!!
415 * The routine works properly only if the
416 * planes are numbered in descending order starting from the
417 * reference plane (ZINI)
418 * -----------------------------------------------------------
419 *
420 *****************************************************************
421
422 SUBROUTINE POSXYZ(AL_P,IFAIL)
423
424 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
425
426 include 'commontracker.f' !tracker general common
427 include 'common_mini_2.f' !common for the tracking procedure
428 c
429 DIMENSION AL_P(5)
430 *
431 DO I=1,nplanes
432 ZV(I)=ZM(I) !
433 ENDDO
434 *
435 * set parameters for GRKUTA
436 *
437 IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
438 IF(AL_P(5).EQ.0) CHARGE=1.
439 VOUT(1)=AL_P(1)
440 VOUT(2)=AL_P(2)
441 VOUT(3)=ZINI ! DBLE(Z0)-DBLE(ZSPEC)
442 VOUT(4)=AL_P(3)*DCOS(AL_P(4))
443 VOUT(5)=AL_P(3)*DSIN(AL_P(4))
444 VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
445 IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
446 IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
447 DO I=1,nplanes
448 step=vout(3)-zv(i)
449 10 DO J=1,7
450 VECT(J)=VOUT(J)
451 VECTINI(J)=VOUT(J)
452 ENDDO
453 11 continue
454 CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
455 IF(VOUT(3).GT.VECT(3)) THEN
456 IFAIL=1
457 PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
458 print*,'charge',charge
459 print*,'vect',vect
460 print*,'vout',vout
461 print*,'step',step
462 RETURN
463 ENDIF
464 Z=VOUT(3)
465 IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100
466 IF(Z.GT.ZM(I)+TOLL) GOTO 10
467 IF(Z.LE.ZM(I)-TOLL) THEN
468 STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
469 DO J=1,7
470 VECT(J)=VECTINI(J)
471 ENDDO
472 GOTO 11
473 ENDIF
474
475 * -----------------------------------------------
476 * evaluate track coordinates
477 100 XV(I)=VOUT(1)
478 YV(I)=VOUT(2)
479 ZV(I)=VOUT(3)
480 AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
481 AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
482 * -----------------------------------------------
483
484 ENDDO
485
486 RETURN
487 END
488
489
490
491
492
493 * **********************************************************
494 * Some initialization routines
495 * **********************************************************
496
497 * ----------------------------------------------------------
498 * Routine to initialize COMMON/TRACK/
499 *
500 subroutine track_init
501
502 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
503
504 include 'commontracker.f' !tracker general common
505 include 'common_mini_2.f' !common for the tracking procedure
506 include 'common_mech.f'
507
508 do i=1,5
509 AL(i) = 0.
510 enddo
511
512 do ip=1,NPLANES
513 ZM(IP) = fitz(nplanes-ip+1) !init to mech. position
514 XM(IP) = -100. !0.
515 YM(IP) = -100. !0.
516 XM_A(IP) = -100. !0.
517 YM_A(IP) = -100. !0.
518 c ZM_A(IP) = 0
519 XM_B(IP) = -100. !0.
520 YM_B(IP) = -100. !0.
521 c ZM_B(IP) = 0
522 RESX(IP) = 1000. !3.d-4
523 RESY(IP) = 1000. !12.d-4
524 XGOOD(IP) = 0
525 YGOOD(IP) = 0
526 enddo
527
528 return
529 end

  ViewVC Help
Powered by ViewVC 1.1.23