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

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

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

revision 1.4 by pam-fi, Thu Sep 28 14:04:40 2006 UTC revision 1.15 by pam-fi, Tue Nov 21 17:13:31 2006 UTC
# Line 12  Line 12 
12        subroutine track_finding(iflag)        subroutine track_finding(iflag)
13    
14        include 'commontracker.f'        include 'commontracker.f'
15          include 'level1.f'
16        include 'common_momanhough.f'        include 'common_momanhough.f'
17        include 'common_mech.f'        include 'common_mech.f'
18        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
19        include 'common_mini_2.f'        include 'common_mini_2.f'
20        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
23        include 'momanhough_init.f'  c      include 'momanhough_init.f'
24                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
25  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
26  *     STEP 1  *     STEP 1
27  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 47  c      common/dbg/DEBUG Line 44  c      common/dbg/DEBUG
44  c      iflag=0  c      iflag=0
45        call cl_to_couples(iflag)        call cl_to_couples(iflag)
46        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
47           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
48        endif        endif
49                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     selezione di tracce pulite per diagnostica  
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 c$$$         if(DEBUG)then  
 c$$$            do ip=1,nplanes  
 c$$$               if(ncp_plane(ip).ne.1)good2=.false.  
 c$$$            enddo  
 c$$$c            if(good2.eq.0)goto 100!next event  
 c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event  
 c$$$         endif  
           
   
   
