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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.27

  ViewVC Help
Powered by ViewVC 1.1.23