/[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.3 by pam-fi, Thu Oct 26 16:22:37 2006 UTC revision 1.16 by pam-fi, Mon May 14 11:03:06 2007 UTC
# Line 22  Line 22 
22  c      logical DEBUG  c      logical DEBUG
23  c      common/dbg/DEBUG  c      common/dbg/DEBUG
24    
25        parameter (inf=1.e8)      !just a huge number...        parameter (dinf=1.d15)      !just a huge number...
26  c------------------------------------------------------------------------  c------------------------------------------------------------------------
27  c     variables used in the tracking procedure (mini and its subroutines)  c     variables used in the tracking procedure (mini and its subroutines)
28  c  c
# Line 41  c      DATA XGOOD,YGOOD/nplanes*1.,nplan Line 41  c      DATA XGOOD,YGOOD/nplanes*1.,nplan
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/inf,inf,inf,inf,0.25e2/      !limits on alpha vector components  c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
45  c      DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !"  c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
46        DATA ALMAX/inf,inf,1.,inf,0.25e2/      !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
47        DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !"        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
   
48    
49        DIMENSION DAL(5)                    !increment of vector alfa        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--------
53          REAL*8 AVRESX,AVRESY
54    c     elena--------
55    
56        INTEGER IFLAG                    INTEGER IFLAG            
57  c--------------------------------------------------------  c--------------------------------------------------------
58  c     IFLAG =1 ---- chi2 derivatives computed by using  c     IFLAG =1 ---- chi2 derivatives computed by using
# Line 62  c     method 2 is faster!! Line 65  c     method 2 is faster!!
65  c--------------------------------------------------------  c--------------------------------------------------------
66        DATA IFLAG/2/          DATA IFLAG/2/  
67                    
68        LOGICAL TRKDEBUG  c      LOGICAL TRKDEBUG,TRKVERBOSE
69        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
70          LOGICAL TRKDEBUG,TRKVERBOSE
71          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
72    
73        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
74           TRKDEBUG = .TRUE.           TRKVERBOSE = .TRUE.
75             TRKDEBUG   = .FALSE.
76          ELSEIF(IPRINT.EQ.2)THEN
77             TRKVERBOSE = .TRUE.
78             TRKDEBUG   = .TRUE.
79        ELSE        ELSE
80           TRKDEBUG = .FALSE.           TRKVERBOSE = .FALSE.
81             TRKDEBUG   = .FALSE.
82        ENDIF                ENDIF        
83    
84  *     ----------------------------------------------------------  *     ----------------------------------------------------------
85    *     evaluate average spatial resolution
86    *     ----------------------------------------------------------
87          AVRESX = RESXAV
88          AVRESY = RESYAV
89          DO IP=1,6
90             IF( XGOOD(IP).EQ.1 )THEN
91                NX=NX+1
92                AVRESX=AVRESX+RESX(IP)
93             ENDIF
94             IF(NX.NE.0)AVRESX=AVRESX/NX
95             IF( YGOOD(IP).EQ.1 )THEN
96                NY=NY+1
97                AVRESY=AVRESY+RESY(IP)
98             ENDIF
99             IF(NX.NE.0)AVRESY=AVRESY/NY        
100          ENDDO
101    
102    *     ----------------------------------------------------------
103  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
104  *      *    
105  *     ----------------------------------------------------------  *     ----------------------------------------------------------
106        FACT=10.                  !scale factor to define  *     changed in order to evaluate energy-dependent
107                                  !tolerance on alfa  *     tolerances on all 5 parameters
108        ALTOL(1)=RESX(1)/FACT     !al(1) = x  cPP      FACT=1.0e10                  !scale factor to define tolerance on alfa
109        ALTOL(2)=RESY(1)/FACT     !al(2) = y  c     deflection error (see PDG)
110        ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta)        DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
111       $     +RESY(1)**2)/44.51/FACT        DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
112        ALTOL(4)=ALTOL(3)         !al(4)=phi  c$$$      ALTOL(1) = AVRESX/FACT     !al(1) = x
113    c$$$      ALTOL(2) = AVRESY/FACT     !al(2) = y
114    c$$$      ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta)
115    c$$$     $     +AVRESY**2)/44.51/FACT
116    c$$$      ALTOL(4) = ALTOL(3)         !al(4)=phi
117  c     deflection error (see PDG)  c     deflection error (see PDG)
118        DELETA1=0.01*RESX(1)/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))  c$$$      DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
119        DELETA2=0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
120  *     ----------------------------------------------------------  *     ----------------------------------------------------------
121  *      *    
122        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
123        JFAIL=0                   !error flag        JFAIL=0                   !error flag
124          CHI2=0
125    
126          if(TRKDEBUG) print*,'guess: ',al
127          if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
128    
129  *      *    
130  *     -----------------------  *     -----------------------
131  *     START MINIMIZATION LOOP  *     START MINIMIZATION LOOP
# Line 99  c     deflection error (see PDG) Line 136  c     deflection error (see PDG)
136        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
137           IFAIL=1           IFAIL=1
138           CHI2=-9999.           CHI2=-9999.
139           if(TRKDEBUG)           if(TRKVERBOSE)
140       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
141           RETURN           RETURN
142        ENDIF        ENDIF
143                
144  *      CHI2_P=CHI2        COST=1e-5
         
 c      print*,'@@@@@ ',istep,' - ',al  
   
       COST=1e-7  
