/[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.6 by pam-fi, Mon Oct 2 13:12:48 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                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
25  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
26  *     STEP 1  *     STEP 1
27  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 47  c      common/dbg/DEBUG Line 44  c      common/dbg/DEBUG
44  c      iflag=0  c      iflag=0
45        call cl_to_couples(iflag)        call cl_to_couples(iflag)
46        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
47           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
48        endif        endif
49                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     selezione di tracce pulite per diagnostica  
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 c$$$         if(DEBUG)then  
 c$$$            do ip=1,nplanes  
 c$$$               if(ncp_plane(ip).ne.1)good2=.false.  
 c$$$            enddo  
 c$$$c            if(good2.eq.0)goto 100!next event  
 c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event  
 c$$$         endif  
           
   
   
50  *-----------------------------------------------------  *-----------------------------------------------------
51  *-----------------------------------------------------  *-----------------------------------------------------
52  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 94  c$$$         endif Line 78  c$$$         endif
78  c      iflag=0  c      iflag=0
79        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
80        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
81           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
82        endif        endif
83                
84                
# Line 123  c      iflag=0 Line 107  c      iflag=0
107  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    *     count number of hit planes
111          planehit=0                
112          do np=1,nplanes          
113            if(ncp_plane(np).ne.0)then  
114              planehit=planehit+1  
115            endif                  
116          enddo                    
117          if(planehit.lt.3) goto 880 ! exit              
118          
119          nptxz_min=x_min_start              
120          nplxz_min=x_min_start              
121                
122          nptyz_min=y_min_start              
123          nplyz_min=y_min_start              
124    
125  c      iflag=0        cutdistyz=cutystart      
126          cutdistxz=cutxstart      
127    
128     878  continue
129        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
130        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
131           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
132        endif        endif
133  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
134            if(cutdistyz.lt.maxcuty/2)then
135              cutdistyz=cutdistyz+cutystep
136            else
137              cutdistyz=cutdistyz+(3*cutystep)
138            endif
139            goto 878                
140          endif                    
141    
142          if(planehit.eq.3) goto 881          
143          
144     879  continue  
145        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
146        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
147           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
148        endif        endif
149                                    
150          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
151            cutdistxz=cutdistxz+cutxstep
152            goto 879                
153          endif                    
154    
155        
156     881  continue  
157    *     if there is at least three planes on the Y view decreases cuts on X view
158          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
159         $     nplxz_min.ne.y_min_start)then
160            nptxz_min=x_min_step    
161            nplxz_min=x_min_start-x_min_step              
162            goto 879                
163          endif                    
164            
165   880  return   880  return
166        end        end
167    
# Line 144  c      iflag=0 Line 171  c      iflag=0
171        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
172    
173        include 'commontracker.f'        include 'commontracker.f'
174          include 'level1.f'
175        include 'common_momanhough.f'        include 'common_momanhough.f'
176        include 'common_mech.f'        include 'common_mech.f'
177        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
178        include 'common_mini_2.f'        include 'common_mini_2.f'
179        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
180        include 'level2.f'        include 'level2.f'
181    
182        include 'momanhough_init.f'  c      include 'momanhough_init.f'
183                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
184        logical FIMAGE            !        logical FIMAGE            !
185          real*8 AL_GUESS(5)
186    
187  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
188  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 198  c         iflag=0 Line 223  c         iflag=0
223           ibest=0                !best track among candidates           ibest=0                !best track among candidates
224           iimage=0               !track image           iimage=0               !track image
225  *     ------------- select the best track -------------  *     ------------- select the best track -------------
226    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
286                AL_GUESS(i)=AL(i)
287             enddo
288    
289           do i=1,5           do i=1,5
290              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
291           enddo           enddo
292            
293           IDCAND = icand         !fitted track-candidate           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 361  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 607  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  c      logical DEBUG  c      logical DEBUG
508  c      common/dbg/DEBUG  c      common/dbg/DEBUG
# Line 667  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 731  c     $     print*,PFAx,icx,angx,stripx, Line 621  c     $     print*,PFAx,icx,angx,stripx,
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 1284  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 1321  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 1349  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 1378  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 1660  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'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1293    
1294  *     output flag  *     output flag
1295  *     --------------  *     --------------
# Line 1691  c      common/dbg/DEBUG Line 1314  c      common/dbg/DEBUG
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
# Line 1746  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
# Line 1794  c     endif Line 1444  c     endif
1444  *     charge correlation  *     charge correlation
1445  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1446    
 *     -------------------------------------------------------------  
 *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<  
 *     -------------------------------------------------------------  
