/[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.6 by pam-fi, Tue Nov 14 16:21:09 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 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/100/        !maximum number of steps in the chi^2 minimization        DATA ISTEPMAX/120/        !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/inf,inf,inf,inf,0.25e2/      !limits on alpha vector components        DATA ALMAX/dinf,dinf,dinf,dinf,dinf/ !limits on alpha vector components
45  c      DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !"        DATA ALMIN/-dinf,-dinf,-dinf,-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                
141  *      CHI2_P=CHI2        COST=1e-9
142                costfac=1.1
143  c      print*,'@@@@@ ',istep,' - ',al        flagstep=0
144    
145     7    continue
146    
147           IF(costfac.le.1.) THEN
148            IFAIL=1
149            PRINT *,'=== WARNING ===> no matrix inversion '
150            RETURN
151           ENDIF
152    
       COST=1e-7  
153        DO I=1,5        DO I=1,5
154           DO J=1,5           DO J=1,5
155              CHI2DD(I,J)=CHI2DD(I,J)*COST              CHI2DD(I,J)=CHI2DD(I,J)*COST
# Line 122  ',istep,' - ',al Line 163  ',istep,' - ',al
163  *     track fitting with FREE deflection  *     track fitting with FREE deflection
164  *------------------------------------------------------------*  *------------------------------------------------------------*
165           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
166           IF(IFA.NE.0) THEN      !not positive-defined                 IF(IFA.NE.0.or.jfa.ne.0) THEN      !not positive-defined      
167              if(TRKDEBUG)then             if(ifa.eq.-1)then
168                 if(TRKVERBOSE)then
169                 PRINT *,                 PRINT *,
170       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
171       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
172       $              ,DET       $              ,DET
173              endif               endif
174              IF(CHI2.EQ.0) CHI2=-9999.               IF(CHI2.EQ.0) CHI2=-9999.
175              IF(CHI2.GT.0) CHI2=-CHI2               IF(CHI2.GT.0) CHI2=-CHI2
176              IFAIL=1               IFAIL=1
177              RETURN                           RETURN  
178               endif
179               if(jfa.eq.-1)then
180                 if(flagstep.eq.-1.) costfac=(costfac-1)/2+1
181                 cost=cost*costfact
182                 flagstep=1.
183                 goto 7
184               elseif(jfa.eq.1)then
185                 if(flagstep.eq.1.) costfac=(costfac-1)/2+1
186                 cost=cost/costfac
187                 flagstep=-1.
188                 goto 7
189               endif
190           ENDIF           ENDIF
191           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion               CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
192  *     *******************************************  *     *******************************************
193  *     find new value of AL-pha                  *     find new value of AL-pha                
194  *                                                *     *******************************************                                        
195           DO I=1,5                                         DO I=1,5                              
196              DAL(I)=0.                        DAL(I)=0.          
197              DO J=1,5                          DO J=1,5            
# Line 161  ',istep,' - ',al Line 215  ',istep,' - ',al
215           ENDDO           ENDDO
216           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
217           IF(IFA.NE.0) THEN           IF(IFA.NE.0) THEN
218              if(TRKDEBUG)then              if(TRKVERBOSE)then
219                 PRINT *,                 PRINT *,
220       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
221       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 173  ',istep,' - ',al Line 227  ',istep,' - ',al
227              RETURN                          RETURN            
228           ENDIF           ENDIF
229           CALL DSFINV(4,CHI2DD_R,4)           CALL DSFINV(4,CHI2DD_R,4)
230    *     *******************************************
231    *     find new value of AL-pha                
232    *     *******************************************                                        
233           DO I=1,4           DO I=1,4
234              DAL(I)=0.              DAL(I)=0.
235              DO J=1,4              DO J=1,4
# Line 185  ',istep,' - ',al Line 242  ',istep,' - ',al
242              AL(I)=AL(I)+DAL(I)              AL(I)=AL(I)+DAL(I)
243           ENDDO           ENDDO
244        ENDIF        ENDIF
245    
246          if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
247    
248  *------------------------------------------------------------*  *------------------------------------------------------------*
249  *     ----------------------------------------------------   *  *     ----------------------------------------------------   *
250  *------------------------------------------------------------*  *------------------------------------------------------------*
 *     *******************************************  
251  *     check parameter bounds:  *     check parameter bounds:
252    *------------------------------------------------------------*
253        DO I=1,5        DO I=1,5
254           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
255              if(TRKDEBUG)then              if(TRKVERBOSE)then
256                 PRINT*,' *** WARNING in mini ***  '                 PRINT*,' *** WARNING in mini ***  '
257                 PRINT*,'MINI_2 ==> AL(',I,') out of range'                 PRINT*,'MINI_2 ==> AL(',I,') out of range'
258                 PRINT*,'         value: ',AL(I),                 PRINT*,'         value: ',AL(I),
# Line 205  ',istep,' - ',al Line 265  ',istep,' - ',al
265              RETURN              RETURN
266           ENDIF           ENDIF
267        ENDDO              ENDDO      
268    *------------------------------------------------------------*
269  *     check number of steps:  *     check number of steps:
270        IF(ISTEP.gt.ISTEPMAX) then  *------------------------------------------------------------*
271          IF(ISTEP.ge.ISTEPMAX) then
272           IFAIL=1           IFAIL=1
273           if(TRKDEBUG)           if(TRKVERBOSE)
274       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
275       $        ISTEPMAX       $        ISTEPMAX
276           goto 11           goto 11
277        endif        endif
278    *------------------------------------------------------------*
279  *     ---------------------------------------------  *     ---------------------------------------------
280  *     evaluate deflection tolerance on the basis of  *     evaluate deflection tolerance on the basis of
281  *     estimated deflection  *     estimated deflection
282  *     ---------------------------------------------  *     ---------------------------------------------
283        ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  *------------------------------------------------------------*
284    c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
285          ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
286          ALTOL(1) = ALTOL(5)/DELETA1
287          ALTOL(2) = ALTOL(1)
288          ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
289          ALTOL(4) = ALTOL(3)        
290          
291  *---- check tolerances:  *---- check tolerances:
292    c$$$      DO I=1,5
293    c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
294    c$$$      ENDDO
295    c$$$      print*,'chi2 -- ',DCHI2
296    
297        DO I=1,5        DO I=1,5
298           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
299        ENDDO        ENDDO
# Line 229  ',istep,' - ',al Line 303  ',istep,' - ',al
303        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
304        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
305           IFAIL=1           IFAIL=1
306           if(TRKDEBUG)THEN           if(TRKVERBOSE)THEN
307              CHI2=-9999.              CHI2=-9999.
308           if(TRKDEBUG)           if(TRKVERBOSE)
309       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
310           ENDIF           ENDIF
311           RETURN           RETURN
# Line 246  ',istep,' - ',al Line 320  ',istep,' - ',al
320        IF(PFIXED.EQ.0.) THEN        IF(PFIXED.EQ.0.) THEN
321           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
322           IF(IFA.NE.0) THEN      !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
323              if(TRKDEBUG)then              if(TRKVERBOSE)then
324                 PRINT *,                 PRINT *,
325       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
326       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 273  ',istep,' - ',al Line 347  ',istep,' - ',al
347           ENDDO           ENDDO
348           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
349           IF(IFA.NE.0) THEN           IF(IFA.NE.0) THEN
350              if(TRKDEBUG)then              if(TRKVERBOSE)then
351                 PRINT *,                 PRINT *,
352       $              '*** ERROR in mini ***'//       $              '*** ERROR in mini ***'//
353       $              'on matrix inversion (not pos-def)'       $              'on matrix inversion (not pos-def)'
# Line 306  ',istep,' - ',al Line 380  ',istep,' - ',al
380        if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***        if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
381        if(ndof.le.0.) then        if(ndof.le.0.) then
382           ndof = 1           ndof = 1
383           if(TRKDEBUG)           if(TRKVERBOSE)
384       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
385        endif        endif
386    
387          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
388    
389  *     ------------------------------------  *     ------------------------------------
390  *     Reduced chi^2  *     Reduced chi^2
391        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
392    
393    c      print*,'mini2: chi2 ',chi2
394    
395   11   CONTINUE         11   CONTINUE      
396    
397        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
# Line 343  ',istep,' - ',al Line 422  ',istep,' - ',al
422       $     ,XV0(nplanes),YV0(nplanes)       $     ,XV0(nplanes),YV0(nplanes)
423        DIMENSION AL_P(5)        DIMENSION AL_P(5)
424    
425        LOGICAL TRKDEBUG  c      LOGICAL TRKVERBOSE
426        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKVERBOSE
427          LOGICAL TRKDEBUG,TRKVERBOSE
428          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
429  *      *    
430  *     chi^2 computation  *     chi^2 computation
431  *      *    
# Line 354  ',istep,' - ',al Line 435  ',istep,' - ',al
435        JFAIL=0                   !error flag        JFAIL=0                   !error flag
436        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
437        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
438           IF(TRKDEBUG)           IF(TRKVERBOSE)
439       $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'       $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
440           IFAIL=1           IFAIL=1
441           RETURN           RETURN
# Line 408  c$$$      ENDDO Line 489  c$$$      ENDDO
489       +        +((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
490       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
491        ENDDO        ENDDO
492    c      print*,'CHISQ ',chi2
493  *     ------------------------------------------------  *     ------------------------------------------------
494  *      *    
495  *     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 508  c$$$      ENDDO
508              JFAIL=0              JFAIL=0
509              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
510              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
511                 IF(TRKDEBUG)                 IF(TRKVERBOSE)
512  *23456789012345678901234567890123456789012345678901234567890123456789012  *23456789012345678901234567890123456789012345678901234567890123456789012
513       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
514                 IFAIL=1                 IFAIL=1
# Line 440  c$$$      ENDDO Line 522  c$$$      ENDDO
522              JFAIL=0              JFAIL=0
523              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
524              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
525                 IF(TRKDEBUG)                 IF(TRKVERBOSE)
526       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'       $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
527                 IFAIL=1                 IFAIL=1
528                 RETURN                 RETURN
# Line 472  c$$$      ENDDO Line 554  c$$$      ENDDO
554                            
555              COSTHE=DSQRT(1.-AL(3)**2)              COSTHE=DSQRT(1.-AL(3)**2)
556              IF(COSTHE.EQ.0.) THEN              IF(COSTHE.EQ.0.) THEN
557                 IF(TRKDEBUG)PRINT *,'=== WARNING ===> COSTHE=0'                 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
558                 IFAIL=1                 IFAIL=1
559                 RETURN                 RETURN
560              ENDIF              ENDIF
# Line 559  c$$$      ENDDO Line 641  c$$$      ENDDO
641        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
642        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
643    
644        LOGICAL TRKDEBUG  c      LOGICAL TRKVERBOSE
645        COMMON/TRKD/TRKDEBUG  c      COMMON/TRKD/TRKVERBOSE
646          LOGICAL TRKDEBUG,TRKVERBOSE
647          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
648  c        c      
649        DIMENSION AL_P(5)        DIMENSION AL_P(5)
650  *      *    
# Line 580  c       Line 664  c      
664        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
665        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))
666        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
667    
668    c$$$      print*,'POSXY ',vout
669    
670        DO I=1,nplanes        DO I=1,nplanes
671           step=vout(3)-zv(i)           step=vout(3)-zv(i)
672   10      DO J=1,7   10      DO J=1,7
# Line 590  c       Line 677  c      
677           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
678           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
679              IFAIL=1              IFAIL=1
680              if(TRKDEBUG)              if(TRKVERBOSE)
681       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
682              if(.TRUE.)print*,'charge',charge  c$$$            if(.TRUE.)print*,'charge',charge
683              if(.TRUE.)print*,'vect',vect  c$$$            if(.TRUE.)print*,'vect',vect
684              if(.TRUE.)print*,'vout',vout  c$$$            if(.TRUE.)print*,'vout',vout
685              if(.TRUE.)print*,'step',step  c$$$            if(.TRUE.)print*,'step',step
686                if(TRKVERBOSE)print*,'charge',charge
687                if(TRKVERBOSE)print*,'vect',vect
688                if(TRKVERBOSE)print*,'vout',vout
689                if(TRKVERBOSE)print*,'step',step
690              RETURN              RETURN
691           ENDIF           ENDIF
692           Z=VOUT(3)           Z=VOUT(3)
# Line 664  c         ZM_B(IP) = 0 Line 755  c         ZM_B(IP) = 0
755    
756        return        return
757        end        end
758    
759    
760    ***************************************************
761    *                                                 *
762    *                                                 *
763    *                                                 *
764    *                                                 *
765    *                                                 *
766    *                                                 *
767    **************************************************
768    
769          subroutine guess()
770    
771    c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
772    
773          include 'commontracker.f' !tracker general common
774          include 'common_mini_2.f' !common for the tracking procedure
775          
776          REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)
777          REAL*4 CHI,XC,ZC,RADIUS
778    *     ----------------------------------------
779    *     Y view
780    *     ----------------------------------------
781    *     ----------------------------------------
782    *     initial guess with a straigth line
783    *     ----------------------------------------
784          SZZ=0.                    
785          SZY=0.
786          SSY=0.
787          SZ=0.
788          S1=0.
789          DO I=1,nplanes
790             IF(YGOOD(I).EQ.1)THEN
791                YY = YM(I)
792                IF(XGOOD(I).EQ.0)THEN
793                   YY = (YM_A(I) + YM_B(I))/2
794                ENDIF
795                SZZ=SZZ+ZM(I)*ZM(I)
796                SZY=SZY+ZM(I)*YY
797                SSY=SSY+YY
798                SZ=SZ+ZM(I)
799                S1=S1+1.
800             ENDIF
801          ENDDO
802          DET=SZZ*S1-SZ*SZ
803          AY=(SZY*S1-SZ*SSY)/DET
804          BY=(SZZ*SSY-SZY*SZ)/DET
805          Y0 = AY*ZINI+BY
806    *     ----------------------------------------
807    *     X view
808    *     ----------------------------------------
809    *     ----------------------------------------
810    *     1) initial guess with a circle
811    *     ----------------------------------------
812          NP=0
813          DO I=1,nplanes
814             IF(XGOOD(I).EQ.1)THEN
815                XX = XM(I)
816                IF(YGOOD(I).EQ.0)THEN
817                   XX = (XM_A(I) + XM_B(I))/2
818                ENDIF
819                NP=NP+1
820                XP(NP)=XX
821                ZP(NP)=ZM(I)
822             ENDIF
823          ENDDO
824          CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
825    c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI
826          IF(IFLAG.NE.0)GOTO 10 !straigth fit
827          ARG = RADIUS**2-(ZINI-ZC)**2
828          IF(ARG.LT.0)GOTO 10       !straigth fit
829          DC = SQRT(ARG)      
830          IF(XC.GT.0)DC=-DC
831          X0=XC+DC
832          AX = -(ZINI-ZC)/DC
833          DEF=100./(RADIUS*0.3*0.43)
834          IF(XC.GT.0)DEF=-DEF
835          GOTO 20                   !guess is ok
836    
837    *     ----------------------------------------
838    *     2) initial guess with a straigth line
839    *     - if circle does not intersect reference plane
840    *     - if bad chi**2
841    *     ----------------------------------------
842     10   CONTINUE
843          SZZ=0.                  
844          SZX=0.
845          SSX=0.
846          SZ=0.
847          S1=0.
848          DO I=1,nplanes
849             IF(XGOOD(I).EQ.1)THEN
850                XX = XM(I)
851                IF(YGOOD(I).EQ.0)THEN
852                   XX = (XM_A(I) + XM_B(I))/2
853                ENDIF
854                SZZ=SZZ+ZM(I)*ZM(I)
855                SZX=SZX+ZM(I)*XX
856                SSX=SSX+XX
857                SZ=SZ+ZM(I)
858                S1=S1+1.
859             ENDIF
860          ENDDO
861          DET=SZZ*S1-SZ*SZ
862          AX=(SZX*S1-SZ*SSX)/DET
863          BX=(SZZ*SSX-SZX*SZ)/DET
864          DEF = 0
865          X0  = AX*ZINI+BX
866    
867     20   CONTINUE
868    *     ----------------------------------------
869    *     guess
870    *     ----------------------------------------
871    
872          AL(1) = X0
873          AL(2) = Y0
874          tath  = sqrt(AY**2+AX**2)
875          AL(3) = tath/sqrt(1+tath**2)
876          IF(AX.NE.0)THEN
877             AL(4)= atan(AY/AX)
878          ELSE
879             AL(4) = acos(-1.)/2
880             IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
881          ENDIF
882          IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
883          AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs
884          AL(5) = DEF
885    
886    c      print*,' guess: ',(al(i),i=1,5)
887    
888          end
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23