/[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.2 by pam-fi, Tue May 30 16:30:37 2006 UTC revision 1.9 by pam-fi, Tue Nov 21 17:12:26 2006 UTC
# Line 12  Line 12 
12  ************************************************************************    ************************************************************************  
13    
14    
15        SUBROUTINE MINI_2(ISTEP,IFAIL)        SUBROUTINE MINI2(ISTEP,IFAIL,IPRINT)
16                
17        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
18    
# 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 31  c     the plane ordering is reversed in Line 31  c     the plane ordering is reversed in
31  c     ordering, but they maintain their Z coordinates. so plane number 1 is  c     ordering, but they maintain their Z coordinates. so plane number 1 is
32  c     the first one that a particle meets, and its Z coordinate is > 0  c     the first one that a particle meets, and its Z coordinate is > 0
33  c------------------------------------------------------------------------  c------------------------------------------------------------------------
34        DATA ZINI/23.5/           !z coordinate of the reference plane        DATA ZINI/23.5/  !!! ***PP*** to be changed !z coordinate of the reference plane
35    
36        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/100/        !maximum number of steps in the chi^2 minimization
# Line 41  c--------------------------------------- Line 41  c---------------------------------------
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,dinf,dinf,dinf/ !limits on alpha vector components
47        DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !"        DATA ALMIN/-dinf,-dinf,-dinf,-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
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 59  c Line 63  c
63  c     NB: the two metods gives equivalent results BUT  c     NB: the two metods gives equivalent results BUT
64  c     method 2 is faster!!  c     method 2 is faster!!
65  c--------------------------------------------------------  c--------------------------------------------------------
66        DATA IFLAG/2/                    DATA IFLAG/2/  
67            
68    c      LOGICAL TRKDEBUG,TRKVERBOSE
69    c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
70          LOGICAL TRKDEBUG,TRKVERBOSE
71          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
72    
73          IF(IPRINT.EQ.1) THEN
74             TRKVERBOSE = .TRUE.
75             TRKDEBUG   = .FALSE.
76          ELSEIF(IPRINT.EQ.2)THEN
77             TRKVERBOSE = .TRUE.
78             TRKDEBUG   = .TRUE.
79          ELSE
80             TRKVERBOSE = .FALSE.
81             TRKDEBUG   = .FALSE.
82          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        FACT=100.                  !scale factor to define tolerance on alfa
       ALTOL(2)=RESY(1)/FACT     !al(2) = y  
       ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta)  
      $     +RESY(1)**2)/44.51/FACT  
       ALTOL(4)=ALTOL(3)         !al(4)=phi  
109  c     deflection error (see PDG)  c     deflection error (see PDG)
110        DELETA1=0.01*RESX(1)/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.))
111        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)
112    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)
118    c$$$      DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
119    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        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  
125        IF(JFAIL.NE.0) THEN        if(TRKDEBUG) print*,'guess: ',al
126           IFAIL=1        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
127           if(DEBUG)  
      $        PRINT *,'mini_2 ===> error on CHISQ computation !!!'  
          RETURN  
       ENDIF  
128  *      *    
129  *     -----------------------  *     -----------------------
130  *     START MINIMIZATION LOOP  *     START MINIMIZATION LOOP
131  *     -----------------------  *     -----------------------
132   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!
       CHI2_P=CHI2  
         
 c      print*,'@@@@@ ',istep,' - ',al  
