/[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.7 by pam-fi, Wed Nov 15 14:54:06 2006 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        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
45  c      DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !"        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
       DATA ALMAX/inf,inf,1.,inf,0.25e2/      !limits on alpha vector components  
       DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !"  
   
46    
47        DIMENSION DAL(5)                    !increment of vector alfa        DIMENSION DAL(5)                    !increment of vector alfa
48        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
49    
50    c     elena--------
51          REAL*8 AVRESX,AVRESY
52    c     elena--------
53    
54        INTEGER IFLAG                    INTEGER IFLAG            
55  c--------------------------------------------------------  c--------------------------------------------------------
56  c     IFLAG =1 ---- chi2 derivatives computed by using  c     IFLAG =1 ---- chi2 derivatives computed by using
# Line 62  c     method 2 is faster!! Line 63  c     method 2 is faster!!
63  c--------------------------------------------------------  c--------------------------------------------------------
64        DATA IFLAG/2/          DATA IFLAG/2/  
65                    
66        LOGICAL TRKDEBUG  c      LOGICAL TRKDEBUG,TRKVERBOSE
67        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
68          LOGICAL TRKDEBUG,TRKVERBOSE
69          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
70    
71        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
72           TRKDEBUG = .TRUE.           TRKVERBOSE = .TRUE.
73             TRKDEBUG   = .FALSE.
74          ELSEIF(IPRINT.EQ.2)THEN
75             TRKVERBOSE = .TRUE.
76             TRKDEBUG   = .TRUE.
77        ELSE        ELSE
78           TRKDEBUG = .FALSE.           TRKVERBOSE = .FALSE.
79             TRKDEBUG   = .FALSE.
80        ENDIF                ENDIF        
81    
82  *     ----------------------------------------------------------  *     ----------------------------------------------------------
83    *     evaluate average spatial resolution
84    *     ----------------------------------------------------------
85          AVRESX = RESXAV
86          AVRESY = RESYAV
87          DO IP=1,6
88             IF( XGOOD(IP).EQ.1 )THEN
89                NX=NX+1
90                AVRESX=AVRESX+RESX(IP)
91             ENDIF
92             IF(NX.NE.0)AVRESX=AVRESX/NX
93             IF( YGOOD(IP).EQ.1 )THEN
94                NY=NY+1
95                AVRESY=AVRESY+RESY(IP)
96             ENDIF
97             IF(NX.NE.0)AVRESY=AVRESY/NY        
98          ENDDO
99    
100    *     ----------------------------------------------------------
101  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
102  *      *    
103  *     ----------------------------------------------------------  *     ----------------------------------------------------------
104        FACT=10.                  !scale factor to define  *     changed in order to evaluate energy-dependent
105                                  !tolerance on alfa  *     tolerances on all 5 parameters
106        ALTOL(1)=RESX(1)/FACT     !al(1) = x        FACT=100.                  !scale factor to define tolerance on alfa
107        ALTOL(2)=RESY(1)/FACT     !al(2) = y  c     deflection error (see PDG)
108        ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta)        DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
109       $     +RESY(1)**2)/44.51/FACT        DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
110        ALTOL(4)=ALTOL(3)         !al(4)=phi  c$$$      ALTOL(1) = AVRESX/FACT     !al(1) = x
111    c$$$      ALTOL(2) = AVRESY/FACT     !al(2) = y
112    c$$$      ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta)
113    c$$$     $     +AVRESY**2)/44.51/FACT
114    c$$$      ALTOL(4) = ALTOL(3)         !al(4)=phi
115  c     deflection error (see PDG)  c     deflection error (see PDG)
116        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.))
117        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)
118  *     ----------------------------------------------------------  *     ----------------------------------------------------------
119  *      *    
120        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
121        JFAIL=0                   !error flag        JFAIL=0                   !error flag
122    
123          if(TRKDEBUG) print*,'guess: ',al
124          if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
125    
126  *      *    
127  *     -----------------------  *     -----------------------
128  *     START MINIMIZATION LOOP  *     START MINIMIZATION LOOP
# Line 99  c     deflection error (see PDG) Line 133  c     deflection error (see PDG)
133        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
134           IFAIL=1           IFAIL=1
135           CHI2=-9999.           CHI2=-9999.
136           if(TRKDEBUG)           if(TRKVERBOSE)
137       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
138           RETURN           RETURN
139        ENDIF        ENDIF
140                
 *      CHI2_P=CHI2  
         
 c      print*,'@@@@@ ',istep,' - ',al  
   
