/[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.1 by mocchiut, Fri May 19 13:15:54 2006 UTC revision 1.14 by pam-fi, Tue Nov 21 14:00:40 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                
       logical DEBUG  
       common/dbg/DEBUG  
   
25  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
26  *     STEP 1  *     STEP 1
27  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 47  Line 44 
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                
       logical DEBUG  
       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    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           rchi2best=1000000000.           rchi2best=1000000000.
237             ndofbest=0             !(1)
238           do i=1,ntracks           do i=1,ntracks
239              if(RCHI2_STORE(i).lt.rchi2best.and.             if(RCHI2_STORE(i).lt.rchi2best.and.
240       $         RCHI2_STORE(i).gt.0)then       $          RCHI2_STORE(i).gt.0)then
241                 ndof=0             !(1)
242                 do ii=1,nplanes    !(1)
243                   ndof=ndof        !(1)
244         $              +int(xgood_store(ii,i)) !(1)
245         $              +int(ygood_store(ii,i)) !(1)
246                 enddo              !(1)
247                 if(ndof.ge.ndofbest)then !(1)
248                 ibest=i                 ibest=i
249                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
250              endif                 ndofbest=ndof    !(1)
251                 endif              !(1)
252               endif
253           enddo           enddo
254           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
255  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
# Line 236  c         iflag=0 Line 281  c         iflag=0
281  *     **********************************************************  *     **********************************************************
282  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
283  *     **********************************************************  *     **********************************************************
284             call guess()
285           do i=1,5           do i=1,5
286              AL(i)=dble(AL_STORE(i,icand))              AL_GUESS(i)=AL(i)
287           enddo           enddo
288    
289             do i=1,5
290                AL(i)=dble(AL_STORE(i,icand))            
291             enddo
292            
293             IDCAND = icand         !fitted track-candidate
294           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
295           jstep=0                !# minimization steps           jstep=0                !# minimization steps
296    
297           call mini_2(jstep,ifail)           iprint=0
298    c         if(DEBUG)iprint=1
299             if(VERBOSE)iprint=1
300             call mini2(jstep,ifail,iprint)
301           if(ifail.ne.0) then           if(ifail.ne.0) then
302              if(DEBUG)then              if(.true.)then
303                 print *,                 print *,
304       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
305       $              ,iev       $              ,iev
306    
307    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
308    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
309    c$$$               print*,'result:   ',(al(i),i=1,5)
310    c$$$               print*,'xgood ',xgood
311    c$$$               print*,'ygood ',ygood
312    c$$$               print*,'----------------------------------------------'
313              endif              endif
314              chi2=-chi2  c            chi2=-chi2
315           endif           endif
316                    
317           if(DEBUG)then           if(DEBUG)then
# Line 310  c         print*,'++++++++++ iimage,fima Line 372  c         print*,'++++++++++ iimage,fima
372  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
373    
374           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
375              if(DEBUG)              if(verbose)
376       $           print*,       $           print*,
377       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
378       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 360  c     $        rchi2best.lt.15..and. Line 422  c     $        rchi2best.lt.15..and.
422        end        end
423    
424    
   
   
 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$$$  
   
425                
426  ************************************************************  ************************************************************
427  ************************************************************  ************************************************************
# Line 606  c                (implemented variable r Line 496  c                (implemented variable r
496  c*****************************************************  c*****************************************************
497                
498        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
499        include 'level1.f'        include 'level1.f'
500          include 'calib.f'
501    c      include 'level1.f'
502        include 'common_align.f'        include 'common_align.f'
503        include 'common_mech.f'        include 'common_mech.f'
504        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
505        include 'common_resxy.f'  c      include 'common_resxy.f'
506    
507        logical DEBUG  c      logical DEBUG
508        common/dbg/DEBUG  c      common/dbg/DEBUG
509    
510        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
511        integer sensor        integer sensor
# Line 666  c      double precision xi_B,yi_B,zi_B Line 557  c      double precision xi_B,yi_B,zi_B
557              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
558           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
559  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
560  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
561  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
562              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
563              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
564              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
565       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
566              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
567           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
568              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
569              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
570              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
571       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
572              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
573           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
574              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
575              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
576              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
577       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
578              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
579           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
580              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
581              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
582              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
583       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 701  c            resxPAM = resxPAM*fbad_cog( Line 592  c            resxPAM = resxPAM*fbad_cog(
592           endif           endif
593    
594        endif        endif
595    c      if(icy.eq.0.and.icx.ne.0)
596    c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'
597                
598  *     -----------------  *     -----------------
599  *     CLUSTER Y  *     CLUSTER Y
# Line 728  c            resxPAM = resxPAM*fbad_cog( Line 621  c            resxPAM = resxPAM*fbad_cog(
621              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
622           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
623  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
624  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
625  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
626              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
627              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
628              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
629              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
630       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
631           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
632              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
633              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
634              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
635       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
636           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
637              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
638              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
639              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
640       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
641           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
642              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
643              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
644  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
645              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 772  C======================================= Line 665  C=======================================
665  c------------------------------------------------------------------------  c------------------------------------------------------------------------
666  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
667  c------------------------------------------------------------------------  c------------------------------------------------------------------------
668             if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
669         $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
670                print*,'xyz_PAM (couple):',
671         $          ' WARNING: false X strip: strip ',stripx
672             endif
673           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
674           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
675           zi = 0.           zi = 0.
# Line 858  C======================================= Line 756  C=======================================
756              nldy = nldx              nldy = nldx
757              viewy = viewx - 1              viewy = viewx - 1
758    
759    c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
760    c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
761                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
762         $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
763                   print*,'xyz_PAM (X-singlet):',
764         $             ' WARNING: false X strip: strip ',stripx
765                endif
766              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
767    
768              xi_A = xi              xi_A = xi
# Line 1136  c$$$         print*,' resolution ',resxP Line 1041  c$$$         print*,' resolution ',resxP
1041  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1042  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1043  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1044                   if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1045         $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1046    c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1047                      print*,'whichsensor: ',
1048         $                ' WARNING: false X strip: strip ',stripx
1049                   endif
1050                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1051                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
1052                 zi = 0.                 zi = 0.
# Line 1263  c     $              ,iv,xvv(iv),yvv(iv) Line 1174  c     $              ,iv,xvv(iv),yvv(iv)
1174  *     it returns the plane number  *     it returns the plane number
1175  *      *    
1176        include 'commontracker.f'        include 'commontracker.f'
1177          include 'level1.f'
1178  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1179        include 'common_momanhough.f'        include 'common_momanhough.f'
1180                
# Line 1300  c      include 'common_analysis.f' Line 1212  c      include 'common_analysis.f'
1212  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1213  *      *    
1214        include 'commontracker.f'        include 'commontracker.f'
1215          include 'level1.f'
1216  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1217        include 'common_momanhough.f'        include 'common_momanhough.f'
1218                
# Line 1328  c      include 'common_analysis.f' Line 1241  c      include 'common_analysis.f'
1241  *     positive if sensor =2  *     positive if sensor =2
1242  *  *
1243        include 'commontracker.f'        include 'commontracker.f'
1244          include 'level1.f'
1245  c      include 'calib.f'  c      include 'calib.f'
1246  c      include 'level1.f'  c      include 'level1.f'
1247  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1357  c      include 'common_analysis.f' Line 1271  c      include 'common_analysis.f'
1271  *************************************************************************  *************************************************************************
1272  *************************************************************************  *************************************************************************
1273  *************************************************************************  *************************************************************************
 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  
1274                
1275    
1276  ***************************************************  ***************************************************
# Line 1639  c$$$      end Line 1285  c$$$      end
1285        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1286    
1287        include 'commontracker.f'        include 'commontracker.f'
1288          include 'level1.f'
1289        include 'common_momanhough.f'        include 'common_momanhough.f'
1290        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1291        include 'calib.f'        include 'calib.f'
1292        include 'level1.f'  c      include 'level1.f'
   
       logical DEBUG  
       common/dbg/DEBUG  
1293    
1294  *     output flag  *     output flag
1295  *     --------------  *     --------------
# Line 1654  c$$$      end Line 1298  c$$$      end
1298  *     --------------  *     --------------
1299        integer iflag        integer iflag
1300    
1301        integer badseed,badcl        integer badseed,badclx,badcly
1302    
1303  *     init variables  *     init variables
1304        ncp_tot=0        ncp_tot=0
# Line 1670  c$$$      end Line 1314  c$$$      end
1314           ncls(ip)=0           ncls(ip)=0
1315        enddo        enddo
1316        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1317           cl_single(icl)=1           cl_single(icl) = 1
1318           cl_good(icl)=0           cl_good(icl)   = 0
1319          enddo
1320          do iv=1,nviews
1321             ncl_view(iv)  = 0
1322             mask_view(iv) = 0      !all included
1323        enddo        enddo
1324                
1325    *     count number of cluster per view
1326          do icl=1,nclstr1
1327             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1328          enddo
1329    *     mask views with too many clusters
1330          do iv=1,nviews
1331             if( ncl_view(iv).gt. nclusterlimit)then
1332                mask_view(iv) = 1
1333                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1334         $           ,nclusterlimit,' on view ', iv,' --> masked!'
1335             endif
1336          enddo
1337    
1338    
1339  *     start association  *     start association
1340        ncouples=0        ncouples=0
1341        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1342           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1343                    
1344  *     ----------------------------------------------------  *     ----------------------------------------------------
1345    *     jump masked views (X VIEW)
1346    *     ----------------------------------------------------
1347             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1348    *     ----------------------------------------------------
1349  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1350    *     ----------------------------------------------------
1351           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
1352              cl_single(icx)=0              cl_single(icx)=0
1353              goto 10              goto 10
1354           endif           endif
1355    *     ----------------------------------------------------
1356  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1357    *     ----------------------------------------------------
1358           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1359           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1360           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1693  c$$$      end Line 1362  c$$$      end
1362           else           else
1363              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1364           endif           endif
1365           badcl=badseed           badclx=badseed
1366           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1367              ibad=1              ibad=1
1368              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1703  c$$$      end Line 1372  c$$$      end
1372       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1373       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1374              endif              endif
1375              badcl=badcl*ibad              badclx=badclx*ibad
1376           enddo           enddo
1377    *     ----------------------------------------------------
1378    *     >>> eliminato il taglio sulle BAD <<<
1379    *     ----------------------------------------------------
1380  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1381  c     cl_single(icx)=0  c     cl_single(icx)=0
1382  c     goto 10  c     goto 10
# Line 1719  c     endif Line 1391  c     endif
1391              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1392                            
1393  *     ----------------------------------------------------  *     ----------------------------------------------------
1394    *     jump masked views (Y VIEW)
1395    *     ----------------------------------------------------
1396                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1397    
1398    *     ----------------------------------------------------
1399  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1400    *     ----------------------------------------------------
1401              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
1402                 cl_single(icy)=0                 cl_single(icy)=0
1403                 goto 20                 goto 20
1404              endif              endif
1405    *     ----------------------------------------------------
1406  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1407    *     ----------------------------------------------------
1408              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1409              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1410              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 1732  c     endif Line 1412  c     endif
1412              else              else
1413                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1414              endif              endif
1415              badcl=badseed              badcly=badseed
1416              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1417                 ibad=1                 ibad=1
1418                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1741  c     endif Line 1421  c     endif
1421       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1422       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1423       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1424                 badcl=badcl*ibad                 badcly=badcly*ibad
1425              enddo              enddo
1426    *     ----------------------------------------------------
1427    *     >>> eliminato il taglio sulle BAD <<<
1428    *     ----------------------------------------------------
1429  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1430  c     cl_single(icy)=0  c     cl_single(icy)=0
1431  c     goto 20  c     goto 20
1432  c     endif  c     endif
1433  *     ----------------------------------------------------  *     ----------------------------------------------------
1434                            
               
1435              cl_good(icy)=1                                cl_good(icy)=1                  
1436              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
1437              nldy=nld(MAXS(icy),VIEW(icy))              nldy=nld(MAXS(icy),VIEW(icy))
# Line 1760  c     endif Line 1442  c     endif
1442  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1443              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1444  *     charge correlation  *     charge correlation
1445                 ddd=(dedx(icy)  *     (modified to be applied only below saturation... obviously)
1446       $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
1447                 ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1448                 cut=chcut*sch(nplx,nldx)       $              .and.
1449                 if(abs(ddd).gt.cut)goto 20 !charge not consistent       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1450                       $              .and.
1451                       $              (badclx.eq.1.and.badcly.eq.1)
1452  *     ------------------> COUPLE <------------------       $              .and.
1453  *     check to do not overflow vector dimentions       $              .true.)then
1454                 if(ncp_plane(nplx).gt.ncouplemax)then  
1455                    if(DEBUG)print*,                    ddd=(dedx(icy)
1456       $                    ' ** warning ** number of identified'//       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1457       $                    ' couples on plane ',nplx,                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1458       $                    ' exceeds vector dimention'//  
1459       $                    ' ( ',ncouplemax,' )'  c                  cut = chcut * sch(nplx,nldx)
1460  c     good2=.false.  
1461  c     goto 880   !fill ntp and go to next event                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1462                    iflag=1       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1463                    return                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1464                 endif                    cut = chcut * (16 + sss/50.)
1465                  
1466                 if(ncp_plane(nplx).eq.ncouplemax)then                    if(abs(ddd).gt.cut)then
1467                    if(DEBUG)print*,                       goto 20    !charge not consistent
1468       $                 '** warning ** number of identified '//                    endif
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
1469                 endif                 endif
1470                                
1471    *     ------------------> COUPLE <------------------
1472                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1473                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1474                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1475                 cl_single(icx)=0                 cl_single(icx)=0
1476                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
 *     ----------------------------------------------  
   
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
1477    
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
       logical DEBUG  
       common/dbg/DEBUG  
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     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  
                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  
                 
1478                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).eq.ncouplemax)then
1479                    if(DEBUG)print*,                    if(verbose)print*,
1480       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1481       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1482       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1483       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1484  c     good2=.false.                    mask_view(nviewx(nplx)) = 2
1485  c     goto 880   !fill ntp and go to next event                                        mask_view(nviewy(nply)) = 2
                   iflag=1  
                   return  
1486                 endif                 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                                
1487  *     ----------------------------------------------  *     ----------------------------------------------
1488    
1489                endif                              
1490    
1491   20         continue   20         continue
1492           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1493                    
# Line 2045  c     goto 880   !fill ntp and go to nex Line 1512  c     goto 880   !fill ntp and go to nex
1512        endif        endif
1513                
1514        do ip=1,6        do ip=1,6
1515           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1516        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  
1517                
1518        return        return
1519        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  
1520                
1521  ***************************************************  ***************************************************
1522  *                                                 *  *                                                 *
# Line 2283  c$$$      end Line 1528  c$$$      end
1528  **************************************************  **************************************************
1529    
1530        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1531    
1532        include 'commontracker.f'        include 'commontracker.f'
1533          include 'level1.f'
1534        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1535        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1536        include 'common_mini_2.f'        include 'common_mini_2.f'
1537        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1538    
       logical DEBUG  
       common/dbg/DEBUG  
1539    
1540  *     output flag  *     output flag
1541  *     --------------  *     --------------
# Line 2329  c      double precision xm3,ym3,zm3 Line 1568  c      double precision xm3,ym3,zm3
1568  *     -----------------------------  *     -----------------------------
1569    
1570    
1571    *     --------------------------------------------
1572    *     put a limit to the maximum number of couples
1573    *     per plane, in order to apply hough transform
1574    *     (couples recovered during track refinement)
1575    *     --------------------------------------------
1576          do ip=1,nplanes
1577             if(ncp_plane(ip).gt.ncouplelimit)then
1578                mask_view(nviewx(ip)) = 8
1579                mask_view(nviewy(ip)) = 8
1580             endif
1581          enddo
1582    
1583    
1584        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1585        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1586                
1587        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1588           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1589                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1590             do is1=1,2             !loop on sensors - COPPIA 1            
1591              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1592                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1593                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
# Line 2345  c               call xyz_PAM(icx1,icy1,i Line 1597  c               call xyz_PAM(icx1,icy1,i
1597                 ym1=yPAM                 ym1=yPAM
1598                 zm1=zPAM                                   zm1=zPAM                  
1599  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1600    
1601                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1602                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1603         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1604                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1605                                            
1606                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2364  c     $                       (icx2,icy2 Line 1619  c     $                       (icx2,icy2
1619  *     (2 couples needed)  *     (2 couples needed)
1620  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1621                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1622                             if(DEBUG)print*,                             if(verbose)print*,
1623       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1624       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1625       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1626  c                           good2=.false.  c                           good2=.false.
1627  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1628                               do iv=1,12
1629                                  mask_view(iv) = 3
1630                               enddo
1631                             iflag=1                             iflag=1
1632                             return                             return
1633                          endif                          endif
# Line 2403  c$$$ Line 1661  c$$$
1661  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1662    
1663    
1664                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1665    
1666                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1667                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1668         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1669                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1670                                                                
1671                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2434  c     $                                 Line 1695  c     $                                
1695  *     (3 couples needed)  *     (3 couples needed)
1696  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1697                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
1698                                      if(DEBUG)print*,                                      if(verbose)print*,
1699       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
1700       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
1701       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
1702  c                                    good2=.false.  c                                    good2=.false.
1703  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
1704                                        do iv=1,nviews
1705                                           mask_view(iv) = 4
1706                                        enddo
1707                                      iflag=1                                      iflag=1
1708                                      return                                      return
1709                                   endif                                   endif
# Line 2478  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1742  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1742                                endif                                endif
1743                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1744                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1745     30                     continue
1746                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1747   30                  continue   31                  continue
1748                                            
1749   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1750                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1751     20            continue
1752              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1753                            
1754           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1755        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1756     10   continue
1757        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1758                
1759        if(DEBUG)then        if(DEBUG)then
# Line 2514  c     goto 880               !ntp fill Line 1781  c     goto 880               !ntp fill
1781        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
1782    
1783        include 'commontracker.f'        include 'commontracker.f'
1784          include 'level1.f'
1785        include 'common_momanhough.f'        include 'common_momanhough.f'
1786        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1787    
       logical DEBUG  
       common/dbg/DEBUG  
1788    
1789  *     output flag  *     output flag
1790  *     --------------  *     --------------
# Line 2550  c     goto 880               !ntp fill Line 1816  c     goto 880               !ntp fill
1816        distance=0        distance=0
1817        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
1818        npt_tot=0        npt_tot=0
1819          nloop=0                  
1820     90   continue                  
1821        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
1822           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
1823                            
# Line 2653  c     print*,'*   idbref,idb2 ',idbref,i Line 1921  c     print*,'*   idbref,idb2 ',idbref,i
1921              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
1922           enddo           enddo
1923  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
1924           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
1925           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
1926           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
1927                    
# Line 2661  c     print*,'>>>> ',ncpused,npt,nplused Line 1929  c     print*,'>>>> ',ncpused,npt,nplused
1929  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
1930    
1931           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
1932              if(DEBUG)print*,              if(verbose)print*,
1933       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
1934       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
1935       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
1936  c               good2=.false.  c               good2=.false.
1937  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
1938                do iv=1,nviews
1939                   mask_view(iv) = 5
1940                enddo
1941              iflag=1              iflag=1
1942              return              return
1943           endif           endif
# Line 2704  c$$$     $           ,(db_cloud(iii),iii Line 1975  c$$$     $           ,(db_cloud(iii),iii
1975        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
1976                
1977                
1978          if(nloop.lt.nstepy)then      
1979            cutdistyz = cutdistyz+cutystep
1980            nloop     = nloop+1          
1981            goto 90                
1982          endif                    
1983          
1984        if(DEBUG)then        if(DEBUG)then
1985           print*,'---------------------- '           print*,'---------------------- '
1986           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2730  c$$$     $           ,(db_cloud(iii),iii Line 2007  c$$$     $           ,(db_cloud(iii),iii
2007        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2008    
2009        include 'commontracker.f'        include 'commontracker.f'
2010          include 'level1.f'
2011        include 'common_momanhough.f'        include 'common_momanhough.f'
2012        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2013    
       logical DEBUG  
       common/dbg/DEBUG  
2014    
2015  *     output flag  *     output flag
2016  *     --------------  *     --------------
# Line 2765  c$$$     $           ,(db_cloud(iii),iii Line 2041  c$$$     $           ,(db_cloud(iii),iii
2041        distance=0        distance=0
2042        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2043        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2044          nloop=0                  
2045     91   continue                  
2046        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2047           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
2048  c     print*,'--------------'  c     print*,'--------------'
# Line 2866  c     print*,'check cp_used' Line 2144  c     print*,'check cp_used'
2144           do ip=1,nplanes           do ip=1,nplanes
2145              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2146           enddo           enddo
2147           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2148           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2149           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2150                    
2151  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2152  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2153           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2154              if(DEBUG)print*,              if(verbose)print*,
2155       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2156       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2157       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2158  c     good2=.false.  c     good2=.false.
2159  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2160                do iv=1,nviews
2161                   mask_view(iv) = 6
2162                enddo
2163              iflag=1              iflag=1
2164              return              return
2165           endif           endif
# Line 2914  c$$$     $           ,(tr_cloud(iii),iii Line 2195  c$$$     $           ,(tr_cloud(iii),iii
2195  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2196  22288    continue  22288    continue
2197        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2198          
2199           if(nloop.lt.nstepx)then      
2200             cutdistxz=cutdistxz+cutxstep
2201             nloop=nloop+1          
2202             goto 91                
2203           endif                    
2204          
2205        if(DEBUG)then        if(DEBUG)then
2206           print*,'---------------------- '           print*,'---------------------- '
2207           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2941  c     02/02/2006 modified by Elena Vannu Line 2228  c     02/02/2006 modified by Elena Vannu
2228  c*****************************************************  c*****************************************************
2229    
2230        include 'commontracker.f'        include 'commontracker.f'
2231          include 'level1.f'
2232        include 'common_momanhough.f'        include 'common_momanhough.f'
2233        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2234        include 'common_mini_2.f'        include 'common_mini_2.f'
2235        include 'common_mech.f'        include 'common_mech.f'
2236        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2237    
       logical DEBUG  
       common/dbg/DEBUG  
2238    
2239  *     output flag  *     output flag
2240  *     --------------  *     --------------
# Line 2964  c*************************************** Line 2250  c***************************************
2250  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2251  *     list of matching couples in the combination  *     list of matching couples in the combination
2252  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2253        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2254        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2255  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2256        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2257  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2258  *     variables for track fitting  *     variables for track fitting
2259        double precision AL_INI(5)        double precision AL_INI(5)
2260        double precision tath  c      double precision tath
2261  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2262  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2263    
# Line 3064  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2350  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2350  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2351                            
2352  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2353              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2354              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2355              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2356       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2357              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2358              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2359              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2360                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2361  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2362              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2363                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2364  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2365  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2366                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2367  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2368  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2369                c$$$            ELSE
2370    c$$$               AL_INI(4) = acos(-1.)/2
2371    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2372    c$$$            ENDIF
2373    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2374    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2375    c$$$            
2376    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2377    c$$$            
2378    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2379                            
2380              if(DEBUG)then              if(DEBUG)then
2381                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2382                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3148  c     $                                 Line 2444  c     $                                
2444  *     **********************************************************  *     **********************************************************
2445  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2446  *     **********************************************************  *     **********************************************************
2447    cccc  scommentare se si usa al_ini della nuvola
2448    c$$$                              do i=1,5
2449    c$$$                                 AL(i)=AL_INI(i)
2450    c$$$                              enddo
2451                                  call guess()
2452                                do i=1,5                                do i=1,5
2453                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2454                                enddo                                enddo
2455                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2456                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2457                                call mini_2(jstep,ifail)                                iprint=0
2458    c                              if(DEBUG)iprint=1
2459                                  if(DEBUG)iprint=1
2460                                  call mini2(jstep,ifail,iprint)
2461                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2462                                   if(DEBUG)then                                   if(DEBUG)then
2463                                      print *,                                      print *,
2464       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2465       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2466                                        print*,'initial guess: '
2467    
2468                                        print*,'AL_INI(1) = ',AL_INI(1)
2469                                        print*,'AL_INI(2) = ',AL_INI(2)
2470                                        print*,'AL_INI(3) = ',AL_INI(3)
2471                                        print*,'AL_INI(4) = ',AL_INI(4)
2472                                        print*,'AL_INI(5) = ',AL_INI(5)
2473                                   endif                                   endif
2474                                   chi2=-chi2  c                                 chi2=-chi2
2475                                endif                                endif
2476  *     **********************************************************  *     **********************************************************
2477  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3173  c     $                                 Line 2484  c     $                                
2484  *     --------------------------  *     --------------------------
2485                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2486                                                                    
2487                                   if(DEBUG)print*,                                   if(verbose)print*,
2488       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2489       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2490       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2491  c                                 good2=.false.  c                                 good2=.false.
2492  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2493                                     do iv=1,nviews
2494                                        mask_view(iv) = 7
2495                                     enddo
2496                                   iflag=1                                   iflag=1
2497                                   return                                   return
2498                                endif                                endif
# Line 3273  c$$$  rchi2=chi2/dble(ndof) Line 2587  c$$$  rchi2=chi2/dble(ndof)
2587  c******************************************************  c******************************************************
2588  cccccc 06/10/2005 modified by elena vannuccini ---> (1)  cccccc 06/10/2005 modified by elena vannuccini ---> (1)
2589  cccccc 31/01/2006 modified by elena vannuccini ---> (2)  cccccc 31/01/2006 modified by elena vannuccini ---> (2)
2590    cccccc 12/08/2006 modified by elena vannucicni ---> (3)
2591  c******************************************************  c******************************************************
2592    
2593        include 'commontracker.f'        include 'commontracker.f'
2594          include 'level1.f'
2595        include 'common_momanhough.f'        include 'common_momanhough.f'
2596        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2597        include 'common_mini_2.f'        include 'common_mini_2.f'
2598        include 'common_mech.f'        include 'common_mech.f'
2599        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2600        include 'level1.f'  c      include 'level1.f'
2601        include 'calib.f'        include 'calib.f'
2602    
       logical DEBUG  
       common/dbg/DEBUG  
2603    
2604  *     flag to chose PFA  *     flag to chose PFA
2605        character*10 PFA        character*10 PFA
# Line 3299  c*************************************** Line 2613  c***************************************
2613        call track_init        call track_init
2614        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2615    
2616    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2617  *     -------------------------------------------------  *     -------------------------------------------------
2618  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2619  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2620  *     using improved PFAs  *     using improved PFAs
2621  *     -------------------------------------------------  *     -------------------------------------------------
2622    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2623           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2624       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2625                            
# Line 3338  c            dedxtrk(nplanes-ip+1) = (de Line 2654  c            dedxtrk(nplanes-ip+1) = (de
2654              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)
2655              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)
2656                            
2657    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2658  *     -------------------------------------------------  *     -------------------------------------------------
2659  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
2660  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
2661  *     -------------------------------------------------  *     -------------------------------------------------
2662    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2663           else                             else                  
2664                                
2665              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3383  c            if(DEBUG)print*,'>>>> try t Line 2701  c            if(DEBUG)print*,'>>>> try t
2701                 icx=clx(ip,icp)                 icx=clx(ip,icp)
2702                 icy=cly(ip,icp)                 icy=cly(ip,icp)
2703                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
2704       $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
2705       $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
2706         $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)
2707         $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
2708       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
2709  *            *          
2710                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3394  c     $              'ETA2','ETA2', Line 2714  c     $              'ETA2','ETA2',
2714       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest))
2715                                
2716                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2717                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2718                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2719                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2720       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3456  c            if(DEBUG)print*,'>>>> try t Line 2777  c            if(DEBUG)print*,'>>>> try t
2777                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
2778  *                                                !jump to the next couple  *                                                !jump to the next couple
2779  *----- try cluster x -----------------------------------------------  *----- try cluster x -----------------------------------------------
2780                 if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
2781                   if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
2782  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
2783                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
2784  c     $              'ETA2','ETA2',  c     $              'ETA2','ETA2',
2785       $              PFA,PFA,       $              PFA,PFA,
2786       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
2787                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2788  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2789                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
2790       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2791                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3483  c                  dedxmm = dedx(icx) !( Line 2805  c                  dedxmm = dedx(icx) !(
2805                 endif                                   endif                  
2806  11881          continue  11881          continue
2807  *----- try cluster y -----------------------------------------------  *----- try cluster y -----------------------------------------------
2808                 if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
2809                   if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
2810  *                                              !jump to the next couple  *                                              !jump to the next couple
2811                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
2812  c     $              'ETA2','ETA2',  c     $              'ETA2','ETA2',
2813       $              PFA,PFA,       $              PFA,PFA,
2814       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest))
2815                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2816                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2817                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
2818       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2819                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3512  c                 dedxmm = dedx(icy)  !( Line 2836  c                 dedxmm = dedx(icy)  !(
2836  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------    
2837              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
2838                 icl=cls(ip,ic)                 icl=cls(ip,ic)
2839                 if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
2840                   if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
2841       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
2842       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
2843                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
# Line 3528  c     $                 'ETA2','ETA2', Line 2853  c     $                 'ETA2','ETA2',
2853                 endif                 endif
2854    
2855                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2856                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2857                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
2858       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
2859                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3557  c                  dedxmm = dedx(icl)   Line 2883  c                  dedxmm = dedx(icl)  
2883                                
2884                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
2885  *              ----------------------------  *              ----------------------------
2886    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2887                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
2888                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
2889                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
2890                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
2891       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
2892         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2893         $                 ,', norm.dist.= ',distmin
2894         $                 ,', cut ',clinc,' )'
2895                 else                 else
2896                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
2897                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
2898                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
2899       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
2900         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2901         $                 ,', norm.dist.= ', distmin
2902         $                 ,', cut ',clinc,' )'
2903                 endif                 endif
2904    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2905  *              ----------------------------  *              ----------------------------
2906                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
2907                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3596  c              dedxtrk(nplanes-ip+1) = d Line 2930  c              dedxtrk(nplanes-ip+1) = d
2930  *                                                 *  *                                                 *
2931  *                                                 *  *                                                 *
2932  **************************************************  **************************************************
2933    cccccc 12/08/2006 modified by elena ---> (1)
2934    *
2935        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
2936    
2937        include 'commontracker.f'        include 'commontracker.f'
2938          include 'level1.f'
2939        include 'common_momanhough.f'        include 'common_momanhough.f'
2940        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2941          include 'level2.f'        !(1)
2942  c      include 'calib.f'  c      include 'calib.f'
2943  c      include 'level1.f'  c      include 'level1.f'
2944    
       logical DEBUG  
       common/dbg/DEBUG  
2945    
2946    
2947        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3617  c      include 'level1.f' Line 2952  c      include 'level1.f'
2952              if(id.ne.0)then              if(id.ne.0)then
2953                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
2954                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
2955                 cl_used(iclx)=1  !tag used clusters  c               cl_used(iclx)=1  !tag used clusters
2956                 cl_used(icly)=1  !tag used clusters  c               cl_used(icly)=1  !tag used clusters
2957                   cl_used(iclx)=ntrk  !tag used clusters !(1)
2958                   cl_used(icly)=ntrk  !tag used clusters !(1)
2959              elseif(icl.ne.0)then              elseif(icl.ne.0)then
2960                 cl_used(icl)=1   !tag used clusters  c               cl_used(icl)=1   !tag used clusters
2961                   cl_used(icl)=ntrk   !tag used clusters !1)
2962              endif              endif
2963                            
2964  c               if(DEBUG)then  c               if(DEBUG)then
# Line 3676  c               endif Line 3014  c               endif
3014    
3015    
3016    
 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$$$  
   
3017    
3018    
3019  *     ****************************************************  *     ****************************************************
3020    
3021        subroutine init_level2        subroutine init_level2
3022    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3023        include 'commontracker.f'        include 'commontracker.f'
3024          include 'level1.f'
3025        include 'common_momanhough.f'        include 'common_momanhough.f'
3026        include 'level2.f'        include 'level2.f'
3027        include 'level1.f'  c      include 'level1.f'
   
3028    
3029          do i=1,nviews
3030             good2(i)=good1(i)
3031          enddo
3032    
       good2 = 0!.false.  
 c$$$      nev2 = nev1  
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
3033    
3034        NTRK = 0        NTRK = 0
3035        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3036           IMAGE(IT)=0           IMAGE(IT)=0
3037           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
          BdL(IT) = 0.  
3038           do ip=1,nplanes           do ip=1,nplanes
3039              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3040              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3826  c*************************************** Line 3043  c***************************************
3043              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3044              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3045              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3046              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3047              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3048  c******************************************************              CLTRX(IP,IT) = 0
3049                CLTRY(IP,IT) = 0
3050           enddo           enddo
3051           do ipa=1,5           do ipa=1,5
3052              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3839  c*************************************** Line 3055  c***************************************
3055              enddo                                enddo                  
3056           enddo                             enddo                  
3057        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3058        nclsx=0        nclsx=0
3059        nclsy=0              nclsy=0      
3060        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3061          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3062          xs(1,ip)=0          xs(1,ip)=0
3063          xs(2,ip)=0          xs(2,ip)=0
3064          sgnlxs(ip)=0          sgnlxs(ip)=0
3065          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3066          ys(1,ip)=0          ys(1,ip)=0
3067          ys(2,ip)=0          ys(2,ip)=0
3068          sgnlys(ip)=0          sgnlys(ip)=0
3069        enddo        enddo
 c*******************************************************  
3070        end        end
3071    
3072    
# Line 3872  c*************************************** Line 3081  c***************************************
3081  ************************************************************  ************************************************************
3082    
3083    
3084          subroutine init_hough
3085    
3086          include 'commontracker.f'
3087          include 'level1.f'
3088          include 'common_momanhough.f'
3089          include 'common_hough.f'
3090          include 'level2.f'
3091    
3092          ntrpt_nt=0
3093          ndblt_nt=0
3094          NCLOUDS_XZ_nt=0
3095          NCLOUDS_YZ_nt=0
3096          do idb=1,ndblt_max_nt
3097             db_cloud_nt(idb)=0
3098             alfayz1_nt(idb)=0      
3099             alfayz2_nt(idb)=0      
3100          enddo
3101          do itr=1,ntrpt_max_nt
3102             tr_cloud_nt(itr)=0
3103             alfaxz1_nt(itr)=0      
3104             alfaxz2_nt(itr)=0      
3105             alfaxz3_nt(itr)=0      
3106          enddo
3107          do idb=1,ncloyz_max      
3108            ptcloud_yz_nt(idb)=0    
3109            alfayz1_av_nt(idb)=0    
3110            alfayz2_av_nt(idb)=0    
3111          enddo                    
3112          do itr=1,ncloxz_max      
3113            ptcloud_xz_nt(itr)=0    
3114            alfaxz1_av_nt(itr)=0    
3115            alfaxz2_av_nt(itr)=0    
3116            alfaxz3_av_nt(itr)=0    
3117          enddo                    
3118    
3119          ntrpt=0                  
3120          ndblt=0                  
3121          NCLOUDS_XZ=0              
3122          NCLOUDS_YZ=0              
3123          do idb=1,ndblt_max        
3124            db_cloud(idb)=0        
3125            cpyz1(idb)=0            
3126            cpyz2(idb)=0            
3127            alfayz1(idb)=0          
3128            alfayz2(idb)=0          
3129          enddo                    
3130          do itr=1,ntrpt_max        
3131            tr_cloud(itr)=0        
3132            cpxz1(itr)=0            
3133            cpxz2(itr)=0            
3134            cpxz3(itr)=0            
3135            alfaxz1(itr)=0          
3136            alfaxz2(itr)=0          
3137            alfaxz3(itr)=0          
3138          enddo                    
3139          do idb=1,ncloyz_max      
3140            ptcloud_yz(idb)=0      
3141            alfayz1_av(idb)=0      
3142            alfayz2_av(idb)=0      
3143            do idbb=1,ncouplemaxtot
3144              cpcloud_yz(idb,idbb)=0
3145            enddo                  
3146          enddo                    
3147          do itr=1,ncloxz_max      
3148            ptcloud_xz(itr)=0      
3149            alfaxz1_av(itr)=0      
3150            alfaxz2_av(itr)=0      
3151            alfaxz3_av(itr)=0      
3152            do itrr=1,ncouplemaxtot
3153              cpcloud_xz(itr,itrr)=0
3154            enddo                  
3155          enddo                    
3156          end
3157    ************************************************************
3158    *
3159    *
3160    *
3161    *
3162    *
3163    *
3164    *
3165    ************************************************************
3166    
3167    
3168        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3169    
3170  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3882  c*************************************** Line 3175  c***************************************
3175    
3176            
3177        include 'commontracker.f'        include 'commontracker.f'
3178    c      include 'level1.f'
3179          include 'level1.f'
3180          include 'common_momanhough.f'
3181        include 'level2.f'        include 'level2.f'
3182        include 'common_mini_2.f'        include 'common_mini_2.f'
3183        real sinth,phi,pig        !(4)        real sinth,phi,pig      
3184        pig=acos(-1.)        pig=acos(-1.)
3185    
       good2=1!.true.  
3186        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3187          nstep_nt(ntr)       = nstep
3188    
3189          phi   = al(4)          
3190          sinth = al(3)            
3191          if(sinth.lt.0)then      
3192             sinth = -sinth        
3193             phi = phi + pig      
3194          endif                    
3195          npig = aint(phi/(2*pig))
3196          phi = phi - npig*2*pig  
3197          if(phi.lt.0)            
3198         $     phi = phi + 2*pig  
3199          al(4) = phi              
3200          al(3) = sinth            
3201    
       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)  
 *****************************************************  
3202        do i=1,5        do i=1,5
3203           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3204           do j=1,5           do j=1,5
3205              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3206           enddo           enddo
 c     print*,al_nt(i,ntr)  
3207        enddo        enddo
3208                
3209        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3924  c     print*,al_nt(i,ntr) Line 3219  c     print*,al_nt(i,ntr)
3219           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3220           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3221           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))
 c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3222           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