50  *-----------------------------------------------------  *-----------------------------------------------------
51  *-----------------------------------------------------  *-----------------------------------------------------
52  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 94  c$$$         endif Line 78  c$$$         endif
78  c      iflag=0  c      iflag=0
79        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
80        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
81           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
82        endif        endif
83                
84                
# Line 123  c      iflag=0 Line 107  c      iflag=0
107  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    *     count number of hit planes
111          planehit=0                
112          do np=1,nplanes          
113            if(ncp_plane(np).ne.0)then  
114              planehit=planehit+1  
115            endif                  
116          enddo                    
117          if(planehit.lt.3) goto 880 ! exit              
118          
119          nptxz_min=x_min_start              
120          nplxz_min=x_min_start              
121                
122          nptyz_min=y_min_start              
123          nplyz_min=y_min_start              
124    
125  c      iflag=0        cutdistyz=cutystart      
126          cutdistxz=cutxstart      
127    
128     878  continue
129        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
130        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
131           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
132        endif        endif
133  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
134            if(cutdistyz.lt.maxcuty/2)then
135              cutdistyz=cutdistyz+cutystep
136            else
137              cutdistyz=cutdistyz+(3*cutystep)
138            endif
139            goto 878                
140          endif                    
141    
142          if(planehit.eq.3) goto 881          
143          
144     879  continue  
145        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
146        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
147           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
148        endif        endif
149                                    
150          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
151            cutdistxz=cutdistxz+cutxstep
152            goto 879                
153          endif                    
154    
155        
156     881  continue  
157    *     if there is at least three planes on the Y view decreases cuts on X view
158          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
159         $     nplxz_min.ne.y_min_start)then
160            nptxz_min=x_min_step    
161            nplxz_min=x_min_start-x_min_step              
162            goto 879                
163          endif                    
164            
165   880  return   880  return
166        end        end
167    
# Line 144  c      iflag=0 Line 171  c      iflag=0
171        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
172    
173        include 'commontracker.f'        include 'commontracker.f'
174          include 'level1.f'
175        include 'common_momanhough.f'        include 'common_momanhough.f'
176        include 'common_mech.f'        include 'common_mech.f'
177        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
178        include 'common_mini_2.f'        include 'common_mini_2.f'
179        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
180        include 'level2.f'        include 'level2.f'
181    
182        include 'momanhough_init.f'  c      include 'momanhough_init.f'
183                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
184        logical FIMAGE            !        logical FIMAGE            !
185          real*8 AL_GUESS(5)
186    
187  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
188  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 198  c         iflag=0 Line 223  c         iflag=0
223           ibest=0                !best track among candidates           ibest=0                !best track among candidates
224           iimage=0               !track image           iimage=0               !track image
225  *     ------------- select the best track -------------  *     ------------- select the best track -------------
226           rchi2best=1000000000.  c$$$         rchi2best=1000000000.
227    c$$$         do i=1,ntracks
228    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
229    c$$$     $         RCHI2_STORE(i).gt.0)then
230    c$$$               ibest=i
231    c$$$               rchi2best=RCHI2_STORE(i)
232    c$$$            endif
233    c$$$         enddo
234    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
235    
236    *     -------------------------------------------------------
237    *     order track-candidates according to:
238    *     1st) decreasing n.points
239    *     2nd) increasing chi**2
240    *     -------------------------------------------------------
241             rchi2best=1000000000.
242             ndofbest=0             !(1)
243           do i=1,ntracks           do i=1,ntracks
244              if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0               !(1)
245       $         RCHI2_STORE(i).gt.0)then             do ii=1,nplanes      !(1)
246                 ndof=ndof          !(1)
247         $            +int(xgood_store(ii,i)) !(1)
248         $            +int(ygood_store(ii,i)) !(1)
249               enddo                !(1)
250               if(ndof.gt.ndofbest)then !(1)
251                 ibest=i
252                 rchi2best=RCHI2_STORE(i)
253                 ndofbest=ndof      !(1)
254               elseif(ndof.eq.ndofbest)then !(1)
255                 if(RCHI2_STORE(i).lt.rchi2best.and.
256         $            RCHI2_STORE(i).gt.0)then
257                 ibest=i                 ibest=i
258                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
259              endif                 ndofbest=ndof    !(1)
260           enddo               endif              !(1)
261               endif
262             enddo
263    
264    c$$$         rchi2best=1000000000.
265    c$$$         ndofbest=0             !(1)
266    c$$$         do i=1,ntracks
267    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
268    c$$$     $          RCHI2_STORE(i).gt.0)then
269    c$$$             ndof=0             !(1)
270    c$$$             do ii=1,nplanes    !(1)
271    c$$$               ndof=ndof        !(1)
272    c$$$     $              +int(xgood_store(ii,i)) !(1)
273    c$$$     $              +int(ygood_store(ii,i)) !(1)
274    c$$$             enddo              !(1)
275    c$$$             if(ndof.ge.ndofbest)then !(1)
276    c$$$               ibest=i
277    c$$$               rchi2best=RCHI2_STORE(i)
278    c$$$               ndofbest=ndof    !(1)
279    c$$$             endif              !(1)
280    c$$$           endif
281    c$$$         enddo
282    
283           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
284  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
285  *     The best track candidate (ibest) is selected and a new fitting is performed.  *     The best track candidate (ibest) is selected and a new fitting is performed.
# Line 236  c         iflag=0 Line 310  c         iflag=0
310  *     **********************************************************  *     **********************************************************
311  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
312  *     **********************************************************  *     **********************************************************
313             call guess()
314             do i=1,5
315                AL_GUESS(i)=AL(i)
316             enddo
317    
318           do i=1,5           do i=1,5
319              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
320           enddo           enddo
321            
322           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
323           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
324           jstep=0                !# minimization steps           jstep=0                !# minimization steps
325    
326           call mini_2(jstep,ifail)           iprint=0
327    c         if(DEBUG)iprint=1
328             if(VERBOSE)iprint=1
329             call mini2(jstep,ifail,iprint)
330           if(ifail.ne.0) then           if(ifail.ne.0) then
331              if(DEBUG)then              if(.true.)then
332                 print *,                 print *,
333       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
334       $              ,iev       $              ,iev
335    
336    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
337    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
338    c$$$               print*,'result:   ',(al(i),i=1,5)
339    c$$$               print*,'xgood ',xgood
340    c$$$               print*,'ygood ',ygood
341    c$$$               print*,'----------------------------------------------'
342              endif              endif
343              chi2=-chi2  c            chi2=-chi2
344           endif           endif
345                    
346           if(DEBUG)then           if(DEBUG)then
# Line 311  c         print*,'++++++++++ iimage,fima Line 401  c         print*,'++++++++++ iimage,fima
401  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
402    
403           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
404              if(DEBUG)              if(verbose)
405       $           print*,       $           print*,
406       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
407       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 361  c     $        rchi2best.lt.15..and. Line 451  c     $        rchi2best.lt.15..and.
451        end        end
452    
453    
   
   
 c$$$************************************************************  
 c$$$  
 c$$$      subroutine readmipparam  
 c$$$              
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('trk-LADDER',i1,'-mip.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READMIPPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do iv=1,nviews  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            mip(int(pip),ilad)  
 c$$$c            print*,ilad,iv,pip,mip(int(pip),ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readchargeparam  
 c$$$        
 c$$$        
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('charge-l',i1,'.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READCHARGEPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do ip=1,nplanes  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)          
 c$$$c            print*,ilad,ip,pip,kch(ip,ilad),  
 c$$$c     $           cch(ip,ilad),sch(ip,ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readetaparam  
 c$$$*     -----------------------------------------  
 c$$$*     read eta2,3,4 calibration parameters  
 c$$$*     and fill variables:  
 c$$$*  
 c$$$*     eta2(netabin,nladders_view,nviews)  
 c$$$*     eta3(2*netabin,nladders_view,nviews)  
 c$$$*     eta4(2*netabin,nladders_view,nviews)  
 c$$$*  
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*40 fname_binning  
 c$$$      character*40 fname_param  
 c$$$c      character*120 cmd1  
 c$$$c      character*120 cmd2  
 c$$$  
 c$$$  
 c$$$******retrieve ANGULAR BINNING info  
 c$$$      fname_binning='binning.dat'  
 c$$$      print *,'Opening file: ',fname_binning  
 c$$$      open(10,  
 c$$$     $     FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning))  
 c$$$     $     ,STATUS='UNKNOWN'  
 c$$$     $     ,IOSTAT=iostat  
 c$$$     $     )  
 c$$$      if(iostat.ne.0)then  
 c$$$         print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$         return  
 c$$$      endif  
 c$$$      print*,'---- ANGULAR BINNING ----'  
 c$$$      print*,'Bin   -   angL   -   angR'  
 c$$$ 101  format(i2,'       ',f6.2,'     ',f6.2)  
 c$$$      do ibin=1,nangmax  
 c$$$         read(10,*  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )xnn,angL(ibin),angR(ibin)  
 c$$$         if(iostat.ne.0)goto 1000  
 c$$$         write(*,101)int(xnn),angL(ibin),angR(ibin)  
 c$$$      enddo          
 c$$$ 1000 nangbin=int(xnn)  
 c$$$      close(10)  
 c$$$      print*,'-------------------------'  
 c$$$        
 c$$$  
 c$$$  
 c$$$      do ieta=2,4               !loop on eta 2,3,4          
 c$$$******retrieve correction parameters  
 c$$$ 200     format(' Opening eta',i1,' files...')  
 c$$$         write(*,200)ieta  
 c$$$  
 c$$$ 201     format('eta',i1,'-bin',i1,'-l',i1,'.dat')  
 c$$$ 202     format('eta',i1,'-bin',i2,'-l',i1,'.dat')  
 c$$$         do iang=1,nangbin  
 c$$$            do ilad=1,nladders_view  
 c$$$               if(iang.lt.10)write(fname_param,201)ieta,iang,ilad  
 c$$$               if(iang.ge.10)write(fname_param,202)ieta,iang,ilad  
 c$$$c               print *,'Opening file: ',fname_param  
 c$$$               open(10,  
 c$$$     $             FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $              ,STATUS='UNKNOWN'  
 c$$$     $              ,IOSTAT=iostat  
 c$$$     $              )  
 c$$$               if(iostat.ne.0)then  
 c$$$                  print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$                  return  
 c$$$               endif  
 c$$$               do ival=1,netavalmax  
 c$$$                  if(ieta.eq.2)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta2(ival,iang),  
 c$$$     $                 (feta2(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.3)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta3(ival,iang),  
 c$$$     $                 (feta3(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.4)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta4(ival,iang),  
 c$$$     $                 (feta4(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(iostat.ne.0)then  
 c$$$                     netaval=ival-1  
 c$$$c$$$                     if(eta2(1,iang).ne.-eta2(netaval,iang))  
 c$$$c$$$     $                    print*,'**** ERROR on parameters !!! ****'  
 c$$$                     goto 2000  
 c$$$                  endif  
 c$$$               enddo  
 c$$$ 2000          close(10)  
 c$$$*               print*,'... done'  
 c$$$            enddo  
 c$$$         enddo  
 c$$$  
 c$$$      enddo                     !end loop on eta 2,3,4  
 c$$$  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
454                
455  ************************************************************  ************************************************************
456  ************************************************************  ************************************************************
# Line 607  c                (implemented variable r Line 525  c                (implemented variable r
525  c*****************************************************  c*****************************************************
526                
527        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
528        include 'level1.f'        include 'level1.f'
529          include 'calib.f'
530    c      include 'level1.f'
531        include 'common_align.f'        include 'common_align.f'
532        include 'common_mech.f'        include 'common_mech.f'
533        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
534        include 'common_resxy.f'  c      include 'common_resxy.f'
535    
536  c      logical DEBUG  c      logical DEBUG
537  c      common/dbg/DEBUG  c      common/dbg/DEBUG
# Line 667  c      double precision xi_B,yi_B,zi_B Line 586  c      double precision xi_B,yi_B,zi_B
586              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
587           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
588  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
589  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
590  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
591              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
592              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
593              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
594       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
595              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
596           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
597              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
598              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
599              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
600       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
601              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
602           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
603              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
604              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
605              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
606       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
607              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
608           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
609              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
610              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
611              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
612       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 731  c     $     print*,PFAx,icx,angx,stripx, Line 650  c     $     print*,PFAx,icx,angx,stripx,
650              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
651           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
652  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
653  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
654  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
655              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
656              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
657              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
658              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
659       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
660           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
661              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
662              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
663              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
664       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
665           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
666              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
667              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
668              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
669       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
670           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
671              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
672              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
673  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
674              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 1284  c     $              ,iv,xvv(iv),yvv(iv) Line 1203  c     $              ,iv,xvv(iv),yvv(iv)
1203  *     it returns the plane number  *     it returns the plane number
1204  *      *    
1205        include 'commontracker.f'        include 'commontracker.f'
1206          include 'level1.f'
1207  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1208        include 'common_momanhough.f'        include 'common_momanhough.f'
1209                
# Line 1321  c      include 'common_analysis.f' Line 1241  c      include 'common_analysis.f'
1241  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1242  *      *    
1243        include 'commontracker.f'        include 'commontracker.f'
1244          include 'level1.f'
1245  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1246        include 'common_momanhough.f'        include 'common_momanhough.f'
1247                
# Line 1349  c      include 'common_analysis.f' Line 1270  c      include 'common_analysis.f'
1270  *     positive if sensor =2  *     positive if sensor =2
1271  *  *
1272        include 'commontracker.f'        include 'commontracker.f'
1273          include 'level1.f'
1274  c      include 'calib.f'  c      include 'calib.f'
1275  c      include 'level1.f'  c      include 'level1.f'
1276  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1378  c      include 'common_analysis.f' Line 1300  c      include 'common_analysis.f'
1300  *************************************************************************  *************************************************************************
1301  *************************************************************************  *************************************************************************
1302  *************************************************************************  *************************************************************************
 c$$$      subroutine book_debug  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$  
 c$$$      character*35 block1,block2,block3!,block4  
 c$$$     $     ,block5!,block6        
 c$$$        
 c$$$* * * * * * * * * * * * * * * * * * * * * * * *  
 c$$$*     HOUGH TRANSFORM PARAMETERS  
 c$$$        
 c$$$      call HBOOK2(1003  
 c$$$     $     ,'y vs tg thyz'  
 c$$$     $     ,300,-1.,1.         !x          
 c$$$     $     ,3000,-70.,70.,0.)   !y  
 c$$$  
 c$$$      call HBOOK1(1004  
 c$$$     $     ,'Dy'  
 c$$$     $     ,100,0.,2.,0.)    
 c$$$  
 c$$$      call HBOOK1(1005  
 c$$$     $     ,'D thyz'  
 c$$$     $     ,100,0.,.05,0.)    
 c$$$  
 c$$$  
 c$$$  
 c$$$*     DEBUG ntuple:  
 c$$$      call HBNT(ntp_level2+1,'LEVEL2',' ')  
 c$$$        
 c$$$      call HBNAME(ntp_level2+1,'EVENT',good2_nt,  
 c$$$     $     'GOOD2:L,NEV2:I')  
 c$$$  
 c$$$ 411  format('NDBLT:I::[0,',I5,']')  
 c$$$      write(block1,411) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,  
 c$$$     $     block1//'  
 c$$$     $     ,ALFAYZ1(NDBLT):R  
 c$$$     $     ,ALFAYZ2(NDBLT):R  
 c$$$     $     ,DB_CLOUD(NDBLT):I  
 c$$$     $     ')    
 c$$$  
 c$$$ 412  format('NTRPT:I::[0,',I5,']')  
 c$$$      write(block2,412) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,  
 c$$$     $     block2//'  
 c$$$     $     ,ALFAXZ1(NTRPT):R  
 c$$$     $     ,ALFAXZ2(NTRPT):R    
 c$$$     $     ,ALFAXZ3(NTRPT):R  
 c$$$     $     ,TR_CLOUD(NTRPT):I  
 c$$$     $     ')  
 c$$$        
 c$$$    
 c$$$ 413  format('NCLOUDS_YZ:I::[0,',I4,']')  
 c$$$c$$$ 414  format('DB_CLOUD(',I4,'):I')  
 c$$$      write(block3,413) ncloyz_max  
 c$$$c$$$      write(block4,414) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,  
 c$$$     $     block3//'  
 c$$$     $     ,ALFAYZ1_AV(NCLOUDS_YZ):R  
 c$$$     $     ,ALFAYZ2_AV(NCLOUDS_YZ):R  
 c$$$     $     ,PTCLOUD_YZ(NCLOUDS_YZ):I'  
 c$$$c$$$     $     ,'//block4  
 c$$$     $     )  
 c$$$  
 c$$$ 415  format('NCLOUDS_XZ:I::[0,',I4,']')  
 c$$$c$$$ 416  format('TR_CLOUD(',I5,'):I')  
 c$$$      write(block5,415) ncloxz_max  
 c$$$c$$$      write(block6,416) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,  
 c$$$     $     block5//'  
 c$$$     $     ,ALFAXZ1_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ2_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ3_AV(NCLOUDS_XZ):R  
 c$$$     $     ,PTCLOUD_XZ(NCLOUDS_XZ):I'  
 c$$$c$$$     $     ,'//block6  
 c$$$     $     )  
 c$$$  
 c$$$        
 c$$$      return  
 c$$$      end  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 c$$$      subroutine book_level2  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2)  
 c$$$c*****************************************************  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      character*35 block1,block2  
 c$$$  
 c$$$c      print*,'__________ booking LEVEL2 n-tuple __________'  
 c$$$  
 c$$$*     LEVEL1 ntuple:  
 c$$$      call HBNT(ntp_level2,'LEVEL2',' ')  
 c$$$        
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I')  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$c     $     ,WHIC_CALIB:I,SWCODE:I')  
 c$$$      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$     $     ,WHICH_CALIB:I,SWCODE:I,CRC(12):L')  
 c$$$c*********************************************************  
 c$$$    
 c$$$  
 c$$$c$$$# ifndef TEST2003  
 c$$$c$$$  
 c$$$c$$$      call HBNAME(ntp_level2,'CPU',pkt_type  
 c$$$c$$$     $     ,'PKT_TYPE:I::[0,50]  
 c$$$c$$$     $     ,PKT_NUM:I  
 c$$$c$$$     $     ,OBT:I'//  
 c$$$c$$$c********************************************************  
 c$$$c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c$$$c     $     ,WHICH_CALIB:I::[0,50]')  
 c$$$c$$$     $     ',CPU_CRC:L')  
 c$$$c$$$c********************************************************  
 c$$$c$$$  
 c$$$c$$$# endif  
 c$$$  
 c$$$ 417  format('NTRK:I::[0,',I4,']')  
 c$$$ 418  format(',IMAGE(NTRK):I::[0,',I4,']')  
 c$$$      write(block1,417)NTRKMAX  
 c$$$      write(block2,418)NTRKMAX  
 c$$$      call HBNAME(ntp_level2,'TRACKS',NTRK,  
 c$$$     $     block1//  
 c$$$     $     block2//'  
 c$$$     $     ,XM(6,NTRK):R  
 c$$$     $     ,YM(6,NTRK):R  
 c$$$     $     ,ZM(6,NTRK):R  
 c$$$     $     ,RESX(6,NTRK):R  
 c$$$     $     ,RESY(6,NTRK):R  
 c$$$     $     ,AL(5,NTRK):R  
 c$$$     $     ,COVAL(5,5,NTRK):R  
 c$$$     $     ,CHI2(NTRK):R  
 c$$$     $     ,XGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,YGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,XV(6,NTRK):R  
 c$$$     $     ,YV(6,NTRK):R  
 c$$$     $     ,ZV(6,NTRK):R  
 c$$$     $     ,AXV(6,NTRK):R  
 c$$$     $     ,AYV(6,NTRK):R'//  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     ,DEDXP(6,NTRK):R'//  
 c$$$c     $     ')  
 c$$$     $     ',DEDX_X(6,NTRK):R  
 c$$$     $     ,DEDX_Y(6,NTRK):R'//  
 c$$$c****************************************************  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c     $     ,CRC(12):L  
 c$$$     $     ',BdL(NTRK):R'  
 c$$$     $     )  
 c$$$c****************************************************  
 c$$$  
 c$$$    
 c$$$      call HBNAME(ntp_level2,'SINGLETX',nclsx,  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     'NCLSX(6):I,NCLSY(6):I')  
 c$$$     $     'NCLSX:I::[0,500],PLANEX(NCLSX):I  
 c$$$     $     ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2)  
 c$$$c    $     ,XS(NCLSX):R,SGNLXS(NCLSX):R')   !(2)  
 c$$$      call HBNAME(ntp_level2,'SINGLETY',nclsy,  
 c$$$     $     'NCLSY:I::[0,500],PLANEY(NCLSY):I  
 c$$$     $     ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2)  
 c$$$c    $     ,YS(NCLSY):R,SGNLYS(NCLSY):R')   !(2)  
 c$$$      return  
 c$$$      end  
   
 c$$$      subroutine fill_level2_clouds  
 c$$$c*****************************************************  
 c$$$c     29/11/2005 created by elena vannuccini  
 c$$$c*****************************************************  
 c$$$  
 c$$$*     -------------------------------------------------------  
 c$$$*     This routine fills the  variables related to the hough  
 c$$$*     transform, for the debig n-tuple  
 c$$$*     -------------------------------------------------------  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      good2_nt=.true.!good2  
 c$$$c      nev2_nt=nev2  
 c$$$        
 c$$$      if(.false.  
 c$$$     $     .or.ntrpt.gt.ntrpt_max_nt  
 c$$$     $     .or.ndblt.gt.ndblt_max_nt  
 c$$$     $     .or.NCLOUDS_XZ.gt.ncloxz_max  
 c$$$     $     .or.NCLOUDS_yZ.gt.ncloyz_max  
 c$$$     $     )then  
 c$$$         good2_nt=.false.  
 c$$$         ntrpt_nt=0  
 c$$$         ndblt_nt=0  
 c$$$         NCLOUDS_XZ_nt=0  
 c$$$         NCLOUDS_YZ_nt=0  
 c$$$      else  
 c$$$         ndblt_nt=ndblt  
 c$$$         ntrpt_nt=ntrpt  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do id=1,ndblt  
 c$$$               alfayz1_nt(id)=alfayz1(id) !Y0  
 c$$$               alfayz2_nt(id)=alfayz2(id) !tg theta-yz  
 c$$$c               db_cloud_nt(id)=db_cloud(id)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do it=1,ntrpt  
 c$$$               alfaxz1_nt(it)=alfaxz1(it) !X0  
 c$$$               alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz  
 c$$$               alfaxz3_nt(it)=alfaxz3(it) !1/r  
 c$$$c               tr_cloud_nt(it)=tr_cloud(it)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         nclouds_yz_nt=nclouds_yz  
 c$$$         nclouds_xz_nt=nclouds_xz  
 c$$$         if(nclouds_yz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do iyz=1,nclouds_yz  
 c$$$               ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)  
 c$$$               alfayz1_av_nt(iyz)=alfayz1_av(iyz)  
 c$$$               alfayz2_av_nt(iyz)=alfayz2_av(iyz)  
 c$$$               nnn=nnn+ptcloud_yz(iyz)  
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               db_cloud_nt(ipt)=db_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_yz '  
 c$$$c     $           ,(ptcloud_yz(i),i=1,nclouds_yz)  
 c$$$c            print*,'#### ntupla #### db_cloud '  
 c$$$c     $           ,(db_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$         if(nclouds_xz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do ixz=1,nclouds_xz  
 c$$$               ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)  
 c$$$               alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)  
 c$$$               alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)  
 c$$$               alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)  
 c$$$               nnn=nnn+ptcloud_xz(ixz)                
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               tr_cloud_nt(ipt)=tr_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_xz '  
 c$$$c     $           ,(ptcloud_xz(i),i=1,nclouds_xz)  
 c$$$c            print*,'#### ntupla #### tr_cloud '  
 c$$$c     $           ,(tr_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$      endif  
 c$$$      end  
1303                
1304    
1305  ***************************************************  ***************************************************
# Line 1660  c$$$      end Line 1314  c$$$      end
1314        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1315    
1316        include 'commontracker.f'        include 'commontracker.f'
1317          include 'level1.f'
1318        include 'common_momanhough.f'        include 'common_momanhough.f'
1319        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1320        include 'calib.f'        include 'calib.f'
1321        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1322    
1323  *     output flag  *     output flag
1324  *     --------------  *     --------------
# Line 1675  c      common/dbg/DEBUG Line 1327  c      common/dbg/DEBUG
1327  *     --------------  *     --------------
1328        integer iflag        integer iflag
1329    
1330        integer badseed,badcl        integer badseed,badclx,badcly
1331    
1332  *     init variables  *     init variables
1333        ncp_tot=0        ncp_tot=0
# Line 1691  c      common/dbg/DEBUG Line 1343  c      common/dbg/DEBUG
1343           ncls(ip)=0           ncls(ip)=0
1344        enddo        enddo
1345        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1346           cl_single(icl)=1           cl_single(icl) = 1
1347           cl_good(icl)=0           cl_good(icl)   = 0
1348          enddo
1349          do iv=1,nviews
1350             ncl_view(iv)  = 0
1351             mask_view(iv) = 0      !all included
1352        enddo        enddo
1353                
1354    *     count number of cluster per view
1355          do icl=1,nclstr1
1356             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1357          enddo
1358    *     mask views with too many clusters
1359          do iv=1,nviews
1360             if( ncl_view(iv).gt. nclusterlimit)then
1361                mask_view(iv) = 1
1362                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1363         $           ,nclusterlimit,' on view ', iv,' --> masked!'
1364             endif
1365          enddo
1366    
1367    
1368  *     start association  *     start association
1369        ncouples=0        ncouples=0
1370        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1371           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1372                    
1373  *     ----------------------------------------------------  *     ----------------------------------------------------
1374    *     jump masked views (X VIEW)
1375    *     ----------------------------------------------------
1376             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1377    *     ----------------------------------------------------
1378  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1379  *     ----------------------------------------------------  *     ----------------------------------------------------
1380           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
# Line 1717  c      common/dbg/DEBUG Line 1391  c      common/dbg/DEBUG
1391           else           else
1392              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1393           endif           endif
1394           badcl=badseed           badclx=badseed
1395           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1396              ibad=1              ibad=1
1397              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1727  c      common/dbg/DEBUG Line 1401  c      common/dbg/DEBUG
1401       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1402       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1403              endif              endif
1404              badcl=badcl*ibad              badclx=badclx*ibad
1405           enddo           enddo
1406  *     ----------------------------------------------------  *     ----------------------------------------------------
1407  *     >>> eliminato il taglio sulle BAD <<<  *     >>> eliminato il taglio sulle BAD <<<
# Line 1746  c     endif Line 1420  c     endif
1420              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1421                            
1422  *     ----------------------------------------------------  *     ----------------------------------------------------
1423    *     jump masked views (Y VIEW)
1424    *     ----------------------------------------------------
1425                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1426    
1427    *     ----------------------------------------------------
1428  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1429  *     ----------------------------------------------------  *     ----------------------------------------------------
1430              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
# Line 1762  c     endif Line 1441  c     endif
1441              else              else
1442                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1443              endif              endif
1444              badcl=badseed              badcly=badseed
1445              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1446                 ibad=1                 ibad=1
1447                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1771  c     endif Line 1450  c     endif
1450       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1451       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1452       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1453                 badcl=badcl*ibad                 badcly=badcly*ibad
1454              enddo              enddo
1455  *     ----------------------------------------------------  *     ----------------------------------------------------
1456  *     >>> eliminato il taglio sulle BAD <<<  *     >>> eliminato il taglio sulle BAD <<<
# Line 1794  c     endif Line 1473  c     endif
1473  *     charge correlation  *     charge correlation
1474  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1475    
1476  *     -------------------------------------------------------------                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1477  *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<       $              .and.
1478  *     -------------------------------------------------------------       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1479  c$$$               if(dedx(icy).lt.chsaty.or.dedx(icx).lt.chsatx)then       $              .and.
1480  c$$$                  ddd=(dedx(icy)       $              (badclx.eq.1.and.badcly.eq.1)
1481  c$$$     $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $              .and.
1482  c$$$                  ddd=ddd/sqrt(kch(nplx,nldx)**2+1)       $              .true.)then
1483  c$$$                  cut=chcut*sch(nplx,nldx)  
1484  c$$$                  if(abs(ddd).gt.cut)goto 20 !charge not consistent                    ddd=(dedx(icy)
1485  c$$$               endif       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1486                                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1487  *     ------------------> COUPLE <------------------  
1488  *     check to do not overflow vector dimentions  c                  cut = chcut * sch(nplx,nldx)
1489                 if(ncp_plane(nplx).gt.ncouplemax)then  
1490                    if(DEBUG)print*,                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1491       $                    ' ** warning ** number of identified'//       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1492       $                    ' couples on plane ',nplx,                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1493       $                    ' exceeds vector dimention'//                    cut = chcut * (16 + sss/50.)
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
                 
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
1494    
1495   20         continue                    if(abs(ddd).gt.cut)then
1496           enddo                  !end loop on clusters(Y)                       goto 20    !charge not consistent
1497                              endif
1498   10      continue                 endif
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
1499    
 *     init variables  
       ncp_tot=0  
       do ip=1,nplanes  
          do ico=1,ncouplemax  
             clx(ip,ico)=0  
             cly(ip,ico)=0  
          enddo  
          ncp_plane(ip)=0  
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
 *     start association  
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
 *     ----------------------------------------------------  
 *     cut on charge (X VIEW)  
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     cut BAD (X VIEW)              
          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
          ifirst=INDSTART(icx)  
          if(icx.ne.nclstr1) then  
             ilast=INDSTART(icx+1)-1  
          else  
             ilast=TOTCLLENGTH  
          endif  
          badcl=badseed  
          do igood=-ngoodstr,ngoodstr  
             ibad=1  
             if((INDMAX(icx)+igood).gt.ifirst.and.  
      $           (INDMAX(icx)+igood).lt.ilast.and.  
      $           .true.)then  
                ibad=BAD(VIEW(icx),  
      $              nvk(MAXS(icx)+igood),  
      $              nst(MAXS(icx)+igood))  
             endif  
             badcl=badcl*ibad  
          enddo  
          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
             goto 10             !<<<<<<<<<<<<<< BAD cut  
          endif                  !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
           
          cl_good(icx)=1  
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
 *     ----------------------------------------------------  
 *     cut on charge (Y VIEW)  
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     cut BAD (Y VIEW)              
             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
             ifirst=INDSTART(icy)  
             if(icy.ne.nclstr1) then  
                ilast=INDSTART(icy+1)-1  
             else  
                ilast=TOTCLLENGTH  
             endif  
             badcl=badseed  
             do igood=-ngoodstr,ngoodstr  
                ibad=1  
                if((INDMAX(icy)+igood).gt.ifirst.and.  
      $              (INDMAX(icy)+igood).lt.ilast.and.  
      $              .true.)  
      $              ibad=BAD(VIEW(icy),  
      $              nvk(MAXS(icy)+igood),  
      $              nst(MAXS(icy)+igood))  
                badcl=badcl*ibad  
             enddo  
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
               
               
             cl_good(icy)=1                    
             nply=npl(VIEW(icy))  
             nldy=nld(MAXS(icy),VIEW(icy))  
               
 *     ----------------------------------------------  
 *     CONDITION TO FORM A COUPLE  
 *     ----------------------------------------------  
 *     geometrical consistency (same plane and ladder)  
             if(nply.eq.nplx.and.nldy.eq.nldx)then  
 *     charge correlation  
 *     ===========================================================  
 *     this version of the subroutine is used for the calibration  
 *     thus charge-correlation selection is obviously removed  
 *     ===========================================================  
 c$$$               ddd=(dedx(icy)  
 c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 *     ===========================================================  
                 
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
1500                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1501                    if(DEBUG)print*,                    if(verbose)print*,
1502       $                    ' ** warning ** number of identified'//       $                 '** warning ** number of identified '//
1503       $                    ' couples on plane ',nplx,       $                 'couples on plane ',nplx,
1504       $                    ' exceeds vector dimention'//       $                 'exceeds vector dimention '
1505       $                    ' ( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1506  c     good2=.false.                    mask_view(nviewx(nplx)) = 2
1507  c     goto 880   !fill ntp and go to next event                    mask_view(nviewy(nply)) = 2
1508                    iflag=1                    goto 10
                   return  
1509                 endif                 endif
1510                                
1511  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  *     ------------------> COUPLE <------------------
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
1512                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1513                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1514                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1515                 cl_single(icx)=0                 cl_single(icx)=0
1516                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1517  *     ----------------------------------------------  *     ----------------------------------------------
1518    
1519                endif                              
1520    
1521   20         continue   20         continue
1522           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1523                    
# Line 2083  c$$$               endif Line 1542  c$$$               endif
1542        endif        endif
1543                
1544        do ip=1,6        do ip=1,6
1545           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1546        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1547                
1548        return        return
1549        end        end
   
 c$$$      subroutine cl_to_couples_2(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$      include 'level1.f'  
 c$$$  
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$*         print*,'icx ',icx,badcl  
 c$$$         if(badcl.eq.0)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$*            print*,'icy ',icy,badcl  
 c$$$            if(badcl.eq.0)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$  
 c$$$c$$$*     charge correlation  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(DEBUG)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
1550                
1551  ***************************************************  ***************************************************
1552  *                                                 *  *                                                 *
# Line 2321  c$$$      end Line 1558  c$$$      end
1558  **************************************************  **************************************************
1559    
1560        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1561    
1562        include 'commontracker.f'        include 'commontracker.f'
1563          include 'level1.f'
1564        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1565        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1566        include 'common_mini_2.f'        include 'common_mini_2.f'
1567        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1568    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1569    
1570  *     output flag  *     output flag
1571  *     --------------  *     --------------
# Line 2367  c      double precision xm3,ym3,zm3 Line 1598  c      double precision xm3,ym3,zm3
1598  *     -----------------------------  *     -----------------------------
1599    
1600    
1601    *     --------------------------------------------
1602    *     put a limit to the maximum number of couples
1603    *     per plane, in order to apply hough transform
1604    *     (couples recovered during track refinement)
1605    *     --------------------------------------------
1606          do ip=1,nplanes
1607             if(ncp_plane(ip).gt.ncouplelimit)then
1608                mask_view(nviewx(ip)) = 8
1609                mask_view(nviewy(ip)) = 8
1610             endif
1611          enddo
1612    
1613    
1614        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1615        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1616                
1617        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1618           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1619                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1620             do is1=1,2             !loop on sensors - COPPIA 1            
1621              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1622                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1623                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
# Line 2383  c               call xyz_PAM(icx1,icy1,i Line 1627  c               call xyz_PAM(icx1,icy1,i
1627                 ym1=yPAM                 ym1=yPAM
1628                 zm1=zPAM                                   zm1=zPAM                  
1629  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1630    
1631                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1632                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1633         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1634                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1635                                            
1636                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2402  c     $                       (icx2,icy2 Line 1649  c     $                       (icx2,icy2
1649  *     (2 couples needed)  *     (2 couples needed)
1650  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1651                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1652                             if(DEBUG)print*,                             if(verbose)print*,
1653       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1654       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1655       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1656  c                           good2=.false.  c                           good2=.false.
1657  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1658                               do iv=1,12
1659                                  mask_view(iv) = 3
1660                               enddo
1661                             iflag=1                             iflag=1
1662                             return                             return
1663                          endif                          endif
# Line 2441  c$$$ Line 1691  c$$$
1691  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1692    
1693    
1694                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1695    
1696                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1697                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1698         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1699                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1700                                                                
1701                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2472  c     $                                 Line 1725  c     $                                
1725  *     (3 couples needed)  *     (3 couples needed)
1726  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1727                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
1728                                      if(DEBUG)print*,                                      if(verbose)print*,
1729       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
1730       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
1731       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
1732  c                                    good2=.false.  c                                    good2=.false.
1733  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
1734                                        do iv=1,nviews
1735                                           mask_view(iv) = 4
1736                                        enddo
1737                                      iflag=1                                      iflag=1
1738                                      return                                      return
1739                                   endif                                   endif
# Line 2516  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1772  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1772                                endif                                endif
1773                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1774                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1775     30                     continue
1776                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1777   30                  continue   31                  continue
1778                                            
1779   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1780                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1781     20            continue
1782              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1783                            
1784           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1785        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1786     10   continue
1787        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1788                
1789        if(DEBUG)then        if(DEBUG)then
# Line 2552  c     goto 880               !ntp fill Line 1811  c     goto 880               !ntp fill
1811        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
1812    
1813        include 'commontracker.f'        include 'commontracker.f'
1814          include 'level1.f'
1815        include 'common_momanhough.f'        include 'common_momanhough.f'
1816        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1817    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1818    
1819  *     output flag  *     output flag
1820  *     --------------  *     --------------
# Line 2588  c      common/dbg/DEBUG Line 1846  c      common/dbg/DEBUG
1846        distance=0        distance=0
1847        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
1848        npt_tot=0        npt_tot=0
1849          nloop=0                  
1850     90   continue                  
1851        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
1852           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
1853                            
# Line 2691  c     print*,'*   idbref,idb2 ',idbref,i Line 1951  c     print*,'*   idbref,idb2 ',idbref,i
1951              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
1952           enddo           enddo
1953  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
1954           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
1955           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
1956           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
1957                    
# Line 2699  c     print*,'>>>> ',ncpused,npt,nplused Line 1959  c     print*,'>>>> ',ncpused,npt,nplused
1959  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
1960    
1961           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
1962              if(DEBUG)print*,              if(verbose)print*,
1963       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
1964       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
1965       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
1966  c               good2=.false.  c               good2=.false.
1967  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
1968                do iv=1,nviews
1969                   mask_view(iv) = 5
1970                enddo
1971              iflag=1              iflag=1
1972              return              return
1973           endif           endif
# Line 2742  c$$$     $           ,(db_cloud(iii),iii Line 2005  c$$$     $           ,(db_cloud(iii),iii
2005        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2006                
2007                
2008          if(nloop.lt.nstepy)then      
2009            cutdistyz = cutdistyz+cutystep
2010            nloop     = nloop+1          
2011            goto 90                
2012          endif                    
2013          
2014        if(DEBUG)then        if(DEBUG)then
2015           print*,'---------------------- '           print*,'---------------------- '
2016           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2768  c$$$     $           ,(db_cloud(iii),iii Line 2037  c$$$     $           ,(db_cloud(iii),iii
2037        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2038    
2039        include 'commontracker.f'        include 'commontracker.f'
2040          include 'level1.f'
2041        include 'common_momanhough.f'        include 'common_momanhough.f'
2042        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2043    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2044    
2045  *     output flag  *     output flag
2046  *     --------------  *     --------------
# Line 2803  c      common/dbg/DEBUG Line 2071  c      common/dbg/DEBUG
2071        distance=0        distance=0
2072        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2073        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2074          nloop=0                  
2075     91   continue                  
2076        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2077           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
2078  c     print*,'--------------'  c     print*,'--------------'
# Line 2904  c     print*,'check cp_used' Line 2174  c     print*,'check cp_used'
2174           do ip=1,nplanes           do ip=1,nplanes
2175              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2176           enddo           enddo
2177           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2178           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2179           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2180                    
2181  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2182  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2183           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2184              if(DEBUG)print*,              if(verbose)print*,
2185       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2186       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2187       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2188  c     good2=.false.  c     good2=.false.
2189  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2190                do iv=1,nviews
2191                   mask_view(iv) = 6
2192                enddo
2193              iflag=1              iflag=1
2194              return              return
2195           endif           endif
# Line 2952  c$$$     $           ,(tr_cloud(iii),iii Line 2225  c$$$     $           ,(tr_cloud(iii),iii
2225  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2226  22288    continue  22288    continue
2227        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2228          
2229           if(nloop.lt.nstepx)then      
2230             cutdistxz=cutdistxz+cutxstep
2231             nloop=nloop+1          
2232             goto 91                
2233           endif                    
2234          
2235        if(DEBUG)then        if(DEBUG)then
2236           print*,'---------------------- '           print*,'---------------------- '
2237           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2979  c     02/02/2006 modified by Elena Vannu Line 2258  c     02/02/2006 modified by Elena Vannu
2258  c*****************************************************  c*****************************************************
2259    
2260        include 'commontracker.f'        include 'commontracker.f'
2261          include 'level1.f'
2262        include 'common_momanhough.f'        include 'common_momanhough.f'
2263        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2264        include 'common_mini_2.f'        include 'common_mini_2.f'
2265        include 'common_mech.f'        include 'common_mech.f'
2266        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2267    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2268    
2269  *     output flag  *     output flag
2270  *     --------------  *     --------------
# Line 3002  c      common/dbg/DEBUG Line 2280  c      common/dbg/DEBUG
2280  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2281  *     list of matching couples in the combination  *     list of matching couples in the combination
2282  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2283        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2284        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2285  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2286        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2287  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2288  *     variables for track fitting  *     variables for track fitting
2289        double precision AL_INI(5)        double precision AL_INI(5)
2290        double precision tath  c      double precision tath
2291  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2292  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2293    
# Line 3075  c      real fitz(nplanes)        !z coor Line 2353  c      real fitz(nplanes)        !z coor
2353                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2354              enddo              enddo
2355                            
2356              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2357                if(nplused.lt.nplyz_min)goto 888 !next doublet
2358              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2359                            
2360              if(DEBUG)then              if(DEBUG)then
# Line 3102  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2381  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2381  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2382                            
2383  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2384              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2385              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2386              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2387       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2388              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2389              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2390              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2391                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2392  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2393              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2394                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2395  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2396  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2397                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2398  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2399  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2400                c$$$            ELSE
2401    c$$$               AL_INI(4) = acos(-1.)/2
2402    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2403    c$$$            ENDIF
2404    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2405    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2406    c$$$            
2407    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2408    c$$$            
2409    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2410                            
2411              if(DEBUG)then              if(DEBUG)then
2412                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2413                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3186  c     $                                 Line 2475  c     $                                
2475  *     **********************************************************  *     **********************************************************
2476  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2477  *     **********************************************************  *     **********************************************************
2478    cccc  scommentare se si usa al_ini della nuvola
2479    c$$$                              do i=1,5
2480    c$$$                                 AL(i)=AL_INI(i)
2481    c$$$                              enddo
2482                                  call guess()
2483                                do i=1,5                                do i=1,5
2484                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2485                                enddo                                enddo
2486                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2487                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2488                                call mini_2(jstep,ifail)                                iprint=0
2489    c                              if(DEBUG)iprint=1
2490                                  if(DEBUG)iprint=1
2491                                  call mini2(jstep,ifail,iprint)
2492                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2493                                   if(DEBUG)then                                   if(DEBUG)then
2494                                      print *,                                      print *,
2495       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2496       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2497                                        print*,'initial guess: '
2498    
2499                                        print*,'AL_INI(1) = ',AL_INI(1)
2500                                        print*,'AL_INI(2) = ',AL_INI(2)
2501                                        print*,'AL_INI(3) = ',AL_INI(3)
2502                                        print*,'AL_INI(4) = ',AL_INI(4)
2503                                        print*,'AL_INI(5) = ',AL_INI(5)
2504                                   endif                                   endif
2505                                   chi2=-chi2  c                                 chi2=-chi2
2506                                endif                                endif
2507  *     **********************************************************  *     **********************************************************
2508  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3211  c     $                                 Line 2515  c     $                                
2515  *     --------------------------  *     --------------------------
2516                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2517                                                                    
2518                                   if(DEBUG)print*,                                   if(verbose)print*,
2519       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2520       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2521       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2522  c                                 good2=.false.  c                                 good2=.false.
2523  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2524                                     do iv=1,nviews
2525                                        mask_view(iv) = 7
2526                                     enddo
2527                                   iflag=1                                   iflag=1
2528                                   return                                   return
2529                                endif                                endif
# Line 3315  cccccc 12/08/2006 modified by elena vann Line 2622  cccccc 12/08/2006 modified by elena vann
2622  c******************************************************  c******************************************************
2623    
2624        include 'commontracker.f'        include 'commontracker.f'
2625          include 'level1.f'
2626        include 'common_momanhough.f'        include 'common_momanhough.f'
2627        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2628        include 'common_mini_2.f'        include 'common_mini_2.f'
2629        include 'common_mech.f'        include 'common_mech.f'
2630        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2631        include 'level1.f'  c      include 'level1.f'
2632        include 'calib.f'        include 'calib.f'
2633    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2634    
2635  *     flag to chose PFA  *     flag to chose PFA
2636        character*10 PFA        character*10 PFA
# Line 3338  c      common/dbg/DEBUG Line 2644  c      common/dbg/DEBUG
2644        call track_init        call track_init
2645        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2646    
2647    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2648  *     -------------------------------------------------  *     -------------------------------------------------
2649  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2650  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2651  *     using improved PFAs  *     using improved PFAs
2652  *     -------------------------------------------------  *     -------------------------------------------------
2653    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2654           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2655       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2656                            
# Line 3377  c            dedxtrk(nplanes-ip+1) = (de Line 2685  c            dedxtrk(nplanes-ip+1) = (de
2685              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2686              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2687                            
2688    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2689  *     -------------------------------------------------  *     -------------------------------------------------
2690  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
2691  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
2692  *     -------------------------------------------------  *     -------------------------------------------------
2693    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2694           else                             else                  
2695                                
2696              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3435  c     $              'ETA2','ETA2', Line 2745  c     $              'ETA2','ETA2',
2745       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest))
2746                                
2747                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2748                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2749                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2750                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2751       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3505  c     $              'ETA2','ETA2', Line 2816  c     $              'ETA2','ETA2',
2816       $              PFA,PFA,       $              PFA,PFA,
2817       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
2818                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2819  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2820                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
2821       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2822                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3533  c     $              'ETA2','ETA2', Line 2844  c     $              'ETA2','ETA2',
2844       $              PFA,PFA,       $              PFA,PFA,
2845       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest))
2846                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2847                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2848                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
2849       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2850                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3572  c     $                 'ETA2','ETA2', Line 2884  c     $                 'ETA2','ETA2',
2884                 endif                 endif
2885    
2886                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2887                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2888                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
2889       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
2890                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3601  c                  dedxmm = dedx(icl)   Line 2914  c                  dedxmm = dedx(icl)  
2914                                
2915                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
2916  *              ----------------------------  *              ----------------------------
2917    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2918                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
2919                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
2920                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
2921                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
2922       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
2923         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2924         $                 ,', norm.dist.= ',distmin
2925         $                 ,', cut ',clinc,' )'
2926                 else                 else
2927                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
2928                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
2929                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
2930       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
2931         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2932         $                 ,', norm.dist.= ', distmin
2933         $                 ,', cut ',clinc,' )'
2934                 endif                 endif
2935    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2936  *              ----------------------------  *              ----------------------------
2937                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
2938                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3645  cccccc 12/08/2006 modified by elena ---> Line 2966  cccccc 12/08/2006 modified by elena --->
2966        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
2967    
2968        include 'commontracker.f'        include 'commontracker.f'
2969          include 'level1.f'
2970        include 'common_momanhough.f'        include 'common_momanhough.f'
2971        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2972        include 'level2.f'        !(1)        include 'level2.f'        !(1)
2973  c      include 'calib.f'  c      include 'calib.f'
2974  c      include 'level1.f'  c      include 'level1.f'
2975    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2976    
2977    
2978        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3725  c               endif Line 3045  c               endif
3045    
3046    
3047    
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      real function fbad_cog(ncog,ic)  
 c$$$  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$*     --> signal of the central strip  
 c$$$      sc = CLSIGNAL(INDMAX(ic)) !center  
 c$$$  
 c$$$*     signal of adjacent strips  
 c$$$*     --> left  
 c$$$      sl1 = 0                  !left 1  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
 c$$$  
 c$$$      sl2 = 0                  !left 2  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
 c$$$  
 c$$$*     --> right  
 c$$$      sr1 = 0                  !right 1  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
 c$$$  
 c$$$      sr2 = 0                  !right 2  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
 c$$$  
 c$$$  
 c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  
 c$$$         f  = 4.  
 c$$$         si = 8.4  
 c$$$      else                              !X-view  
 c$$$         f  = 6.  
 c$$$         si = 3.9  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = 1.  
 c$$$      f0 = 1  
 c$$$      f1 = 1  
 c$$$      f2 = 1  
 c$$$      f3 = 1    
 c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then  
 c$$$          
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$          
 c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then  
 c$$$  
 c$$$  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sr1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = sqrt(fbad_cog)  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
3048    
3049    
3050  *     ****************************************************  *     ****************************************************
3051    
3052        subroutine init_level2        subroutine init_level2
3053    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3054        include 'commontracker.f'        include 'commontracker.f'
3055          include 'level1.f'
3056        include 'common_momanhough.f'        include 'common_momanhough.f'
3057        include 'level2.f'        include 'level2.f'
3058        include 'level1.f'  c      include 'level1.f'
3059    
3060        do i=1,nviews        do i=1,nviews
3061           good2(i)=good1(i)           good2(i)=good1(i)
3062        enddo        enddo
3063    
 c      good2 = 0!.false.  
 c$$$      nev2 = nev1  
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
3064    
3065        NTRK = 0        NTRK = 0
3066        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3067           IMAGE(IT)=0           IMAGE(IT)=0
3068           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3069           do ip=1,nplanes           do ip=1,nplanes
3070              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3071              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3877  c         BdL(IT) = 0. Line 3074  c         BdL(IT) = 0.
3074              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3075              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3076              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3077              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3078              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3079              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3080              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3081           enddo           enddo
# Line 3893  cccccc 17/8/2006 modified by elena Line 3086  cccccc 17/8/2006 modified by elena
3086              enddo                                enddo                  
3087           enddo                             enddo                  
3088        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3089        nclsx=0        nclsx=0
3090        nclsy=0              nclsy=0      
3091        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3092          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3093          xs(1,ip)=0          xs(1,ip)=0
3094          xs(2,ip)=0          xs(2,ip)=0
3095          sgnlxs(ip)=0          sgnlxs(ip)=0
3096          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3097          ys(1,ip)=0          ys(1,ip)=0
3098          ys(2,ip)=0          ys(2,ip)=0
3099          sgnlys(ip)=0          sgnlys(ip)=0
3100        enddo        enddo
 c*******************************************************  
3101        end        end
3102    
3103    
# Line 3926  c*************************************** Line 3112  c***************************************
3112  ************************************************************  ************************************************************
3113    
3114    
3115          subroutine init_hough
3116    
3117          include 'commontracker.f'
3118          include 'level1.f'
3119          include 'common_momanhough.f'
3120          include 'common_hough.f'
3121          include 'level2.f'
3122    
3123          ntrpt_nt=0
3124          ndblt_nt=0
3125          NCLOUDS_XZ_nt=0
3126          NCLOUDS_YZ_nt=0
3127          do idb=1,ndblt_max_nt
3128             db_cloud_nt(idb)=0
3129             alfayz1_nt(idb)=0      
3130             alfayz2_nt(idb)=0      
3131          enddo
3132          do itr=1,ntrpt_max_nt
3133             tr_cloud_nt(itr)=0
3134             alfaxz1_nt(itr)=0      
3135             alfaxz2_nt(itr)=0      
3136             alfaxz3_nt(itr)=0      
3137          enddo
3138          do idb=1,ncloyz_max      
3139            ptcloud_yz_nt(idb)=0    
3140            alfayz1_av_nt(idb)=0    
3141            alfayz2_av_nt(idb)=0    
3142          enddo                    
3143          do itr=1,ncloxz_max      
3144            ptcloud_xz_nt(itr)=0    
3145            alfaxz1_av_nt(itr)=0    
3146            alfaxz2_av_nt(itr)=0    
3147            alfaxz3_av_nt(itr)=0    
3148          enddo                    
3149    
3150          ntrpt=0                  
3151          ndblt=0                  
3152          NCLOUDS_XZ=0              
3153          NCLOUDS_YZ=0              
3154          do idb=1,ndblt_max        
3155            db_cloud(idb)=0        
3156            cpyz1(idb)=0            
3157            cpyz2(idb)=0            
3158            alfayz1(idb)=0          
3159            alfayz2(idb)=0          
3160          enddo                    
3161          do itr=1,ntrpt_max        
3162            tr_cloud(itr)=0        
3163            cpxz1(itr)=0            
3164            cpxz2(itr)=0            
3165            cpxz3(itr)=0            
3166            alfaxz1(itr)=0          
3167            alfaxz2(itr)=0          
3168            alfaxz3(itr)=0          
3169          enddo                    
3170          do idb=1,ncloyz_max      
3171            ptcloud_yz(idb)=0      
3172            alfayz1_av(idb)=0      
3173            alfayz2_av(idb)=0      
3174            do idbb=1,ncouplemaxtot
3175              cpcloud_yz(idb,idbb)=0
3176            enddo                  
3177          enddo                    
3178          do itr=1,ncloxz_max      
3179            ptcloud_xz(itr)=0      
3180            alfaxz1_av(itr)=0      
3181            alfaxz2_av(itr)=0      
3182            alfaxz3_av(itr)=0      
3183            do itrr=1,ncouplemaxtot
3184              cpcloud_xz(itr,itrr)=0
3185            enddo                  
3186          enddo                    
3187          end
3188    ************************************************************
3189    *
3190    *
3191    *
3192    *
3193    *
3194    *
3195    *
3196    ************************************************************
3197    
3198    
3199        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3200    
3201  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3936  c*************************************** Line 3206  c***************************************
3206    
3207            
3208        include 'commontracker.f'        include 'commontracker.f'
3209    c      include 'level1.f'
3210        include 'level1.f'        include 'level1.f'
3211          include 'common_momanhough.f'
3212        include 'level2.f'        include 'level2.f'
3213        include 'common_mini_2.f'        include 'common_mini_2.f'
3214        include 'common_momanhough.f'        real sinth,phi,pig      
       real sinth,phi,pig        !(4)  
3215        pig=acos(-1.)        pig=acos(-1.)
3216    
 c      good2=1!.true.  
3217        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3218        nstep_nt(ntr)       = 0!nstep        nstep_nt(ntr)       = nstep
3219    
3220          phi   = al(4)          
3221          sinth = al(3)            
3222          if(sinth.lt.0)then      
3223             sinth = -sinth        
3224             phi = phi + pig      
3225          endif                    
3226          npig = aint(phi/(2*pig))
3227          phi = phi - npig*2*pig  
3228          if(phi.lt.0)            
3229         $     phi = phi + 2*pig  
3230          al(4) = phi              
3231          al(3) = sinth            
3232    
       phi   = al(4)             !(4)  
       sinth = al(3)             !(4)  
       if(sinth.lt.0)then        !(4)  
          sinth = -sinth         !(4)  
          phi = phi + pig        !(4)  
       endif                     !(4)  
       npig = aint(phi/(2*pig))  !(4)  
       phi = phi - npig*2*pig    !(4)  
       if(phi.lt.0)              !(4)  
      $     phi = phi + 2*pig    !(4)  
       al(4) = phi               !(4)  
       al(3) = sinth             !(4)  
 *****************************************************  
3233        do i=1,5        do i=1,5
3234           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3235           do j=1,5           do j=1,5
3236              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3237           enddo           enddo
 c     print*,al_nt(i,ntr)  
3238        enddo        enddo
3239                
3240        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3981  c     print*,al_nt(i,ntr) Line 3250  c     print*,al_nt(i,ntr)
3250           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3251           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3252           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))
 c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3253           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
3254           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)  
3255        
# Line 3998  c            print*,ip,' ',cltrx(ip,ntr) Line 3266  c            print*,ip,' ',cltrx(ip,ntr)
3266           endif                     endif          
3267    
3268        enddo        enddo
 c      call CalcBdL(100,xxxx,IFAIL)  
 c      if(ifps(xxxx).eq.1)BdL(ntr) = xxxx  
 c$$$      print*,'xgood(ip,ntr) ',(xgood_nt(ip,ntr),ip=1,6)  
 c$$$      print*,'ygood(ip,ntr) ',(ygood_nt(ip,ntr),ip=1,6)  
 c$$$      print*,'dedx_x(ip,ntr) ',(dedx_x(ip,ntr),ip=1,6)  
 c$$$      print*,'dedx_y(ip,ntr) ',(dedx_y(ip,ntr),ip=1,6)  
3269    
3270    
3271        end        end
3272    
3273        subroutine fill_level2_siglets        subroutine fill_level2_siglets
 c*****************************************************  
 c     07/10/2005 created by elena vannuccini  
 c     31/01/2006 modified by elena vannuccini  
 *     to convert adc to mip  --> (2)  
 c*****************************************************  
3274    
3275  *     -------------------------------------------------------  *     -------------------------------------------------------
3276  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 4022  c*************************************** Line 3279  c***************************************
3279  *     -------------------------------------------------------  *     -------------------------------------------------------
3280    
3281        include 'commontracker.f'        include 'commontracker.f'
3282        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
3283        include 'calib.f'        include 'calib.f'
3284          include 'level1.f'
3285        include 'common_momanhough.f'        include 'common_momanhough.f'
3286          include 'level2.f'
3287        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3288    
3289  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
# Line 4033  c      good2=1!.true. Line 3291  c      good2=1!.true.
3291        nclsx = 0        nclsx = 0
3292        nclsy = 0        nclsy = 0
3293    
3294          do iv = 1,nviews
3295             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3296          enddo
3297    
3298        do icl=1,nclstr1        do icl=1,nclstr1
3299           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
3300              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
# Line 4076  c      print*,icl,cl_used(icl),cl_good(i Line 3338  c      print*,icl,cl_used(icl),cl_good(i
3338        enddo        enddo
3339        end        end
3340    
3341    ***************************************************
3342    *                                                 *
3343    *                                                 *
3344    *                                                 *
3345    *                                                 *
3346    *                                                 *
3347    *                                                 *
3348    **************************************************
3349    
3350          subroutine fill_hough
3351    
3352    *     -------------------------------------------------------
3353    *     This routine fills the  variables related to the hough
3354    *     transform, for the debig n-tuple
3355    *     -------------------------------------------------------
3356    
3357          include 'commontracker.f'
3358          include 'level1.f'
3359          include 'common_momanhough.f'
3360          include 'common_hough.f'
3361          include 'level2.f'
3362    
3363          if(.false.
3364         $     .or.ntrpt.gt.ntrpt_max_nt
3365         $     .or.ndblt.gt.ndblt_max_nt
3366         $     .or.NCLOUDS_XZ.gt.ncloxz_max
3367         $     .or.NCLOUDS_yZ.gt.ncloyz_max
3368         $     )then
3369             ntrpt_nt=0
3370             ndblt_nt=0
3371             NCLOUDS_XZ_nt=0
3372             NCLOUDS_YZ_nt=0
3373          else
3374             ndblt_nt=ndblt
3375             ntrpt_nt=ntrpt
3376             if(ndblt.ne.0)then
3377                do id=1,ndblt
3378                   alfayz1_nt(id)=alfayz1(id) !Y0
3379                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3380                enddo
3381             endif
3382             if(ndblt.ne.0)then
3383                do it=1,ntrpt
3384                   alfaxz1_nt(it)=alfaxz1(it) !X0
3385                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3386                   alfaxz3_nt(it)=alfaxz3(it) !1/r
3387                enddo
3388             endif
3389             nclouds_yz_nt=nclouds_yz
3390             nclouds_xz_nt=nclouds_xz
3391             if(nclouds_yz.ne.0)then
3392                nnn=0
3393                do iyz=1,nclouds_yz
3394                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3395                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3396                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3397                   nnn=nnn+ptcloud_yz(iyz)
3398                enddo
3399                do ipt=1,nnn
3400                   db_cloud_nt(ipt)=db_cloud(ipt)
3401                 enddo
3402             endif
3403             if(nclouds_xz.ne.0)then
3404                nnn=0
3405                do ixz=1,nclouds_xz
3406                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
3407                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
3408                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
3409                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
3410                   nnn=nnn+ptcloud_xz(ixz)              
3411                enddo
3412                do ipt=1,nnn
3413                  tr_cloud_nt(ipt)=tr_cloud(ipt)
3414                 enddo
3415             endif
3416          endif
3417          end
3418          

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.15

  ViewVC Help
Powered by ViewVC 1.1.23