/[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.5 by pam-fi, Wed Nov 8 16:42:28 2006 UTC revision 1.22 by bongi, Thu Nov 20 15:06:27 2008 UTC
# Line 36  c--------------------------------------- Line 36  c---------------------------------------
36  c      DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking  c      DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking
37    
38        DATA STEPAL/5*1.d-7/      !alpha vector step        DATA STEPAL/5*1.d-7/      !alpha vector step
39        DATA ISTEPMAX/120/        !maximum number of steps in the chi^2 minimization        DATA ISTEPMAX/100/        !maximum number of steps in the chi^2 minimization
40        DATA TOLL/1.d-8/          !tolerance in reaching the next plane during        DATA TOLL/1.d-8/          !tolerance in reaching the next plane during
41  *                               !the tracking procedure  *                               !the tracking procedure
42        DATA STEPMAX/100./        !maximum number of steps in the trackin gprocess        DATA STEPMAX/100./        !maximum number of steps in the trackin gprocess
43    
44    c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
45    c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
46        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
47        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
48    
49        DIMENSION DAL(5)                    !increment of vector alfa  c$$$      DIMENSION DAL(5)                    !increment of vector alfa
50        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
51    
52  c     elena--------  c     elena--------
# Line 65  c--------------------------------------- Line 67  c---------------------------------------
67                    
68  c      LOGICAL TRKDEBUG,TRKVERBOSE  c      LOGICAL TRKDEBUG,TRKVERBOSE
69  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
70        LOGICAL TRKDEBUG,TRKVERBOSE        LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT
71        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
72    
73          DIMENSION AL0(5)
74          LOGICAL SUCCESS_NEW,SUCCESS_OLD
75    
76    c$$$      PRINT*,'==========' ! TEST
77    c$$$      PRINT*,'START MINI' ! TEST
78    c$$$      PRINT*,'==========' ! TEST
79    
80    *
81    *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
82    *
83          STUDENT = .false.
84          FIRSTSTEPS = .true.
85          FIRSTSTUDENT = .true.
86          IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.
87    
88        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
89           TRKVERBOSE = .TRUE.           TRKVERBOSE = .TRUE.
90           TRKDEBUG   = .FALSE.           TRKDEBUG   = .FALSE.
# Line 84  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 101  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
101  *     ----------------------------------------------------------  *     ----------------------------------------------------------
102        AVRESX = RESXAV        AVRESX = RESXAV
103        AVRESY = RESYAV        AVRESY = RESYAV
104          NX = 0.0
105          NY = 0.0
106        DO IP=1,6        DO IP=1,6
107           IF( XGOOD(IP).EQ.1 )THEN           IF( XGOOD(IP).EQ.1 )THEN
108              NX=NX+1              NX=NX+1.0
109              AVRESX=AVRESX+RESX(IP)              AVRESX=AVRESX+RESX(IP)
110           ENDIF           ENDIF
          IF(NX.NE.0)AVRESX=AVRESX/NX  
111           IF( YGOOD(IP).EQ.1 )THEN           IF( YGOOD(IP).EQ.1 )THEN
112              NY=NY+1              NY=NY+1.0
113              AVRESY=AVRESY+RESY(IP)              AVRESY=AVRESY+RESY(IP)
114           ENDIF           ENDIF
          IF(NX.NE.0)AVRESY=AVRESY/NY          
115        ENDDO        ENDDO
116          IF(NX.NE.0.0)AVRESX=AVRESX/NX
117          IF(NY.NE.0.0)AVRESY=AVRESY/NY        
118    
119  *     ----------------------------------------------------------  *     ----------------------------------------------------------
120  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
# Line 103  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 122  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
122  *     ----------------------------------------------------------  *     ----------------------------------------------------------
123  *     changed in order to evaluate energy-dependent  *     changed in order to evaluate energy-dependent
124  *     tolerances on all 5 parameters  *     tolerances on all 5 parameters
125        FACT=100.                  !scale factor to define tolerance on alfa  cPP      FACT=1.0e10                  !scale factor to define tolerance on alfa
126  c     deflection error (see PDG)  c     deflection error (see PDG)
127        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.))
128        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 119  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 138  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
138  *      *    
139        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
140        JFAIL=0                   !error flag        JFAIL=0                   !error flag
141          CHI2=0
142    
143        if(TRKDEBUG) print*,'guess: ',al        if(TRKDEBUG) print*,'guess: ',al
144        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
# Line 129  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 149  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
149  *     -----------------------  *     -----------------------
150   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!
151    
152        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  * -------------------------------
153        IF(JFAIL.NE.0) THEN  * **** Chi2+gaussian minimization
154           IFAIL=1  * -------------------------------
155           CHI2=-9999.  
156           if(TRKVERBOSE)        IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN
157       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'  
158           RETURN           IF(ISTEP.GE.3) FIRSTSTEPS = .false.
159        ENDIF  
160                   CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
161        COST=1e-7           IF(JFAIL.NE.0) THEN
162        DO I=1,5              IFAIL=1
163           DO J=1,5              CHI2=-9999.
164              CHI2DD(I,J)=CHI2DD(I,J)*COST              if(TRKVERBOSE)
165         $           PRINT *,'*** ERROR in mini *** wrong CHISQ'
166                RETURN
167             ENDIF
168            
169    c      COST=1e-5
170             COST=1.
171             DO I=1,5
172                IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
173             ENDDO
174             DO I=1,5
175                DO J=1,5
176                   CHI2DD(I,J)=CHI2DD(I,J)*COST
177                ENDDO
178    c$$$            CHI2D(I)=CHI2D(I)*COST
179           ENDDO           ENDDO
          CHI2D(I)=CHI2D(I)*COST  
       ENDDO        
