/[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.16 by pam-fi, Mon Aug 20 16:07:16 2007 UTC revision 1.17 by pam-fi, Fri Dec 5 08:26:47 2008 UTC
# Line 31  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 40  c      parameter (inf=1.e8)      !just a Line 45  c      parameter (inf=1.e8)      !just a
45  ************************************************************  ************************************************************
46  ************************************************************  ************************************************************
47  ************************************************************  ************************************************************
 c$$$      TRACKMODE = 0  
 c$$$      FACT = 100.  
 c$$$      ISTEPMIN = 3  
48    
49        call idtoc(pfaid,PFA)        call idtoc(pfaid,PFA)
50                
 c$$$      PFA='COG4'  
 c$$$      if(pfaid.eq.0)PFA='ETA'  
 c$$$      if(pfaid.eq.2)PFA='ETA2'  
 c$$$      if(pfaid.eq.3)PFA='ETA3'  
 c$$$      if(pfaid.eq.4)PFA='ETA4'  
 c$$$      if(pfaid.eq.10)PFA='COG'  
 c$$$      if(pfaid.eq.11)PFA='COG1'  
 c$$$      if(pfaid.eq.12)PFA='COG2'  
 c$$$      if(pfaid.eq.13)PFA='COG3'  
 c$$$      if(pfaid.eq.14)PFA='COG4'  
 ***********************************************************  
   
 c      if(DEBUG.EQ.1)PRINT*,'P.F.A. --> ',fin,PFA  
51    
52        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
53           print*,'----------------------------------'           print*,'----------------------------------'
# Line 70  c      if(DEBUG.EQ.1)PRINT*,'P.F.A. --> Line 59  c      if(DEBUG.EQ.1)PRINT*,'P.F.A. -->
59                    
60        endif        endif
61    
62  *------------------------------------------------------        RECOVER_SINGLETS = .false.
63        call init_level2        RECOVER_NUCLEI   = .false.
64        call init_hough        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  c$$$      if(DEBUG.EQ.1)then        dedx_y_min = dedx_y_min_mip
73  c$$$         print*,'----------------------------------'  *------------------------------------------------------
74  c$$$         print*,iev,'   ** ',nev2        call init_level2
75  c$$$      endif        call init_hough      
76    *------------------------------------------------------
77    
78     888  continue
79    
80        
81  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
82  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
83    
# Line 100  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        call fill_hough  *     -------------------------------------
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             RECOVER_SINGLETS=.true.
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 119  c$$$      endif Line 304  c$$$      endif
304        call fill_level2_siglets        call fill_level2_siglets
305                
306        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
               
307           print*,''           print*,''
308           print*,'DONE!'           print*,'DONE!'
309           print*,''           print*,''
# Line 134  c$$$      endif Line 318  c$$$      endif
318        do iv = 1,nviews        do iv = 1,nviews
319           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
320        enddo        enddo
 c$$$      if(ngood.ne.0)print*,'* WARNING * Event '  
 c$$$     $     ,':LEVEL2 event status: '  
 c$$$     $     ,(good2(i),i=1,nviews)  
321                
322   8800 continue   8800 continue
323                
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.23