/[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.23 by pam-ts, Wed Jun 4 07:57:04 2014 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.
 *------------------------------------------------------  
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
         
 c$$$      if(DEBUG.EQ.1)then  
 c$$$         print*,'----------------------------------'  
 c$$$         print*,iev,'   ** ',nev2  
 c$$$      endif  
77    
78          call init_level2
79          call init_hough      
80    *------------------------------------------------------
81    
82     888  continue
83    
84        
85  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
86  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
87    
# Line 100  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    *
136    *     EV may 2014 ==> SKIP NUCLEI PATCH
137    *     (done differently...)
138    *
139    *     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@
140          goto 7575
141    
142    
143    *     ////////////////////////////////////////////////
144    *     PATCH to recover nuclei
145    *     ////////////////////////////////////////////////
146          itrk = ntrk               !check the last track first
147          isimage = 0
148          BAD_NUCLEUS=.false.
149     150  continue
150    
151          dedxav  = 0.
152          npt     = 0
153          dedxmin = 1000000
154          if(itrk.gt.0)then
155             do ip=1,NPLANES
156    cc            print*,ip,' ** ',dedx_x(ip,itrk),dedx_y(ip,itrk)
157                if(xgood_nt(ip,itrk).eq.1)then
158                   dedxav = dedxav + abs(dedx_x(ip,itrk))
159                   if(abs(dedx_x(ip,itrk)).lt.dedxmin)
160         $              dedxmin=abs(dedx_x(ip,itrk))
161                   npt = npt + 1
162                endif        
163                if(ygood_nt(ip,itrk).eq.1)then
164                   dedxav = dedxav + abs(dedx_y(ip,itrk))
165                   if(abs(dedx_y(ip,itrk)).lt.dedxmin)
166         $              dedxmin=abs(dedx_y(ip,itrk))
167                   npt = npt + 1
168                endif            
169             enddo
170             if(npt.gt.0)dedxav = dedxav/npt
171             if(  dedxav.gt.dedx_min_nuclei )then
172                if( (dedxav-dedxmin)/dedxav.gt.ddedx_min_nuclei )
173         $           BAD_NUCLEUS=.true.
174             endif
175          endif
176          if(DEBUG.eq.1)then
177    cc      if(.true.)then
178             print*,''
179             print*,'-----------------------------------'
180             print*,'itrk          ',itrk
181             print*,'dedx(average) ',dedxav,' cut ',dedx_min_nuclei
182             print*,'dedx(minimum) ',dedxmin
183             if(dedxav.gt.0)print*,'(av-min)/av',(dedxav-dedxmin)/dedxav
184         $        ,' cut ',ddedx_min_nuclei
185             print*,'BAD_NUCLEUS ? ',BAD_NUCLEUS
186          endif
187    
188          itrk=0
189          if(ntrk.gt.0)itrk = image(ntrk)              
190          if(itrk.gt.0.and.isimage.eq.0)then      !hence check the image
191             isimage=1        
192             goto 150
193          endif
194    
195    
196          if(
197         $     .not.RECOVER_NUCLEI.and. !recover only once
198         $     (ntrk.eq.0.or.BAD_NUCLEUS).and. !no tracks found or bad nucleus
199         $     .true.)then
200    
201             RECOVER_NUCLEI   = .true.
202             dedx_x_min = dedx_x_min_nuclei
203             dedx_y_min = dedx_y_min_nuclei
204    cc         if(.true.)
205             if(DEBUG.EQ.1)
206         $     print*,' -----> ((( TRY TO RECOVER NUCLEI )))'
207    *     ----------------------------------------------
208    *     in case of re-tracking, un-tag used clusters
209    *     ----------------------------------------------
210             if(ntrk.ne.0)then
211                do ip=1,NPLANES
212                   if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
213                   if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
214                enddo
215                ntrk = ntrk-1       !decrement track
216                if(isimage.eq.1)then
217                   do ip=1,NPLANES
218                      if(cltrx(ip,ntrk).ne.0)cl_used(cltrx(ip,ntrk)) = 0 !un-tag used clusters
219                      if(cltry(ip,ntrk).ne.0)cl_used(cltry(ip,ntrk)) = 0 !un-tag used clusters
220                   enddo
221                   ntrk = ntrk-1    !decrement track image
222                endif
223             endif
224             goto 888               !back to track finding
225            
226          endif
227    *     ----------------------------------------------
228    *     back to default parameters for the second-track search
229    *     ----------------------------------------------
230          if(RECOVER_NUCLEI)then
231             dedx_x_min = dedx_x_min_mip
232             dedx_y_min = dedx_y_min_mip
233             RECOVER_NUCLEI   = .false.    
234          endif
235    
236    *     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@
237     7575 continue
238    *     @~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@~@
239    
240          if(ntrk.eq.0)goto 880
241    
242    *     /////////////////////////////////////////////////////
243    *     PATCH to perform a more efficient second-track search
244    *     /////////////////////////////////////////////////////
245    *     -------------------------------------------------
246    *     n. hit planes and number of clusters
247    *     -------------------------------------------------
248          do ip=1,NPLANES
249             hitplanex(ip)=0
250             hitplaney(ip)=0
251             nclusters=0
252          enddo
253          do icl=1,nclstr1
254             ip = npl(VIEW(icl))
255             if(cl_good(icl).eq.0.or.cl_used(icl).eq.1)goto 9
256             if(mod(VIEW(icl),2).eq.1)hitplaney(ip)=hitplaney(ip)+1 !YY        
257             if(mod(VIEW(icl),2).eq.0)hitplanex(ip)=hitplanex(ip)+1 !XX        
258             nclusters=nclusters+1
259     9       continue
260          enddo
261          nhitplaney=0
262          nhitplanex=0
263          do ip=1,NPLANES
264             if(hitplaney(ip).gt.0)nhitplaney=nhitplaney+1
265             if(hitplanex(ip).gt.0)nhitplanex=nhitplanex+1
266          enddo      
267    *     -------------------------------------------------
268    *     chi2 cut
269    *     -------------------------------------------------
270          p0 = 1.111588e+00
271          p1 = 1.707656e+00
272          p2 = 1.489693e-01
273          chi2m025 = p0
274         $     + abs(al_nt(5,ntrk))*p1
275         $     + al_nt(5,ntrk)*al_nt(5,ntrk)*p2
276          chi2m025_0 = p0 + 0.07*p1 + 0.07*0.07*p2
277    cc      chi2m = (chi2m025-chi2m025_0+12**0.25)**4. !95%
278          chi2m = (chi2m025-chi2m025_0+20.**0.25)**4. !???
279    *     -------------------------------------------------
280    *     image track
281    *     -------------------------------------------------
282          ntrk_best = ntrk
283          chi2_best = chi2_nt(ntrk)
284          if(image(ntrk).ne.0)
285         $     chi2_best = min(chi2_nt(ntrk),chi2_nt(image(ntrk)))
286          if(DEBUG.eq.1)then
287    cc      if(.true.)then
288             print*,''
289             print*,'-----------------------------------'
290             print*,'chi2_nt(',ntrk_best,') = ',chi2_best,' < ',chi2m,' ?'
291             print*,'nhitplanex',nhitplanex
292             print*,'nhitplaney',nhitplaney
293             print*,'nclusters',nclusters
294          endif
295          if(
296         $     ntrk.gt.ntrk_old.and. !one track must have been already found
297         $     chi2_best.gt.0 .and.     !this track must not be garbage
298         $     chi2_best.lt.chi2m .and. !"
299         $     nhitplaney.ge.3.and. !a minimum number of hit planes is required
300         $     nhitplanex.ge.3.and. !"
301         $     nclusters.le.20.and. !limit on number of clusters to save time
302         $     .true.)then
303    
304             SECOND_SEARCH=.true.
305    c         RECOVER_SINGLETS=.true.      !OPTIMIZATION  
306             ntrk_old=ntrk
307    
308             if(DEBUG.EQ.1)
309    cc         if(.true.)
310         $        print*,'***** NEW SEARCH ****'
311             goto 888               !back to track finding
312    
313          endif
314    
315   880  continue   880  continue
316    
317  *     **********************************************************  *     **********************************************************
# Line 119  c$$$      endif Line 321  c$$$      endif
321        call fill_level2_siglets        call fill_level2_siglets
322                
323        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
               
324           print*,''           print*,''
325           print*,'DONE!'           print*,'DONE!'
326           print*,''           print*,''
# Line 134  c$$$      endif Line 335  c$$$      endif
335        do iv = 1,nviews        do iv = 1,nviews
336           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
337        enddo        enddo
 c$$$      if(ngood.ne.0)print*,'* WARNING * Event '  
 c$$$     $     ,':LEVEL2 event status: '  
 c$$$     $     ,(good2(i),i=1,nviews)  
338                
339   8800 continue  c 8800 continue
340          continue
341                
342                
343  *     ///////////////////////////////////////////////  *     ///////////////////////////////////////////////
344  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\  *     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
345    
346   100  continue  c 100  continue
347          continue
348                
349        return        return
350        end        end
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23