180    
181        IF(PFIXED.EQ.0.) THEN           IF(PFIXED.EQ.0.) THEN
182    
183  *------------------------------------------------------------*  *------------------------------------------------------------*
184  *     track fitting with FREE deflection  *     track fitting with FREE deflection
185  *------------------------------------------------------------*  *------------------------------------------------------------*
186           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant              CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
187           IF(IFA.NE.0) THEN      !not positive-defined                    IF(IFA.NE.0) THEN   !not positive-defined      
188              if(TRKVERBOSE)then                 if(TRKVERBOSE)then
189                 PRINT *,                    PRINT *,
190       $              '*** ERROR in mini ***'//       $                 '*** ERROR in mini ***'//
191       $              'on matrix inversion (not pos-def)'       $                 'on matrix inversion (not pos-def)'
192       $              ,DET       $                 ,DET
193              endif                 endif
194              IF(CHI2.EQ.0) CHI2=-9999.                 IF(CHI2.EQ.0) CHI2=-9999.
195              IF(CHI2.GT.0) CHI2=-CHI2                 IF(CHI2.GT.0) CHI2=-CHI2
196              IFAIL=1                 IFAIL=1
197              RETURN                             RETURN            
198           ENDIF              ENDIF
199           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion                  CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
200  *     *******************************************  *     *******************************************
201  *     find new value of AL-pha                  *     find new value of AL-pha                
202  *     *******************************************                                          *     *******************************************                                        
203           DO I=1,5                                            DO I=1,5                              
204              DAL(I)=0.                           DAL(I)=0.          
205              DO J=1,5                             DO J=1,5            
206                 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)                    DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)  *COST
207                 COV(I,J)=2.*COST*CHI2DD(I,J)                    COV(I,J)=2.*COST*CHI2DD(I,J)
208              ENDDO                               ENDDO  
209           ENDDO                                ENDDO                  
210           DO I=1,5                            DO I=1,5              
211              AL(I)=AL(I)+DAL(I)                   AL(I)=AL(I)+DAL(I)  
212           ENDDO                                ENDDO                  
213  *------------------------------------------------------------*  *------------------------------------------------------------*
214  *     track fitting with FIXED deflection  *     track fitting with FIXED deflection
215  *------------------------------------------------------------*  *------------------------------------------------------------*
216        ELSE           ELSE
217           AL(5)=1./PFIXED              AL(5)=1./PFIXED
218           DO I=1,4              DO I=1,4
219              CHI2D_R(I)=CHI2D(I)                 CHI2D_R(I)=CHI2D(I)
220              DO J=1,4                 DO J=1,4
221                 CHI2DD_R(I,J)=CHI2DD(I,J)                    CHI2DD_R(I,J)=CHI2DD(I,J)
222                   ENDDO
223              ENDDO              ENDDO
224           ENDDO              CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
225           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)              IF(IFA.NE.0) THEN
226           IF(IFA.NE.0) THEN                 if(TRKVERBOSE)then
227              if(TRKVERBOSE)then                    PRINT *,
228                 PRINT *,       $                 '*** ERROR in mini ***'//
229       $              '*** ERROR in mini ***'//       $                 'on matrix inversion (not pos-def)'
230       $              'on matrix inversion (not pos-def)'       $                 ,DET
231       $              ,DET                 endif
232              endif                 IF(CHI2.EQ.0) CHI2=-9999.
233              IF(CHI2.EQ.0) CHI2=-9999.                 IF(CHI2.GT.0) CHI2=-CHI2
234              IF(CHI2.GT.0) CHI2=-CHI2                 IFAIL=1
235              IFAIL=1                 RETURN            
236              RETURN                          ENDIF
237           ENDIF              CALL DSFINV(4,CHI2DD_R,4)
          CALL DSFINV(4,CHI2DD_R,4)  
