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