141        COST=1e-7        COST=1e-7
142        DO I=1,5        DO I=1,5
143           DO J=1,5           DO J=1,5
# Line 123  ',istep,' - ',al Line 153  ',istep,' - ',al
153  *------------------------------------------------------------*  *------------------------------------------------------------*
154           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
155           IF(IFA.NE.0) THEN      !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
156              if(TRKDEBUG)then              if(TRKVERBOSE)then
157                 PRINT *,                 PRINT *,
158       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
159       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 137  ',istep,' - ',al Line 167  ',istep,' - ',al
167           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion               CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
168  *     *******************************************  *     *******************************************
169  *     find new value of AL-pha                  *     find new value of AL-pha                
170  *                                                *     *******************************************                                        
171           DO I=1,5                                         DO I=1,5                              
172              DAL(I)=0.                        DAL(I)=0.          
173              DO J=1,5                          DO J=1,5            
# Line 161  ',istep,' - ',al Line 191  ',istep,' - ',al
191           ENDDO           ENDDO
192           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
193           IF(IFA.NE.0) THEN           IF(IFA.NE.0) THEN
194              if(TRKDEBUG)then              if(TRKVERBOSE)then
195                 PRINT *,                 PRINT *,
196       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
197       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 173  ',istep,' - ',al Line 203  ',istep,' - ',al
203              RETURN                          RETURN            
204           ENDIF           ENDIF
205           CALL DSFINV(4,CHI2DD_R,4)           CALL DSFINV(4,CHI2DD_R,4)
206    *     *******************************************
207    *     find new value of AL-pha                
208    *     *******************************************                                        
209           DO I=1,4           DO I=1,4
210              DAL(I)=0.              DAL(I)=0.
211              DO J=1,4              DO J=1,4
# Line 185  ',istep,' - ',al Line 218  ',istep,' - ',al
218              AL(I)=AL(I)+DAL(I)              AL(I)=AL(I)+DAL(I)
219           ENDDO           ENDDO
220        ENDIF        ENDIF
221    
222          if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
223    
224  *------------------------------------------------------------*  *------------------------------------------------------------*
225  *     ----------------------------------------------------   *  *     ----------------------------------------------------   *
226  *------------------------------------------------------------*  *------------------------------------------------------------*
 *     *******************************************  