133    
134        cost=1e-7        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
135          IF(JFAIL.NE.0) THEN
136             IFAIL=1
137             CHI2=-9999.
138             if(TRKVERBOSE)
139         $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
140             RETURN
141          ENDIF
142          
143          COST=1e-7
144        DO I=1,5        DO I=1,5
145           DO J=1,5           DO J=1,5
146              CHI2DD(I,J)=CHI2DD(I,J)*COST              CHI2DD(I,J)=CHI2DD(I,J)*COST
147           ENDDO           ENDDO
148           CHI2D(I)=CHI2D(I)*COST           CHI2D(I)=CHI2D(I)*COST
149        ENDDO              ENDDO      
150    
151          IF(PFIXED.EQ.0.) THEN
152    
153  *------------------------------------------------------------*  *------------------------------------------------------------*
154  *     track fitting with FREE deflection  *     track fitting with FREE deflection
155  *------------------------------------------------------------*  *------------------------------------------------------------*
156        CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
157        IF(IFA.NE.0) THEN         !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
158           if(DEBUG)then              if(TRKVERBOSE)then
159              PRINT *,                 PRINT *,
160       $     'MINI_HOUGH ==> '//       $              '*** ERROR in mini ***'//
161       $     '** ERROR ** on matrix inversion (not positive-defined)!!!'       $              'on matrix inversion (not pos-def)'
162       $     ,DET       $              ,DET
163           endif              endif
164           IFAIL=1              IF(CHI2.EQ.0) CHI2=-9999.
165           RETURN                          IF(CHI2.GT.0) CHI2=-CHI2
166        ENDIF              IFAIL=1
167        CALL DSFINV(5,CHI2DD,5)   !CHI2DD matrix inversion                  RETURN            
168             ENDIF
169             CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
170  *     *******************************************  *     *******************************************
171  *     find new value of AL-pha                 !*  *     find new value of AL-pha                
172  *                                              !*  *     *******************************************                                        
173        DO I=1,5                                 !*           DO I=1,5                              
174           DAL(I)=0.                             !*              DAL(I)=0.          
175           DO J=1,5                              !*              DO J=1,5            
176              DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) !*                 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)
177              COV(I,J)=CHI2DD(I,J)               !*                 COV(I,J)=2.*COST*CHI2DD(I,J)
178           ENDDO                                 !*              ENDDO              
179        ENDDO                                    !*           ENDDO                  
180        DO I=1,5                                 !*           DO I=1,5              
181           AL(I)=AL(I)+DAL(I)                    !*              AL(I)=AL(I)+DAL(I)  
182        ENDDO                                    !*           ENDDO                  
183    *------------------------------------------------------------*
184    *     track fitting with FIXED deflection
185    *------------------------------------------------------------*
186          ELSE
187             AL(5)=1./PFIXED
188             DO I=1,4
189                CHI2D_R(I)=CHI2D(I)
190                DO J=1,4
191                   CHI2DD_R(I,J)=CHI2DD(I,J)
192                ENDDO
193             ENDDO
194             CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
195             IF(IFA.NE.0) THEN
196                if(TRKVERBOSE)then
197                   PRINT *,
198         $              '*** ERROR in mini ***'//
199         $              'on matrix inversion (not pos-def)'
200         $              ,DET
201                endif
202                IF(CHI2.EQ.0) CHI2=-9999.
203                IF(CHI2.GT.0) CHI2=-CHI2
204                IFAIL=1
205                RETURN            
206             ENDIF
207             CALL DSFINV(4,CHI2DD_R,4)
208  *     *******************************************  *     *******************************************
209    *     find new value of AL-pha                
210    *     *******************************************                                        
211             DO I=1,4
212                DAL(I)=0.
213                DO J=1,4
214                   DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)
215                   COV(I,J)=2.*COST*CHI2DD_R(I,J)
216                ENDDO
217             ENDDO
218             DAL(5)=0.
219             DO I=1,4
220                AL(I)=AL(I)+DAL(I)
221             ENDDO
222          ENDIF
223    
224          if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
225    
226    *------------------------------------------------------------*
227    *     ----------------------------------------------------   *
228    *------------------------------------------------------------*
229  *     check parameter bounds:  *     check parameter bounds:
230    *------------------------------------------------------------*
231        DO I=1,5        DO I=1,5
232           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
233              if(DEBUG)then              if(TRKVERBOSE)then
234                 PRINT*,' **WARNING**  '                 PRINT*,' *** WARNING in mini ***  '
235                 PRINT*,'MINI_2 ==> AL(',I,') out of range'                 PRINT*,'MINI_2 ==> AL(',I,') out of range'
236                 PRINT*,'         value: ',AL(I),                 PRINT*,'         value: ',AL(I),
237       $              '  limits: ',ALMIN(I),ALMAX(I)                     $              '  limits: ',ALMIN(I),ALMAX(I)              
238                 print*,'istep ',istep                 print*,'istep ',istep
239              endif              endif
240                IF(CHI2.EQ.0) CHI2=-9999.
241                IF(CHI2.GT.0) CHI2=-CHI2
242              IFAIL=1              IFAIL=1
243              RETURN              RETURN
244           ENDIF           ENDIF
245        ENDDO              ENDDO      
246  *     new estimate of chi^2:  *------------------------------------------------------------*
       JFAIL=0                   !error flag  
       CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  
       IF(JFAIL.NE.0) THEN  
          IFAIL=1  
          if(DEBUG)  
      $        PRINT *,'mini_2:  ===> error on CHISQ computation !!!'  
          RETURN  
       ENDIF  