145        DO I=1,5        DO I=1,5
146           DO J=1,5           DO J=1,5
147              CHI2DD(I,J)=CHI2DD(I,J)*COST              CHI2DD(I,J)=CHI2DD(I,J)*COST
# Line 123  ',istep,' - ',al Line 156  ',istep,' - ',al
156  *------------------------------------------------------------*  *------------------------------------------------------------*
157           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
158           IF(IFA.NE.0) THEN      !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
159              if(TRKDEBUG)then              if(TRKVERBOSE)then
160                 PRINT *,                 PRINT *,
161       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
162       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 137  ',istep,' - ',al Line 170  ',istep,' - ',al
170           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion               CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
171  *     *******************************************  *     *******************************************
172  *     find new value of AL-pha                  *     find new value of AL-pha                
173  *                                                *     *******************************************                                        
174           DO I=1,5                                         DO I=1,5                              
175              DAL(I)=0.                        DAL(I)=0.          
176              DO J=1,5                          DO J=1,5            
# Line 161  ',istep,' - ',al Line 194  ',istep,' - ',al
194           ENDDO           ENDDO
195           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
196           IF(IFA.NE.0) THEN           IF(IFA.NE.0) THEN
197              if(TRKDEBUG)then              if(TRKVERBOSE)then
198                 PRINT *,                 PRINT *,
199       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
200       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 173  ',istep,' - ',al Line 206  ',istep,' - ',al
206              RETURN                          RETURN            
207           ENDIF           ENDIF
208           CALL DSFINV(4,CHI2DD_R,4)           CALL DSFINV(4,CHI2DD_R,4)
209    *     *******************************************
210    *     find new value of AL-pha                
211    *     *******************************************                                        
212           DO I=1,4           DO I=1,4
213              DAL(I)=0.              DAL(I)=0.
214              DO J=1,4              DO J=1,4
# Line 185  ',istep,' - ',al Line 221  ',istep,' - ',al
221              AL(I)=AL(I)+DAL(I)              AL(I)=AL(I)+DAL(I)
222           ENDDO           ENDDO
223        ENDIF        ENDIF
224    
225          if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
226    
227  *------------------------------------------------------------*  *------------------------------------------------------------*
228  *     ----------------------------------------------------   *  *     ----------------------------------------------------   *
229  *------------------------------------------------------------*  *------------------------------------------------------------*
 *     *******************************************  