238  *     *******************************************  *     *******************************************
239  *     find new value of AL-pha                  *     find new value of AL-pha                
240  *     *******************************************                                          *     *******************************************                                        
241           DO I=1,4              DO I=1,4
242              DAL(I)=0.                 DAL(I)=0.
243              DO J=1,4                 DO J=1,4
244                 DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)                    DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)  *COST
245                 COV(I,J)=2.*COST*CHI2DD_R(I,J)                    COV(I,J)=2.*COST*CHI2DD_R(I,J)
246                   ENDDO
247                ENDDO
248                DAL(5)=0.
249                DO I=1,4
250                   AL(I)=AL(I)+DAL(I)
251                ENDDO
252             ENDIF
253          
254             if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
255    
256    c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
257    c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
258    
259    
260          ENDIF
261    
262    * -------------------------------
263    * **** Likelihood+Student minimization
264    * -------------------------------
265    
266          IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN
267    
268             IF(FIRSTSTUDENT) THEN
269                FIRSTSTUDENT = .false.
270                ISTEP = 1
271             ENDIF
272    
273             CALL CHISQSTT(1,JFAIL)
274             DO I=1,5                              
275                DAL(I)=0.          
276                DO J=1,5            
277                   DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)
278                ENDDO  
279             ENDDO                  
280            
281             DO I=1,5
282                DO j=1,5
283                   COV(I,J) = 2.*CHI2DD(I,J)
284              ENDDO              ENDDO
285           ENDDO           ENDDO
286           DAL(5)=0.          
287           DO I=1,4           CHI2TOLL = 1.E-3
288              AL(I)=AL(I)+DAL(I)           ALPHA = 3.0
289             BETA = -0.4
290             E=1.
291             EA=1.
292             EB=1.
293             EC=1.
294             FA=1.
295             FB=1.
296             FC=1.
297             SUCCESS_OLD = .FALSE.
298             SUCCESS_NEW = .FALSE.
299    
300             CALL CHISQSTT(0,JFAIL)
301    c$$$         PRINT*,CHI2
302             CHI2_NEW = CHI2
303             FC = CHI2
304             EC = 0.
305    
306             ICOUNT = 0
307     100     CONTINUE
308             ICOUNT = ICOUNT+1
309    
310             DO I=1,5
311                AL0(I)=AL(I)
312             ENDDO
313             DO I=1,5
314                AL(I)=AL(I)+E*DAL(I)
315             ENDDO
316             CALL CHISQSTT(0,JFAIL)
317             CHI2_OLD = CHI2_NEW
318             CHI2_NEW = CHI2
319             FA = FB
320             FB = FC
321             FC = CHI2
322             EA = EB
323             EB = EC
324             EC = E
325    
326    c$$$         PRINT*,E,CHI2_NEW
327    
328             IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success
329                IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101
330                SUCCESS_OLD = SUCCESS_NEW
331                SUCCESS_NEW = .TRUE.
332                E = E*ALPHA
333             ELSE                   ! failure
334                SUCCESS_OLD = SUCCESS_NEW
335                SUCCESS_NEW = .FALSE.
336                CHI2_NEW = CHI2_OLD
337                DO I=1,5
338                   AL(I)=AL0(I)
339                ENDDO
340                IF(SUCCESS_OLD) THEN
341                   DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA)
342                   IF(DENOM.NE.0.) THEN
343                      E = EB - 0.5*( (EB-EA)**2*(FB-FC)
344         $                 - (EB-EC)**2*(FB-FA) ) / DENOM
345                   ELSE
346                      E = BETA*E
347                   ENDIF
348                ELSE
349                   E = BETA*E
350                ENDIF
351    c$$$            E = BETA*E
352             ENDIF
353             IF(ICOUNT.GT.20) GOTO 101
354             GOTO 100
355    
356     101     CONTINUE
357    
358             DO I=1,5
359                DAL(I)=E*DAL(I)
360           ENDDO           ENDDO
361    
362    c$$$         print*,' '
363    c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
364    c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
365    c$$$         print*,'==== CHI2 ===='
366    c$$$         print*,chi2
367    c$$$         print*,'==== CHI2d ===='
368    c$$$         print*,(chi2d(i),i=1,5)
369    c$$$         print*,'==== CHI2dd ===='
370    c$$$         do j=1,5
371    c$$$            print*,(chi2dd(j,i),i=1,5)
372    c$$$         enddo
373    c$$$         print*,'================'
374    c$$$         print*,' '
375            
376    *========= FIN QUI =============
377    
378        ENDIF        ENDIF
379    
380        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)          
381    
382    
383    
384  *------------------------------------------------------------*  *------------------------------------------------------------*
385  *     ----------------------------------------------------   *  *     ----------------------------------------------------   *
# Line 245  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 405  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
405  *     check number of steps:  *     check number of steps:
406  *------------------------------------------------------------*  *------------------------------------------------------------*
407        IF(ISTEP.ge.ISTEPMAX) then        IF(ISTEP.ge.ISTEPMAX) then
408           IFAIL=1  c$$$         IFAIL=1
409           if(TRKVERBOSE)  c$$$         if(TRKVERBOSE)
410       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',  c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
411       $        ISTEPMAX  c$$$     $        ISTEPMAX
412           goto 11           goto 11
413        endif        endif
414  *------------------------------------------------------------*  *------------------------------------------------------------*
# Line 264  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE Line 424  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE
424        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
425        ALTOL(4) = ALTOL(3)                ALTOL(4) = ALTOL(3)        
426                
427    c$$$      print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!
428    
429  *---- check tolerances:  *---- check tolerances:
430  c$$$      DO I=1,5  c$$$      DO I=1,5
431  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
432  c$$$      ENDDO  c$$$      ENDDO
433  c$$$      print*,'chi2 -- ',DCHI2  c$$$      print*,'chi2 -- ',DCHI2
434    
435          IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
436        DO I=1,5        DO I=1,5
437           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
438        ENDDO        ENDDO
439    
440  *     new estimate of chi^2:  *****************************
441        JFAIL=0                   !error flag  *     final estimate of chi^2
442        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  *****************************
443        IF(JFAIL.NE.0) THEN  
444           IFAIL=1  * -------------------------------
445           if(TRKVERBOSE)THEN  * **** Chi2+gaussian minimization
446              CHI2=-9999.  * -------------------------------
447           if(TRKVERBOSE)  
448       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'        IF(.NOT.STUDENT) THEN
449           ENDIF  
450           RETURN           JFAIL=0                !error flag
451        ENDIF           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
452        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  
453              IFAIL=1              IFAIL=1
454              RETURN                          if(TRKVERBOSE)THEN
455                   CHI2=-9999.
456                   if(TRKVERBOSE)
457         $              PRINT *,'*** ERROR in mini *** wrong CHISQ'
458                ENDIF
459                RETURN
460           ENDIF           ENDIF
461           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion      c      COST=1e-7
462           DO I=1,5                                         COST=1.
463              DAL(I)=0.                     DO I=1,5
464              DO J=1,5                          IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
465                 COV(I,J)=2.*COST*CHI2DD(I,J)           ENDDO
466              ENDDO                         DO I=1,5
467           ENDDO                                DO J=1,5
468        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)  
469              ENDDO              ENDDO
470           ENDDO           ENDDO
471           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           IF(PFIXED.EQ.0.) THEN
472           IF(IFA.NE.0) THEN              CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
473              if(TRKVERBOSE)then              IF(IFA.NE.0) THEN   !not positive-defined      
474                 PRINT *,                 if(TRKVERBOSE)then
475       $              '*** ERROR in mini ***'//                    PRINT *,
476       $              'on matrix inversion (not pos-def)'       $                 '*** ERROR in mini ***'//
477       $              ,DET       $                 'on matrix inversion (not pos-def)'
478              endif       $                 ,DET
479              IF(CHI2.EQ.0) CHI2=-9999.                 endif
480              IF(CHI2.GT.0) CHI2=-CHI2                 IF(CHI2.EQ.0) CHI2=-9999.
481              IFAIL=1                 IF(CHI2.GT.0) CHI2=-CHI2
482              RETURN                             IFAIL=1
483                   RETURN            
484                ENDIF
485                CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
486                DO I=1,5                              
487    c$$$               DAL(I)=0.          
488                   DO J=1,5            
489                      COV(I,J)=2.*COST*CHI2DD(I,J)
490                   ENDDO              
491                ENDDO                  
492             ELSE
493                DO I=1,4
494                   CHI2D_R(I)=CHI2D(I)
495                   DO J=1,4
496                      CHI2DD_R(I,J)=CHI2DD(I,J)
497                   ENDDO
498                ENDDO
499                CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
500                IF(IFA.NE.0) THEN
501                   if(TRKVERBOSE)then
502                      PRINT *,
503         $                 '*** ERROR in mini ***'//
504         $                 'on matrix inversion (not pos-def)'
505         $                 ,DET
506                   endif
507                   IF(CHI2.EQ.0) CHI2=-9999.
508                   IF(CHI2.GT.0) CHI2=-CHI2
509                   IFAIL=1
510                   RETURN            
511                ENDIF
512                CALL DSFINV(4,CHI2DD_R,4)
513                DO I=1,4
514    c$$$               DAL(I)=0.
515                   DO J=1,4
516                      COV(I,J)=2.*COST*CHI2DD_R(I,J)
517                   ENDDO
518                ENDDO
519           ENDIF           ENDIF
520           CALL DSFINV(4,CHI2DD_R,4)  
521           DO I=1,4        ENDIF
522              DAL(I)=0.  
523              DO J=1,4  * -------------------------------
524                 COV(I,J)=2.*COST*CHI2DD_R(I,J)  * **** Likelihood+student minimization
525    * -------------------------------
526    
527          IF(STUDENT) THEN
528             CALL CHISQSTT(1,JFAIL)
529             DO I=1,5
530                DO j=1,5
531                   COV(I,J) = 2.*CHI2DD(I,J)
532              ENDDO              ENDDO
533           ENDDO           ENDDO
534        ENDIF        ENDIF
535    
536  *****************************  *****************************
537    
538  *     ------------------------------------  *     ------------------------------------
# Line 360  c$$$      print*,'chi2 -- ',DCHI2 Line 551  c$$$      print*,'chi2 -- ',DCHI2
551       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
552        endif        endif
553    
       if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)  
   