3223           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)                     dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)  
3224      
3225             id  = CP_STORE(ip,IDCAND)
3226             icl = CLS_STORE(ip,IDCAND)
3227             if(id.ne.0)then
3228                cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3229                cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3230    c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3231             elseif(icl.ne.0)then
3232                if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl
3233                if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl
3234    c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3235             endif          
3236    
3237        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)  
3238    
3239    
3240        end        end
3241    
3242        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*****************************************************  
3243    
3244  *     -------------------------------------------------------  *     -------------------------------------------------------
3245  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3952  c*************************************** Line 3248  c***************************************
3248  *     -------------------------------------------------------  *     -------------------------------------------------------
3249    
3250        include 'commontracker.f'        include 'commontracker.f'
3251        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
3252        include 'calib.f'        include 'calib.f'
3253          include 'level1.f'
3254        include 'common_momanhough.f'        include 'common_momanhough.f'
3255          include 'level2.f'
3256        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3257    
3258  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
3259        good2=1!.true.  c      good2=1!.true.
3260        nclsx = 0        nclsx = 0
3261        nclsy = 0        nclsy = 0
3262    
3263          do iv = 1,nviews
3264             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3265          enddo
3266    
3267        do icl=1,nclstr1        do icl=1,nclstr1
3268           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
3269              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
# Line 3970  c*************************************** Line 3271  c***************************************
3271                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3272                 planex(nclsx) = ip                 planex(nclsx) = ip
3273                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3274                   clsx(nclsx)   = icl
3275                 do is=1,2                 do is=1,2
3276  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3277                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
# Line 3984  c$$$               print*,'xs(2,nclsx)   Line 3286  c$$$               print*,'xs(2,nclsx)  
3286                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3287                 planey(nclsy) = ip                 planey(nclsy) = ip
3288                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3289                   clsy(nclsy)   = icl
3290                 do is=1,2                 do is=1,2
3291  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3292                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
# Line 3997  c$$$               print*,'ys(2,nclsy)   Line 3300  c$$$               print*,'ys(2,nclsy)  
3300              endif              endif
3301           endif           endif
3302  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
3303    
3304    ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3305             whichtrack(icl) = cl_used(icl)
3306    
3307        enddo        enddo
3308        end        end
3309    
3310    ***************************************************
3311    *                                                 *
3312    *                                                 *
3313    *                                                 *
3314    *                                                 *
3315    *                                                 *
3316    *                                                 *
3317    **************************************************
3318    
3319          subroutine fill_hough
3320    
3321    *     -------------------------------------------------------
3322    *     This routine fills the  variables related to the hough
3323    *     transform, for the debig n-tuple
3324    *     -------------------------------------------------------
3325    
3326          include 'commontracker.f'
3327          include 'level1.f'
3328          include 'common_momanhough.f'
3329          include 'common_hough.f'
3330          include 'level2.f'
3331    
3332          if(.false.
3333         $     .or.ntrpt.gt.ntrpt_max_nt
3334         $     .or.ndblt.gt.ndblt_max_nt
3335         $     .or.NCLOUDS_XZ.gt.ncloxz_max
3336         $     .or.NCLOUDS_yZ.gt.ncloyz_max
3337         $     )then
3338             ntrpt_nt=0
3339             ndblt_nt=0
3340             NCLOUDS_XZ_nt=0
3341             NCLOUDS_YZ_nt=0
3342          else
3343             ndblt_nt=ndblt
3344             ntrpt_nt=ntrpt
3345             if(ndblt.ne.0)then
3346                do id=1,ndblt
3347                   alfayz1_nt(id)=alfayz1(id) !Y0
3348                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3349                enddo
3350             endif
3351             if(ndblt.ne.0)then
3352                do it=1,ntrpt
3353                   alfaxz1_nt(it)=alfaxz1(it) !X0
3354                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3355                   alfaxz3_nt(it)=alfaxz3(it) !1/r
3356                enddo
3357             endif
3358             nclouds_yz_nt=nclouds_yz
3359             nclouds_xz_nt=nclouds_xz
3360             if(nclouds_yz.ne.0)then
3361                nnn=0
3362                do iyz=1,nclouds_yz
3363                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3364                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3365                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3366                   nnn=nnn+ptcloud_yz(iyz)
3367                enddo
3368                do ipt=1,nnn
3369                   db_cloud_nt(ipt)=db_cloud(ipt)
3370                 enddo
3371             endif
3372             if(nclouds_xz.ne.0)then
3373                nnn=0
3374                do ixz=1,nclouds_xz
3375                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
3376                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
3377                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
3378                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
3379                   nnn=nnn+ptcloud_xz(ixz)              
3380                enddo
3381                do ipt=1,nnn
3382                  tr_cloud_nt(ipt)=tr_cloud(ipt)
3383                 enddo
3384             endif
3385          endif
3386          end
3387          
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.23