230  *     check parameter bounds:  *     check parameter bounds:
231    *------------------------------------------------------------*
232        DO I=1,5        DO I=1,5
233           IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN           IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN
234              if(TRKDEBUG)then              if(TRKVERBOSE)then
235                 PRINT*,' *** WARNING in mini ***  '                 PRINT*,' *** WARNING in mini ***  '
236                 PRINT*,'MINI_2 ==> AL(',I,') out of range'                 PRINT*,'MINI_2 ==> AL(',I,') out of range'
237                 PRINT*,'         value: ',AL(I),                 PRINT*,'         value: ',AL(I),
# Line 205  ',istep,' - ',al Line 244  ',istep,' - ',al
244              RETURN              RETURN
245           ENDIF           ENDIF
246        ENDDO              ENDDO      
247    *------------------------------------------------------------*
248  *     check number of steps:  *     check number of steps:
249        IF(ISTEP.gt.ISTEPMAX) then  *------------------------------------------------------------*
250           IFAIL=1        IF(ISTEP.ge.ISTEPMAX) then
251           if(TRKDEBUG)  c$$$         IFAIL=1
252       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',  c$$$         if(TRKVERBOSE)
253       $        ISTEPMAX  c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
254    c$$$     $        ISTEPMAX
255           goto 11           goto 11
256        endif        endif
257    *------------------------------------------------------------*
258  *     ---------------------------------------------  *     ---------------------------------------------
259  *     evaluate deflection tolerance on the basis of  *     evaluate deflection tolerance on the basis of
260  *     estimated deflection  *     estimated deflection
261  *     ---------------------------------------------  *     ---------------------------------------------
262        ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  *------------------------------------------------------------*
263    c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
264          ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
265          ALTOL(1) = ALTOL(5)/DELETA1
266          ALTOL(2) = ALTOL(1)
267          ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
268          ALTOL(4) = ALTOL(3)        
269          
270    c$$$      print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!
271    
272  *---- check tolerances:  *---- check tolerances:
273    c$$$      DO I=1,5
274    c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
275    c$$$      ENDDO
276    c$$$      print*,'chi2 -- ',DCHI2
277    
278          IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
279        DO I=1,5        DO I=1,5
280           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
281        ENDDO        ENDDO
# Line 229  ',istep,' - ',al Line 285  ',istep,' - ',al
285        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
286        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
287           IFAIL=1           IFAIL=1
288           if(TRKDEBUG)THEN           if(TRKVERBOSE)THEN
289              CHI2=-9999.              CHI2=-9999.
290           if(TRKDEBUG)           if(TRKVERBOSE)
291       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
292           ENDIF           ENDIF
293           RETURN           RETURN
# Line 246  ',istep,' - ',al Line 302  ',istep,' - ',al
302        IF(PFIXED.EQ.0.) THEN        IF(PFIXED.EQ.0.) THEN
303           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
304           IF(IFA.NE.0) THEN      !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
305              if(TRKDEBUG)then              if(TRKVERBOSE)then
306                 PRINT *,                 PRINT *,
307       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
308       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 273  ',istep,' - ',al Line 329  ',istep,' - ',al
329           ENDDO           ENDDO
330           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
331           IF(IFA.NE.0) THEN           IF(IFA.NE.0) THEN
332              if(TRKDEBUG)then              if(TRKVERBOSE)then
333                 PRINT *,                 PRINT *,
334       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
335       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 306  ',istep,' - ',al Line 362  ',istep,' - ',al
362        if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***        if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
363        if(ndof.le.0.) then        if(ndof.le.0.) then
364           ndof = 1           ndof = 1
365           if(TRKDEBUG)           if(TRKVERBOSE)
366       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
367        endif        endif
368    
369  *     ------------------------------------  *     ------------------------------------
370  *     Reduced chi^2  *     Reduced chi^2
371        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
372    
373    c      print*,'mini2: chi2 ',chi2
374    
375   11   CONTINUE         11   CONTINUE      
376    
377          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
378    
379        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
380    
381    c$$$      print*,'>>>>> NSTEP = ',NSTEP
382    
383        RETURN        RETURN
384        END                END        
385    
# Line 343  ',istep,' - ',al Line 406  ',istep,' - ',al
406       $     ,XV0(nplanes),YV0(nplanes)       $     ,XV0(nplanes),YV0(nplanes)
407        DIMENSION AL_P(5)        DIMENSION AL_P(5)
408    
409        LOGICAL TRKDEBUG  c      LOGICAL TRKVERBOSE
410        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKVERBOSE
411          LOGICAL TRKDEBUG,TRKVERBOSE
412          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
413  *      *    
414  *     chi^2 computation  *     chi^2 computation
415  *      *    
# Line 354  ',istep,' - ',al Line 419  ',istep,' - ',al
419        JFAIL=0                   !error flag        JFAIL=0                   !error flag
420        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
421        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
422           IF(TRKDEBUG)           IF(TRKVERBOSE)
423       $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'       $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
424           IFAIL=1           IFAIL=1
425           RETURN           RETURN
# Line 380  c$$$      ENDDO Line 445  c$$$      ENDDO
445  *     measured position of the cluster.  *     measured position of the cluster.
446  *     ---------------------------------------------------------  *     ---------------------------------------------------------
447        CHI2=0.        CHI2=0.
448        DO I=1,nplanes        DO I=1,nplanes        
449           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
450              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))
451              ALFA = XM_A(I) - BETA * YM_A(I)              ALFA = XM_A(I) - BETA * YM_A(I)
# Line 407  c$$$      ENDDO Line 472  c$$$      ENDDO
472       +                                       *( XGOOD(I)*(1-YGOOD(I)) )       +                                       *( XGOOD(I)*(1-YGOOD(I)) )
473       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2       +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
474       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
475    c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
476    c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
477    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
478    c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
479    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
480    c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
481    c$$$         print*,XV(I),XM(I),XGOOD(I)
482    c$$$         print*,YV(I),YM(I),YGOOD(I)
483        ENDDO        ENDDO
484    c$$$      print*,'CHISQ ',chi2
485  *     ------------------------------------------------  *     ------------------------------------------------
486  *      *    
487  *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)      *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)    
# Line 426  c$$$      ENDDO Line 500  c$$$      ENDDO
500              JFAIL=0              JFAIL=0
501              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
502              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
503                 IF(TRKDEBUG)                 IF(TRKVERBOSE)
504  *23456789012345678901234567890123456789012345678901234567890123456789012  *23456789012345678901234567890123456789012345678901234567890123456789012
505       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
506                 IFAIL=1                 IFAIL=1
# Line 440  c$$$      ENDDO Line 514  c$$$      ENDDO
514              JFAIL=0              JFAIL=0
515              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
516              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
517                 IF(TRKDEBUG)                 IF(TRKVERBOSE)
518       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
519                 IFAIL=1                 IFAIL=1
520                 RETURN                 RETURN
# Line 472  c$$$      ENDDO Line 546  c$$$      ENDDO
546                            
547              COSTHE=DSQRT(1.-AL(3)**2)              COSTHE=DSQRT(1.-AL(3)**2)
548              IF(COSTHE.EQ.0.) THEN              IF(COSTHE.EQ.0.) THEN
549                 IF(TRKDEBUG)PRINT *,'=== WARNING ===> COSTHE=0'                 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
550                 IFAIL=1                 IFAIL=1
551                 RETURN                 RETURN
552              ENDIF              ENDIF
# Line 559  c$$$      ENDDO Line 633  c$$$      ENDDO
633        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
634        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
635    
636        LOGICAL TRKDEBUG  c      LOGICAL TRKVERBOSE
637        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKVERBOSE
638          LOGICAL TRKDEBUG,TRKVERBOSE
639          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
640  c        c      
641        DIMENSION AL_P(5)        DIMENSION AL_P(5)
642  *      *    
643        DO I=1,nplanes  cpp      DO I=1,nplanes
644           ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
645        ENDDO        cpp      ENDDO      
646  *      *    
647  *     set parameters for GRKUTA  *     set parameters for GRKUTA
648  *      *    
# Line 580  c       Line 656  c      
656        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
657        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))
658        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
659    
660    c$$$      print*,'POSXY (prima) ',vout
661    
662        DO I=1,nplanes        DO I=1,nplanes
663           step=vout(3)-zv(i)  cpp         step=vout(3)-zv(i)
664             step=vout(3)-zm(i)
665   10      DO J=1,7   10      DO J=1,7
666              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
667              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
# Line 590  c       Line 670  c      
670           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
671           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
672              IFAIL=1              IFAIL=1
673              if(TRKDEBUG)              if(TRKVERBOSE)
674       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
675              if(.TRUE.)print*,'charge',charge  c$$$            if(.TRUE.)print*,'charge',charge
676              if(.TRUE.)print*,'vect',vect  c$$$            if(.TRUE.)print*,'vect',vect
677              if(.TRUE.)print*,'vout',vout  c$$$            if(.TRUE.)print*,'vout',vout
678              if(.TRUE.)print*,'step',step  c$$$            if(.TRUE.)print*,'step',step
679                if(TRKVERBOSE)print*,'charge',charge
680                if(TRKVERBOSE)print*,'vect',vect
681                if(TRKVERBOSE)print*,'vout',vout
682                if(TRKVERBOSE)print*,'step',step
683              RETURN              RETURN
684           ENDIF           ENDIF
685           Z=VOUT(3)           Z=VOUT(3)
# Line 609  c       Line 693  c      
693              GOTO 11              GOTO 11
694           ENDIF           ENDIF
695    
696    
697  *     -----------------------------------------------  *     -----------------------------------------------
698  *        evaluate track coordinates  *        evaluate track coordinates
699   100     XV(I)=VOUT(1)   100     XV(I)=VOUT(1)
# Line 618  c       Line 703  c      
703           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
704  *     -----------------------------------------------  *     -----------------------------------------------
705    
706             IF(TRACKMODE.EQ.1) THEN
707    *     -----------------------------------------------
708    *        change of energy by bremsstrahlung for electrons
709                VOUT(7) = VOUT(7) * 0.997 !0.9968
710    *     -----------------------------------------------
711             ENDIF
712    
713        ENDDO        ENDDO
714    
715    c$$$      print*,'POSXY (dopo) ',vout
716    
717    
718        RETURN        RETURN
719        END        END
720    
# Line 660  c         ZM_B(IP) = 0 Line 755  c         ZM_B(IP) = 0
755           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
756           XGOOD(IP) = 0           XGOOD(IP) = 0
757           YGOOD(IP) = 0           YGOOD(IP) = 0
758             DEDXTRK_X(IP) = 0
759             DEDXTRK_Y(IP) = 0
760             AXV(IP) = 0
761             AYV(IP) = 0
762             XV(IP) = -100
763             YV(IP) = -100
764        enddo        enddo
765    
766        return        return
767        end        end
768    
769    
770    ***************************************************
771    *                                                 *
772    *                                                 *
773    *                                                 *
774    *                                                 *
775    *                                                 *
776    *                                                 *
777    **************************************************
778    
779          subroutine guess()
780    
781    c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
782    
783          include 'commontracker.f' !tracker general common
784          include 'common_mini_2.f' !common for the tracking procedure
785          
786          REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)
787          REAL*4 CHI,XC,ZC,RADIUS
788    *     ----------------------------------------
789    *     Y view
790    *     ----------------------------------------
791    *     ----------------------------------------
792    *     initial guess with a straigth line
793    *     ----------------------------------------
794          SZZ=0.                    
795          SZY=0.
796          SSY=0.
797          SZ=0.
798          S1=0.
799          DO I=1,nplanes
800             IF(YGOOD(I).EQ.1)THEN
801                YY = YM(I)
802                IF(XGOOD(I).EQ.0)THEN
803                   YY = (YM_A(I) + YM_B(I))/2
804                ENDIF
805                SZZ=SZZ+ZM(I)*ZM(I)
806                SZY=SZY+ZM(I)*YY
807                SSY=SSY+YY
808                SZ=SZ+ZM(I)
809                S1=S1+1.
810             ENDIF
811          ENDDO
812          DET=SZZ*S1-SZ*SZ
813          AY=(SZY*S1-SZ*SSY)/DET
814          BY=(SZZ*SSY-SZY*SZ)/DET
815          Y0 = AY*ZINI+BY
816    *     ----------------------------------------
817    *     X view
818    *     ----------------------------------------
819    *     ----------------------------------------
820    *     1) initial guess with a circle
821    *     ----------------------------------------
822          NP=0
823          DO I=1,nplanes
824             IF(XGOOD(I).EQ.1)THEN
825                XX = XM(I)
826                IF(YGOOD(I).EQ.0)THEN
827                   XX = (XM_A(I) + XM_B(I))/2
828                ENDIF
829                NP=NP+1
830                XP(NP)=XX
831                ZP(NP)=ZM(I)
832             ENDIF
833          ENDDO
834          IFLAG=0                   !no debug mode
835          CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
836    
837    c$$$      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
838    c$$$      print*,' XP ',(xp(i),i=1,np)
839    c$$$      print*,' ZP ',(zp(i),i=1,np)
840    c$$$      print*,' AP ',(ap(i),i=1,np)
841    c$$$      print*,' XP ',(rp(i),i=1,np)
842    
843          IF(IFLAG.NE.0)GOTO 10 !straigth fit
844    c      if(CHI.gt.100)GOTO 10 !straigth fit
845          ARG = RADIUS**2-(ZINI-ZC)**2
846          IF(ARG.LT.0)GOTO 10       !straigth fit
847          DC = SQRT(ARG)      
848          IF(XC.GT.0)DC=-DC
849          X0=XC+DC
850          AX = -(ZINI-ZC)/DC
851          DEF=100./(RADIUS*0.3*0.43)
852          IF(XC.GT.0)DEF=-DEF
853          
854    
855    
856          IF(ABS(X0).GT.30)THEN
857    c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
858    c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
859             GOTO 10       !straigth fit
860          ENDIF
861          GOTO 20                   !guess is ok
862    
863    *     ----------------------------------------
864    *     2) initial guess with a straigth line
865    *     - if circle does not intersect reference plane
866    *     - if bad chi**2
867    *     ----------------------------------------
868     10   CONTINUE
869          SZZ=0.                  
870          SZX=0.
871          SSX=0.
872          SZ=0.
873          S1=0.
874          DO I=1,nplanes
875             IF(XGOOD(I).EQ.1)THEN
876                XX = XM(I)
877                IF(YGOOD(I).EQ.0)THEN
878                   XX = (XM_A(I) + XM_B(I))/2
879                ENDIF
880                SZZ=SZZ+ZM(I)*ZM(I)
881                SZX=SZX+ZM(I)*XX
882                SSX=SSX+XX
883                SZ=SZ+ZM(I)
884                S1=S1+1.
885             ENDIF
886          ENDDO
887          DET=SZZ*S1-SZ*SZ
888          AX=(SZX*S1-SZ*SSX)/DET
889          BX=(SZZ*SSX-SZX*SZ)/DET
890          DEF = 0
891          X0  = AX*ZINI+BX
892    
893     20   CONTINUE
894    *     ----------------------------------------
895    *     guess
896    *     ----------------------------------------
897    
898          AL(1) = X0
899          AL(2) = Y0
900          tath  = sqrt(AY**2+AX**2)
901          AL(3) = tath/sqrt(1+tath**2)
902    c$$$      IF(AX.NE.0)THEN
903    c$$$         AL(4)= atan(AY/AX)
904    c$$$      ELSE
905    c$$$         AL(4) = acos(-1.)/2
906    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
907    c$$$  ENDIF
908    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
909    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
910    
911    c$$$      AL(4) =  0.
912    c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN
913    c$$$         AL(4)= atan(AY/AX)
914    c$$$      ELSEIF(AY.EQ.0)THEN
915    c$$$         AL(4) = 0.
916    c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)        
917    c$$$      ELSEIF(AX.EQ.0)THEN
918    c$$$         AL(4) = acos(-1.)/2
919    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
920    c$$$      ENDIF
921    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
922    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
923    
924    c$$$      AL(4)=0.
925    c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN
926    c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
927    c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
928    c$$$      ENDIF
929    
930          AL(4)=0.
931          IF( AX.NE.0.OR.AY.NE.0. ) THEN
932             AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
933             IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
934             IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
935          ENDIF
936          IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
937          IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
938    
939          AL(5) = DEF
940    
941    c      print*,' guess: ',(al(i),i=1,5)
942    
943          end
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23