/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/analysisflight.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/analysisflight.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by pam-fi, Fri Oct 27 16:08:19 2006 UTC revision 1.20 by mocchiut, Tue Aug 4 14:01:35 2009 UTC
# Line 17  Line 17 
17        include 'calib.f'        include 'calib.f'
18        include 'level2.f'        include 'level2.f'
19    
20    *     input flag
21    *
22    c      integer fin
23    
24  *     flag to chose PFA  *     flag to chose PFA
25        character*10 PFA        character*10 PFA
26        common/FINALPFA/PFA        common/FINALPFA/PFA
# Line 27  c      parameter (inf=1.e8)      !just a Line 31  c      parameter (inf=1.e8)      !just a
31        external npl        external npl
32        external acoordsi,coordsi,nld,coord,dcoord        external acoordsi,coordsi,nld,coord,dcoord
33    
34          integer hitplanex(6)
35          integer hitplaney(6)
36    
37          logical BAD_NUCLEUS
38    
39  ************************************************************  ************************************************************
40  ************************************************************  ************************************************************
41  ************************************************************  ************************************************************
# Line 36  c      parameter (inf=1.e8)      !just a Line 45  c      parameter (inf=1.e8)      !just a
45  ************************************************************  ************************************************************
46  ************************************************************  ************************************************************
47  ************************************************************  ************************************************************
       PFA='ETA'  
48    
49        if(DEBUG)then        call idtoc(pfaid,PFA)
50          
51    
52          if(DEBUG.EQ.1)then
53           print*,'----------------------------------'           print*,'----------------------------------'
54  c         print*,'START... ',good1,nclstr1,nclstrmax_level2           print*,'Settings: '
55             print*,'PFA                  ',pfaid,' ---> ',PFA
56             print*,'tracking mode        ',trackmode
57             print*,'fit-tolerance factor ',fact
58             print*,'minimum n.step       ',istepmin
59            
60        endif        endif
61    
62  *------------------------------------------------------        RECOVER_SINGLETS = .false.
63        call init_level2        RECOVER_NUCLEI   = .false.
64  *------------------------------------------------------        SECOND_SEARCH    = .false.
         
 *------------------------------------------------------  
 *     cut on maximum number of clusters  
 *------------------------------------------------------  
 c$$$      if(nclstr1.gt.nclstrmax_level2)then  
 c$$$         goto 8800              !fill nt-uple and go to next event  
 c$$$      endif  
65                
66          ntrk_old = 0
67    
68        do i=1,nclstr1        do i=1,nclstr1
69           cl_used(i)=0           !init mask of clusters associated to a track           cl_used(i)=0           !init mask of clusters associated to a track
70        enddo        enddo
71                dedx_x_min = dedx_x_min_mip
72        if(DEBUG)then        dedx_y_min = dedx_y_min_mip
73           print*,'----------------------------------'  *------------------------------------------------------
74           print*,iev,'   ** ',nev2        call init_level2
75        endif        call init_hough      
76    *------------------------------------------------------
77    
78     888  continue
79    
80        
81  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
82  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
83    
# Line 72  c$$$      endif Line 87  c$$$      endif
87           goto 880               !fill ntp and go to next event               goto 880               !fill ntp and go to next event    
88        endif        endif
89    
90    *     ----------------------------------------------------
91    *     fill hough output variables only for the first searc
92    *     ----------------------------------------------------
93          if(.not.SECOND_SEARCH)call fill_hough
94    
95        iflag=0        iflag=0
96        call track_fitting(iflag)        call track_fitting(iflag)
       if(iflag.eq.1)then        !bad event  
          goto 880               !fill ntp and go to next event      
       endif  