554  *     ------------------------------------  *     ------------------------------------
555  *     Reduced chi^2  *     Reduced chi^2
556        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
# Line 370  c      print*,'mini2: chi2 ',chi2 Line 559  c      print*,'mini2: chi2 ',chi2
559    
560   11   CONTINUE         11   CONTINUE      
561    
562          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
563    
564        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
565    
566    c$$$      print*,'>>>>> NSTEP = ',NSTEP
567    
568        RETURN        RETURN
569        END                END        
570    
# Line 437  c$$$      ENDDO Line 630  c$$$      ENDDO
630  *     measured position of the cluster.  *     measured position of the cluster.
631  *     ---------------------------------------------------------  *     ---------------------------------------------------------
632        CHI2=0.        CHI2=0.
633    
634          DO I=1,nplanes                                          
635             IF( XGOOD(I).NE.YGOOD(I) ) THEN ! singlet
636                IF(XGOOD(I).EQ.1) THEN
637                   Z =
638         $         ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) -
639         $           (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/
640         $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -
641         $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )
642                   ZM(I) = Z
643                   ZV(I) = Z
644                   XV(I) = XV_A(I)+(XV_B(I)-XV_A(I))*
645         $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))
646                   Y =
647         $         ( (ZM_A(I)*YM_B(I)-YM_A(I)*ZM_B(I))*(YV_A(I)-YV_B(I)) -
648         $           (ZV_A(I)*YV_B(I)-YV_A(I)*ZV_B(I))*(YM_A(I)-YM_B(I)) )/
649         $         ( (ZM_A(I)-ZM_B(I))*(YV_A(I)-YV_B(I)) -
650         $           (ZV_A(I)-ZV_B(I))*(YM_A(I)-YM_B(I)) )
651                   YM(I) = Y
652                   YV(I) = Y
653                   XM(I) = XM_A(I)+(XM_B(I)-XM_A(I))*
654         $                 (Y-YM_A(I))/(YM_B(I)-YM_A(I))
655    
656                   CHI2=CHI2+(XV(I)-XM(I))**2/RESX(I)**2
657    
658                ENDIF
659                IF(YGOOD(I).EQ.1) THEN
660                   Z =
661         $         ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(ZV_A(I)-ZV_B(I)) -
662         $           (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(ZM_A(I)-ZM_B(I)) )/
663         $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -
664         $           (ZV_A(I)-ZV_B(I))*(XM_A(I)-XM_B(I)) )
665                   ZM(I) = Z
666                   ZV(I) = Z
667                   YV(I) = YV_A(I)+(YV_B(I)-YV_A(I))*
668         $                 (Z-ZV_A(I))/(ZV_B(I)-ZV_A(I))
669                   X =
670         $         ( (ZM_A(I)*XM_B(I)-XM_A(I)*ZM_B(I))*(XV_A(I)-XV_B(I)) -
671         $           (ZV_A(I)*XV_B(I)-XV_A(I)*ZV_B(I))*(XM_A(I)-XM_B(I)) )/
672         $         ( (ZM_A(I)-ZM_B(I))*(XV_A(I)-XV_B(I)) -
673         $           (ZV_A(I)-ZV_B(I))*(XM_A(I)-XM_B(I)) )
674                   XM(I) = X
675                   XV(I) = X
676                   YM(I) = YM_A(I)+(YM_B(I)-YM_A(I))*
677         $                 (X-XM_A(I))/(XM_B(I)-XM_A(I))
678    
679                   CHI2=CHI2+(YV(I)-YM(I))**2/RESY(I)**2
680    
681                ENDIF
682             ELSEIF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.1)THEN !Y-cl
683                CHI2=CHI2
684         +           +(XV(I)-XM(I))**2/RESX(i)**2
685         +           +(YV(I)-YM(I))**2/RESY(i)**2
686             ENDIF
687          ENDDO
688        DO I=1,nplanes        DO I=1,nplanes
689           IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl           XV0(I)=XV(I)
690              BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))           YV0(I)=YV(I)
             ALFA = XM_A(I) - BETA * YM_A(I)  
             YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)  
             if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))  
      $           YM(I)=dmin1(YM_A(I),YM_B(I))  
             if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))  
      $           YM(I)=dmax1(YM_A(I),YM_B(I))  
             XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates  
          ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl  
             BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))  
             ALFA = YM_A(I) - BETA * XM_A(I)  
             XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)  
             if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))  
      $           XM(I)=dmin1(XM_A(I),XM_B(I))  
             if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))  
      $           XM(I)=dmax1(XM_A(I),XM_B(I))  
             YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates  
          ENDIF  
          CHI2=CHI2  
      +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )  
      +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )  
      +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2  
      +                                       *( XGOOD(I)*(1-YGOOD(I)) )  
      +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2  
      +                                       *( (1-XGOOD(I))*YGOOD(I) )  