247  *     check number of steps:  *     check number of steps:
248        IF(ISTEP.gt.ISTEPMAX) then  *------------------------------------------------------------*
249           IFAIL=1        IF(ISTEP.ge.ISTEPMAX) then
250           if(DEBUG)  c$$$         IFAIL=1
251       $        PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX  c$$$         if(TRKVERBOSE)
252    c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
253    c$$$     $        ISTEPMAX
254           goto 11           goto 11
255        endif        endif
256    *------------------------------------------------------------*
257  *     ---------------------------------------------  *     ---------------------------------------------
258  *     evaluate deflection tolerance on the basis of  *     evaluate deflection tolerance on the basis of
259  *     estimated deflection  *     estimated deflection
260  *     ---------------------------------------------  *     ---------------------------------------------
261        ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT  *------------------------------------------------------------*
262    c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
263          ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
264          ALTOL(1) = ALTOL(5)/DELETA1
265          ALTOL(2) = ALTOL(1)
266          ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
267          ALTOL(4) = ALTOL(3)        
268          
269  *---- check tolerances:  *---- check tolerances:
270    c$$$      DO I=1,5
271    c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
272    c$$$      ENDDO
273    c$$$      print*,'chi2 -- ',DCHI2
274    
275        DO I=1,5        DO I=1,5
276           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
277        ENDDO        ENDDO
278    
279    *     new estimate of chi^2:
280          JFAIL=0                   !error flag
281          CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
282          IF(JFAIL.NE.0) THEN
283             IFAIL=1
284             if(TRKVERBOSE)THEN
285                CHI2=-9999.
286             if(TRKVERBOSE)
287         $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
288             ENDIF
289             RETURN
290          ENDIF
291          COST=1e-7
292          DO I=1,5
293             DO J=1,5
294                CHI2DD(I,J)=CHI2DD(I,J)*COST
295             ENDDO
296             CHI2D(I)=CHI2D(I)*COST
297          ENDDO      
298          IF(PFIXED.EQ.0.) THEN
299             CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
300             IF(IFA.NE.0) THEN      !not positive-defined      
301                if(TRKVERBOSE)then
302                   PRINT *,
303         $              '*** ERROR in mini ***'//
304         $              'on matrix inversion (not pos-def)'
305         $              ,DET
306                endif
307                IF(CHI2.EQ.0) CHI2=-9999.
308                IF(CHI2.GT.0) CHI2=-CHI2
309                IFAIL=1
310                RETURN            
311             ENDIF
312             CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
313             DO I=1,5                              
314                DAL(I)=0.          
315                DO J=1,5            
316                   COV(I,J)=2.*COST*CHI2DD(I,J)
317                ENDDO              
318             ENDDO                  
319          ELSE
320             DO I=1,4
321                CHI2D_R(I)=CHI2D(I)
322                DO J=1,4
323                   CHI2DD_R(I,J)=CHI2DD(I,J)
324                ENDDO
325             ENDDO
326             CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
327             IF(IFA.NE.0) THEN
328                if(TRKVERBOSE)then
329                   PRINT *,
330         $              '*** ERROR in mini ***'//
331         $              'on matrix inversion (not pos-def)'
332         $              ,DET
333                endif
334                IF(CHI2.EQ.0) CHI2=-9999.
335                IF(CHI2.GT.0) CHI2=-CHI2
336                IFAIL=1
337                RETURN            
338             ENDIF
339             CALL DSFINV(4,CHI2DD_R,4)
340             DO I=1,4
341                DAL(I)=0.
342                DO J=1,4
343                   COV(I,J)=2.*COST*CHI2DD_R(I,J)
344                ENDDO
345             ENDDO
346          ENDIF
347    *****************************
348    
349  *     ------------------------------------  *     ------------------------------------
350  *     Number of Degree Of Freedom  *     Number of Degree Of Freedom
# Line 180  ',istep,' - ',al Line 354  ',istep,' - ',al
354       $        +int(xgood(ip))       $        +int(xgood(ip))
355       $        +int(ygood(ip))       $        +int(ygood(ip))
356        enddo        enddo
357        ndof=ndof-5                                  if(pfixed.eq.0.) ndof=ndof-5 ! ***PP***                          
358          if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
359          if(ndof.le.0.) then
360             ndof = 1
361             if(TRKVERBOSE)
362         $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
363          endif
364    
365          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
366    
367  *     ------------------------------------  *     ------------------------------------
368  *     Reduced chi^2  *     Reduced chi^2
369        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
370    
371          c      print*,'mini2: chi2 ',chi2
  11   CONTINUE        