227  *     check parameter bounds:  *     check parameter bounds:
228    *------------------------------------------------------------*
229        DO I=1,5        DO I=1,5
230           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
231              if(TRKDEBUG)then              if(TRKVERBOSE)then
232                 PRINT*,' *** WARNING in mini ***  '                 PRINT*,' *** WARNING in mini ***  '
233                 PRINT*,'MINI_2 ==> AL(',I,') out of range'                 PRINT*,'MINI_2 ==> AL(',I,') out of range'
234                 PRINT*,'         value: ',AL(I),                 PRINT*,'         value: ',AL(I),
# Line 205  ',istep,' - ',al Line 241  ',istep,' - ',al
241              RETURN              RETURN
242           ENDIF           ENDIF
243        ENDDO              ENDDO      
244    *------------------------------------------------------------*
245  *     check number of steps:  *     check number of steps:
246        IF(ISTEP.gt.ISTEPMAX) then  *------------------------------------------------------------*
247           IFAIL=1        IF(ISTEP.ge.ISTEPMAX) then
248           if(TRKDEBUG)  c$$$         IFAIL=1
249       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',  c$$$         if(TRKVERBOSE)
250       $        ISTEPMAX  c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
251    c$$$     $        ISTEPMAX
252           goto 11           goto 11
253        endif        endif
254    *------------------------------------------------------------*
255  *     ---------------------------------------------  *     ---------------------------------------------
256  *     evaluate deflection tolerance on the basis of  *     evaluate deflection tolerance on the basis of
257  *     estimated deflection  *     estimated deflection
258  *     ---------------------------------------------  *     ---------------------------------------------
259        ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  *------------------------------------------------------------*
260    c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
261          ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
262          ALTOL(1) = ALTOL(5)/DELETA1
263          ALTOL(2) = ALTOL(1)
264          ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
265          ALTOL(4) = ALTOL(3)        
266          
267  *---- check tolerances:  *---- check tolerances:
268    c$$$      DO I=1,5
269    c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
270    c$$$      ENDDO
271    c$$$      print*,'chi2 -- ',DCHI2
272    
273        DO I=1,5        DO I=1,5
274           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
275        ENDDO        ENDDO
# Line 229  ',istep,' - ',al Line 279  ',istep,' - ',al
279        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
280        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
281           IFAIL=1           IFAIL=1
282           if(TRKDEBUG)THEN           if(TRKVERBOSE)THEN
283              CHI2=-9999.              CHI2=-9999.
284           if(TRKDEBUG)           if(TRKVERBOSE)
285       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
286           ENDIF           ENDIF
287           RETURN           RETURN
# Line 246  ',istep,' - ',al Line 296  ',istep,' - ',al
296        IF(PFIXED.EQ.0.) THEN        IF(PFIXED.EQ.0.) THEN
297           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
298           IF(IFA.NE.0) THEN      !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
299              if(TRKDEBUG)then              if(TRKVERBOSE)then
300                 PRINT *,                 PRINT *,
301       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
302       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 273  ',istep,' - ',al Line 323  ',istep,' - ',al
323           ENDDO           ENDDO
324           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
325           IF(IFA.NE.0) THEN           IF(IFA.NE.0) THEN
326              if(TRKDEBUG)then              if(TRKVERBOSE)then
327                 PRINT *,                 PRINT *,
328       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
329       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 306  ',istep,' - ',al Line 356  ',istep,' - ',al
356        if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***        if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
357        if(ndof.le.0.) then        if(ndof.le.0.) then
358           ndof = 1           ndof = 1
359           if(TRKDEBUG)           if(TRKVERBOSE)
360       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
361        endif        endif
362    
363          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
364    
365  *     ------------------------------------  *     ------------------------------------
366  *     Reduced chi^2  *     Reduced chi^2
367        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
368    
369    c      print*,'mini2: chi2 ',chi2
370    
371   11   CONTINUE         11   CONTINUE      
372    
373        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
# Line 343  ',istep,' - ',al Line 398  ',istep,' - ',al
398       $     ,XV0(nplanes),YV0(nplanes)       $     ,XV0(nplanes),YV0(nplanes)
399        DIMENSION AL_P(5)        DIMENSION AL_P(5)
400    
401        LOGICAL TRKDEBUG  c      LOGICAL TRKVERBOSE
402        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKVERBOSE
403          LOGICAL TRKDEBUG,TRKVERBOSE
404          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
405  *      *    
406  *     chi^2 computation  *     chi^2 computation
407  *      *    
# Line 354  ',istep,' - ',al Line 411  ',istep,' - ',al
411        JFAIL=0                   !error flag        JFAIL=0                   !error flag
412        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
413        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
414           IF(TRKDEBUG)           IF(TRKVERBOSE)
415       $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'       $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
416           IFAIL=1           IFAIL=1
417           RETURN           RETURN
# Line 408  c$$$      ENDDO Line 465  c$$$      ENDDO
465       +        +((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
466       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
467        ENDDO        ENDDO
468    c      print*,'CHISQ ',chi2
469  *     ------------------------------------------------  *     ------------------------------------------------
470  *      *    
471  *     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 484  c$$$      ENDDO
484              JFAIL=0              JFAIL=0
485              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
486              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
487                 IF(TRKDEBUG)                 IF(TRKVERBOSE)
488  *23456789012345678901234567890123456789012345678901234567890123456789012  *23456789012345678901234567890123456789012345678901234567890123456789012
489       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
490                 IFAIL=1                 IFAIL=1
# Line 440  c$$$      ENDDO Line 498  c$$$      ENDDO
498              JFAIL=0              JFAIL=0
499              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
500              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
501                 IF(TRKDEBUG)                 IF(TRKVERBOSE)
502       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
503                 IFAIL=1                 IFAIL=1
504                 RETURN                 RETURN
# Line 472  c$$$      ENDDO Line 530  c$$$      ENDDO
530                            
531              COSTHE=DSQRT(1.-AL(3)**2)              COSTHE=DSQRT(1.-AL(3)**2)
532              IF(COSTHE.EQ.0.) THEN              IF(COSTHE.EQ.0.) THEN
533                 IF(TRKDEBUG)PRINT *,'=== WARNING ===> COSTHE=0'                 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
534                 IFAIL=1                 IFAIL=1
535                 RETURN                 RETURN
536              ENDIF              ENDIF
# Line 559  c$$$      ENDDO Line 617  c$$$      ENDDO
617        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
618        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
619    
620        LOGICAL TRKDEBUG  c      LOGICAL TRKVERBOSE
621        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKVERBOSE
622          LOGICAL TRKDEBUG,TRKVERBOSE
623          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
624  c        c      
625        DIMENSION AL_P(5)        DIMENSION AL_P(5)
626  *      *    
# Line 580  c       Line 640  c      
640        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
641        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))
642        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
643    
644    c$$$      print*,'POSXY ',vout
645    
646        DO I=1,nplanes        DO I=1,nplanes
647           step=vout(3)-zv(i)           step=vout(3)-zv(i)
648   10      DO J=1,7   10      DO J=1,7
# Line 590  c       Line 653  c      
653           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
654           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
655              IFAIL=1              IFAIL=1
656              if(TRKDEBUG)              if(TRKVERBOSE)
657       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
658              if(.TRUE.)print*,'charge',charge  c$$$            if(.TRUE.)print*,'charge',charge
659              if(.TRUE.)print*,'vect',vect  c$$$            if(.TRUE.)print*,'vect',vect
660              if(.TRUE.)print*,'vout',vout  c$$$            if(.TRUE.)print*,'vout',vout
661              if(.TRUE.)print*,'step',step  c$$$            if(.TRUE.)print*,'step',step
662                if(TRKVERBOSE)print*,'charge',charge
663                if(TRKVERBOSE)print*,'vect',vect
664                if(TRKVERBOSE)print*,'vout',vout
665                if(TRKVERBOSE)print*,'step',step
666              RETURN              RETURN
667           ENDIF           ENDIF
668           Z=VOUT(3)           Z=VOUT(3)
# Line 664  c         ZM_B(IP) = 0 Line 731  c         ZM_B(IP) = 0
731    
732        return        return
733        end        end
734    
735    
736    ***************************************************
737    *                                                 *
738    *                                                 *
739    *                                                 *
740    *                                                 *
741    *                                                 *
742    *                                                 *
743    **************************************************
744    
745          subroutine guess()
746    
747    c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
748    
749          include 'commontracker.f' !tracker general common
750          include 'common_mini_2.f' !common for the tracking procedure
751          
752          REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)
753          REAL*4 CHI,XC,ZC,RADIUS
754    *     ----------------------------------------
755    *     Y view
756    *     ----------------------------------------
757    *     ----------------------------------------
758    *     initial guess with a straigth line
759    *     ----------------------------------------
760          SZZ=0.                    
761          SZY=0.
762          SSY=0.
763          SZ=0.
764          S1=0.
765          DO I=1,nplanes
766             IF(YGOOD(I).EQ.1)THEN
767                YY = YM(I)
768                IF(XGOOD(I).EQ.0)THEN
769                   YY = (YM_A(I) + YM_B(I))/2
770                ENDIF
771                SZZ=SZZ+ZM(I)*ZM(I)
772                SZY=SZY+ZM(I)*YY
773                SSY=SSY+YY
774                SZ=SZ+ZM(I)
775                S1=S1+1.
776             ENDIF
777          ENDDO
778          DET=SZZ*S1-SZ*SZ
779          AY=(SZY*S1-SZ*SSY)/DET
780          BY=(SZZ*SSY-SZY*SZ)/DET
781          Y0 = AY*ZINI+BY
782    *     ----------------------------------------
783    *     X view
784    *     ----------------------------------------
785    *     ----------------------------------------
786    *     1) initial guess with a circle
787    *     ----------------------------------------
788          NP=0
789          DO I=1,nplanes
790             IF(XGOOD(I).EQ.1)THEN
791                XX = XM(I)
792                IF(YGOOD(I).EQ.0)THEN
793                   XX = (XM_A(I) + XM_B(I))/2
794                ENDIF
795                NP=NP+1
796                XP(NP)=XX
797                ZP(NP)=ZM(I)
798             ENDIF
799          ENDDO
800          CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
801    c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI
802          IF(IFLAG.NE.0)GOTO 10 !straigth fit
803          ARG = RADIUS**2-(ZINI-ZC)**2
804          IF(ARG.LT.0)GOTO 10       !straigth fit
805          DC = SQRT(ARG)      
806          IF(XC.GT.0)DC=-DC
807          X0=XC+DC
808          AX = -(ZINI-ZC)/DC
809          DEF=100./(RADIUS*0.3*0.43)
810          IF(XC.GT.0)DEF=-DEF
811          GOTO 20                   !guess is ok
812    
813    *     ----------------------------------------
814    *     2) initial guess with a straigth line
815    *     - if circle does not intersect reference plane
816    *     - if bad chi**2
817    *     ----------------------------------------
818     10   CONTINUE
819          SZZ=0.                  
820          SZX=0.
821          SSX=0.
822          SZ=0.
823          S1=0.
824          DO I=1,nplanes
825             IF(XGOOD(I).EQ.1)THEN
826                XX = XM(I)
827                IF(YGOOD(I).EQ.0)THEN
828                   XX = (XM_A(I) + XM_B(I))/2
829                ENDIF
830                SZZ=SZZ+ZM(I)*ZM(I)
831                SZX=SZX+ZM(I)*XX
832                SSX=SSX+XX
833                SZ=SZ+ZM(I)
834                S1=S1+1.
835             ENDIF
836          ENDDO
837          DET=SZZ*S1-SZ*SZ
838          AX=(SZX*S1-SZ*SSX)/DET
839          BX=(SZZ*SSX-SZX*SZ)/DET
840          DEF = 0
841          X0  = AX*ZINI+BX
842    
843     20   CONTINUE
844    *     ----------------------------------------
845    *     guess
846    *     ----------------------------------------
847    
848          AL(1) = X0
849          AL(2) = Y0
850          tath  = sqrt(AY**2+AX**2)
851          AL(3) = tath/sqrt(1+tath**2)
852          IF(AX.NE.0)THEN
853             AL(4)= atan(AY/AX)
854          ELSE
855             AL(4) = acos(-1.)/2
856             IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
857          ENDIF
858          IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
859          AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs
860          AL(5) = DEF
861    
862    c      print*,' guess: ',(al(i),i=1,5)
863    
864          end
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23