/[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.12 by pam-fi, Thu Mar 15 12:17:10 2007 UTC revision 1.22 by mocchiut, Mon Feb 15 20:47:11 2010 UTC
# Line 6  Line 6 
6  *     - create LEVEL2 ntuple  *     - create LEVEL2 ntuple
7  *  *
8  *************************************************************************  *************************************************************************
9        subroutine analysisflight(fin)        subroutine analysisflight
10    
11        include 'commontracker.f'        include 'commontracker.f'
12        include 'level1.f'        include 'level1.f'
# Line 19  Line 19 
19    
20  *     input flag  *     input flag
21  *  *
22  *     0 - PFA = ETA  c      integer fin
 *     1 - PFA = COG4  
       integer fin  
23    
24  *     flag to chose PFA  *     flag to chose PFA
25        character*10 PFA        character*10 PFA
# Line 33  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 42  c      parameter (inf=1.e8)      !just a Line 45  c      parameter (inf=1.e8)      !just a
45  ************************************************************  ************************************************************
46  ************************************************************  ************************************************************
47  ************************************************************  ************************************************************
       PFA='ETA'  
       TRACKMODE = 0  
       if(fin.eq.0)PFA='ETA'  
       if(fin.eq.2)PFA='ETA2'  
       if(fin.eq.3)PFA='ETA3'  
       if(fin.eq.4)PFA='ETA4'  
       if(fin.eq.10)PFA='COG'  
       if(fin.eq.11)PFA='COG1'  
       if(fin.eq.12)PFA='COG2'  
       if(fin.eq.13)PFA='COG3'  
       if(fin.eq.14)PFA='COG4'  
   
 c$$$      PFA='ETA'  
 c$$$      DEBUG=.true.  
 c$$$      DEBUG=.false.  
48    
49        if(DEBUG)PRINT*,'P.F.A. --> ',fin,PFA        call idtoc(pfaid,PFA)
50          
51    
52        if(DEBUG)then        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        call init_hough        SECOND_SEARCH    = .false.
 *------------------------------------------------------  
65                
66          ntrk_old = 0
67    
68          dedx_x_min = dedx_x_min_mip !mip search
69          dedx_y_min = dedx_y_min_mip !
70    
71    c 887  continue ! Emiliano, fix compilation warning (gcc 4.1), Label 887 defined but not used
72    
73  *------------------------------------------------------  *------------------------------------------------------
 *     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  
         
74        do i=1,nclstr1        do i=1,nclstr1
75           cl_used(i)=0           !init mask of clusters associated to a track           cl_used(i)=0           !init mask of clusters associated to a track
76        enddo        enddo
         
       if(DEBUG)then  
          print*,'----------------------------------'  
          print*,iev,'   ** ',nev2  
       endif  
77    
78          call init_level2
79          call init_hough      
80    *------------------------------------------------------
81    
82     888  continue
83    
84        
85  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
86  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
87    
# Line 95  c$$$      endif Line 91  c$$$      endif
91           goto 880               !fill ntp and go to next event               goto 880               !fill ntp and go to next event    
92        endif        endif
93    
94        call fill_hough  *     ----------------------------------------------------
95    *     fill hough output variables only for the first searc
96    *     ----------------------------------------------------
97          if(.not.SECOND_SEARCH)call fill_hough
98    
99        iflag=0        iflag=0
100        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  
101                
102        
103    *     ////////////////////////////////////////////////
104    *     PATCH to recover events with less than 3 couples
105    *     ////////////////////////////////////////////////    
106    *     if no tracks have been found
107    *     and n.planes with at least one couple (planehit)
108    *     is less than 3 try to recover applying the hough
109    *     transform to singlet
110          nplanehit=0
111          do ip=1,NPLANES
112             if(ncp_plane(ip).gt.0)nplanehit=nplanehit+1
113          enddo
114          if(DEBUG.eq.1)then
115    cc      if(.true.)then
116             print*,''
117             print*,'-----------------------------------'
118             print*,'n.tracks     ',ntrk
119             print*,'n.hit-planes ',nplanehit
120          endif
121          if(
122         $     .not.RECOVER_SINGLETS.and. !recover only once
123         $     ntrk.eq.0            .and. !no tracks has been found
124         $     nplanehit.le.3       .and. !if less than 3 hit planes (n.real couples)
125         $     .true.)then
126             RECOVER_SINGLETS=.true.
127    cc         if(.true.)
128             if(DEBUG.EQ.1)
129         $        print*,' -----> ((( TRY TO RECOVER SINGLETS )))'
130             goto 888               !back to track finding
131          endif
132    
133    
134    *     ////////////////////////////////////////////////
135    *     PATCH to recover nuclei
136    *     ////////////////////////////////////////////////
137          itrk = ntrk               !check the last track first
138          isimage = 0
139          BAD_NUCLEUS=.false.
140     150  continue
141    
142          dedxav  = 0.
143          npt     = 0
144          dedxmin = 1000000
145          if(itrk.gt.0)then
146             do ip=1,NPLANES
147    cc            print*,ip,' ** ',dedx_x(ip,itrk),dedx_y(ip,itrk)
148                if(xgood_nt(ip,itrk).eq.1)then
149                   dedxav = dedxav + abs(dedx_x(ip,itrk))
150                   if(abs(dedx_x(ip,itrk)).lt.dedxmin)
151         $              dedxmin=abs(dedx_x(ip,itrk))
152                   npt = npt + 1
153                endif        
154                if(ygood_nt(ip,itrk).eq.1)then
155                   dedxav = dedxav + abs(dedx_y(ip,itrk))
156                   if(abs(dedx_y(ip,itrk)).lt.dedxmin)
157         $              dedxmin=abs(dedx_y(ip,itrk))
158                   npt = npt + 1
159                endif            
160             enddo
161             if(npt.gt.0)dedxav = dedxav/npt
162             if(  dedxav.gt.dedx_min_nuclei )then
163                if( (dedxav-dedxmin)/dedxav.gt.ddedx_min_nuclei )
164         $           BAD_NUCLEUS=.true.
165             endif
166          endif
167          if(DEBUG.eq.1)then
168    cc      if(.true.)then
169             print*,''
170             print*,'-----------------------------------'
171             print*,'itrk          ',itrk
172             print*,'dedx(average) ',dedxav,' cut ',dedx_min_nuclei
173             print*,'dedx(minimum) ',dedxmin
174             if(dedxav.gt.0)print*,'(av-min)/av',(dedxav-dedxmin)/dedxav
175         $        ,' cut ',ddedx_min_nuclei
176             print*,'BAD_NUCLEUS ? ',BAD_NUCLEUS
177          endif
178    
179          itrk=0
180          if(ntrk.gt.0)itrk = image(ntrk)              
181          if(itrk.gt.0.and.isimage.eq.0)then      !hence check the image
182             isimage=1        
183             goto 150
184          endif
185    
186    
187          if(
188         $     .not.RECOVER_NUCLEI.and. !recover only once
189         $     (ntrk.eq.0.or.BAD_NUCLEUS).and. !no tracks found or bad nucleus
190         $     .true.)then
191    
192             RECOVER_NUCLEI   = .true.
193             dedx_x_min = dedx_x_min_nuclei
194             dedx_y_min = dedx_y_min_nuclei
195    cc         if(.true.)
196             if(DEBUG.EQ.1)
197         $     print*,' -----> ((( TRY TO RECOVER NUCLEI )))'
198    *     ----------------------------------------------
199    *     in case of re-tracking, un-tag used clusters
200    *     ----------------------------------------------
201             if(ntrk.ne.0)then
202                do ip=1,NPLANES
203                   if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
204                   if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
205                enddo
206                ntrk = ntrk-1       !decrement track
207                if(isimage.eq.1)then
208                   do ip=1,NPLANES
209                      if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
210                      if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
211                   enddo
212                   ntrk = ntrk-1    !decrement track image
213                endif
214             endif
215             goto 888               !back to track finding
216            
217          endif
218    *     ----------------------------------------------
219    *     back to default parameters for the second-track search
220    *     ----------------------------------------------
221          if(RECOVER_NUCLEI)then
222             dedx_x_min = dedx_x_min_mip
223             dedx_y_min = dedx_y_min_mip
224             RECOVER_NUCLEI   = .false.    
225          endif
226    
227          if(ntrk.eq.0)goto 880
228    
229    *     /////////////////////////////////////////////////////
230    *     PATCH to perform a more efficient second-track search
231    *     /////////////////////////////////////////////////////
232    *     -------------------------------------------------
233    *     n. hit planes and number of clusters
234    *     -------------------------------------------------
235          do ip=1,NPLANES
236             hitplanex(ip)=0
237             hitplaney(ip)=0
238             nclusters=0
239          enddo
240          do icl=1,nclstr1
241             ip = npl(VIEW(icl))
242             if(cl_good(icl).eq.0.or.cl_used(icl).eq.1)goto 9
243             if(mod(VIEW(icl),2).eq.1)hitplaney(ip)=hitplaney(ip)+1 !YY        
244             if(mod(VIEW(icl),2).eq.0)hitplanex(ip)=hitplanex(ip)+1 !XX        
245             nclusters=nclusters+1
246     9       continue
247          enddo
248          nhitplaney=0
249          nhitplanex=0
250          do ip=1,NPLANES
251             if(hitplaney(ip).gt.0)nhitplaney=nhitplaney+1
252             if(hitplanex(ip).gt.0)nhitplanex=nhitplanex+1
253          enddo      
254    *     -------------------------------------------------
255    *     chi2 cut
256    *     -------------------------------------------------
257          p0 = 1.111588e+00
258          p1 = 1.707656e+00
259          p2 = 1.489693e-01
260          chi2m025 = p0
261         $     + abs(al_nt(5,ntrk))*p1
262         $     + al_nt(5,ntrk)*al_nt(5,ntrk)*p2
263          chi2m025_0 = p0 + 0.07*p1 + 0.07*0.07*p2
264    cc      chi2m = (chi2m025-chi2m025_0+12**0.25)**4. !95%
265          chi2m = (chi2m025-chi2m025_0+20.**0.25)**4. !???
266    *     -------------------------------------------------
267    *     image track
268    *     -------------------------------------------------
269          ntrk_best = ntrk
270          chi2_best = chi2_nt(ntrk)
271          if(image(ntrk).ne.0)
272         $     chi2_best = min(chi2_nt(ntrk),chi2_nt(image(ntrk)))
273          if(DEBUG.eq.1)then
274    cc      if(.true.)then
275             print*,''
276             print*,'-----------------------------------'
277             print*,'chi2_nt(',ntrk_best,') = ',chi2_best,' < ',chi2m,' ?'
278             print*,'nhitplanex',nhitplanex
279             print*,'nhitplaney',nhitplaney
280             print*,'nclusters',nclusters
281          endif
282          if(
283         $     ntrk.gt.ntrk_old.and. !one track must have been already found
284         $     chi2_best.gt.0 .and.     !this track must not be garbage
285         $     chi2_best.lt.chi2m .and. !"
286         $     nhitplaney.ge.3.and. !a minimum number of hit planes is required
287         $     nhitplanex.ge.3.and. !"
288         $     nclusters.le.20.and. !limit on number of clusters to save time
289         $     .true.)then
290    
291             SECOND_SEARCH=.true.
292    c         RECOVER_SINGLETS=.true.      !OPTIMIZATION  
293             ntrk_old=ntrk
294    
295             if(DEBUG.EQ.1)
296    cc         if(.true.)
297         $        print*,'***** NEW SEARCH ****'
298             goto 888               !back to track finding
299    
300          endif
301    
302   880  continue   880  continue
303    
# Line 113  c$$$      endif Line 307  c$$$      endif
307    
308        call fill_level2_siglets        call fill_level2_siglets
309                
310        if(DEBUG)then        if(DEBUG.EQ.1)then
               
311           print*,''           print*,''
312           print*,'DONE!'           print*,'DONE!'
313           print*,''           print*,''
# Line 129  c$$$      endif Line 322  c$$$      endif
322        do iv = 1,nviews        do iv = 1,nviews
323           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
324        enddo        enddo
 c$$$      if(ngood.ne.0)print*,'* WARNING * Event '  
 c$$$     $     ,':LEVEL2 event status: '  
 c$$$     $     ,(good2(i),i=1,nviews)  
325                
326   8800 continue  c 8800 continue
327          continue
328                
329                
330  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
331  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
332    
333   100  continue  c 100  continue
334          continue
335                
336        return        return
337        end        end
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23