97                
98        
99    *     ////////////////////////////////////////////////
100    *     PATCH to recover events with less than 3 couples
101    *     ////////////////////////////////////////////////    
102    *     if no tracks have been found
103    *     and n.planes with at least one couple (planehit)
104    *     is less than 3 try to recover applying the hough
105    *     transform to singlet
106          nplanehit=0
107          do ip=1,NPLANES
108             if(ncp_plane(ip).gt.0)nplanehit=nplanehit+1
109          enddo
110          if(DEBUG.eq.1)then
111    cc      if(.true.)then
112             print*,''
113             print*,'-----------------------------------'
114             print*,'n.tracks     ',ntrk
115             print*,'n.hit-planes ',nplanehit
116          endif
117          if(
118         $     .not.RECOVER_SINGLETS.and. !recover only once
119         $     ntrk.eq.0            .and. !no tracks has been found
120         $     nplanehit.le.3       .and. !if less than 3 hit planes (n.real couples)
121         $     .true.)then
122             RECOVER_SINGLETS=.true.
123    cc         if(.true.)
124             if(DEBUG.EQ.1)
125         $        print*,' -----> ((( TRY TO RECOVER SINGLETS )))'
126             goto 888               !back to track finding
127          endif
128    
129    
130    *     ////////////////////////////////////////////////
131    *     PATCH to recover nuclei
132    *     ////////////////////////////////////////////////
133          itrk = ntrk               !check the last track first
134          isimage = 0
135          BAD_NUCLEUS=.false.
136     150  continue
137    
138          dedxav  = 0.
139          npt     = 0
140          dedxmin = 1000000
141          if(itrk.gt.0)then
142             do ip=1,NPLANES
143    cc            print*,ip,' ** ',dedx_x(ip,itrk),dedx_y(ip,itrk)
144                if(xgood_nt(ip,itrk).eq.1)then
145                   dedxav = dedxav + abs(dedx_x(ip,itrk))
146                   if(abs(dedx_x(ip,itrk)).lt.dedxmin)
147         $              dedxmin=abs(dedx_x(ip,itrk))
148                   npt = npt + 1
149                endif        
150                if(ygood_nt(ip,itrk).eq.1)then
151                   dedxav = dedxav + abs(dedx_y(ip,itrk))
152                   if(abs(dedx_y(ip,itrk)).lt.dedxmin)
153         $              dedxmin=abs(dedx_y(ip,itrk))
154                   npt = npt + 1
155                endif            
156             enddo
157             if(npt.gt.0)dedxav = dedxav/npt
158             if(  dedxav.gt.dedx_min_nuclei )then
159                if( (dedxav-dedxmin)/dedxav.gt.ddedx_min_nuclei )
160         $           BAD_NUCLEUS=.true.
161             endif
162          endif
163          if(DEBUG.eq.1)then
164    cc      if(.true.)then
165             print*,''
166             print*,'-----------------------------------'
167             print*,'itrk          ',itrk
168             print*,'dedx(average) ',dedxav,' cut ',dedx_min_nuclei
169             print*,'dedx(minimum) ',dedxmin
170             if(dedxav.gt.0)print*,'(av-min)/av',(dedxav-dedxmin)/dedxav
171         $        ,' cut ',ddedx_min_nuclei
172             print*,'BAD_NUCLEUS ? ',BAD_NUCLEUS
173          endif
174    
175          itrk=0
176          if(ntrk.gt.0)itrk = image(ntrk)              
177          if(itrk.gt.0.and.isimage.eq.0)then      !hence check the image
178             isimage=1        
179             goto 150
180          endif
181    
182    
183          if(
184         $     .not.RECOVER_NUCLEI.and. !recover only once
185         $     (ntrk.eq.0.or.BAD_NUCLEUS).and. !no tracks found or bad nucleus
186         $     .true.)then
187    
188             RECOVER_NUCLEI   = .true.
189             dedx_x_min = dedx_x_min_nuclei
190             dedx_y_min = dedx_y_min_nuclei
191    cc         if(.true.)
192             if(DEBUG.EQ.1)
193         $     print*,' -----> ((( TRY TO RECOVER NUCLEI )))'
194    *     ----------------------------------------------
195    *     in case of re-tracking, un-tag used clusters
196    *     ----------------------------------------------
197             if(ntrk.ne.0)then
198                do ip=1,NPLANES
199                   if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
200                   if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
201                enddo
202                ntrk = ntrk-1       !decrement track
203                if(isimage.eq.1)then
204                   do ip=1,NPLANES
205                      if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
206                      if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
207                   enddo
208                   ntrk = ntrk-1    !decrement track image
209                endif
210             endif
211             goto 888               !back to track finding
212            
213          endif
214    *     ----------------------------------------------
215    *     back to default parameters for the second-track search
216    *     ----------------------------------------------
217          if(RECOVER_NUCLEI)then
218             dedx_x_min = dedx_x_min_mip
219             dedx_y_min = dedx_y_min_mip
220             RECOVER_NUCLEI   = .false.    
221          endif
222    
223          if(ntrk.eq.0)goto 880
224    
225    *     /////////////////////////////////////////////////////
226    *     PATCH to perform a more efficient second-track search
227    *     /////////////////////////////////////////////////////
228    *     -------------------------------------------------
229    *     n. hit planes and number of clusters
230    *     -------------------------------------------------
231          do ip=1,NPLANES
232             hitplanex(ip)=0
233             hitplaney(ip)=0
234             nclusters=0
235          enddo
236          do icl=1,nclstr1
237             ip = npl(VIEW(icl))
238             if(cl_good(icl).eq.0.or.cl_used(icl).eq.1)goto 9
239             if(mod(VIEW(icl),2).eq.1)hitplaney(ip)=hitplaney(ip)+1 !YY        
240             if(mod(VIEW(icl),2).eq.0)hitplanex(ip)=hitplanex(ip)+1 !XX        
241             nclusters=nclusters+1
242     9       continue
243          enddo
244          nhitplaney=0
245          nhitplanex=0
246          do ip=1,NPLANES
247             if(hitplaney(ip).gt.0)nhitplaney=nhitplaney+1
248             if(hitplanex(ip).gt.0)nhitplanex=nhitplanex+1
249          enddo      
250    *     -------------------------------------------------
251    *     chi2 cut
252    *     -------------------------------------------------
253          p0 = 1.111588e+00
254          p1 = 1.707656e+00
255          p2 = 1.489693e-01
256          chi2m025 = p0
257         $     + abs(al_nt(5,ntrk))*p1
258         $     + al_nt(5,ntrk)*al_nt(5,ntrk)*p2
259          chi2m025_0 = p0 + 0.07*p1 + 0.07*0.07*p2
260    cc      chi2m = (chi2m025-chi2m025_0+12**0.25)**4. !95%
261          chi2m = (chi2m025-chi2m025_0+20.**0.25)**4. !???
262    *     -------------------------------------------------
263    *     image track
264    *     -------------------------------------------------
265          ntrk_best = ntrk
266          chi2_best = chi2_nt(ntrk)
267          if(image(ntrk).ne.0)
268         $     chi2_best = min(chi2_nt(ntrk),chi2_nt(image(ntrk)))
269          if(DEBUG.eq.1)then
270    cc      if(.true.)then
271             print*,''
272             print*,'-----------------------------------'
273             print*,'chi2_nt(',ntrk_best,') = ',chi2_best,' < ',chi2m,' ?'
274             print*,'nhitplanex',nhitplanex
275             print*,'nhitplaney',nhitplaney
276             print*,'nclusters',nclusters
277          endif
278          if(
279         $     ntrk.gt.ntrk_old.and. !one track must have been already found
280         $     chi2_best.gt.0 .and.     !this track must not be garbage
281         $     chi2_best.lt.chi2m .and. !"
282         $     nhitplaney.ge.3.and. !a minimum number of hit planes is required
283         $     nhitplanex.ge.3.and. !"
284         $     nclusters.le.20.and. !limit on number of clusters to save time
285         $     .true.)then
286    
287             SECOND_SEARCH=.true.
288    c         RECOVER_SINGLETS=.true.      !OPTIMIZATION  
289             ntrk_old=ntrk
290    
291             if(DEBUG.EQ.1)
292    cc         if(.true.)
293         $        print*,'***** NEW SEARCH ****'
294             goto 888               !back to track finding
295    
296          endif
297    
298   880  continue   880  continue
299    
# Line 88  c$$$      endif Line 303  c$$$      endif
303    
304        call fill_level2_siglets        call fill_level2_siglets
305                
306        if(DEBUG)then        if(DEBUG.EQ.1)then
               
307           print*,''           print*,''
308           print*,'DONE!'           print*,'DONE!'
309           print*,''           print*,''
# Line 99  c$$$      endif Line 313  c$$$      endif
313           print*,''           print*,''
314           print*,''           print*,''
315        endif        endif
316    
317          ngood = 0
318          do iv = 1,nviews
319             ngood = ngood + good1(iv)
320          enddo
321                
322   8800 continue  c 8800 continue
323          continue
324                
325                
326  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
327  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
328    
329   100  continue  c 100  continue
330          continue
331                
332        return        return
333        end        end
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23