/[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.3 by pam-fi, Tue Sep 5 12:52:20 2006 UTC revision 1.16 by pam-fi, Thu Nov 30 17:04:27 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             if(DEBUG)iprint=2
330             call mini2(jstep,ifail,iprint)
331           if(ifail.ne.0) then           if(ifail.ne.0) then
332              if(DEBUG)then              if(VERBOSE)then
333                 print *,                 print *,
334       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
335       $              ,iev       $              ,iev
336    
337    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
338    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
339    c$$$               print*,'result:   ',(al(i),i=1,5)
340    c$$$               print*,'xgood ',xgood
341    c$$$               print*,'ygood ',ygood
342    c$$$               print*,'----------------------------------------------'
343              endif              endif
344              chi2=-chi2  c            chi2=-chi2
345           endif           endif
346                    
347           if(DEBUG)then           if(DEBUG)then
# Line 311  c         print*,'++++++++++ iimage,fima Line 402  c         print*,'++++++++++ iimage,fima
402  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
403    
404           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
405              if(DEBUG)              if(verbose)
406       $           print*,       $           print*,
407       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
408       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 361  c     $        rchi2best.lt.15..and. Line 452  c     $        rchi2best.lt.15..and.
452        end        end
453    
454    
   
   
 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$$$  
   
455                
456  ************************************************************  ************************************************************
457  ************************************************************  ************************************************************
# Line 607  c                (implemented variable r Line 526  c                (implemented variable r
526  c*****************************************************  c*****************************************************
527                
528        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
529        include 'level1.f'        include 'level1.f'
530          include 'calib.f'
531    c      include 'level1.f'
532        include 'common_align.f'        include 'common_align.f'
533        include 'common_mech.f'        include 'common_mech.f'
534        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
535        include 'common_resxy.f'  c      include 'common_resxy.f'
536    
537  c      logical DEBUG  c      logical DEBUG
538  c      common/dbg/DEBUG  c      common/dbg/DEBUG
# Line 667  c      double precision xi_B,yi_B,zi_B Line 587  c      double precision xi_B,yi_B,zi_B
587              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
588           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
589  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
590  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
591  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
592              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
593              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
594              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
595       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
596              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
597           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
598              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
599              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
600              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
601       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
602              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
603           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
604              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
605              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
606              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
607       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
608              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
609           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
610              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
611              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
612              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
613       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 702  c            resxPAM = resxPAM*fbad_cog( Line 622  c            resxPAM = resxPAM*fbad_cog(
622           endif           endif
623    
624        endif        endif
625    c      if(icy.eq.0.and.icx.ne.0)
626    c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'
627                
628  *     -----------------  *     -----------------
629  *     CLUSTER Y  *     CLUSTER Y
# Line 729  c            resxPAM = resxPAM*fbad_cog( Line 651  c            resxPAM = resxPAM*fbad_cog(
651              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
652           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
653  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
654  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
655  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
656              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
657              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
658              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
659              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
660       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
661           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
662              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
663              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
664              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
665       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
666           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
667              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
668              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
669              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
670       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
671           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
672              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
673              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
674  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
675              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 773  C======================================= Line 695  C=======================================
695  c------------------------------------------------------------------------  c------------------------------------------------------------------------
696  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
697  c------------------------------------------------------------------------  c------------------------------------------------------------------------
698             if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
699         $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
700                print*,'xyz_PAM (couple):',
701         $          ' WARNING: false X strip: strip ',stripx
702             endif
703           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
704           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
705           zi = 0.           zi = 0.
# Line 859  C======================================= Line 786  C=======================================
786              nldy = nldx              nldy = nldx
787              viewy = viewx - 1              viewy = viewx - 1
788    
789    c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
790    c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
791                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
792         $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
793                   print*,'xyz_PAM (X-singlet):',
794         $             ' WARNING: false X strip: strip ',stripx
795                endif
796              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
797    
798              xi_A = xi              xi_A = xi
# Line 1137  c$$$         print*,' resolution ',resxP Line 1071  c$$$         print*,' resolution ',resxP
1071  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1072  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1073  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1074                   if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1075         $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1076    c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1077                      print*,'whichsensor: ',
1078         $                ' WARNING: false X strip: strip ',stripx
1079                   endif
1080                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1081                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
1082                 zi = 0.                 zi = 0.
# Line 1264  c     $              ,iv,xvv(iv),yvv(iv) Line 1204  c     $              ,iv,xvv(iv),yvv(iv)
1204  *     it returns the plane number  *     it returns the plane number
1205  *      *    
1206        include 'commontracker.f'        include 'commontracker.f'
1207          include 'level1.f'
1208  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1209        include 'common_momanhough.f'        include 'common_momanhough.f'
1210                
# Line 1301  c      include 'common_analysis.f' Line 1242  c      include 'common_analysis.f'
1242  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1243  *      *    
1244        include 'commontracker.f'        include 'commontracker.f'
1245          include 'level1.f'
1246  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1247        include 'common_momanhough.f'        include 'common_momanhough.f'
1248                
# Line 1329  c      include 'common_analysis.f' Line 1271  c      include 'common_analysis.f'
1271  *     positive if sensor =2  *     positive if sensor =2
1272  *  *
1273        include 'commontracker.f'        include 'commontracker.f'
1274          include 'level1.f'
1275  c      include 'calib.f'  c      include 'calib.f'
1276  c      include 'level1.f'  c      include 'level1.f'
1277  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1358  c      include 'common_analysis.f' Line 1301  c      include 'common_analysis.f'
1301  *************************************************************************  *************************************************************************
1302  *************************************************************************  *************************************************************************
1303  *************************************************************************  *************************************************************************
 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  
1304                
1305    
1306  ***************************************************  ***************************************************
# Line 1640  c$$$      end Line 1315  c$$$      end
1315        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1316    
1317        include 'commontracker.f'        include 'commontracker.f'
1318          include 'level1.f'
1319        include 'common_momanhough.f'        include 'common_momanhough.f'
1320        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1321        include 'calib.f'        include 'calib.f'
1322        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1323    
1324  *     output flag  *     output flag
1325  *     --------------  *     --------------
# Line 1655  c      common/dbg/DEBUG Line 1328  c      common/dbg/DEBUG
1328  *     --------------  *     --------------
1329        integer iflag        integer iflag
1330    
1331        integer badseed,badcl        integer badseed,badclx,badcly
1332    
1333  *     init variables  *     init variables
1334        ncp_tot=0        ncp_tot=0
# Line 1671  c      common/dbg/DEBUG Line 1344  c      common/dbg/DEBUG
1344           ncls(ip)=0           ncls(ip)=0
1345        enddo        enddo
1346        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1347           cl_single(icl)=1           cl_single(icl) = 1
1348           cl_good(icl)=0           cl_good(icl)   = 0
1349          enddo
1350          do iv=1,nviews
1351             ncl_view(iv)  = 0
1352             mask_view(iv) = 0      !all included
1353        enddo        enddo
1354                
1355    *     count number of cluster per view
1356          do icl=1,nclstr1
1357             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1358          enddo
1359    *     mask views with too many clusters
1360          do iv=1,nviews
1361             if( ncl_view(iv).gt. nclusterlimit)then
1362                mask_view(iv) = 1
1363                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1364         $           ,nclusterlimit,' on view ', iv,' --> masked!'
1365             endif
1366          enddo
1367    
1368    
1369  *     start association  *     start association
1370        ncouples=0        ncouples=0
1371        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1372           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1373                    
1374  *     ----------------------------------------------------  *     ----------------------------------------------------
1375    *     jump masked views (X VIEW)
1376    *     ----------------------------------------------------
1377             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1378    *     ----------------------------------------------------
1379  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1380    *     ----------------------------------------------------
1381           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
1382              cl_single(icx)=0              cl_single(icx)=0
1383              goto 10              goto 10
1384           endif           endif
1385    *     ----------------------------------------------------
1386  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1387    *     ----------------------------------------------------
1388           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1389           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1390           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1694  c      common/dbg/DEBUG Line 1392  c      common/dbg/DEBUG
1392           else           else
1393              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1394           endif           endif
1395           badcl=badseed           badclx=badseed
1396           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1397              ibad=1              ibad=1
1398              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1704  c      common/dbg/DEBUG Line 1402  c      common/dbg/DEBUG
1402       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1403       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1404              endif              endif
1405              badcl=badcl*ibad              badclx=badclx*ibad
1406           enddo           enddo
1407    *     ----------------------------------------------------
1408    *     >>> eliminato il taglio sulle BAD <<<
1409    *     ----------------------------------------------------
1410  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1411  c     cl_single(icx)=0  c     cl_single(icx)=0
1412  c     goto 10  c     goto 10
# Line 1720  c     endif Line 1421  c     endif
1421              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1422                            
1423  *     ----------------------------------------------------  *     ----------------------------------------------------
1424    *     jump masked views (Y VIEW)
1425    *     ----------------------------------------------------
1426                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1427    
1428    *     ----------------------------------------------------
1429  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1430    *     ----------------------------------------------------
1431              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
1432                 cl_single(icy)=0                 cl_single(icy)=0
1433                 goto 20                 goto 20
1434              endif              endif
1435    *     ----------------------------------------------------
1436  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1437    *     ----------------------------------------------------
1438              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1439              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1440              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 1733  c     endif Line 1442  c     endif
1442              else              else
1443                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1444              endif              endif
1445              badcl=badseed              badcly=badseed
1446              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1447                 ibad=1                 ibad=1
1448                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1742  c     endif Line 1451  c     endif
1451       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1452       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1453       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1454                 badcl=badcl*ibad                 badcly=badcly*ibad
1455              enddo              enddo
1456    *     ----------------------------------------------------
1457    *     >>> eliminato il taglio sulle BAD <<<
1458    *     ----------------------------------------------------
1459  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1460  c     cl_single(icy)=0  c     cl_single(icy)=0
1461  c     goto 20  c     goto 20
1462  c     endif  c     endif
1463  *     ----------------------------------------------------  *     ----------------------------------------------------
1464                            
               
1465              cl_good(icy)=1                                cl_good(icy)=1                  
1466              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
1467              nldy=nld(MAXS(icy),VIEW(icy))              nldy=nld(MAXS(icy),VIEW(icy))
# Line 1761  c     endif Line 1472  c     endif
1472  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1473              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1474  *     charge correlation  *     charge correlation
1475                 ddd=(dedx(icy)  *     (modified to be applied only below saturation... obviously)
      $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
                ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
                cut=chcut*sch(nplx,nldx)  
                if(abs(ddd).gt.cut)goto 20 !charge not consistent  
                 
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
                if(ncp_plane(nplx).gt.ncouplemax)then  
                   if(DEBUG)print*,  
      $                    ' ** warning ** number of identified'//  
      $                    ' couples on plane ',nplx,  
      $                    ' exceeds vector dimention'//  
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(DEBUG)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                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                                
 *     ----------------------------------------------  
1476    
1477   20         continue                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1478           enddo                  !end loop on clusters(Y)       $              .and.
1479                 $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1480   10      continue       $              .and.
1481        enddo                     !end loop on clusters(X)       $              (badclx.eq.1.and.badcly.eq.1)
1482               $              .and.
1483               $              .true.)then
1484        do icl=1,nclstr1  
1485           if(cl_single(icl).eq.1)then                    ddd=(dedx(icy)
1486              ip=npl(VIEW(icl))       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1487              ncls(ip)=ncls(ip)+1                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1488              cls(ip,ncls(ip))=icl  
1489           endif  c                  cut = chcut * sch(nplx,nldx)
1490        enddo  
1491                            sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1492               $                 -kch(nplx,nldx)*cch(nplx,nldx))
1493        if(DEBUG)then                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1494           print*,'clusters  ',nclstr1                    cut = chcut * (16 + sss/50.)
          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  
1495    
1496  *     output flag                    if(abs(ddd).gt.cut)then
1497  *     --------------                       goto 20    !charge not consistent
1498  *     0 = good event                    endif
1499  *     1 = bad event                 endif
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
1500    
 *     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  
1501                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1502                    if(DEBUG)print*,                    if(verbose)print*,
      $                    ' ** warning ** number of identified'//  
      $                    ' couples on plane ',nplx,  
      $                    ' exceeds vector dimention'//  
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(DEBUG)print*,  
1503       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1504       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1505       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1506       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1507  c     good2=.false.                    mask_view(nviewx(nplx)) = 2
1508  c     goto 880   !fill ntp and go to next event                                        mask_view(nviewy(nply)) = 2
1509                    iflag=1                    goto 10
                   return  
1510                 endif                 endif
1511                                
1512    *     ------------------> COUPLE <------------------
1513                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1514                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1515                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1516                 cl_single(icx)=0                 cl_single(icx)=0
1517                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1518  *     ----------------------------------------------  *     ----------------------------------------------
1519    
1520                endif                              
1521    
1522   20         continue   20         continue
1523           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1524                    
# Line 2046  c     goto 880   !fill ntp and go to nex Line 1543  c     goto 880   !fill ntp and go to nex
1543        endif        endif
1544                
1545        do ip=1,6        do ip=1,6
1546           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1547        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  
1548                
1549        return        return
1550        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  
1551                
1552  ***************************************************  ***************************************************
1553  *                                                 *  *                                                 *
# Line 2284  c$$$      end Line 1559  c$$$      end
1559  **************************************************  **************************************************
1560    
1561        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1562    
1563        include 'commontracker.f'        include 'commontracker.f'
1564          include 'level1.f'
1565        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1566        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1567        include 'common_mini_2.f'        include 'common_mini_2.f'
1568        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1569    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1570    
1571  *     output flag  *     output flag
1572  *     --------------  *     --------------
# Line 2330  c      double precision xm3,ym3,zm3 Line 1599  c      double precision xm3,ym3,zm3
1599  *     -----------------------------  *     -----------------------------
1600    
1601    
1602    *     --------------------------------------------
1603    *     put a limit to the maximum number of couples
1604    *     per plane, in order to apply hough transform
1605    *     (couples recovered during track refinement)
1606    *     --------------------------------------------
1607          do ip=1,nplanes
1608             if(ncp_plane(ip).gt.ncouplelimit)then
1609                mask_view(nviewx(ip)) = 8
1610                mask_view(nviewy(ip)) = 8
1611             endif
1612          enddo
1613    
1614    
1615        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1616        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1617                
1618        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1619           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1620                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1621             do is1=1,2             !loop on sensors - COPPIA 1            
1622              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1623                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1624                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
# Line 2346  c               call xyz_PAM(icx1,icy1,i Line 1628  c               call xyz_PAM(icx1,icy1,i
1628                 ym1=yPAM                 ym1=yPAM
1629                 zm1=zPAM                                   zm1=zPAM                  
1630  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1631    
1632                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1633                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1634         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1635                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1636                                            
1637                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2365  c     $                       (icx2,icy2 Line 1650  c     $                       (icx2,icy2
1650  *     (2 couples needed)  *     (2 couples needed)
1651  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1652                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1653                             if(DEBUG)print*,                             if(verbose)print*,
1654       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1655       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1656       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1657  c                           good2=.false.  c                           good2=.false.
1658  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1659                               do iv=1,12
1660                                  mask_view(iv) = 3
1661                               enddo
1662                             iflag=1                             iflag=1
1663                             return                             return
1664                          endif                          endif
# Line 2404  c$$$ Line 1692  c$$$
1692  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1693    
1694    
1695                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1696    
1697                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1698                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1699         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1700                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1701                                                                
1702                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2435  c     $                                 Line 1726  c     $                                
1726  *     (3 couples needed)  *     (3 couples needed)
1727  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1728                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
1729                                      if(DEBUG)print*,                                      if(verbose)print*,
1730       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
1731       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
1732       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
1733  c                                    good2=.false.  c                                    good2=.false.
1734  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
1735                                        do iv=1,nviews
1736                                           mask_view(iv) = 4
1737                                        enddo
1738                                      iflag=1                                      iflag=1
1739                                      return                                      return
1740                                   endif                                   endif
# Line 2479  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1773  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1773                                endif                                endif
1774                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1775                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1776     30                     continue
1777                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1778   30                  continue   31                  continue
1779                                            
1780   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1781                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1782     20            continue
1783              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1784                            
1785           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1786        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1787     10   continue
1788        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1789                
1790        if(DEBUG)then        if(DEBUG)then
# Line 2515  c     goto 880               !ntp fill Line 1812  c     goto 880               !ntp fill
1812        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
1813    
1814        include 'commontracker.f'        include 'commontracker.f'
1815          include 'level1.f'
1816        include 'common_momanhough.f'        include 'common_momanhough.f'
1817        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1818    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1819    
1820  *     output flag  *     output flag
1821  *     --------------  *     --------------
# Line 2551  c      common/dbg/DEBUG Line 1847  c      common/dbg/DEBUG
1847        distance=0        distance=0
1848        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
1849        npt_tot=0        npt_tot=0
1850          nloop=0                  
1851     90   continue                  
1852        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
1853           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
1854                            
# Line 2654  c     print*,'*   idbref,idb2 ',idbref,i Line 1952  c     print*,'*   idbref,idb2 ',idbref,i
1952              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
1953           enddo           enddo
1954  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
1955           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
1956           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
1957           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
1958                    
# Line 2662  c     print*,'>>>> ',ncpused,npt,nplused Line 1960  c     print*,'>>>> ',ncpused,npt,nplused
1960  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
1961    
1962           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
1963              if(DEBUG)print*,              if(verbose)print*,
1964       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
1965       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
1966       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
1967  c               good2=.false.  c               good2=.false.
1968  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
1969                do iv=1,nviews
1970                   mask_view(iv) = 5
1971                enddo
1972              iflag=1              iflag=1
1973              return              return
1974           endif           endif
# Line 2705  c$$$     $           ,(db_cloud(iii),iii Line 2006  c$$$     $           ,(db_cloud(iii),iii
2006        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2007                
2008                
2009          if(nloop.lt.nstepy)then      
2010            cutdistyz = cutdistyz+cutystep
2011            nloop     = nloop+1          
2012            goto 90                
2013          endif                    
2014          
2015        if(DEBUG)then        if(DEBUG)then
2016           print*,'---------------------- '           print*,'---------------------- '
2017           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2731  c$$$     $           ,(db_cloud(iii),iii Line 2038  c$$$     $           ,(db_cloud(iii),iii
2038        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2039    
2040        include 'commontracker.f'        include 'commontracker.f'
2041          include 'level1.f'
2042        include 'common_momanhough.f'        include 'common_momanhough.f'
2043        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2044    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2045    
2046  *     output flag  *     output flag
2047  *     --------------  *     --------------
# Line 2766  c      common/dbg/DEBUG Line 2072  c      common/dbg/DEBUG
2072        distance=0        distance=0
2073        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2074        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2075          nloop=0                  
2076     91   continue                  
2077        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2078           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
2079  c     print*,'--------------'  c     print*,'--------------'
# Line 2867  c     print*,'check cp_used' Line 2175  c     print*,'check cp_used'
2175           do ip=1,nplanes           do ip=1,nplanes
2176              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2177           enddo           enddo
2178           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2179           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2180           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2181                    
2182  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2183  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2184           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2185              if(DEBUG)print*,              if(verbose)print*,
2186       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2187       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2188       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2189  c     good2=.false.  c     good2=.false.
2190  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2191                do iv=1,nviews
2192                   mask_view(iv) = 6
2193                enddo
2194              iflag=1              iflag=1
2195              return              return
2196           endif           endif
# Line 2915  c$$$     $           ,(tr_cloud(iii),iii Line 2226  c$$$     $           ,(tr_cloud(iii),iii
2226  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2227  22288    continue  22288    continue
2228        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2229          
2230           if(nloop.lt.nstepx)then      
2231             cutdistxz=cutdistxz+cutxstep
2232             nloop=nloop+1          
2233             goto 91                
2234           endif                    
2235          
2236        if(DEBUG)then        if(DEBUG)then
2237           print*,'---------------------- '           print*,'---------------------- '
2238           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2942  c     02/02/2006 modified by Elena Vannu Line 2259  c     02/02/2006 modified by Elena Vannu
2259  c*****************************************************  c*****************************************************
2260    
2261        include 'commontracker.f'        include 'commontracker.f'
2262          include 'level1.f'
2263        include 'common_momanhough.f'        include 'common_momanhough.f'
2264        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2265        include 'common_mini_2.f'        include 'common_mini_2.f'
2266        include 'common_mech.f'        include 'common_mech.f'
2267        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2268    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2269    
2270  *     output flag  *     output flag
2271  *     --------------  *     --------------
# Line 2965  c      common/dbg/DEBUG Line 2281  c      common/dbg/DEBUG
2281  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2282  *     list of matching couples in the combination  *     list of matching couples in the combination
2283  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2284        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2285        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2286  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2287        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2288  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2289  *     variables for track fitting  *     variables for track fitting
2290        double precision AL_INI(5)        double precision AL_INI(5)
2291        double precision tath  c      double precision tath
2292  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2293  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2294    
# Line 3038  c      real fitz(nplanes)        !z coor Line 2354  c      real fitz(nplanes)        !z coor
2354                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2355              enddo              enddo
2356                            
2357              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2358                if(nplused.lt.nplyz_min)goto 888 !next doublet
2359              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2360                            
2361              if(DEBUG)then              if(DEBUG)then
# Line 3065  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2382  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2382  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2383                            
2384  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2385              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2386              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2387              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2388       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2389              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2390              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2391              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2392                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2393  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2394              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2395                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2396  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2397  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2398                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2399  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2400  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2401                c$$$            ELSE
2402    c$$$               AL_INI(4) = acos(-1.)/2
2403    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2404    c$$$            ENDIF
2405    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2406    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2407    c$$$            
2408    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2409    c$$$            
2410    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2411                            
2412              if(DEBUG)then              if(DEBUG)then
2413                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2414                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3149  c     $                                 Line 2476  c     $                                
2476  *     **********************************************************  *     **********************************************************
2477  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2478  *     **********************************************************  *     **********************************************************
2479    cccc  scommentare se si usa al_ini della nuvola
2480    c$$$                              do i=1,5
2481    c$$$                                 AL(i)=AL_INI(i)
2482    c$$$                              enddo
2483                                  call guess()
2484                                do i=1,5                                do i=1,5
2485                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2486                                enddo                                enddo
2487                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2488                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2489                                call mini_2(jstep,ifail)                                iprint=0
2490    c                              if(DEBUG)iprint=1
2491                                  if(DEBUG)iprint=2
2492                                  call mini2(jstep,ifail,iprint)
2493                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2494                                   if(DEBUG)then                                   if(DEBUG)then
2495                                      print *,                                      print *,
2496       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2497       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2498                                        print*,'initial guess: '
2499    
2500                                        print*,'AL_INI(1) = ',AL_INI(1)
2501                                        print*,'AL_INI(2) = ',AL_INI(2)
2502                                        print*,'AL_INI(3) = ',AL_INI(3)
2503                                        print*,'AL_INI(4) = ',AL_INI(4)
2504                                        print*,'AL_INI(5) = ',AL_INI(5)
2505                                   endif                                   endif
2506                                   chi2=-chi2  c                                 chi2=-chi2
2507                                endif                                endif
2508  *     **********************************************************  *     **********************************************************
2509  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3174  c     $                                 Line 2516  c     $                                
2516  *     --------------------------  *     --------------------------
2517                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2518                                                                    
2519                                   if(DEBUG)print*,                                   if(verbose)print*,
2520       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2521       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2522       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2523  c                                 good2=.false.  c                                 good2=.false.
2524  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2525                                     do iv=1,nviews
2526                                        mask_view(iv) = 7
2527                                     enddo
2528                                   iflag=1                                   iflag=1
2529                                   return                                   return
2530                                endif                                endif
# Line 3278  cccccc 12/08/2006 modified by elena vann Line 2623  cccccc 12/08/2006 modified by elena vann
2623  c******************************************************  c******************************************************
2624    
2625        include 'commontracker.f'        include 'commontracker.f'
2626          include 'level1.f'
2627        include 'common_momanhough.f'        include 'common_momanhough.f'
2628        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2629        include 'common_mini_2.f'        include 'common_mini_2.f'
2630        include 'common_mech.f'        include 'common_mech.f'
2631        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2632        include 'level1.f'  c      include 'level1.f'
2633        include 'calib.f'        include 'calib.f'
2634    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2635    
2636  *     flag to chose PFA  *     flag to chose PFA
2637        character*10 PFA        character*10 PFA
# Line 3301  c      common/dbg/DEBUG Line 2645  c      common/dbg/DEBUG
2645        call track_init        call track_init
2646        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2647    
2648    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2649  *     -------------------------------------------------  *     -------------------------------------------------
2650  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2651  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2652  *     using improved PFAs  *     using improved PFAs
2653  *     -------------------------------------------------  *     -------------------------------------------------
2654    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2655           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2656       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2657                            
# Line 3340  c            dedxtrk(nplanes-ip+1) = (de Line 2686  c            dedxtrk(nplanes-ip+1) = (de
2686              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)
2687              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)
2688                            
2689    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2690  *     -------------------------------------------------  *     -------------------------------------------------
2691  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
2692  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
2693  *     -------------------------------------------------  *     -------------------------------------------------
2694    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2695           else                             else                  
2696                                
2697              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3398  c     $              'ETA2','ETA2', Line 2746  c     $              'ETA2','ETA2',
2746       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest))
2747                                
2748                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2749                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2750                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2751                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2752       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3468  c     $              'ETA2','ETA2', Line 2817  c     $              'ETA2','ETA2',
2817       $              PFA,PFA,       $              PFA,PFA,
2818       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
2819                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2820  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2821                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
2822       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2823                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3496  c     $              'ETA2','ETA2', Line 2845  c     $              'ETA2','ETA2',
2845       $              PFA,PFA,       $              PFA,PFA,
2846       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest))
2847                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2848                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2849                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
2850       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2851                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3535  c     $                 'ETA2','ETA2', Line 2885  c     $                 'ETA2','ETA2',
2885                 endif                 endif
2886    
2887                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2888                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2889                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
2890       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
2891                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3564  c                  dedxmm = dedx(icl)   Line 2915  c                  dedxmm = dedx(icl)  
2915                                
2916                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
2917  *              ----------------------------  *              ----------------------------
2918    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2919                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
2920                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
2921                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
2922                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
2923       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
2924         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2925         $                 ,', norm.dist.= ',distmin
2926         $                 ,', cut ',clinc,' )'
2927                 else                 else
2928                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
2929                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
2930                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
2931       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
2932         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2933         $                 ,', norm.dist.= ', distmin
2934         $                 ,', cut ',clinc,' )'
2935                 endif                 endif
2936    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2937  *              ----------------------------  *              ----------------------------
2938                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
2939                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3608  cccccc 12/08/2006 modified by elena ---> Line 2967  cccccc 12/08/2006 modified by elena --->
2967        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
2968    
2969        include 'commontracker.f'        include 'commontracker.f'
2970          include 'level1.f'
2971        include 'common_momanhough.f'        include 'common_momanhough.f'
2972        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2973        include 'level2.f'        !(1)        include 'level2.f'        !(1)
2974  c      include 'calib.f'  c      include 'calib.f'
2975  c      include 'level1.f'  c      include 'level1.f'
2976    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2977    
2978    
2979        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3688  c               endif Line 3046  c               endif
3046    
3047    
3048    
 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$$$  
   
3049    
3050    
3051  *     ****************************************************  *     ****************************************************
3052    
3053        subroutine init_level2        subroutine init_level2
3054    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3055        include 'commontracker.f'        include 'commontracker.f'
3056          include 'level1.f'
3057        include 'common_momanhough.f'        include 'common_momanhough.f'
3058        include 'level2.f'        include 'level2.f'
3059        include 'level1.f'  c      include 'level1.f'
   
   
3060    
3061        good2 = 0!.false.        do i=1,nviews
3062  c$$$      nev2 = nev1           good2(i)=good1(i)
3063          enddo
3064    
 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*****************************************************  
3065    
3066        NTRK = 0        NTRK = 0
3067        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3068           IMAGE(IT)=0           IMAGE(IT)=0
3069           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3070           do ip=1,nplanes           do ip=1,nplanes
3071              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3072              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3838  c         BdL(IT) = 0. Line 3075  c         BdL(IT) = 0.
3075              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3076              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3077              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3078              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3079              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3080              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3081              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3082           enddo           enddo
# Line 3854  cccccc 17/8/2006 modified by elena Line 3087  cccccc 17/8/2006 modified by elena
3087              enddo                                enddo                  
3088           enddo                             enddo                  
3089        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3090        nclsx=0        nclsx=0
3091        nclsy=0              nclsy=0      
3092        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3093          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3094          xs(1,ip)=0          xs(1,ip)=0
3095          xs(2,ip)=0          xs(2,ip)=0
3096          sgnlxs(ip)=0          sgnlxs(ip)=0
3097          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3098          ys(1,ip)=0          ys(1,ip)=0
3099          ys(2,ip)=0          ys(2,ip)=0
3100          sgnlys(ip)=0          sgnlys(ip)=0
3101        enddo        enddo
 c*******************************************************  
3102        end        end
3103    
3104    
# Line 3887  c*************************************** Line 3113  c***************************************
3113  ************************************************************  ************************************************************
3114    
3115    
3116          subroutine init_hough
3117    
3118          include 'commontracker.f'
3119          include 'level1.f'
3120          include 'common_momanhough.f'
3121          include 'common_hough.f'
3122          include 'level2.f'
3123    
3124          ntrpt_nt=0
3125          ndblt_nt=0
3126          NCLOUDS_XZ_nt=0
3127          NCLOUDS_YZ_nt=0
3128          do idb=1,ndblt_max_nt
3129             db_cloud_nt(idb)=0
3130             alfayz1_nt(idb)=0      
3131             alfayz2_nt(idb)=0      
3132          enddo
3133          do itr=1,ntrpt_max_nt
3134             tr_cloud_nt(itr)=0
3135             alfaxz1_nt(itr)=0      
3136             alfaxz2_nt(itr)=0      
3137             alfaxz3_nt(itr)=0      
3138          enddo
3139          do idb=1,ncloyz_max      
3140            ptcloud_yz_nt(idb)=0    
3141            alfayz1_av_nt(idb)=0    
3142            alfayz2_av_nt(idb)=0    
3143          enddo                    
3144          do itr=1,ncloxz_max      
3145            ptcloud_xz_nt(itr)=0    
3146            alfaxz1_av_nt(itr)=0    
3147            alfaxz2_av_nt(itr)=0    
3148            alfaxz3_av_nt(itr)=0    
3149          enddo                    
3150    
3151          ntrpt=0                  
3152          ndblt=0                  
3153          NCLOUDS_XZ=0              
3154          NCLOUDS_YZ=0              
3155          do idb=1,ndblt_max        
3156            db_cloud(idb)=0        
3157            cpyz1(idb)=0            
3158            cpyz2(idb)=0            
3159            alfayz1(idb)=0          
3160            alfayz2(idb)=0          
3161          enddo                    
3162          do itr=1,ntrpt_max        
3163            tr_cloud(itr)=0        
3164            cpxz1(itr)=0            
3165            cpxz2(itr)=0            
3166            cpxz3(itr)=0            
3167            alfaxz1(itr)=0          
3168            alfaxz2(itr)=0          
3169            alfaxz3(itr)=0          
3170          enddo                    
3171          do idb=1,ncloyz_max      
3172            ptcloud_yz(idb)=0      
3173            alfayz1_av(idb)=0      
3174            alfayz2_av(idb)=0      
3175            do idbb=1,ncouplemaxtot
3176              cpcloud_yz(idb,idbb)=0
3177            enddo                  
3178          enddo                    
3179          do itr=1,ncloxz_max      
3180            ptcloud_xz(itr)=0      
3181            alfaxz1_av(itr)=0      
3182            alfaxz2_av(itr)=0      
3183            alfaxz3_av(itr)=0      
3184            do itrr=1,ncouplemaxtot
3185              cpcloud_xz(itr,itrr)=0
3186            enddo                  
3187          enddo                    
3188          end
3189    ************************************************************
3190    *
3191    *
3192    *
3193    *
3194    *
3195    *
3196    *
3197    ************************************************************
3198    
3199    
3200        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3201    
3202  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3897  c*************************************** Line 3207  c***************************************
3207    
3208            
3209        include 'commontracker.f'        include 'commontracker.f'
3210    c      include 'level1.f'
3211        include 'level1.f'        include 'level1.f'
3212          include 'common_momanhough.f'
3213        include 'level2.f'        include 'level2.f'
3214        include 'common_mini_2.f'        include 'common_mini_2.f'
3215        include 'common_momanhough.f'        real sinth,phi,pig      
       real sinth,phi,pig        !(4)  
3216        pig=acos(-1.)        pig=acos(-1.)
3217    
       good2=1!.true.  
3218        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3219          nstep_nt(ntr)       = nstep
3220    
3221          phi   = al(4)          
3222          sinth = al(3)            
3223          if(sinth.lt.0)then      
3224             sinth = -sinth        
3225             phi = phi + pig      
3226          endif                    
3227          npig = aint(phi/(2*pig))
3228          phi = phi - npig*2*pig  
3229          if(phi.lt.0)            
3230         $     phi = phi + 2*pig  
3231          al(4) = phi              
3232          al(3) = sinth            
3233    
       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)  
 *****************************************************  