691        ENDDO        ENDDO
692  c      print*,'CHISQ ',chi2  
693    c$$$      DO I=1,nplanes        
694    c$$$         IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
695    c$$$            BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
696    c$$$            ALFA = XM_A(I) - BETA * YM_A(I)
697    c$$$            YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
698    c$$$            if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
699    c$$$     $           YM(I)=dmin1(YM_A(I),YM_B(I))
700    c$$$            if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
701    c$$$     $           YM(I)=dmax1(YM_A(I),YM_B(I))
702    c$$$            XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
703    c$$$         ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
704    c$$$            BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
705    c$$$            ALFA = YM_A(I) - BETA * XM_A(I)
706    c$$$            XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
707    c$$$            if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
708    c$$$     $           XM(I)=dmin1(XM_A(I),XM_B(I))
709    c$$$            if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
710    c$$$     $           XM(I)=dmax1(XM_A(I),XM_B(I))
711    c$$$            YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
712    c$$$         ENDIF
713    c$$$         CHI2=CHI2
714    c$$$     +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
715    c$$$     +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
716    c$$$     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
717    c$$$     +                                       *( XGOOD(I)*(1-YGOOD(I)) )
718    c$$$     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
719    c$$$     +                                       *( (1-XGOOD(I))*YGOOD(I) )
720    c$$$      ENDDO
721    
722    c$$$      print*,'CHISQ ',chi2
723  *     ------------------------------------------------  *     ------------------------------------------------
724  *      *    
725  *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)      *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)    
# Line 585  c      print*,'CHISQ ',chi2 Line 839  c      print*,'CHISQ ',chi2
839        RETURN        RETURN
840        END        END
841    
842    ******************************************************************************    
843    *    
844    *     routine to compute Likelihodd+Student and its derivatives
845    *    
846    *     (modified in respect to the previous one in order to include
847    *     single clusters. In this case the residual is evaluated by
848    *     calculating the distance between the track intersection and the
849    *     segment AB associated to the single cluster)
850    *
851    ******************************************************************************
852    
853          SUBROUTINE CHISQSTT(IFLAG,JFAIL)
854    
855          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
856    
857          include 'commontracker.f' !tracker general common
858          include 'common_mini_2.f' !common for the tracking procedure
859    
860          LOGICAL TRKDEBUG,TRKVERBOSE
861          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
862    
863          DIMENSION AL_P(5)
864          DIMENSION VECTEMP(5)
865    c$$$      DIMENSION U(5) ! BFGS
866    
867          DO I=1,5
868             AL_P(I)=AL(I)
869          ENDDO
870          JFAIL=0                   !error flag
871          CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
872          IF(JFAIL.NE.0) THEN
873             IF(TRKVERBOSE)
874         $        PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!'
875             IFAIL=1
876             RETURN
877          ENDIF
878    
879          DO I=1,nplanes
880             DXDAL(I,1)=1.
881             DYDAL(I,1)=0.
882             DXDAL(I,2)=0.
883             DYDAL(I,2)=1.
884             COSTHE=DSQRT(1.-AL(3)**2)
885             IF(COSTHE.EQ.0.) THEN
886                IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
887                IFAIL=1
888                RETURN
889             ENDIF
890             DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
891             DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
892             DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
893             DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
894             IF(AL(5).NE.0.) THEN      
895                DXDAL(I,5)=
896         +           (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
897         +           *DCOS(AL(4))))/AL(5)
898                DYDAL(I,5)=
899         +           (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
900         +           *DSIN(AL(4))))/AL(5)
901             ELSE
902                DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
903                DYDAL(I,5)=0.
904             ENDIF
905          ENDDO
906    
907          IF(IFLAG.EQ.0) THEN ! function calulation
908             CHI2=0.
909             DO I=1,nplanes
910                IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
911                   BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
912                   ALFA = XM_A(I) - BETA * YM_A(I)
913                   YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
914                   if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
915         $              YM(I)=dmin1(YM_A(I),YM_B(I))
916                   if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
917         $              YM(I)=dmax1(YM_A(I),YM_B(I))
918                   XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
919                ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
920                   BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
921                   ALFA = YM_A(I) - BETA * XM_A(I)
922                   XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
923                   if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
924         $              XM(I)=dmin1(XM_A(I),XM_B(I))
925                   if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
926         $              XM(I)=dmax1(XM_A(I),XM_B(I))
927                   YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
928                ENDIF
929                TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/
930         $           (TAILX(I)*RESX(I)**2) )
931                TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/
932         $           (TAILY(I)*RESY(I)**2) )
933                CHI2=CHI2
934         $           +(TAILX(I)+1.0)*TERMX   *( XGOOD(I) )
935         $           +(TAILY(I)+1.0)*TERMY   *( YGOOD(I) )
936             ENDDO
937          ENDIF
938    
939          IF(IFLAG.EQ.1) THEN ! derivative calulation
940             DO I=1,5
941                CHI2DOLD(I)=CHI2D(I)
942             ENDDO
943             DO J=1,5
944                CHI2D(J)=0.      
945                DO I=1,nplanes
946                   CHI2D(J)=CHI2D(J)
947         $              +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/
948         $              (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)*
949         $              DXDAL(I,J)   *XGOOD(I)
950         $              +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/
951         $              (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)*
952         $              DYDAL(I,J)   *YGOOD(I)
953                ENDDO
954             ENDDO
955             DO K=1,5
956                VECTEMP(K)=0.
957                DO M=1,5
958                   VECTEMP(K) = VECTEMP(K) +
959         $              COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M))
960                ENDDO
961             ENDDO
962             DOWN1 = 0.
963             DO K=1,5
964                DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K))
965             ENDDO
966             IF(DOWN1.EQ.0.) THEN
967                PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0'
968                IFAIL=1
969                RETURN
970             ENDIF
971             DOWN2 = 0.
972             DO K=1,5
973                DO M=1,5
974                   DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K)
975                ENDDO
976             ENDDO
977             IF(DOWN2.EQ.0.) THEN
978                PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0'
979                IFAIL=1
980                RETURN
981             ENDIF
982    c$$$         DO K=1,5 ! BFGS
983    c$$$            U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2
984    c$$$         ENDDO
985             DO I=1,5
986                DO J=1,5
987                   CHI2DD(I,J) = COV(I,J)/2.
988         $              +DAL(I)*DAL(J)/DOWN1
989         $              -VECTEMP(I)*VECTEMP(J)/DOWN2
990    c$$$     $              +DOWN2*U(I)*U(J) ! BFGS
991                ENDDO
992             ENDDO
993          ENDIF
994    
995          RETURN
996          END
997          
998  *****************************************************************  *****************************************************************
999  *  *
1000  *     Routine to compute the track intersection points  *     Routine to compute the track intersection points
# Line 610  c      print*,'CHISQ ',chi2 Line 1019  c      print*,'CHISQ ',chi2
1019  *      *    
1020  *****************************************************************  *****************************************************************
1021    
1022    cPPP --- new --- (with singlets in 3D)
1023        SUBROUTINE POSXYZ(AL_P,IFAIL)        SUBROUTINE POSXYZ(AL_P,IFAIL)
1024    
1025        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# Line 623  c      COMMON/TRKD/TRKVERBOSE Line 1033  c      COMMON/TRKD/TRKVERBOSE
1033        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
1034  c        c      
1035        DIMENSION AL_P(5)        DIMENSION AL_P(5)
1036          LOGICAL SINGLET,SINGLET_FIRST,ZDEGENER
1037  *      *    
1038        DO I=1,nplanes  cpp      DO I=1,nplanes
1039           ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
1040        ENDDO        cpp      ENDDO      
1041  *      *    
1042  *     set parameters for GRKUTA  *     set parameters for GRKUTA
1043  *      *    
# Line 641  c       Line 1052  c      
1052        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))
1053        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
1054    
1055  c$$$      print*,'POSXY ',vout  c$$$      print*,'POSXY (prima) ',vout
1056    
1057        DO I=1,nplanes        DO I=1,nplanes
1058           step=vout(3)-zv(i)           IF(XGOOD(I).EQ.YGOOD(I)) SINGLET = .false.
1059             IF(XGOOD(I).NE.YGOOD(I)) SINGLET = .true.
1060             ZNEXT = ZM(I)
1061             IF(SINGLET) THEN
1062                IF(ZM_A(I).GT.ZM_B(I)+TOLL) THEN
1063                   ZNEXT = ZM_A(I)
1064                   ZNEXT2 = ZM_B(I)
1065                   SINGLET_FIRST = .true.
1066                   ZDEGENER = .false.
1067                ELSEIF(ZM_B(I).GT.ZM_A(I)+TOLL) THEN
1068                   ZNEXT = ZM_B(I)
1069                   ZNEXT2 = ZM_A(I)
1070                   SINGLET_FIRST = .true.
1071                   ZDEGENER = .false.
1072                ELSE
1073                   ZNEXT = 0.5*(ZM_A(I)+ZM_B(I))
1074                   SINGLET_FIRST = .false.
1075                   ZDEGENER = .true.
1076                ENDIF
1077             ENDIF
1078    c$$$         IF(SINGLET) PRINT*,'SINGLET!!!'
1079   10      DO J=1,7   10      DO J=1,7
1080              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
1081              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
1082           ENDDO           ENDDO
1083    
1084             IF(VOUT(6).GE.0.) THEN
1085                IFAIL=1
1086                if(TRKVERBOSE)      
1087         $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
1088                RETURN
1089             ENDIF
1090    cPP         step=(zm(i)-vect(3))/VOUT(6)
1091             step=(ZNEXT-vect(3))/VOUT(6)
1092   11      continue   11      continue
1093           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
1094    c$$$         ipass = ipass + 1 ! TEST
1095    c$$$         PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST
1096           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
1097              IFAIL=1              IFAIL=1
1098              if(TRKVERBOSE)              if(TRKVERBOSE)
# Line 663  c$$$            if(.TRUE.)print*,'step', Line 1105  c$$$            if(.TRUE.)print*,'step',
1105              if(TRKVERBOSE)print*,'vect',vect              if(TRKVERBOSE)print*,'vect',vect
1106              if(TRKVERBOSE)print*,'vout',vout              if(TRKVERBOSE)print*,'vout',vout
1107              if(TRKVERBOSE)print*,'step',step              if(TRKVERBOSE)print*,'step',step
1108                if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
1109              RETURN              RETURN
1110           ENDIF           ENDIF
1111           Z=VOUT(3)           Z=VOUT(3)
1112           IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100           IF(Z.LE.ZNEXT+TOLL.AND.Z.GE.ZNEXT-TOLL) GOTO 100
1113           IF(Z.GT.ZM(I)+TOLL) GOTO 10                 IF(Z.GT.ZNEXT+TOLL) GOTO 10
1114           IF(Z.LE.ZM(I)-TOLL) THEN           IF(Z.LE.ZNEXT-TOLL) THEN
1115              STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZNEXT-VECT(3))/(Z-VECT(3))
1116              DO J=1,7              DO J=1,7
1117                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
1118              ENDDO              ENDDO
1119              GOTO 11              GOTO 11
1120           ENDIF           ENDIF
1121    c$$$         IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100
1122    c$$$         IF(Z.GT.ZM(I)+TOLL) GOTO 10      
1123    c$$$         IF(Z.LE.ZM(I)-TOLL) THEN
1124    c$$$            STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
1125    c$$$            DO J=1,7
1126    c$$$               VECT(J)=VECTINI(J)
1127    c$$$            ENDDO
1128    c$$$            GOTO 11
1129    c$$$         ENDIF
1130    
1131    
1132  *     -----------------------------------------------  *     -----------------------------------------------
1133  *        evaluate track coordinates  *        evaluate track coordinates
1134   100     XV(I)=VOUT(1)  
1135           YV(I)=VOUT(2)   100     IF(SINGLET.AND.(.NOT.ZDEGENER).AND.SINGLET_FIRST) THEN
1136           ZV(I)=VOUT(3)              IF(ZM_A(I).GT.ZM_B(I)) THEN
1137           AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)                 XV_A(I) = VOUT(1)
1138           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)                 YV_A(I) = VOUT(2)
1139                   ZV_A(I) = VOUT(3)
1140                ELSE
1141                   XV_B(I) = VOUT(1)
1142                   YV_B(I) = VOUT(2)
1143                   ZV_B(I) = VOUT(3)
1144                ENDIF
1145                AXVUP = DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1146                AYVUP = DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1147                ZNEXT = ZNEXT2
1148                SINGLET_FIRST = .false.
1149                GOTO 10
1150             ENDIF
1151    
1152             IF(SINGLET.AND.(.NOT.ZDEGENER).AND.(.NOT.SINGLET_FIRST)) THEN
1153                IF(ZM_A(I).LT.ZM_B(I)) THEN
1154                   XV_A(I) = VOUT(1)
1155                   YV_A(I) = VOUT(2)
1156                   ZV_A(I) = VOUT(3)
1157                ELSE
1158                   XV_B(I) = VOUT(1)
1159                   YV_B(I) = VOUT(2)
1160                   ZV_B(I) = VOUT(3)
1161                ENDIF
1162                AXV(I)=0.5*(DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)+AXVUP)
1163                AYV(I)=0.5*(DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)+AYVUP)
1164             ENDIF
1165    
1166             IF(SINGLET.AND.ZDEGENER) THEN
1167                XV_A(I) = VOUT(1)
1168                YV_A(I) = VOUT(2)
1169                ZV_A(I) = VOUT(3)+0.1
1170                XV_B(I) = VOUT(1)
1171                YV_B(I) = VOUT(2)
1172                ZV_B(I) = VOUT(3)-0.1
1173                AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1174                AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1175             ENDIF
1176    
1177             IF(.NOT.SINGLET) THEN
1178                XV(I)=VOUT(1)
1179                YV(I)=VOUT(2)
1180                ZV(I)=VOUT(3)
1181                AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
1182                AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1183             ENDIF
1184    
1185    *     -----------------------------------------------
1186            
1187             IF(TRACKMODE.EQ.1) THEN
1188    *     -----------------------------------------------
1189    *        change of energy by bremsstrahlung for electrons
1190                VOUT(7) = VOUT(7) * 0.997 !0.9968
1191  *     -----------------------------------------------  *     -----------------------------------------------
1192             ENDIF
1193    c$$$         PRINT *,'TRACKING -> END' ! TEST
1194    
1195        ENDDO        ENDDO
1196    
1197    c$$$      print*,'POSXY (dopo) ',vout
1198    
1199    
1200        RETURN        RETURN
1201        END        END
1202    
# Line 719  c$$$            if(.TRUE.)print*,'step', Line 1229  c$$$            if(.TRUE.)print*,'step',
1229           YM(IP) = -100.         !0.           YM(IP) = -100.         !0.
1230           XM_A(IP) = -100.         !0.           XM_A(IP) = -100.         !0.
1231           YM_A(IP) = -100.         !0.           YM_A(IP) = -100.         !0.
1232  c         ZM_A(IP) = 0           ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position
1233           XM_B(IP) =  -100.         !0.           XM_B(IP) =  -100.         !0.
1234           YM_B(IP) =  -100.         !0.           YM_B(IP) =  -100.         !0.
1235  c         ZM_B(IP) = 0           ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position
1236           RESX(IP) = 1000.       !3.d-4           RESX(IP) = 1000.       !3.d-4
1237           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
1238           XGOOD(IP) = 0           XGOOD(IP) = 0
1239           YGOOD(IP) = 0           YGOOD(IP) = 0
1240             DEDXTRK_X(IP) = 0
1241             DEDXTRK_Y(IP) = 0
1242             AXV(IP) = 0
1243             AYV(IP) = 0
1244             XV(IP) = -100
1245             YV(IP) = -100
1246        enddo        enddo
1247    
1248        return        return
# Line 797  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1313  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1313              ZP(NP)=ZM(I)              ZP(NP)=ZM(I)
1314           ENDIF           ENDIF
1315        ENDDO        ENDDO
1316          IFLAG=0                   !no debug mode
1317        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
1318  c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI  
1319    c$$$      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
1320    c$$$      print*,' XP ',(xp(i),i=1,np)
1321    c$$$      print*,' ZP ',(zp(i),i=1,np)
1322    c$$$      print*,' AP ',(ap(i),i=1,np)
1323    c$$$      print*,' XP ',(rp(i),i=1,np)
1324    
1325        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
1326    c      if(CHI.gt.100)GOTO 10 !straigth fit
1327        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = RADIUS**2-(ZINI-ZC)**2
1328        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
1329        DC = SQRT(ARG)              DC = SQRT(ARG)      
# Line 808  c      print*,' circle: ',XC,ZC,RADIUS,' Line 1332  c      print*,' circle: ',XC,ZC,RADIUS,'
1332        AX = -(ZINI-ZC)/DC        AX = -(ZINI-ZC)/DC
1333        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
1334        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
1335          
1336    
1337    
1338          IF(ABS(X0).GT.30)THEN
1339    c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
1340    c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
1341             GOTO 10       !straigth fit
1342          ENDIF
1343        GOTO 20                   !guess is ok        GOTO 20                   !guess is ok
1344    
1345  *     ----------------------------------------  *     ----------------------------------------
# Line 849  c      print*,' circle: ',XC,ZC,RADIUS,' Line 1381  c      print*,' circle: ',XC,ZC,RADIUS,'
1381        AL(2) = Y0        AL(2) = Y0
1382        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1383        AL(3) = tath/sqrt(1+tath**2)        AL(3) = tath/sqrt(1+tath**2)
1384        IF(AX.NE.0)THEN  c$$$      IF(AX.NE.0)THEN
1385           AL(4)= atan(AY/AX)  c$$$         AL(4)= atan(AY/AX)
1386        ELSE  c$$$      ELSE
1387           AL(4) = acos(-1.)/2  c$$$         AL(4) = acos(-1.)/2
1388           IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1389    c$$$  ENDIF
1390    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1391    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1392    
1393    c$$$      AL(4) =  0.
1394    c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN
1395    c$$$         AL(4)= atan(AY/AX)
1396    c$$$      ELSEIF(AY.EQ.0)THEN
1397    c$$$         AL(4) = 0.
1398    c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)        
1399    c$$$      ELSEIF(AX.EQ.0)THEN
1400    c$$$         AL(4) = acos(-1.)/2
1401    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1402    c$$$      ENDIF
1403    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1404    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1405    
1406    c$$$      AL(4)=0.
1407    c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN
1408    c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1409    c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
1410    c$$$      ENDIF
1411    
1412          AL(4)=0.
1413          IF( AX.NE.0.OR.AY.NE.0. ) THEN
1414             AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1415             IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
1416             IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
1417        ENDIF        ENDIF
1418        IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)        IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
1419        AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs        IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
1420    
1421        AL(5) = DEF        AL(5) = DEF
1422    
1423  c      print*,' guess: ',(al(i),i=1,5)  c      print*,' guess: ',(al(i),i=1,5)
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23