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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.20 - (hide annotations) (download)
Tue Jan 15 14:22:45 2008 UTC (16 years, 11 months ago) by pam-fi
Branch: MAIN
CVS Tags: v5r00
Changes since 1.19: +26 -6 lines
tracking optimization

1 mocchiut 1.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 pam-fi 1.3 SUBROUTINE MINI2(ISTEP,IFAIL,IPRINT)
16 mocchiut 1.1
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 pam-fi 1.2 c logical DEBUG
23     c common/dbg/DEBUG
24 mocchiut 1.1
25 pam-fi 1.4 parameter (dinf=1.d15) !just a huge number...
26 mocchiut 1.1 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 pam-fi 1.3 DATA ZINI/23.5/ !!! ***PP*** to be changed !z coordinate of the reference plane
35 mocchiut 1.1
36 pam-fi 1.3 c DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking
37 mocchiut 1.1
38     DATA STEPAL/5*1.d-7/ !alpha vector step
39 pam-fi 1.7 DATA ISTEPMAX/100/ !maximum number of steps in the chi^2 minimization
40 mocchiut 1.1 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 pam-fi 1.8 c DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
45     c DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
46 pam-fi 1.10 DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
47     DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
48 mocchiut 1.1
49 pam-fi 1.17 c$$$ DIMENSION DAL(5) !increment of vector alfa
50 pam-fi 1.3 DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2
51 pam-fi 1.4
52     c elena--------
53     REAL*8 AVRESX,AVRESY
54     c elena--------
55    
56 mocchiut 1.1 INTEGER IFLAG
57     c--------------------------------------------------------
58     c IFLAG =1 ---- chi2 derivatives computed by using
59     c incremental ratios and posxyz.f
60     c IFLAG =2 ---- the approximation of Golden is used
61     c (see chisq.f)
62     c
63     c NB: the two metods gives equivalent results BUT
64     c method 2 is faster!!
65     c--------------------------------------------------------
66 pam-fi 1.3 DATA IFLAG/2/
67    
68 pam-fi 1.4 c LOGICAL TRKDEBUG,TRKVERBOSE
69     c COMMON/TRKD/TRKDEBUG,TRKVERBOSE
70 pam-fi 1.18 LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT
71 pam-fi 1.4 COMMON/TRKD/TRKDEBUG,TRKVERBOSE
72 pam-fi 1.3
73 pam-fi 1.17 DIMENSION AL0(5)
74     LOGICAL SUCCESS_NEW,SUCCESS_OLD
75 pam-fi 1.20
76     c$$$ PRINT*,'==========' ! TEST
77     c$$$ PRINT*,'START MINI' ! TEST
78     c$$$ PRINT*,'==========' ! TEST
79    
80 pam-fi 1.17 *
81     * define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
82     *
83     STUDENT = .false.
84 pam-fi 1.18 FIRSTSTEPS = .true.
85     FIRSTSTUDENT = .true.
86 pam-fi 1.17 IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.
87    
88 pam-fi 1.3 IF(IPRINT.EQ.1) THEN
89 pam-fi 1.4 TRKVERBOSE = .TRUE.
90     TRKDEBUG = .FALSE.
91     ELSEIF(IPRINT.EQ.2)THEN
92     TRKVERBOSE = .TRUE.
93     TRKDEBUG = .TRUE.
94 pam-fi 1.3 ELSE
95 pam-fi 1.4 TRKVERBOSE = .FALSE.
96     TRKDEBUG = .FALSE.
97 pam-fi 1.3 ENDIF
98 mocchiut 1.1
99     * ----------------------------------------------------------
100 pam-fi 1.4 * evaluate average spatial resolution
101     * ----------------------------------------------------------
102     AVRESX = RESXAV
103     AVRESY = RESYAV
104 pam-fi 1.20 NX = 0.0
105     NY = 0.0
106 pam-fi 1.4 DO IP=1,6
107     IF( XGOOD(IP).EQ.1 )THEN
108 pam-fi 1.20 NX=NX+1.0
109 pam-fi 1.4 AVRESX=AVRESX+RESX(IP)
110     ENDIF
111     IF( YGOOD(IP).EQ.1 )THEN
112 pam-fi 1.20 NY=NY+1.0
113 pam-fi 1.4 AVRESY=AVRESY+RESY(IP)
114     ENDIF
115     ENDDO
116 pam-fi 1.20 IF(NX.NE.0.0)AVRESX=AVRESX/NX
117     IF(NY.NE.0.0)AVRESY=AVRESY/NY
118 pam-fi 1.4
119     * ----------------------------------------------------------
120 mocchiut 1.1 * define ALTOL(5) ---> tolerances on state vector
121     *
122     * ----------------------------------------------------------
123 pam-fi 1.4 * changed in order to evaluate energy-dependent
124     * tolerances on all 5 parameters
125 pam-fi 1.14 cPP FACT=1.0e10 !scale factor to define tolerance on alfa
126 mocchiut 1.1 c deflection error (see PDG)
127 pam-fi 1.4 DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
128     DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
129     c$$$ ALTOL(1) = AVRESX/FACT !al(1) = x
130     c$$$ ALTOL(2) = AVRESY/FACT !al(2) = y
131     c$$$ ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta)
132     c$$$ $ +AVRESY**2)/44.51/FACT
133     c$$$ ALTOL(4) = ALTOL(3) !al(4)=phi
134     c deflection error (see PDG)
135     c$$$ DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
136     c$$$ DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
137 mocchiut 1.1 * ----------------------------------------------------------
138     *
139     ISTEP=0 !num. steps to minimize chi^2
140     JFAIL=0 !error flag
141 pam-fi 1.12 CHI2=0
142 pam-fi 1.4
143 pam-fi 1.5 if(TRKDEBUG) print*,'guess: ',al
144 pam-fi 1.4 if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
145    
146 pam-fi 1.3 *
147     * -----------------------
148     * START MINIMIZATION LOOP
149     * -----------------------
150     10 ISTEP=ISTEP+1 !<<<<<<<<<<<<<< NEW STEP !!
151    
152 pam-fi 1.17 * -------------------------------
153     * **** Chi2+gaussian minimization
154     * -------------------------------
155    
156 pam-fi 1.19 IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN
157 pam-fi 1.18
158     IF(ISTEP.GE.3) FIRSTSTEPS = .false.
159 pam-fi 1.17
160     CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
161     IF(JFAIL.NE.0) THEN
162     IFAIL=1
163     CHI2=-9999.
164     if(TRKVERBOSE)
165     $ PRINT *,'*** ERROR in mini *** wrong CHISQ'
166     RETURN
167     ENDIF
168    
169     c COST=1e-5
170     COST=1.
171     DO I=1,5
172     IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
173     ENDDO
174     DO I=1,5
175     DO J=1,5
176     CHI2DD(I,J)=CHI2DD(I,J)*COST
177     ENDDO
178     c$$$ CHI2D(I)=CHI2D(I)*COST
179 mocchiut 1.1 ENDDO
180 pam-fi 1.3
181 pam-fi 1.17 IF(PFIXED.EQ.0.) THEN
182 pam-fi 1.3
183 mocchiut 1.1 *------------------------------------------------------------*
184     * track fitting with FREE deflection
185     *------------------------------------------------------------*
186 pam-fi 1.17 CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
187     IF(IFA.NE.0) THEN !not positive-defined
188     if(TRKVERBOSE)then
189     PRINT *,
190     $ '*** ERROR in mini ***'//
191     $ 'on matrix inversion (not pos-def)'
192     $ ,DET
193     endif
194     IF(CHI2.EQ.0) CHI2=-9999.
195     IF(CHI2.GT.0) CHI2=-CHI2
196     IFAIL=1
197     RETURN
198     ENDIF
199     CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion
200 pam-fi 1.3 * *******************************************
201     * find new value of AL-pha
202 pam-fi 1.4 * *******************************************
203 pam-fi 1.17 DO I=1,5
204     DAL(I)=0.
205     DO J=1,5
206     DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) *COST
207     COV(I,J)=2.*COST*CHI2DD(I,J)
208     ENDDO
209     ENDDO
210     DO I=1,5
211     AL(I)=AL(I)+DAL(I)
212     ENDDO
213 pam-fi 1.3 *------------------------------------------------------------*
214     * track fitting with FIXED deflection
215     *------------------------------------------------------------*
216 pam-fi 1.17 ELSE
217     AL(5)=1./PFIXED
218     DO I=1,4
219     CHI2D_R(I)=CHI2D(I)
220     DO J=1,4
221     CHI2DD_R(I,J)=CHI2DD(I,J)
222     ENDDO
223 pam-fi 1.3 ENDDO
224 pam-fi 1.17 CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
225     IF(IFA.NE.0) THEN
226     if(TRKVERBOSE)then
227     PRINT *,
228     $ '*** ERROR in mini ***'//
229     $ 'on matrix inversion (not pos-def)'
230     $ ,DET
231     endif
232     IF(CHI2.EQ.0) CHI2=-9999.
233     IF(CHI2.GT.0) CHI2=-CHI2
234     IFAIL=1
235     RETURN
236     ENDIF
237     CALL DSFINV(4,CHI2DD_R,4)
238 pam-fi 1.4 * *******************************************
239     * find new value of AL-pha
240     * *******************************************
241 pam-fi 1.17 DO I=1,4
242     DAL(I)=0.
243     DO J=1,4
244     DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J) *COST
245     COV(I,J)=2.*COST*CHI2DD_R(I,J)
246     ENDDO
247     ENDDO
248     DAL(5)=0.
249     DO I=1,4
250     AL(I)=AL(I)+DAL(I)
251 pam-fi 1.3 ENDDO
252 pam-fi 1.17 ENDIF
253    
254     if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
255    
256     c$$$ PRINT*,'DAL ',(DAL(K),K=1,5)
257     c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
258    
259    
260     ENDIF
261    
262     * -------------------------------
263     * **** Likelihood+Student minimization
264     * -------------------------------
265    
266 pam-fi 1.18 IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN
267    
268     IF(FIRSTSTUDENT) THEN
269     FIRSTSTUDENT = .false.
270     ISTEP = 1
271     ENDIF
272    
273 pam-fi 1.17 CALL CHISQSTT(1,JFAIL)
274     DO I=1,5
275     DAL(I)=0.
276     DO J=1,5
277     DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)
278     ENDDO
279     ENDDO
280    
281     DO I=1,5
282     DO j=1,5
283     COV(I,J) = 2.*CHI2DD(I,J)
284     ENDDO
285     ENDDO
286    
287     CHI2TOLL = 1.E-3
288     ALPHA = 3.0
289     BETA = -0.4
290     E=1.
291     EA=1.
292     EB=1.
293     EC=1.
294     FA=1.
295     FB=1.
296     FC=1.
297     SUCCESS_OLD = .FALSE.
298     SUCCESS_NEW = .FALSE.
299    
300     CALL CHISQSTT(0,JFAIL)
301     c$$$ PRINT*,CHI2
302     CHI2_NEW = CHI2
303     FC = CHI2
304     EC = 0.
305    
306 pam-fi 1.19 ICOUNT = 0
307 pam-fi 1.17 100 CONTINUE
308 pam-fi 1.19 ICOUNT = ICOUNT+1
309    
310 pam-fi 1.17 DO I=1,5
311     AL0(I)=AL(I)
312     ENDDO
313     DO I=1,5
314     AL(I)=AL(I)+E*DAL(I)
315 pam-fi 1.3 ENDDO
316 pam-fi 1.17 CALL CHISQSTT(0,JFAIL)
317     CHI2_OLD = CHI2_NEW
318     CHI2_NEW = CHI2
319     FA = FB
320     FB = FC
321     FC = CHI2
322     EA = EB
323     EB = EC
324     EC = E
325    
326     c$$$ PRINT*,E,CHI2_NEW
327    
328     IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success
329     IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101
330     SUCCESS_OLD = SUCCESS_NEW
331     SUCCESS_NEW = .TRUE.
332     E = E*ALPHA
333     ELSE ! failure
334     SUCCESS_OLD = SUCCESS_NEW
335     SUCCESS_NEW = .FALSE.
336     CHI2_NEW = CHI2_OLD
337     DO I=1,5
338     AL(I)=AL0(I)
339     ENDDO
340     IF(SUCCESS_OLD) THEN
341     DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA)
342     IF(DENOM.NE.0.) THEN
343     E = EB - 0.5*( (EB-EA)**2*(FB-FC)
344     $ - (EB-EC)**2*(FB-FA) ) / DENOM
345     ELSE
346     E = BETA*E
347     ENDIF
348     ELSE
349     E = BETA*E
350     ENDIF
351     c$$$ E = BETA*E
352     ENDIF
353 pam-fi 1.19 IF(ICOUNT.GT.20) GOTO 101
354 pam-fi 1.17 GOTO 100
355    
356     101 CONTINUE
357    
358     DO I=1,5
359     DAL(I)=E*DAL(I)
360 pam-fi 1.3 ENDDO
361 pam-fi 1.17
362     c$$$ print*,' '
363     c$$$ PRINT*,'DAL ',(DAL(K),K=1,5)
364     c$$$ PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
365     c$$$ print*,'==== CHI2 ===='
366     c$$$ print*,chi2
367     c$$$ print*,'==== CHI2d ===='
368     c$$$ print*,(chi2d(i),i=1,5)
369     c$$$ print*,'==== CHI2dd ===='
370     c$$$ do j=1,5
371     c$$$ print*,(chi2dd(j,i),i=1,5)
372     c$$$ enddo
373     c$$$ print*,'================'
374     c$$$ print*,' '
375    
376     *========= FIN QUI =============
377    
378 mocchiut 1.1 ENDIF
379 pam-fi 1.4
380 pam-fi 1.17
381    
382    
383 pam-fi 1.4
384 pam-fi 1.3 *------------------------------------------------------------*
385     * ---------------------------------------------------- *
386     *------------------------------------------------------------*
387 mocchiut 1.1 * check parameter bounds:
388 pam-fi 1.4 *------------------------------------------------------------*
389 mocchiut 1.1 DO I=1,5
390     IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN
391 pam-fi 1.4 if(TRKVERBOSE)then
392 pam-fi 1.3 PRINT*,' *** WARNING in mini *** '
393 mocchiut 1.1 PRINT*,'MINI_2 ==> AL(',I,') out of range'
394     PRINT*,' value: ',AL(I),
395     $ ' limits: ',ALMIN(I),ALMAX(I)
396     print*,'istep ',istep
397     endif
398 pam-fi 1.3 IF(CHI2.EQ.0) CHI2=-9999.
399     IF(CHI2.GT.0) CHI2=-CHI2
400 mocchiut 1.1 IFAIL=1
401     RETURN
402     ENDIF
403     ENDDO
404 pam-fi 1.4 *------------------------------------------------------------*
405 mocchiut 1.1 * check number of steps:
406 pam-fi 1.4 *------------------------------------------------------------*
407     IF(ISTEP.ge.ISTEPMAX) then
408 pam-fi 1.7 c$$$ IFAIL=1
409     c$$$ if(TRKVERBOSE)
410     c$$$ $ PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
411     c$$$ $ ISTEPMAX
412 mocchiut 1.1 goto 11
413     endif
414 pam-fi 1.4 *------------------------------------------------------------*
415 mocchiut 1.1 * ---------------------------------------------
416     * evaluate deflection tolerance on the basis of
417     * estimated deflection
418     * ---------------------------------------------
419 pam-fi 1.4 *------------------------------------------------------------*
420     c$$$ ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
421     ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
422     ALTOL(1) = ALTOL(5)/DELETA1
423     ALTOL(2) = ALTOL(1)
424     ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
425     ALTOL(4) = ALTOL(3)
426    
427 pam-fi 1.14 c$$$ print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!
428    
429 mocchiut 1.1 *---- check tolerances:
430 pam-fi 1.4 c$$$ DO I=1,5
431     c$$$ if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
432     c$$$ ENDDO
433     c$$$ print*,'chi2 -- ',DCHI2
434    
435 pam-fi 1.14 IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
436 mocchiut 1.1 DO I=1,5
437     IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
438     ENDDO
439    
440 pam-fi 1.17 *****************************
441     * final estimate of chi^2
442     *****************************
443    
444     * -------------------------------
445     * **** Chi2+gaussian minimization
446     * -------------------------------
447    
448     IF(.NOT.STUDENT) THEN
449    
450     JFAIL=0 !error flag
451     CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
452     IF(JFAIL.NE.0) THEN
453     IFAIL=1
454     if(TRKVERBOSE)THEN
455     CHI2=-9999.
456     if(TRKVERBOSE)
457     $ PRINT *,'*** ERROR in mini *** wrong CHISQ'
458     ENDIF
459     RETURN
460 pam-fi 1.3 ENDIF
461 pam-fi 1.17 c COST=1e-7
462     COST=1.
463     DO I=1,5
464     IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
465 pam-fi 1.3 ENDDO
466 pam-fi 1.17 DO I=1,5
467     DO J=1,5
468     CHI2DD(I,J)=CHI2DD(I,J)*COST
469 pam-fi 1.3 ENDDO
470     ENDDO
471 pam-fi 1.17 IF(PFIXED.EQ.0.) THEN
472     CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
473     IF(IFA.NE.0) THEN !not positive-defined
474     if(TRKVERBOSE)then
475     PRINT *,
476     $ '*** ERROR in mini ***'//
477     $ 'on matrix inversion (not pos-def)'
478     $ ,DET
479     endif
480     IF(CHI2.EQ.0) CHI2=-9999.
481     IF(CHI2.GT.0) CHI2=-CHI2
482     IFAIL=1
483     RETURN
484     ENDIF
485     CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion
486     DO I=1,5
487     c$$$ DAL(I)=0.
488     DO J=1,5
489     COV(I,J)=2.*COST*CHI2DD(I,J)
490     ENDDO
491     ENDDO
492     ELSE
493     DO I=1,4
494     CHI2D_R(I)=CHI2D(I)
495     DO J=1,4
496     CHI2DD_R(I,J)=CHI2DD(I,J)
497     ENDDO
498     ENDDO
499     CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
500     IF(IFA.NE.0) THEN
501     if(TRKVERBOSE)then
502     PRINT *,
503     $ '*** ERROR in mini ***'//
504     $ 'on matrix inversion (not pos-def)'
505     $ ,DET
506     endif
507     IF(CHI2.EQ.0) CHI2=-9999.
508     IF(CHI2.GT.0) CHI2=-CHI2
509     IFAIL=1
510     RETURN
511     ENDIF
512     CALL DSFINV(4,CHI2DD_R,4)
513     DO I=1,4
514     c$$$ DAL(I)=0.
515     DO J=1,4
516     COV(I,J)=2.*COST*CHI2DD_R(I,J)
517     ENDDO
518     ENDDO
519 pam-fi 1.3 ENDIF
520 pam-fi 1.17
521     ENDIF
522    
523     * -------------------------------
524     * **** Likelihood+student minimization
525     * -------------------------------
526    
527     IF(STUDENT) THEN
528     CALL CHISQSTT(1,JFAIL)
529     DO I=1,5
530     DO j=1,5
531     COV(I,J) = 2.*CHI2DD(I,J)
532 pam-fi 1.3 ENDDO
533     ENDDO
534     ENDIF
535 pam-fi 1.17
536 pam-fi 1.3 *****************************
537 mocchiut 1.1
538     * ------------------------------------
539     * Number of Degree Of Freedom
540     ndof=0
541     do ip=1,nplanes
542     ndof=ndof
543     $ +int(xgood(ip))
544     $ +int(ygood(ip))
545     enddo
546 pam-fi 1.3 if(pfixed.eq.0.) ndof=ndof-5 ! ***PP***
547     if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
548     if(ndof.le.0.) then
549     ndof = 1
550 pam-fi 1.4 if(TRKVERBOSE)
551 pam-fi 1.3 $ print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
552     endif
553 pam-fi 1.4
554 mocchiut 1.1 * ------------------------------------
555     * Reduced chi^2
556     CHI2 = CHI2/dble(ndof)
557    
558 pam-fi 1.4 c print*,'mini2: chi2 ',chi2
559    
560 mocchiut 1.1 11 CONTINUE
561    
562 pam-fi 1.14 if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
563    
564 pam-fi 1.3 NSTEP=ISTEP ! ***PP***
565 mocchiut 1.1
566 pam-fi 1.14 c$$$ print*,'>>>>> NSTEP = ',NSTEP
567    
568 mocchiut 1.1 RETURN
569     END
570    
571     ******************************************************************************
572     *
573     * routine to compute chi^2 and its derivatives
574     *
575     *
576     * (modified in respect to the previous one in order to include
577     * single clusters. In this case the residual is evaluated by
578     * calculating the distance between the track intersection and the
579     * segment AB associated to the single cluster)
580     *
581     ******************************************************************************
582    
583     SUBROUTINE CHISQ(IFLAG,IFAIL)
584    
585     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
586    
587     include 'commontracker.f' !tracker general common
588     include 'common_mini_2.f' !common for the tracking procedure
589    
590     DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)
591     $ ,XV0(nplanes),YV0(nplanes)
592     DIMENSION AL_P(5)
593 pam-fi 1.3
594 pam-fi 1.4 c LOGICAL TRKVERBOSE
595     c COMMON/TRKD/TRKVERBOSE
596     LOGICAL TRKDEBUG,TRKVERBOSE
597     COMMON/TRKD/TRKDEBUG,TRKVERBOSE
598 mocchiut 1.1 *
599     * chi^2 computation
600     *
601     DO I=1,5
602     AL_P(I)=AL(I)
603     ENDDO
604     JFAIL=0 !error flag
605     CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes
606     IF(JFAIL.NE.0) THEN
607 pam-fi 1.4 IF(TRKVERBOSE)
608 pam-fi 1.3 $ PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
609 mocchiut 1.1 IFAIL=1
610     RETURN
611     ENDIF
612     DO I=1,nplanes
613     XV0(I)=XV(I)
614     YV0(I)=YV(I)
615     ENDDO
616     * ------------------------------------------------
617     c$$$ CHI2=0.
618     c$$$ DO I=1,nplanes
619     c$$$ CHI2=CHI2
620     c$$$ + +(XV(I)-XM(I))**2/RESX(i)**2 *XGOOD(I)*YGOOD(I)
621     c$$$ + +(YV(I)-YM(I))**2/RESY(i)**2 *YGOOD(I)*XGOOD(I)
622     c$$$ ENDDO
623     * ---------------------------------------------------------
624     * For planes with only a X or Y-cl included, instead of
625     * a X-Y couple, the residual for chi^2 calculation is
626     * evaluated by finding the point x-y, along the segment AB,
627     * closest to the track.
628     * The X or Y coordinate, respectivelly for X and Y-cl, is
629     * then assigned to XM or YM, which is then considered the
630     * measured position of the cluster.
631     * ---------------------------------------------------------
632     CHI2=0.
633 pam-fi 1.16 DO I=1,nplanes
634 mocchiut 1.1 IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
635     BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
636     ALFA = XM_A(I) - BETA * YM_A(I)
637     YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
638     if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
639     $ YM(I)=dmin1(YM_A(I),YM_B(I))
640     if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
641     $ YM(I)=dmax1(YM_A(I),YM_B(I))
642     XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
643     ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
644     BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
645     ALFA = YM_A(I) - BETA * XM_A(I)
646     XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
647     if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
648     $ XM(I)=dmin1(XM_A(I),XM_B(I))
649     if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
650     $ XM(I)=dmax1(XM_A(I),XM_B(I))
651     YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
652     ENDIF
653     CHI2=CHI2
654     + +(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) )
655     + +(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) )
656     + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
657     + *( XGOOD(I)*(1-YGOOD(I)) )
658     + +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
659     + *( (1-XGOOD(I))*YGOOD(I) )
660 pam-fi 1.10 c$$$ print*,(XV(I)-XM(I))**2/RESX(i)**2 *( XGOOD(I)*YGOOD(I) )
661     c$$$ print*,(YV(I)-YM(I))**2/RESY(i)**2 *( YGOOD(I)*XGOOD(I) )
662     c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
663     c$$$ + *( XGOOD(I)*(1-YGOOD(I)) )
664     c$$$ print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
665     c$$$ + *( (1-XGOOD(I))*YGOOD(I) )
666     c$$$ print*,XV(I),XM(I),XGOOD(I)
667     c$$$ print*,YV(I),YM(I),YGOOD(I)
668 mocchiut 1.1 ENDDO
669 pam-fi 1.10 c$$$ print*,'CHISQ ',chi2
670 mocchiut 1.1 * ------------------------------------------------
671     *
672     * calculation of derivatives (dX/dAL_fa and dY/dAL_fa)
673     *
674     * //////////////////////////////////////////////////
675     * METHOD 1 -- incremental ratios
676     * //////////////////////////////////////////////////
677    
678     IF(IFLAG.EQ.1) THEN
679    
680     DO J=1,5
681     DO JJ=1,5
682     AL_P(JJ)=AL(JJ)
683     ENDDO
684     AL_P(J)=AL_P(J)+STEPAL(J)/2.
685     JFAIL=0
686     CALL POSXYZ(AL_P,JFAIL)
687     IF(JFAIL.NE.0) THEN
688 pam-fi 1.4 IF(TRKVERBOSE)
689 pam-fi 1.3 *23456789012345678901234567890123456789012345678901234567890123456789012
690     $ PRINT *,'CHISQ ==> error from trk routine POSXYZ'
691 mocchiut 1.1 IFAIL=1
692     RETURN
693     ENDIF
694     DO I=1,nplanes
695     XV2(I)=XV(I)
696     YV2(I)=YV(I)
697     ENDDO
698     AL_P(J)=AL_P(J)-STEPAL(J)
699     JFAIL=0
700     CALL POSXYZ(AL_P,JFAIL)
701     IF(JFAIL.NE.0) THEN
702 pam-fi 1.4 IF(TRKVERBOSE)
703 pam-fi 1.3 $ PRINT *,'CHISQ ==> error from trk routine POSXYZ'
704 mocchiut 1.1 IFAIL=1
705     RETURN
706     ENDIF
707     DO I=1,nplanes
708     XV1(I)=XV(I)
709     YV1(I)=YV(I)
710     ENDDO
711     DO I=1,nplanes
712     DXDAL(I,J)=(XV2(I)-XV1(I))/STEPAL(J)
713     DYDAL(I,J)=(YV2(I)-YV1(I))/STEPAL(J)
714     ENDDO
715     ENDDO
716    
717     ENDIF
718    
719     * //////////////////////////////////////////////////
720     * METHOD 2 -- Bob Golden
721     * //////////////////////////////////////////////////
722    
723     IF(IFLAG.EQ.2) THEN
724    
725     DO I=1,nplanes
726     DXDAL(I,1)=1.
727     DYDAL(I,1)=0.
728    
729     DXDAL(I,2)=0.
730     DYDAL(I,2)=1.
731    
732     COSTHE=DSQRT(1.-AL(3)**2)
733     IF(COSTHE.EQ.0.) THEN
734 pam-fi 1.4 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
735 pam-fi 1.3 IFAIL=1
736     RETURN
737 mocchiut 1.1 ENDIF
738    
739     DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
740     DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
741    
742     DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
743     DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
744    
745     IF(AL(5).NE.0.) THEN
746     DXDAL(I,5)=
747     + (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
748     + *DCOS(AL(4))))/AL(5)
749     DYDAL(I,5)=
750     + (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
751     + *DSIN(AL(4))))/AL(5)
752     ELSE
753     DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
754     DYDAL(I,5)=0.
755     ENDIF
756    
757     ENDDO
758     ENDIF
759     *
760     * 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
761     * >>> CHI2D evaluation
762     *
763     DO J=1,5
764     CHI2D(J)=0.
765     DO I=1,nplanes
766     CHI2D(J)=CHI2D(J)
767     + +2.*(XV0(I)-XM(I))/RESX(i)**2*DXDAL(I,J) *XGOOD(I)
768     + +2.*(YV0(I)-YM(I))/RESY(i)**2*DYDAL(I,J) *YGOOD(I)
769     ENDDO
770     ENDDO
771     *
772     * >>> CHI2DD evaluation
773     *
774     DO I=1,5
775     DO J=1,5
776     CHI2DD(I,J)=0.
777     DO K=1,nplanes
778     CHI2DD(I,J)=CHI2DD(I,J)
779     + +2.*DXDAL(K,I)*DXDAL(K,J)/RESX(k)**2 *XGOOD(K)
780     + +2.*DYDAL(K,I)*DYDAL(K,J)/RESY(k)**2 *YGOOD(K)
781     ENDDO
782     ENDDO
783     ENDDO
784     * 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
785    
786     RETURN
787     END
788    
789 pam-fi 1.17 ******************************************************************************
790     *
791     * routine to compute Likelihodd+Student and its derivatives
792     *
793     * (modified in respect to the previous one in order to include
794     * single clusters. In this case the residual is evaluated by
795     * calculating the distance between the track intersection and the
796     * segment AB associated to the single cluster)
797     *
798     ******************************************************************************
799    
800     SUBROUTINE CHISQSTT(IFLAG,JFAIL)
801    
802     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
803    
804     include 'commontracker.f' !tracker general common
805     include 'common_mini_2.f' !common for the tracking procedure
806    
807     LOGICAL TRKDEBUG,TRKVERBOSE
808     COMMON/TRKD/TRKDEBUG,TRKVERBOSE
809    
810     DIMENSION AL_P(5)
811     DIMENSION VECTEMP(5)
812     c$$$ DIMENSION U(5) ! BFGS
813    
814     DO I=1,5
815     AL_P(I)=AL(I)
816     ENDDO
817     JFAIL=0 !error flag
818     CALL POSXYZ(AL_P,JFAIL) !track intersection with tracking planes
819     IF(JFAIL.NE.0) THEN
820     IF(TRKVERBOSE)
821     $ PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!'
822     IFAIL=1
823     RETURN
824     ENDIF
825    
826     DO I=1,nplanes
827     DXDAL(I,1)=1.
828     DYDAL(I,1)=0.
829     DXDAL(I,2)=0.
830     DYDAL(I,2)=1.
831     COSTHE=DSQRT(1.-AL(3)**2)
832     IF(COSTHE.EQ.0.) THEN
833     IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
834     IFAIL=1
835     RETURN
836     ENDIF
837     DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
838     DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
839     DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
840     DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
841     IF(AL(5).NE.0.) THEN
842     DXDAL(I,5)=
843     + (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
844     + *DCOS(AL(4))))/AL(5)
845     DYDAL(I,5)=
846     + (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
847     + *DSIN(AL(4))))/AL(5)
848     ELSE
849     DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
850     DYDAL(I,5)=0.
851     ENDIF
852     ENDDO
853    
854     IF(IFLAG.EQ.0) THEN ! function calulation
855     CHI2=0.
856     DO I=1,nplanes
857     IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
858     BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
859     ALFA = XM_A(I) - BETA * YM_A(I)
860     YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
861     if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
862     $ YM(I)=dmin1(YM_A(I),YM_B(I))
863     if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
864     $ YM(I)=dmax1(YM_A(I),YM_B(I))
865     XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
866     ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
867     BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
868     ALFA = YM_A(I) - BETA * XM_A(I)
869     XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
870     if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
871     $ XM(I)=dmin1(XM_A(I),XM_B(I))
872     if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
873     $ XM(I)=dmax1(XM_A(I),XM_B(I))
874     YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
875     ENDIF
876     TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/
877     $ (TAILX(I)*RESX(I)**2) )
878     TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/
879     $ (TAILY(I)*RESY(I)**2) )
880     CHI2=CHI2
881     $ +(TAILX(I)+1.0)*TERMX *( XGOOD(I) )
882     $ +(TAILY(I)+1.0)*TERMY *( YGOOD(I) )
883     ENDDO
884     ENDIF
885    
886     IF(IFLAG.EQ.1) THEN ! derivative calulation
887     DO I=1,5
888     CHI2DOLD(I)=CHI2D(I)
889     ENDDO
890     DO J=1,5
891     CHI2D(J)=0.
892     DO I=1,nplanes
893     CHI2D(J)=CHI2D(J)
894     $ +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/
895     $ (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)*
896     $ DXDAL(I,J) *XGOOD(I)
897     $ +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/
898     $ (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)*
899     $ DYDAL(I,J) *YGOOD(I)
900     ENDDO
901     ENDDO
902     DO K=1,5
903     VECTEMP(K)=0.
904     DO M=1,5
905     VECTEMP(K) = VECTEMP(K) +
906     $ COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M))
907     ENDDO
908     ENDDO
909     DOWN1 = 0.
910     DO K=1,5
911     DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K))
912     ENDDO
913     IF(DOWN1.EQ.0.) THEN
914     PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0'
915     IFAIL=1
916     RETURN
917     ENDIF
918     DOWN2 = 0.
919     DO K=1,5
920     DO M=1,5
921     DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K)
922     ENDDO
923     ENDDO
924     IF(DOWN2.EQ.0.) THEN
925     PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0'
926     IFAIL=1
927     RETURN
928     ENDIF
929     c$$$ DO K=1,5 ! BFGS
930     c$$$ U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2
931     c$$$ ENDDO
932     DO I=1,5
933     DO J=1,5
934     CHI2DD(I,J) = COV(I,J)/2.
935     $ +DAL(I)*DAL(J)/DOWN1
936     $ -VECTEMP(I)*VECTEMP(J)/DOWN2
937     c$$$ $ +DOWN2*U(I)*U(J) ! BFGS
938     ENDDO
939     ENDDO
940     ENDIF
941 mocchiut 1.1
942 pam-fi 1.17 RETURN
943     END
944    
945 mocchiut 1.1 *****************************************************************
946     *
947     * Routine to compute the track intersection points
948     * on the tracking-system planes, given the track parameters
949     *
950     * The routine is based on GRKUTA, which computes the
951     * trajectory of a charged particle in a magnetic field
952     * by solving the equatins of motion with Runge-Kuta method.
953     *
954     * Variables that have to be assigned when the subroutine
955     * is called are:
956     *
957     * ZM(1,NPLANES) ----> z coordinates of the planes
958     * AL_P(1,5) ----> track-parameter vector
959     *
960     * -----------------------------------------------------------
961     * NB !!!
962     * The routine works properly only if the
963     * planes are numbered in descending order starting from the
964     * reference plane (ZINI)
965     * -----------------------------------------------------------
966     *
967     *****************************************************************
968    
969     SUBROUTINE POSXYZ(AL_P,IFAIL)
970    
971     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
972    
973     include 'commontracker.f' !tracker general common
974     include 'common_mini_2.f' !common for the tracking procedure
975 pam-fi 1.3
976 pam-fi 1.4 c LOGICAL TRKVERBOSE
977     c COMMON/TRKD/TRKVERBOSE
978     LOGICAL TRKDEBUG,TRKVERBOSE
979     COMMON/TRKD/TRKDEBUG,TRKVERBOSE
980 mocchiut 1.1 c
981     DIMENSION AL_P(5)
982     *
983 pam-fi 1.14 cpp DO I=1,nplanes
984     cpp ZV(I)=ZM(I) !
985     cpp ENDDO
986 mocchiut 1.1 *
987     * set parameters for GRKUTA
988     *
989     IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
990     IF(AL_P(5).EQ.0) CHARGE=1.
991     VOUT(1)=AL_P(1)
992     VOUT(2)=AL_P(2)
993     VOUT(3)=ZINI ! DBLE(Z0)-DBLE(ZSPEC)
994     VOUT(4)=AL_P(3)*DCOS(AL_P(4))
995     VOUT(5)=AL_P(3)*DSIN(AL_P(4))
996     VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
997     IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
998     IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
999 pam-fi 1.5
1000 pam-fi 1.10 c$$$ print*,'POSXY (prima) ',vout
1001 pam-fi 1.5
1002 mocchiut 1.1 DO I=1,nplanes
1003 pam-fi 1.20 c$$$ ipass = 0 ! TEST
1004     c$$$ PRINT *,'TRACKING -> START PLANE: ',I ! TEST
1005     cPPP step=vout(3)-zm(i)
1006     cPP step=(zm(i)-vout(3))/VOUT(6)
1007 mocchiut 1.1 10 DO J=1,7
1008     VECT(J)=VOUT(J)
1009     VECTINI(J)=VOUT(J)
1010     ENDDO
1011 pam-fi 1.20 cPPP step=vect(3)-zm(i)
1012     IF(VOUT(6).GE.0.) THEN
1013     IFAIL=1
1014     if(TRKVERBOSE)
1015     $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1016     RETURN
1017     ENDIF
1018     step=(zm(i)-vect(3))/VOUT(6)
1019 mocchiut 1.1 11 continue
1020     CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1021 pam-fi 1.20 c$$$ ipass = ipass + 1 ! TEST
1022     c$$$ PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST
1023 mocchiut 1.1 IF(VOUT(3).GT.VECT(3)) THEN
1024     IFAIL=1
1025 pam-fi 1.4 if(TRKVERBOSE)
1026 pam-fi 1.2 $ PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1027 pam-fi 1.4 c$$$ if(.TRUE.)print*,'charge',charge
1028     c$$$ if(.TRUE.)print*,'vect',vect
1029     c$$$ if(.TRUE.)print*,'vout',vout
1030     c$$$ if(.TRUE.)print*,'step',step
1031     if(TRKVERBOSE)print*,'charge',charge
1032     if(TRKVERBOSE)print*,'vect',vect
1033     if(TRKVERBOSE)print*,'vout',vout
1034     if(TRKVERBOSE)print*,'step',step
1035 mocchiut 1.1 RETURN
1036     ENDIF
1037     Z=VOUT(3)
1038     IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100
1039     IF(Z.GT.ZM(I)+TOLL) GOTO 10
1040     IF(Z.LE.ZM(I)-TOLL) THEN
1041     STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
1042     DO J=1,7
1043     VECT(J)=VECTINI(J)
1044     ENDDO
1045     GOTO 11
1046     ENDIF
1047    
1048 pam-fi 1.10
1049 mocchiut 1.1 * -----------------------------------------------
1050     * evaluate track coordinates
1051     100 XV(I)=VOUT(1)
1052     YV(I)=VOUT(2)
1053     ZV(I)=VOUT(3)
1054     AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1055     AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1056     * -----------------------------------------------
1057    
1058 pam-fi 1.13 IF(TRACKMODE.EQ.1) THEN
1059     * -----------------------------------------------
1060     * change of energy by bremsstrahlung for electrons
1061     VOUT(7) = VOUT(7) * 0.997 !0.9968
1062     * -----------------------------------------------
1063     ENDIF
1064 pam-fi 1.20 c$$$ PRINT *,'TRACKING -> END' ! TEST
1065 pam-fi 1.13
1066 mocchiut 1.1 ENDDO
1067    
1068 pam-fi 1.10 c$$$ print*,'POSXY (dopo) ',vout
1069    
1070    
1071 mocchiut 1.1 RETURN
1072     END
1073    
1074    
1075    
1076    
1077    
1078     * **********************************************************
1079     * Some initialization routines
1080     * **********************************************************
1081    
1082     * ----------------------------------------------------------
1083     * Routine to initialize COMMON/TRACK/
1084     *
1085     subroutine track_init
1086    
1087     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1088    
1089     include 'commontracker.f' !tracker general common
1090     include 'common_mini_2.f' !common for the tracking procedure
1091     include 'common_mech.f'
1092    
1093     do i=1,5
1094     AL(i) = 0.
1095     enddo
1096    
1097     do ip=1,NPLANES
1098     ZM(IP) = fitz(nplanes-ip+1) !init to mech. position
1099     XM(IP) = -100. !0.
1100     YM(IP) = -100. !0.
1101     XM_A(IP) = -100. !0.
1102     YM_A(IP) = -100. !0.
1103     c ZM_A(IP) = 0
1104     XM_B(IP) = -100. !0.
1105     YM_B(IP) = -100. !0.
1106     c ZM_B(IP) = 0
1107     RESX(IP) = 1000. !3.d-4
1108     RESY(IP) = 1000. !12.d-4
1109     XGOOD(IP) = 0
1110     YGOOD(IP) = 0
1111 pam-fi 1.15 DEDXTRK_X(IP) = 0
1112     DEDXTRK_Y(IP) = 0
1113     AXV(IP) = 0
1114     AYV(IP) = 0
1115     XV(IP) = -100
1116     YV(IP) = -100
1117 mocchiut 1.1 enddo
1118    
1119     return
1120     end
1121 pam-fi 1.4
1122    
1123     ***************************************************
1124     * *
1125     * *
1126     * *
1127     * *
1128     * *
1129     * *
1130     **************************************************
1131    
1132     subroutine guess()
1133    
1134     c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1135    
1136     include 'commontracker.f' !tracker general common
1137     include 'common_mini_2.f' !common for the tracking procedure
1138    
1139     REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)
1140     REAL*4 CHI,XC,ZC,RADIUS
1141     * ----------------------------------------
1142     * Y view
1143     * ----------------------------------------
1144     * ----------------------------------------
1145     * initial guess with a straigth line
1146     * ----------------------------------------
1147     SZZ=0.
1148     SZY=0.
1149     SSY=0.
1150     SZ=0.
1151     S1=0.
1152     DO I=1,nplanes
1153     IF(YGOOD(I).EQ.1)THEN
1154     YY = YM(I)
1155     IF(XGOOD(I).EQ.0)THEN
1156     YY = (YM_A(I) + YM_B(I))/2
1157     ENDIF
1158     SZZ=SZZ+ZM(I)*ZM(I)
1159     SZY=SZY+ZM(I)*YY
1160     SSY=SSY+YY
1161     SZ=SZ+ZM(I)
1162     S1=S1+1.
1163     ENDIF
1164     ENDDO
1165     DET=SZZ*S1-SZ*SZ
1166     AY=(SZY*S1-SZ*SSY)/DET
1167     BY=(SZZ*SSY-SZY*SZ)/DET
1168     Y0 = AY*ZINI+BY
1169     * ----------------------------------------
1170     * X view
1171     * ----------------------------------------
1172     * ----------------------------------------
1173     * 1) initial guess with a circle
1174     * ----------------------------------------
1175     NP=0
1176     DO I=1,nplanes
1177     IF(XGOOD(I).EQ.1)THEN
1178     XX = XM(I)
1179     IF(YGOOD(I).EQ.0)THEN
1180     XX = (XM_A(I) + XM_B(I))/2
1181     ENDIF
1182     NP=NP+1
1183     XP(NP)=XX
1184     ZP(NP)=ZM(I)
1185     ENDIF
1186     ENDDO
1187 pam-fi 1.9 IFLAG=0 !no debug mode
1188 pam-fi 1.4 CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
1189 pam-fi 1.14
1190     c$$$ print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
1191     c$$$ print*,' XP ',(xp(i),i=1,np)
1192     c$$$ print*,' ZP ',(zp(i),i=1,np)
1193     c$$$ print*,' AP ',(ap(i),i=1,np)
1194     c$$$ print*,' XP ',(rp(i),i=1,np)
1195    
1196 pam-fi 1.4 IF(IFLAG.NE.0)GOTO 10 !straigth fit
1197 pam-fi 1.14 c if(CHI.gt.100)GOTO 10 !straigth fit
1198 pam-fi 1.4 ARG = RADIUS**2-(ZINI-ZC)**2
1199     IF(ARG.LT.0)GOTO 10 !straigth fit
1200     DC = SQRT(ARG)
1201     IF(XC.GT.0)DC=-DC
1202     X0=XC+DC
1203     AX = -(ZINI-ZC)/DC
1204     DEF=100./(RADIUS*0.3*0.43)
1205     IF(XC.GT.0)DEF=-DEF
1206 pam-fi 1.8
1207 pam-fi 1.14
1208    
1209 pam-fi 1.8 IF(ABS(X0).GT.30)THEN
1210 pam-fi 1.10 c$$$ PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
1211     c$$$ $ ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
1212 pam-fi 1.8 GOTO 10 !straigth fit
1213     ENDIF
1214 pam-fi 1.4 GOTO 20 !guess is ok
1215    
1216     * ----------------------------------------
1217     * 2) initial guess with a straigth line
1218     * - if circle does not intersect reference plane
1219     * - if bad chi**2
1220     * ----------------------------------------
1221     10 CONTINUE
1222     SZZ=0.
1223     SZX=0.
1224     SSX=0.
1225     SZ=0.
1226     S1=0.
1227     DO I=1,nplanes
1228     IF(XGOOD(I).EQ.1)THEN
1229     XX = XM(I)
1230     IF(YGOOD(I).EQ.0)THEN
1231     XX = (XM_A(I) + XM_B(I))/2
1232     ENDIF
1233     SZZ=SZZ+ZM(I)*ZM(I)
1234     SZX=SZX+ZM(I)*XX
1235     SSX=SSX+XX
1236     SZ=SZ+ZM(I)
1237     S1=S1+1.
1238     ENDIF
1239     ENDDO
1240     DET=SZZ*S1-SZ*SZ
1241     AX=(SZX*S1-SZ*SSX)/DET
1242     BX=(SZZ*SSX-SZX*SZ)/DET
1243     DEF = 0
1244     X0 = AX*ZINI+BX
1245    
1246     20 CONTINUE
1247     * ----------------------------------------
1248     * guess
1249     * ----------------------------------------
1250    
1251     AL(1) = X0
1252     AL(2) = Y0
1253     tath = sqrt(AY**2+AX**2)
1254     AL(3) = tath/sqrt(1+tath**2)
1255 pam-fi 1.10 c$$$ IF(AX.NE.0)THEN
1256     c$$$ AL(4)= atan(AY/AX)
1257     c$$$ ELSE
1258     c$$$ AL(4) = acos(-1.)/2
1259     c$$$ IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1260     c$$$ ENDIF
1261     c$$$ IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1262     c$$$ AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1263    
1264     c$$$ AL(4) = 0.
1265     c$$$ IF(AX.NE.0.AND.AY.NE.0)THEN
1266     c$$$ AL(4)= atan(AY/AX)
1267     c$$$ ELSEIF(AY.EQ.0)THEN
1268     c$$$ AL(4) = 0.
1269     c$$$ IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)
1270     c$$$ ELSEIF(AX.EQ.0)THEN
1271     c$$$ AL(4) = acos(-1.)/2
1272     c$$$ IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1273     c$$$ ENDIF
1274     c$$$ IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1275     c$$$ AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1276    
1277     c$$$ AL(4)=0.
1278     c$$$ IF( AX.NE.0.OR.AY.NE.0. ) THEN
1279     c$$$ AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1280     c$$$ IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
1281     c$$$ ENDIF
1282    
1283     AL(4)=0.
1284     IF( AX.NE.0.OR.AY.NE.0. ) THEN
1285     AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1286     IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
1287     IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
1288 pam-fi 1.4 ENDIF
1289 pam-fi 1.10 IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
1290     IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
1291    
1292 pam-fi 1.4 AL(5) = DEF
1293    
1294     c print*,' guess: ',(al(i),i=1,5)
1295    
1296     end

  ViewVC Help
Powered by ViewVC 1.1.23