372    
373   101  CONTINUE   11   CONTINUE      
374    
375  c      print*,'END MINI'        NSTEP=ISTEP ! ***PP***
376    
377        RETURN        RETURN
378        END                END        
# Line 217  c      print*,'END MINI' Line 399  c      print*,'END MINI'
399        DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)        DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)
400       $     ,XV0(nplanes),YV0(nplanes)       $     ,XV0(nplanes),YV0(nplanes)
401        DIMENSION AL_P(5)        DIMENSION AL_P(5)
402    
403    c      LOGICAL TRKVERBOSE
404    c      COMMON/TRKD/TRKVERBOSE
405          LOGICAL TRKDEBUG,TRKVERBOSE
406          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
407  *      *    
408  *     chi^2 computation  *     chi^2 computation
409  *      *    
# Line 226  c      print*,'END MINI' Line 413  c      print*,'END MINI'
413        JFAIL=0                   !error flag        JFAIL=0                   !error flag
414        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
415        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
416           PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'           IF(TRKVERBOSE)
417         $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
418           IFAIL=1           IFAIL=1
419           RETURN           RETURN
420        ENDIF        ENDIF
# Line 279  c$$$      ENDDO Line 467  c$$$      ENDDO
467       +        +((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
468       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
469        ENDDO        ENDDO
470    c      print*,'CHISQ ',chi2
471  *     ------------------------------------------------  *     ------------------------------------------------
472  *      *    
473  *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)      *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)    
# Line 297  c$$$      ENDDO Line 486  c$$$      ENDDO
486              JFAIL=0              JFAIL=0
487              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
488              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
489                 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'                 IF(TRKVERBOSE)
490    *23456789012345678901234567890123456789012345678901234567890123456789012
491         $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
492                 IFAIL=1                 IFAIL=1
493                 RETURN                 RETURN
494              ENDIF              ENDIF
# Line 309  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                 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'                 IF(TRKVERBOSE)
504         $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
505                 IFAIL=1                 IFAIL=1
506                 RETURN                 RETURN
507              ENDIF              ENDIF
# Line 340  c$$$      ENDDO Line 532  c$$$      ENDDO
532                            
533              COSTHE=DSQRT(1.-AL(3)**2)              COSTHE=DSQRT(1.-AL(3)**2)
534              IF(COSTHE.EQ.0.) THEN              IF(COSTHE.EQ.0.) THEN
535                 PRINT *,'=== WARNING ===> COSTHE=0'                 IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
536                 STOP                 IFAIL=1
537                   RETURN
538              ENDIF              ENDIF
539                            
540              DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3              DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
# Line 425  c$$$      ENDDO Line 618  c$$$      ENDDO
618                
619        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
620        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
621    
622    c      LOGICAL TRKVERBOSE
623    c      COMMON/TRKD/TRKVERBOSE
624          LOGICAL TRKDEBUG,TRKVERBOSE
625          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
626  c        c      
627        DIMENSION AL_P(5)        DIMENSION AL_P(5)
628  *      *    
# Line 444  c       Line 642  c      
642        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
643        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))
644        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
645    
646    c$$$      print*,'POSXY ',vout
647    
648        DO I=1,nplanes        DO I=1,nplanes
649           step=vout(3)-zv(i)           step=vout(3)-zv(i)
650   10      DO J=1,7   10      DO J=1,7
# Line 454  c       Line 655  c      
655           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
656           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
657              IFAIL=1              IFAIL=1
658              if(WARNING)              if(TRKVERBOSE)
659       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
660              if(WARNING)print*,'charge',charge  c$$$            if(.TRUE.)print*,'charge',charge
661              if(WARNING)print*,'vect',vect  c$$$            if(.TRUE.)print*,'vect',vect
662              if(WARNING)print*,'vout',vout  c$$$            if(.TRUE.)print*,'vout',vout
663              if(WARNING)print*,'step',step  c$$$            if(.TRUE.)print*,'step',step
664                if(TRKVERBOSE)print*,'charge',charge
665                if(TRKVERBOSE)print*,'vect',vect
666                if(TRKVERBOSE)print*,'vout',vout
667                if(TRKVERBOSE)print*,'step',step
668              RETURN              RETURN
669           ENDIF           ENDIF
670           Z=VOUT(3)           Z=VOUT(3)
# Line 528  c         ZM_B(IP) = 0 Line 733  c         ZM_B(IP) = 0
733    
734        return        return
735        end        end
736    
737    
738    ***************************************************
739    *                                                 *
740    *                                                 *
741    *                                                 *
742    *                                                 *
743    *                                                 *
744    *                                                 *
745    **************************************************
746    
747          subroutine guess()
748    
749    c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
750    
751          include 'commontracker.f' !tracker general common
752          include 'common_mini_2.f' !common for the tracking procedure
753          
754          REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)
755          REAL*4 CHI,XC,ZC,RADIUS
756    *     ----------------------------------------
757    *     Y view
758    *     ----------------------------------------
759    *     ----------------------------------------
760    *     initial guess with a straigth line
761    *     ----------------------------------------
762          SZZ=0.                    
763          SZY=0.
764          SSY=0.
765          SZ=0.
766          S1=0.
767          DO I=1,nplanes
768             IF(YGOOD(I).EQ.1)THEN
769                YY = YM(I)
770                IF(XGOOD(I).EQ.0)THEN
771                   YY = (YM_A(I) + YM_B(I))/2
772                ENDIF
773                SZZ=SZZ+ZM(I)*ZM(I)
774                SZY=SZY+ZM(I)*YY
775                SSY=SSY+YY
776                SZ=SZ+ZM(I)
777                S1=S1+1.
778             ENDIF
779          ENDDO
780          DET=SZZ*S1-SZ*SZ
781          AY=(SZY*S1-SZ*SSY)/DET
782          BY=(SZZ*SSY-SZY*SZ)/DET
783          Y0 = AY*ZINI+BY
784    *     ----------------------------------------
785    *     X view
786    *     ----------------------------------------
787    *     ----------------------------------------
788    *     1) initial guess with a circle
789    *     ----------------------------------------
790          NP=0
791          DO I=1,nplanes
792             IF(XGOOD(I).EQ.1)THEN
793                XX = XM(I)
794                IF(YGOOD(I).EQ.0)THEN
795                   XX = (XM_A(I) + XM_B(I))/2
796                ENDIF
797                NP=NP+1
798                XP(NP)=XX
799                ZP(NP)=ZM(I)
800             ENDIF
801          ENDDO
802          IFLAG=0                   !no debug mode
803          CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
804    c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
805          IF(IFLAG.NE.0)GOTO 10 !straigth fit
806          if(CHI.gt.100)GOTO 10 !straigth fit
807          ARG = RADIUS**2-(ZINI-ZC)**2
808          IF(ARG.LT.0)GOTO 10       !straigth fit
809          DC = SQRT(ARG)      
810          IF(XC.GT.0)DC=-DC
811          X0=XC+DC
812          AX = -(ZINI-ZC)/DC
813          DEF=100./(RADIUS*0.3*0.43)
814          IF(XC.GT.0)DEF=-DEF
815          
816          IF(ABS(X0).GT.30)THEN
817             PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
818         $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
819             GOTO 10       !straigth fit
820          ENDIF
821          GOTO 20                   !guess is ok
822    
823    *     ----------------------------------------
824    *     2) initial guess with a straigth line
825    *     - if circle does not intersect reference plane
826    *     - if bad chi**2
827    *     ----------------------------------------
828     10   CONTINUE
829          SZZ=0.                  
830          SZX=0.
831          SSX=0.
832          SZ=0.
833          S1=0.
834          DO I=1,nplanes
835             IF(XGOOD(I).EQ.1)THEN
836                XX = XM(I)
837                IF(YGOOD(I).EQ.0)THEN
838                   XX = (XM_A(I) + XM_B(I))/2
839                ENDIF
840                SZZ=SZZ+ZM(I)*ZM(I)
841                SZX=SZX+ZM(I)*XX
842                SSX=SSX+XX
843                SZ=SZ+ZM(I)
844                S1=S1+1.
845             ENDIF
846          ENDDO
847          DET=SZZ*S1-SZ*SZ
848          AX=(SZX*S1-SZ*SSX)/DET
849          BX=(SZZ*SSX-SZX*SZ)/DET
850          DEF = 0
851          X0  = AX*ZINI+BX
852    
853     20   CONTINUE
854    *     ----------------------------------------
855    *     guess
856    *     ----------------------------------------
857    
858          AL(1) = X0
859          AL(2) = Y0
860          tath  = sqrt(AY**2+AX**2)
861          AL(3) = tath/sqrt(1+tath**2)
862          IF(AX.NE.0)THEN
863             AL(4)= atan(AY/AX)
864          ELSE
865             AL(4) = acos(-1.)/2
866             IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
867          ENDIF
868          IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
869          AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs
870          AL(5) = DEF
871    
872    c      print*,' guess: ',(al(i),i=1,5)
873    
874          end

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9

  ViewVC Help
Powered by ViewVC 1.1.23