1447                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1448       $              .and.       $              .and.
1449       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
# Line 1822  c                  cut = chcut * sch(npl Line 1469  c                  cut = chcut * sch(npl
1469                 endif                 endif
1470                                
1471  *     ------------------> COUPLE <------------------  *     ------------------> COUPLE <------------------
 *     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  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
                 
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                                
 *     ----------------------------------------------  
1477    
  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(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     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  
 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  
                 
1478                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).eq.ncouplemax)then
1479                    if(verbose)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 2099  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(verbose)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
   
1520                
1521  ***************************************************  ***************************************************
1522  *                                                 *  *                                                 *
# Line 2128  c     goto 880       !fill ntp and go to Line 1528  c     goto 880       !fill ntp and go to
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    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1539    
1540  *     output flag  *     output flag
1541  *     --------------  *     --------------
# Line 2174  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 2190  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 2215  c     $                       (icx2,icy2 Line 1625  c     $                       (icx2,icy2
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 2248  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 2285  c     $                                 Line 1701  c     $                                
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 2323  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 2359  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    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1788    
1789  *     output flag  *     output flag
1790  *     --------------  *     --------------
# Line 2395  c      common/dbg/DEBUG Line 1816  c      common/dbg/DEBUG
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 2498  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 2512  c     print*,'>>>> ',ncpused,npt,nplused Line 1935  c     print*,'>>>> ',ncpused,npt,nplused
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 2549  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 2575  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    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2014    
2015  *     output flag  *     output flag
2016  *     --------------  *     --------------
# Line 2610  c      common/dbg/DEBUG Line 2041  c      common/dbg/DEBUG
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 2711  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                    
# Line 2724  c     print*,'check cp_used' Line 2157  c     print*,'check cp_used'
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 2759  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 2786  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    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2238    
2239  *     output flag  *     output flag
2240  *     --------------  *     --------------
# Line 2809  c      common/dbg/DEBUG Line 2250  c      common/dbg/DEBUG
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 2909  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 2993  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 3024  c     $                                 Line 2490  c     $                                
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 3122  cccccc 12/08/2006 modified by elena vann Line 2591  cccccc 12/08/2006 modified by elena vann
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    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2603    
2604  *     flag to chose PFA  *     flag to chose PFA
2605        character*10 PFA        character*10 PFA
# Line 3145  c      common/dbg/DEBUG Line 2613  c      common/dbg/DEBUG
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 3184  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 3242  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 3312  c     $              'ETA2','ETA2', Line 2785  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 3340  c     $              'ETA2','ETA2', Line 2813  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 3379  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 3408  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 3452  cccccc 12/08/2006 modified by elena ---> Line 2935  cccccc 12/08/2006 modified by elena --->
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)        include 'level2.f'        !(1)
2942  c      include 'calib.f'  c      include 'calib.f'
2943  c      include 'level1.f'  c      include 'level1.f'
2944    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2945    
2946    
2947        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3532  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        do i=1,nviews
3030           good2(i)=good1(i)           good2(i)=good1(i)
3031        enddo        enddo
3032    
 c      good2 = 0!.false.  
 c$$$      nev2 = nev1  
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
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.
 c         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 3684  c         BdL(IT) = 0. Line 3043  c         BdL(IT) = 0.
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
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3048              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3049              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3050           enddo           enddo
# Line 3700  cccccc 17/8/2006 modified by elena Line 3055  cccccc 17/8/2006 modified by elena
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 3733  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 3743  c*************************************** Line 3175  c***************************************
3175    
3176            
3177        include 'commontracker.f'        include 'commontracker.f'
3178    c      include 'level1.f'
3179        include 'level1.f'        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        include 'common_momanhough.f'        real sinth,phi,pig      
       real sinth,phi,pig        !(4)  
3184        pig=acos(-1.)        pig=acos(-1.)
3185    
 c      good2=1!.true.  
3186        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3187        nstep_nt(ntr)       = 0!nstep        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 3788  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        
# Line 3805  c            print*,ip,' ',cltrx(ip,ntr) Line 3235  c            print*,ip,' ',cltrx(ip,ntr)
3235           endif                     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 3829  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
# Line 3840  c      good2=1!.true. Line 3260  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 3883  c      print*,icl,cl_used(icl),cl_good(i Line 3307  c      print*,icl,cl_used(icl),cl_good(i
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.6  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.23