3234        do i=1,5        do i=1,5
3235           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3236           do j=1,5           do j=1,5
3237              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3238           enddo           enddo
 c     print*,al_nt(i,ntr)  
3239        enddo        enddo
3240                
3241        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3941  c     print*,al_nt(i,ntr) Line 3251  c     print*,al_nt(i,ntr)
3251           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3252           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3253           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))
 c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3254           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
3255           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)  
3256        
# Line 3958  c            print*,ip,' ',cltrx(ip,ntr) Line 3267  c            print*,ip,' ',cltrx(ip,ntr)
3267           endif                     endif          
3268    
3269        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)  
3270    
3271    
3272        end        end
3273    
3274        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*****************************************************  
3275    
3276  *     -------------------------------------------------------  *     -------------------------------------------------------
3277  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3982  c*************************************** Line 3280  c***************************************
3280  *     -------------------------------------------------------  *     -------------------------------------------------------
3281    
3282        include 'commontracker.f'        include 'commontracker.f'
3283        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
3284        include 'calib.f'        include 'calib.f'
3285          include 'level1.f'
3286        include 'common_momanhough.f'        include 'common_momanhough.f'
3287          include 'level2.f'
3288        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3289    
3290  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
3291        good2=1!.true.  c      good2=1!.true.
3292        nclsx = 0        nclsx = 0
3293        nclsy = 0        nclsy = 0
3294    
3295          do iv = 1,nviews
3296             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3297          enddo
3298    
3299        do icl=1,nclstr1        do icl=1,nclstr1
3300           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
3301              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
# Line 4036  c      print*,icl,cl_used(icl),cl_good(i Line 3339  c      print*,icl,cl_used(icl),cl_good(i
3339        enddo        enddo
3340        end        end
3341    
3342    ***************************************************
3343    *                                                 *
3344    *                                                 *
3345    *                                                 *
3346    *                                                 *
3347    *                                                 *
3348    *                                                 *
3349    **************************************************
3350    
3351          subroutine fill_hough
3352    
3353    *     -------------------------------------------------------
3354    *     This routine fills the  variables related to the hough
3355    *     transform, for the debig n-tuple
3356    *     -------------------------------------------------------
3357    
3358          include 'commontracker.f'
3359          include 'level1.f'
3360          include 'common_momanhough.f'
3361          include 'common_hough.f'
3362          include 'level2.f'
3363    
3364          if(.false.
3365         $     .or.ntrpt.gt.ntrpt_max_nt
3366         $     .or.ndblt.gt.ndblt_max_nt
3367         $     .or.NCLOUDS_XZ.gt.ncloxz_max
3368         $     .or.NCLOUDS_yZ.gt.ncloyz_max
3369         $     )then
3370             ntrpt_nt=0
3371             ndblt_nt=0
3372             NCLOUDS_XZ_nt=0
3373             NCLOUDS_YZ_nt=0
3374          else
3375             ndblt_nt=ndblt
3376             ntrpt_nt=ntrpt
3377             if(ndblt.ne.0)then
3378                do id=1,ndblt
3379                   alfayz1_nt(id)=alfayz1(id) !Y0
3380                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3381                enddo
3382             endif
3383             if(ndblt.ne.0)then
3384                do it=1,ntrpt
3385                   alfaxz1_nt(it)=alfaxz1(it) !X0
3386                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3387                   alfaxz3_nt(it)=alfaxz3(it) !1/r
3388                enddo
3389             endif
3390             nclouds_yz_nt=nclouds_yz
3391             nclouds_xz_nt=nclouds_xz
3392             if(nclouds_yz.ne.0)then
3393                nnn=0
3394                do iyz=1,nclouds_yz
3395                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3396                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3397                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3398                   nnn=nnn+ptcloud_yz(iyz)
3399                enddo
3400                do ipt=1,nnn
3401                   db_cloud_nt(ipt)=db_cloud(ipt)
3402                 enddo
3403             endif
3404             if(nclouds_xz.ne.0)then
3405                nnn=0
3406                do ixz=1,nclouds_xz
3407                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
3408                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
3409                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
3410                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
3411                   nnn=nnn+ptcloud_xz(ixz)              
3412                enddo
3413                do ipt=1,nnn
3414                  tr_cloud_nt(ipt)=tr_cloud(ipt)
3415                 enddo
3416             endif
3417          endif
3418          end
3419          
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23