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

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

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

revision 1.1 by mocchiut, Fri May 19 13:15:54 2006 UTC revision 1.11 by pam-fi, Tue Nov 7 15:55:11 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'
21        include 'level1.f'  c      include 'level1.f'
22        include 'level2.f'        include 'level2.f'
23    
24        include 'momanhough_init.f'  c      include 'momanhough_init.f'
25                
       logical DEBUG  
       common/dbg/DEBUG  
   
26  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
27  *     STEP 1  *     STEP 1
28  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 47  Line 45 
45  c      iflag=0  c      iflag=0
46        call cl_to_couples(iflag)        call cl_to_couples(iflag)
47        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
48           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
49        endif        endif
50                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     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  
           
   
   
51  *-----------------------------------------------------  *-----------------------------------------------------
52  *-----------------------------------------------------  *-----------------------------------------------------
53  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 94  c$$$         endif Line 79  c$$$         endif
79  c      iflag=0  c      iflag=0
80        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
81        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
82           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
83        endif        endif
84                
85                
# Line 123  c      iflag=0 Line 108  c      iflag=0
108  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
111    *     count number of hit planes
112          planehit=0                
113          do np=1,nplanes          
114            if(ncp_plane(np).ne.0)then  
115              planehit=planehit+1  
116            endif                  
117          enddo                    
118          if(planehit.lt.3) goto 880 ! exit              
119          
120          nptxz_min=x_min_start              
121          nplxz_min=x_min_start              
122                
123          nptyz_min=y_min_start              
124          nplyz_min=y_min_start              
125    
126  c      iflag=0        cutdistyz=cutystart      
127          cutdistxz=cutxstart      
128    
129     878  continue
130        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
131        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
132           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
133        endif        endif
134  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
135            if(cutdistyz.lt.maxcuty/2)then
136              cutdistyz=cutdistyz+cutystep
137            else
138              cutdistyz=cutdistyz+(3*cutystep)
139            endif
140            goto 878                
141          endif                    
142    
143          if(planehit.eq.3) goto 881          
144          
145     879  continue  
146        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
147        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
148           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
149        endif        endif
150                                    
151          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
152            cutdistxz=cutdistxz+cutxstep
153            goto 879                
154          endif                    
155    
156        
157     881  continue  
158    *     if there is at least three planes on the Y view decreases cuts on X view
159          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
160         $     nplxz_min.ne.y_min_start)then
161            nptxz_min=x_min_step    
162            nplxz_min=x_min_start-x_min_step              
163            goto 879                
164          endif                    
165            
166   880  return   880  return
167        end        end
168    
# Line 144  c      iflag=0 Line 172  c      iflag=0
172        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
173    
174        include 'commontracker.f'        include 'commontracker.f'
175          include 'level1.f'
176        include 'common_momanhough.f'        include 'common_momanhough.f'
177        include 'common_mech.f'        include 'common_mech.f'
178        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
179        include 'common_mini_2.f'        include 'common_mini_2.f'
180        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
181        include 'level2.f'        include 'level2.f'
182    
183        include 'momanhough_init.f'  c      include 'momanhough_init.f'
184                
       logical DEBUG  
       common/dbg/DEBUG  
   
185        logical FIMAGE            !        logical FIMAGE            !
186          real*8 AL_GUESS(5)
187    
188  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
189  *     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 224  c         iflag=0
224           ibest=0                !best track among candidates           ibest=0                !best track among candidates
225           iimage=0               !track image           iimage=0               !track image
226  *     ------------- select the best track -------------  *     ------------- select the best track -------------
227    c$$$         rchi2best=1000000000.
228    c$$$         do i=1,ntracks
229    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
230    c$$$     $         RCHI2_STORE(i).gt.0)then
231    c$$$               ibest=i
232    c$$$               rchi2best=RCHI2_STORE(i)
233    c$$$            endif
234    c$$$         enddo
235    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
236    
237           rchi2best=1000000000.           rchi2best=1000000000.
238             ndofbest=0             !(1)
239           do i=1,ntracks           do i=1,ntracks
240              if(RCHI2_STORE(i).lt.rchi2best.and.             if(RCHI2_STORE(i).lt.rchi2best.and.
241       $         RCHI2_STORE(i).gt.0)then       $          RCHI2_STORE(i).gt.0)then
242                 ndof=0             !(1)
243                 do ii=1,nplanes    !(1)
244                   ndof=ndof        !(1)
245         $              +int(xgood_store(ii,i)) !(1)
246         $              +int(ygood_store(ii,i)) !(1)
247                 enddo              !(1)
248                 if(ndof.ge.ndofbest)then !(1)
249                 ibest=i                 ibest=i
250                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
251              endif                 ndofbest=ndof    !(1)
252                 endif              !(1)
253               endif
254           enddo           enddo
255           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
256  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
# Line 236  c         iflag=0 Line 282  c         iflag=0
282  *     **********************************************************  *     **********************************************************
283  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
284  *     **********************************************************  *     **********************************************************
285             call guess()
286             do i=1,5
287                AL_GUESS(i)=AL(i)
288             enddo
289    
290           do i=1,5           do i=1,5
291              AL(i)=dble(AL_STORE(i,icand))              AL(i)=dble(AL_STORE(i,icand))            
292           enddo           enddo
293            
294             IDCAND = icand         !fitted track-candidate
295           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
296           jstep=0                !# minimization steps           jstep=0                !# minimization steps
297    
298           call mini_2(jstep,ifail)           iprint=0
299    c         if(DEBUG)iprint=1
300             if(.true.)iprint=1
301             call mini2(jstep,ifail,iprint)
302           if(ifail.ne.0) then           if(ifail.ne.0) then
303              if(DEBUG)then              if(.true.)then
304                 print *,                 print *,
305       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
306       $              ,iev       $              ,iev
307    
308                   print*,'guess:   ',(al_guess(i),i=1,5)
309                   print*,'previous: ',(al_store(i,icand),i=1,5)
310                   print*,'result:   ',(al(i),i=1,5)
311                   print*,'xgood ',xgood
312                   print*,'ygood ',ygood
313                   print*,'----------------------------------------------'
314              endif              endif
315              chi2=-chi2  c            chi2=-chi2
316           endif           endif
317                    
318           if(DEBUG)then           if(DEBUG)then
# Line 310  c         print*,'++++++++++ iimage,fima Line 373  c         print*,'++++++++++ iimage,fima
373  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
374    
375           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
376              if(DEBUG)              if(verbose)
377       $           print*,       $           print*,
378       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
379       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 606  c                (implemented variable r Line 669  c                (implemented variable r
669  c*****************************************************  c*****************************************************
670                
671        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
672        include 'level1.f'        include 'level1.f'
673          include 'calib.f'
674    c      include 'level1.f'
675        include 'common_align.f'        include 'common_align.f'
676        include 'common_mech.f'        include 'common_mech.f'
677        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
678        include 'common_resxy.f'  c      include 'common_resxy.f'
679    
680        logical DEBUG  c      logical DEBUG
681        common/dbg/DEBUG  c      common/dbg/DEBUG
682    
683        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
684        integer sensor        integer sensor
# Line 666  c      double precision xi_B,yi_B,zi_B Line 730  c      double precision xi_B,yi_B,zi_B
730              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
731           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
732  c            cog2 = cog(2,icx)  c            cog2 = cog(2,icx)
733  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)              c            etacorr = pfaeta2(cog2,viewx,nldx,angx)            
734  c            stripx = stripx + etacorr  c            stripx = stripx + etacorr
735              stripx = stripx + pfa_eta2(icx,angx)            !(3)              stripx = stripx + pfaeta2(icx,angx)            !(3)
736              resxPAM = risx_eta2(angx)                       !   (4)              resxPAM = risx_eta2(angx)                       !   (4)
737              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
738       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
739              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
740           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                         !(3)
741              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)            !(3)
742              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(angx)                       !   (4)
743              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)
744       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
745              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)
746           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                         !(3)
747              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            !(3)
748              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(angx)                       !   (4)
749              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)
750       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
751              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)
752           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then                          !(3)
753              stripx = stripx + pfa_eta(icx,angx)             !(3)              stripx = stripx + pfaeta(icx,angx)             !(3)
754              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = ris_eta(icx,angx)                     !   (4)
755              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)
756       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
# Line 701  c            resxPAM = resxPAM*fbad_cog( Line 765  c            resxPAM = resxPAM*fbad_cog(
765           endif           endif
766    
767        endif        endif
768    c      if(icy.eq.0.and.icx.ne.0)
769    c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'
770                
771  *     -----------------  *     -----------------
772  *     CLUSTER Y  *     CLUSTER Y
# Line 728  c            resxPAM = resxPAM*fbad_cog( Line 794  c            resxPAM = resxPAM*fbad_cog(
794              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
795           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
796  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
797  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
798  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
799              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
800              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(angy)                       !   (4)
801              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
802              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
803       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
804           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
805              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
806              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
807              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
808       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
809           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
810              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
811              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
812              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
813       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
814           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
815              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
816              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
817  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
818              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 772  C======================================= Line 838  C=======================================
838  c------------------------------------------------------------------------  c------------------------------------------------------------------------
839  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
840  c------------------------------------------------------------------------  c------------------------------------------------------------------------
841             if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
842         $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
843                print*,'xyz_PAM (couple):',
844         $          ' WARNING: false X strip: strip ',stripx
845             endif
846           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
847           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
848           zi = 0.           zi = 0.
# Line 858  C======================================= Line 929  C=======================================
929              nldy = nldx              nldy = nldx
930              viewy = viewx - 1              viewy = viewx - 1
931    
932    c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
933    c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
934                if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
935         $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
936                   print*,'xyz_PAM (X-singlet):',
937         $             ' WARNING: false X strip: strip ',stripx
938                endif
939              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
940    
941              xi_A = xi              xi_A = xi
# Line 1136  c$$$         print*,' resolution ',resxP Line 1214  c$$$         print*,' resolution ',resxP
1214  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1215  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1216  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1217                   if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1218         $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1219    c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1220                      print*,'whichsensor: ',
1221         $                ' WARNING: false X strip: strip ',stripx
1222                   endif
1223                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1224                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
1225                 zi = 0.                 zi = 0.
# Line 1263  c     $              ,iv,xvv(iv),yvv(iv) Line 1347  c     $              ,iv,xvv(iv),yvv(iv)
1347  *     it returns the plane number  *     it returns the plane number
1348  *      *    
1349        include 'commontracker.f'        include 'commontracker.f'
1350          include 'level1.f'
1351  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1352        include 'common_momanhough.f'        include 'common_momanhough.f'
1353                
# Line 1300  c      include 'common_analysis.f' Line 1385  c      include 'common_analysis.f'
1385  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1386  *      *    
1387        include 'commontracker.f'        include 'commontracker.f'
1388          include 'level1.f'
1389  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1390        include 'common_momanhough.f'        include 'common_momanhough.f'
1391                
# Line 1328  c      include 'common_analysis.f' Line 1414  c      include 'common_analysis.f'
1414  *     positive if sensor =2  *     positive if sensor =2
1415  *  *
1416        include 'commontracker.f'        include 'commontracker.f'
1417          include 'level1.f'
1418  c      include 'calib.f'  c      include 'calib.f'
1419  c      include 'level1.f'  c      include 'level1.f'
1420  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1639  c$$$      end Line 1726  c$$$      end
1726        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1727    
1728        include 'commontracker.f'        include 'commontracker.f'
1729          include 'level1.f'
1730        include 'common_momanhough.f'        include 'common_momanhough.f'
1731        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1732        include 'calib.f'        include 'calib.f'
1733        include 'level1.f'  c      include 'level1.f'
   
       logical DEBUG  
       common/dbg/DEBUG  
1734    
1735  *     output flag  *     output flag
1736  *     --------------  *     --------------
# Line 1654  c$$$      end Line 1739  c$$$      end
1739  *     --------------  *     --------------
1740        integer iflag        integer iflag
1741    
1742        integer badseed,badcl        integer badseed,badclx,badcly
1743    
1744  *     init variables  *     init variables
1745        ncp_tot=0        ncp_tot=0
# Line 1670  c$$$      end Line 1755  c$$$      end
1755           ncls(ip)=0           ncls(ip)=0
1756        enddo        enddo
1757        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1758           cl_single(icl)=1           cl_single(icl) = 1
1759           cl_good(icl)=0           cl_good(icl)   = 0
1760          enddo
1761          do iv=1,nviews
1762             ncl_view(iv)  = 0
1763             mask_view(iv) = 0      !all included
1764        enddo        enddo
1765                
1766    *     count number of cluster per view
1767          do icl=1,nclstr1
1768             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1769          enddo
1770    *     mask views with too many clusters
1771          do iv=1,nviews
1772             if( ncl_view(iv).gt. nclustermax)then
1773                mask_view(iv) = 1
1774                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1775         $           ,nclustermax,' on view ', iv,' --> masked!'
1776             endif
1777          enddo
1778    
1779    
1780  *     start association  *     start association
1781        ncouples=0        ncouples=0
1782        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1783           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1784                    
1785  *     ----------------------------------------------------  *     ----------------------------------------------------
1786    *     jump masked views (X VIEW)
1787    *     ----------------------------------------------------
1788             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1789    *     ----------------------------------------------------
1790  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1791    *     ----------------------------------------------------
1792           if(dedx(icx).lt.dedx_x_min)then           if(dedx(icx).lt.dedx_x_min)then
1793              cl_single(icx)=0              cl_single(icx)=0
1794              goto 10              goto 10
1795           endif           endif
1796    *     ----------------------------------------------------
1797  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1798    *     ----------------------------------------------------
1799           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1800           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1801           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1693  c$$$      end Line 1803  c$$$      end
1803           else           else
1804              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1805           endif           endif
1806           badcl=badseed           badclx=badseed
1807           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1808              ibad=1              ibad=1
1809              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1703  c$$$      end Line 1813  c$$$      end
1813       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1814       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1815              endif              endif
1816              badcl=badcl*ibad              badclx=badclx*ibad
1817           enddo           enddo
1818    *     ----------------------------------------------------
1819    *     >>> eliminato il taglio sulle BAD <<<
1820    *     ----------------------------------------------------
1821  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1822  c     cl_single(icx)=0  c     cl_single(icx)=0
1823  c     goto 10  c     goto 10
# Line 1719  c     endif Line 1832  c     endif
1832              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1833                            
1834  *     ----------------------------------------------------  *     ----------------------------------------------------
1835    *     jump masked views (Y VIEW)
1836    *     ----------------------------------------------------
1837                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1838    
1839    *     ----------------------------------------------------
1840  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1841    *     ----------------------------------------------------
1842              if(dedx(icy).lt.dedx_y_min)then              if(dedx(icy).lt.dedx_y_min)then
1843                 cl_single(icy)=0                 cl_single(icy)=0
1844                 goto 20                 goto 20
1845              endif              endif
1846    *     ----------------------------------------------------
1847  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1848    *     ----------------------------------------------------
1849              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1850              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1851              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 1732  c     endif Line 1853  c     endif
1853              else              else
1854                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1855              endif              endif
1856              badcl=badseed              badcly=badseed
1857              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1858                 ibad=1                 ibad=1
1859                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1741  c     endif Line 1862  c     endif
1862       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1863       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1864       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1865                 badcl=badcl*ibad                 badcly=badcly*ibad
1866              enddo              enddo
1867    *     ----------------------------------------------------
1868    *     >>> eliminato il taglio sulle BAD <<<
1869    *     ----------------------------------------------------
1870  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1871  c     cl_single(icy)=0  c     cl_single(icy)=0
1872  c     goto 20  c     goto 20
1873  c     endif  c     endif
1874  *     ----------------------------------------------------  *     ----------------------------------------------------
1875                            
               
1876              cl_good(icy)=1                                cl_good(icy)=1                  
1877              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
1878              nldy=nld(MAXS(icy),VIEW(icy))              nldy=nld(MAXS(icy),VIEW(icy))
# Line 1760  c     endif Line 1883  c     endif
1883  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1884              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1885  *     charge correlation  *     charge correlation
1886                 ddd=(dedx(icy)  *     (modified to be applied only below saturation... obviously)
1887       $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
1888                 ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1889                 cut=chcut*sch(nplx,nldx)       $              .and.
1890                 if(abs(ddd).gt.cut)goto 20 !charge not consistent       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1891                       $              .and.
1892         $              (badclx.eq.1.and.badcly.eq.1)
1893         $              .and.
1894         $              .true.)then
1895    
1896                      ddd=(dedx(icy)
1897         $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1898                      ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1899    
1900    c                  cut = chcut * sch(nplx,nldx)
1901    
1902                      sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1903         $                 -kch(nplx,nldx)*cch(nplx,nldx))
1904                      sss=sss/sqrt(kch(nplx,nldx)**2+1)
1905                      cut = chcut * (16 + sss/50.)
1906    
1907                      if(abs(ddd).gt.cut)then
1908                         goto 20    !charge not consistent
1909                      endif
1910                   endif
1911                                
1912  *     ------------------> COUPLE <------------------  *     ------------------> COUPLE <------------------
1913  *     check to do not overflow vector dimentions  *     check to do not overflow vector dimentions
1914                 if(ncp_plane(nplx).gt.ncouplemax)then  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
1915                    if(DEBUG)print*,  c$$$                  if(DEBUG)print*,
1916       $                    ' ** warning ** number of identified'//  c$$$     $                    ' ** warning ** number of identified'//
1917       $                    ' couples on plane ',nplx,  c$$$     $                    ' couples on plane ',nplx,
1918       $                    ' exceeds vector dimention'//  c$$$     $                    ' exceeds vector dimention'//
1919       $                    ' ( ',ncouplemax,' )'  c$$$     $                    ' ( ',ncouplemax,' )'
1920  c     good2=.false.  c$$$c     good2=.false.
1921  c     goto 880   !fill ntp and go to next event  c$$$c     goto 880   !fill ntp and go to next event
1922                    iflag=1  c$$$                  iflag=1
1923                    return  c$$$                  return
1924                 endif  c$$$               endif
1925                                
1926                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).eq.ncouplemax)then
1927                    if(DEBUG)print*,                    if(verbose)print*,
1928       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1929       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1930       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1931       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' )'
1932  c     good2=.false.  c     good2=.false.
1933  c     goto 880   !fill ntp and go to next event                      c     goto 880   !fill ntp and go to next event
1934                    iflag=1                    mask_view(nviewx(nplx)) = 2
1935                    return                    mask_view(nviewy(nply)) = 2
1936    c                  iflag=1
1937    c                  return
1938                 endif                 endif
1939                                
1940                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
# Line 1825  c     goto 880   !fill ntp and go to nex Line 1969  c     goto 880   !fill ntp and go to nex
1969        endif        endif
1970                
1971        do ip=1,6        do ip=1,6
1972           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1973        enddo        enddo
1974  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
1975                
1976        if(ncp_tot.gt.ncp_max)then  c$$$      if(ncp_tot.gt.ncp_max)then
1977           if(DEBUG)print*,  c$$$         if(verbose)print*,
1978       $           '** warning ** number of identified '//  c$$$     $           '** warning ** number of identified '//
1979       $           'couples exceeds upper limit for Hough tr. '  c$$$     $           'couples exceeds upper limit for Hough tr. '
1980       $           ,'( ',ncp_max,' )'              c$$$     $           ,'( ',ncp_max,' )'            
1981  c            good2=.false.  c$$$         iflag=1
1982  c     goto 880       !fill ntp and go to next event  c$$$         return
1983           iflag=1  c$$$      endif
          return  
       endif  
1984                
1985        return        return
1986        end        end
# Line 1851  c     goto 880       !fill ntp and go to Line 1993  c     goto 880       !fill ntp and go to
1993  *                                                 *  *                                                 *
1994  *                                                 *  *                                                 *
1995  **************************************************  **************************************************
1996        subroutine cl_to_couples_nocharge(iflag)  c$$$      subroutine cl_to_couples_nocharge(iflag)
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
       logical DEBUG  
       common/dbg/DEBUG  
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     init variables  
       ncp_tot=0  
       do ip=1,nplanes  
          do ico=1,ncouplemax  
             clx(ip,ico)=0  
             cly(ip,ico)=0  
          enddo  
          ncp_plane(ip)=0  
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
 *     start association  
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
 *     ----------------------------------------------------  
 *     cut on charge (X VIEW)  
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     cut BAD (X VIEW)              
          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
          ifirst=INDSTART(icx)  
          if(icx.ne.nclstr1) then  
             ilast=INDSTART(icx+1)-1  
          else  
             ilast=TOTCLLENGTH  
          endif  
          badcl=badseed  
          do igood=-ngoodstr,ngoodstr  
             ibad=1  
             if((INDMAX(icx)+igood).gt.ifirst.and.  
      $           (INDMAX(icx)+igood).lt.ilast.and.  
      $           .true.)then  
                ibad=BAD(VIEW(icx),  
      $              nvk(MAXS(icx)+igood),  
      $              nst(MAXS(icx)+igood))  
             endif  
             badcl=badcl*ibad  
          enddo  
          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
             goto 10             !<<<<<<<<<<<<<< BAD cut  
          endif                  !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
           
          cl_good(icx)=1  
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
 *     ----------------------------------------------------  
 *     cut on charge (Y VIEW)  
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     cut BAD (Y VIEW)              
             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
             ifirst=INDSTART(icy)  
             if(icy.ne.nclstr1) then  
                ilast=INDSTART(icy+1)-1  
             else  
                ilast=TOTCLLENGTH  
             endif  
             badcl=badseed  
             do igood=-ngoodstr,ngoodstr  
                ibad=1  
                if((INDMAX(icy)+igood).gt.ifirst.and.  
      $              (INDMAX(icy)+igood).lt.ilast.and.  
      $              .true.)  
      $              ibad=BAD(VIEW(icy),  
      $              nvk(MAXS(icy)+igood),  
      $              nst(MAXS(icy)+igood))  
                badcl=badcl*ibad  
             enddo  
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
               
               
             cl_good(icy)=1                    
             nply=npl(VIEW(icy))  
             nldy=nld(MAXS(icy),VIEW(icy))  
               
 *     ----------------------------------------------  
 *     CONDITION TO FORM A COUPLE  
 *     ----------------------------------------------  
 *     geometrical consistency (same plane and ladder)  
             if(nply.eq.nplx.and.nldy.eq.nldx)then  
 *     charge correlation  
 *     ===========================================================  
 *     this version of the subroutine is used for the calibration  
 *     thus charge-correlation selection is obviously removed  
 *     ===========================================================  
 c$$$               ddd=(dedx(icy)  
 c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 *     ===========================================================  
                 
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
                if(ncp_plane(nplx).gt.ncouplemax)then  
                   if(DEBUG)print*,  
      $                    ' ** warning ** number of identified'//  
      $                    ' couples on plane ',nplx,  
      $                    ' exceeds vector dimention'//  
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(DEBUG)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
                 
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
   
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
   
 c$$$      subroutine cl_to_couples_2(iflag)  
1997  c$$$  c$$$
1998  c$$$      include 'commontracker.f'  c$$$      include 'commontracker.f'
1999    c$$$      include 'level1.f'
2000  c$$$      include 'common_momanhough.f'  c$$$      include 'common_momanhough.f'
2001  c$$$      include 'momanhough_init.f'  c$$$c      include 'momanhough_init.f'
2002  c$$$      include 'calib.f'  c$$$      include 'calib.f'
2003  c$$$      include 'level1.f'  c$$$c      include 'level1.f'
2004  c$$$  c$$$
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
2005  c$$$  c$$$
2006  c$$$*     output flag  c$$$*     output flag
2007  c$$$*     --------------  c$$$*     --------------
# Line 2132  c$$$     $              nst(MAXS(icx)+ig Line 2061  c$$$     $              nst(MAXS(icx)+ig
2061  c$$$            endif  c$$$            endif
2062  c$$$            badcl=badcl*ibad  c$$$            badcl=badcl*ibad
2063  c$$$         enddo  c$$$         enddo
2064  c$$$*         print*,'icx ',icx,badcl  c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut
2065  c$$$         if(badcl.eq.0)then  c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut
2066  c$$$            cl_single(icx)=0  c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut
2067  c$$$            goto 10  c$$$         endif                  !<<<<<<<<<<<<<< BAD cut
 c$$$         endif  
2068  c$$$*     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2069  c$$$          c$$$        
2070  c$$$         cl_good(icx)=1  c$$$         cl_good(icx)=1
# Line 2171  c$$$     $              nvk(MAXS(icy)+ig Line 2099  c$$$     $              nvk(MAXS(icy)+ig
2099  c$$$     $              nst(MAXS(icy)+igood))  c$$$     $              nst(MAXS(icy)+igood))
2100  c$$$               badcl=badcl*ibad  c$$$               badcl=badcl*ibad
2101  c$$$            enddo  c$$$            enddo
2102  c$$$*            print*,'icy ',icy,badcl  c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut
2103  c$$$            if(badcl.eq.0)then  c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut
2104  c$$$               cl_single(icy)=0  c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut
2105  c$$$               goto 20  c$$$            endif               !<<<<<<<<<<<<<< BAD cut
 c$$$            endif  
2106  c$$$*     ----------------------------------------------------  c$$$*     ----------------------------------------------------
2107  c$$$              c$$$            
2108  c$$$              c$$$            
# Line 2188  c$$$*     CONDITION TO FORM A COUPLE Line 2115  c$$$*     CONDITION TO FORM A COUPLE
2115  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
2116  c$$$*     geometrical consistency (same plane and ladder)  c$$$*     geometrical consistency (same plane and ladder)
2117  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then
2118  c$$$  c$$$*     charge correlation
2119  c$$$c$$$*     charge correlation  c$$$*     ===========================================================
2120    c$$$*     this version of the subroutine is used for the calibration
2121    c$$$*     thus charge-correlation selection is obviously removed
2122    c$$$*     ===========================================================
2123  c$$$c$$$               ddd=(dedx(icy)  c$$$c$$$               ddd=(dedx(icy)
2124  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
2125  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
2126  c$$$c$$$               cut=chcut*sch(nplx,nldx)  c$$$c$$$               cut=chcut*sch(nplx,nldx)
2127  c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent
2128    c$$$*     ===========================================================
2129    c$$$              
2130  c$$$                c$$$              
2131  c$$$*     ------------------> COUPLE <------------------  c$$$*     ------------------> COUPLE <------------------
2132  c$$$*     check to do not overflow vector dimentions  c$$$*     check to do not overflow vector dimentions
2133  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  c$$$c$$$               if(ncp_plane(nplx).gt.ncouplemax)then
2134  c$$$                  if(DEBUG)print*,  c$$$c$$$                  if(DEBUG)print*,
2135  c$$$     $                    ' ** warning ** number of identified'//  c$$$c$$$     $                    ' ** warning ** number of identified'//
2136  c$$$     $                    ' couples on plane ',nplx,  c$$$c$$$     $                    ' couples on plane ',nplx,
2137  c$$$     $                    ' exceeds vector dimention'//  c$$$c$$$     $                    ' exceeds vector dimention'//
2138  c$$$     $                    ' ( ',ncouplemax,' )'  c$$$c$$$     $                    ' ( ',ncouplemax,' )'
2139  c$$$c     good2=.false.  c$$$c$$$c     good2=.false.
2140  c$$$c     goto 880   !fill ntp and go to next event  c$$$c$$$c     goto 880   !fill ntp and go to next event
2141  c$$$                  iflag=1  c$$$c$$$                  iflag=1
2142  c$$$                  return  c$$$c$$$                  return
2143  c$$$               endif  c$$$c$$$               endif
2144  c$$$                c$$$              
2145  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then
2146  c$$$                  if(DEBUG)print*,  c$$$                  if(verbose)print*,
2147  c$$$     $                 '** warning ** number of identified '//  c$$$     $                 '** warning ** number of identified '//
2148  c$$$     $                 'couples on plane ',nplx,  c$$$     $                 'couples on plane ',nplx,
2149  c$$$     $                 'exceeds vector dimention '  c$$$     $                 'exceeds vector dimention '
# Line 2227  c$$$               clx(nplx,ncp_plane(np Line 2159  c$$$               clx(nplx,ncp_plane(np
2159  c$$$               cly(nply,ncp_plane(nplx))=icy  c$$$               cly(nply,ncp_plane(nplx))=icy
2160  c$$$               cl_single(icx)=0  c$$$               cl_single(icx)=0
2161  c$$$               cl_single(icy)=0  c$$$               cl_single(icy)=0
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
2162  c$$$            endif                                c$$$            endif                              
2163  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
2164  c$$$  c$$$
# Line 2260  c$$$      enddo Line 2191  c$$$      enddo
2191  c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
2192  c$$$        c$$$      
2193  c$$$      if(ncp_tot.gt.ncp_max)then  c$$$      if(ncp_tot.gt.ncp_max)then
2194  c$$$         if(DEBUG)print*,  c$$$         if(verbose)print*,
2195  c$$$     $           '** warning ** number of identified '//  c$$$     $           '** warning ** number of identified '//
2196  c$$$     $           'couples exceeds upper limit for Hough tr. '  c$$$     $           'couples exceeds upper limit for Hough tr. '
2197  c$$$     $           ,'( ',ncp_max,' )'              c$$$     $           ,'( ',ncp_max,' )'            
# Line 2272  c$$$      endif Line 2203  c$$$      endif
2203  c$$$        c$$$      
2204  c$$$      return  c$$$      return
2205  c$$$      end  c$$$      end
2206    c$$$
2207                
2208  ***************************************************  ***************************************************
2209  *                                                 *  *                                                 *
# Line 2288  c     02/02/2006 modified by Elena Vannu Line 2220  c     02/02/2006 modified by Elena Vannu
2220  c*****************************************************  c*****************************************************
2221    
2222        include 'commontracker.f'        include 'commontracker.f'
2223          include 'level1.f'
2224        include 'common_momanhough.f'        include 'common_momanhough.f'
2225        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2226        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2227        include 'common_mini_2.f'        include 'common_mini_2.f'
2228        include 'calib.f'        include 'calib.f'
2229        include 'level1.f'  c      include 'level1.f'
2230    
       logical DEBUG  
       common/dbg/DEBUG  
2231    
2232  *     output flag  *     output flag
2233  *     --------------  *     --------------
# Line 2364  c     $                       (icx2,icy2 Line 2295  c     $                       (icx2,icy2
2295  *     (2 couples needed)  *     (2 couples needed)
2296  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2297                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2298                             if(DEBUG)print*,                             if(verbose)print*,
2299       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2300       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2301       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2302  c                           good2=.false.  c                           good2=.false.
2303  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2304                               do iv=1,12
2305                                  mask_view(iv) = 3
2306                               enddo
2307                             iflag=1                             iflag=1
2308                             return                             return
2309                          endif                          endif
# Line 2434  c     $                                 Line 2368  c     $                                
2368  *     (3 couples needed)  *     (3 couples needed)
2369  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2370                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2371                                      if(DEBUG)print*,                                      if(verbose)print*,
2372       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2373       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2374       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2375  c                                    good2=.false.  c                                    good2=.false.
2376  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2377                                        do iv=1,nviews
2378                                           mask_view(iv) = 4
2379                                        enddo
2380                                      iflag=1                                      iflag=1
2381                                      return                                      return
2382                                   endif                                   endif
# Line 2514  c     goto 880               !ntp fill Line 2451  c     goto 880               !ntp fill
2451        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2452    
2453        include 'commontracker.f'        include 'commontracker.f'
2454          include 'level1.f'
2455        include 'common_momanhough.f'        include 'common_momanhough.f'
2456        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2457    
       logical DEBUG  
       common/dbg/DEBUG  
2458    
2459  *     output flag  *     output flag
2460  *     --------------  *     --------------
# Line 2550  c     goto 880               !ntp fill Line 2486  c     goto 880               !ntp fill
2486        distance=0        distance=0
2487        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2488        npt_tot=0        npt_tot=0
2489          nloop=0                  
2490     90   continue                  
2491        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2492           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
2493                            
# Line 2653  c     print*,'*   idbref,idb2 ',idbref,i Line 2591  c     print*,'*   idbref,idb2 ',idbref,i
2591              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2592           enddo           enddo
2593  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
2594           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
2595           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
2596           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2597                    
# Line 2661  c     print*,'>>>> ',ncpused,npt,nplused Line 2599  c     print*,'>>>> ',ncpused,npt,nplused
2599  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2600    
2601           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2602              if(DEBUG)print*,              if(verbose)print*,
2603       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2604       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2605       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2606  c               good2=.false.  c               good2=.false.
2607  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2608                do iv=1,nviews
2609                   mask_view(iv) = 5
2610                enddo
2611              iflag=1              iflag=1
2612              return              return
2613           endif           endif
# Line 2704  c$$$     $           ,(db_cloud(iii),iii Line 2645  c$$$     $           ,(db_cloud(iii),iii
2645        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2646                
2647                
2648          if(nloop.lt.nstepy)then      
2649            cutdistyz = cutdistyz+cutystep
2650            nloop     = nloop+1          
2651            goto 90                
2652          endif                    
2653          
2654        if(DEBUG)then        if(DEBUG)then
2655           print*,'---------------------- '           print*,'---------------------- '
2656           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2730  c$$$     $           ,(db_cloud(iii),iii Line 2677  c$$$     $           ,(db_cloud(iii),iii
2677        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2678    
2679        include 'commontracker.f'        include 'commontracker.f'
2680          include 'level1.f'
2681        include 'common_momanhough.f'        include 'common_momanhough.f'
2682        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2683    
       logical DEBUG  
       common/dbg/DEBUG  
2684    
2685  *     output flag  *     output flag
2686  *     --------------  *     --------------
# Line 2765  c$$$     $           ,(db_cloud(iii),iii Line 2711  c$$$     $           ,(db_cloud(iii),iii
2711        distance=0        distance=0
2712        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2713        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2714          nloop=0                  
2715     91   continue                  
2716        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2717           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
2718  c     print*,'--------------'  c     print*,'--------------'
# Line 2866  c     print*,'check cp_used' Line 2814  c     print*,'check cp_used'
2814           do ip=1,nplanes           do ip=1,nplanes
2815              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2816           enddo           enddo
2817           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2818           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2819           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2820                    
2821  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2822  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2823           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2824              if(DEBUG)print*,              if(verbose)print*,
2825       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2826       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2827       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2828  c     good2=.false.  c     good2=.false.
2829  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2830                do iv=1,nviews
2831                   mask_view(iv) = 6
2832                enddo
2833              iflag=1              iflag=1
2834              return              return
2835           endif           endif
# Line 2914  c$$$     $           ,(tr_cloud(iii),iii Line 2865  c$$$     $           ,(tr_cloud(iii),iii
2865  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2866  22288    continue  22288    continue
2867        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2868          
2869           if(nloop.lt.nstepx)then      
2870             cutdistxz=cutdistxz+cutxstep
2871             nloop=nloop+1          
2872             goto 91                
2873           endif                    
2874          
2875        if(DEBUG)then        if(DEBUG)then
2876           print*,'---------------------- '           print*,'---------------------- '
2877           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2941  c     02/02/2006 modified by Elena Vannu Line 2898  c     02/02/2006 modified by Elena Vannu
2898  c*****************************************************  c*****************************************************
2899    
2900        include 'commontracker.f'        include 'commontracker.f'
2901          include 'level1.f'
2902        include 'common_momanhough.f'        include 'common_momanhough.f'
2903        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2904        include 'common_mini_2.f'        include 'common_mini_2.f'
2905        include 'common_mech.f'        include 'common_mech.f'
2906        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2907    
       logical DEBUG  
       common/dbg/DEBUG  
2908    
2909  *     output flag  *     output flag
2910  *     --------------  *     --------------
# Line 2964  c*************************************** Line 2920  c***************************************
2920  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2921  *     list of matching couples in the combination  *     list of matching couples in the combination
2922  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2923        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2924        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2925  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2926        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2927  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2928  *     variables for track fitting  *     variables for track fitting
2929        double precision AL_INI(5)        double precision AL_INI(5)
2930        double precision tath  c      double precision tath
2931  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2932  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2933    
# Line 3064  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 3020  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
3020  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
3021                            
3022  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
3023              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
3024              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
3025              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
3026       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
3027              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
3028              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
3029              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
3030                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
3031  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
3032              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
3033                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
3034  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
3035  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
3036                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
3037  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
3038  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
3039                c$$$            ELSE
3040    c$$$               AL_INI(4) = acos(-1.)/2
3041    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
3042    c$$$            ENDIF
3043    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
3044    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
3045    c$$$            
3046    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
3047    c$$$            
3048    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
3049                            
3050              if(DEBUG)then              if(DEBUG)then
3051                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
3052                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3148  c     $                                 Line 3114  c     $                                
3114  *     **********************************************************  *     **********************************************************
3115  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
3116  *     **********************************************************  *     **********************************************************
3117    cccc  scommentare se si usa al_ini della nuvola
3118    c$$$                              do i=1,5
3119    c$$$                                 AL(i)=AL_INI(i)
3120    c$$$                              enddo
3121                                  call guess()
3122                                do i=1,5                                do i=1,5
3123                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
3124                                enddo                                enddo
3125                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3126                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3127                                call mini_2(jstep,ifail)                                iprint=0
3128    c                              if(DEBUG)iprint=1
3129                                  if(DEBUG)iprint=1
3130                                  call mini2(jstep,ifail,iprint)
3131                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3132                                   if(DEBUG)then                                   if(DEBUG)then
3133                                      print *,                                      print *,
3134       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3135       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
3136                                        print*,'initial guess: '
3137    
3138                                        print*,'AL_INI(1) = ',AL_INI(1)
3139                                        print*,'AL_INI(2) = ',AL_INI(2)
3140                                        print*,'AL_INI(3) = ',AL_INI(3)
3141                                        print*,'AL_INI(4) = ',AL_INI(4)
3142                                        print*,'AL_INI(5) = ',AL_INI(5)
3143                                   endif                                   endif
3144                                   chi2=-chi2  c                                 chi2=-chi2
3145                                endif                                endif
3146  *     **********************************************************  *     **********************************************************
3147  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3173  c     $                                 Line 3154  c     $                                
3154  *     --------------------------  *     --------------------------
3155                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3156                                                                    
3157                                   if(DEBUG)print*,                                   if(verbose)print*,
3158       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3159       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3160       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3161  c                                 good2=.false.  c                                 good2=.false.
3162  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3163                                     do iv=1,nviews
3164                                        mask_view(iv) = 7
3165                                     enddo
3166                                   iflag=1                                   iflag=1
3167                                   return                                   return
3168                                endif                                endif
# Line 3273  c$$$  rchi2=chi2/dble(ndof) Line 3257  c$$$  rchi2=chi2/dble(ndof)
3257  c******************************************************  c******************************************************
3258  cccccc 06/10/2005 modified by elena vannuccini ---> (1)  cccccc 06/10/2005 modified by elena vannuccini ---> (1)
3259  cccccc 31/01/2006 modified by elena vannuccini ---> (2)  cccccc 31/01/2006 modified by elena vannuccini ---> (2)
3260    cccccc 12/08/2006 modified by elena vannucicni ---> (3)
3261  c******************************************************  c******************************************************
3262    
3263        include 'commontracker.f'        include 'commontracker.f'
3264          include 'level1.f'
3265        include 'common_momanhough.f'        include 'common_momanhough.f'
3266        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3267        include 'common_mini_2.f'        include 'common_mini_2.f'
3268        include 'common_mech.f'        include 'common_mech.f'
3269        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3270        include 'level1.f'  c      include 'level1.f'
3271        include 'calib.f'        include 'calib.f'
3272    
       logical DEBUG  
       common/dbg/DEBUG  
3273    
3274  *     flag to chose PFA  *     flag to chose PFA
3275        character*10 PFA        character*10 PFA
# Line 3299  c*************************************** Line 3283  c***************************************
3283        call track_init        call track_init
3284        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3285    
3286    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3287  *     -------------------------------------------------  *     -------------------------------------------------
3288  *     If the plane has been already included, it just  *     If the plane has been already included, it just
3289  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
3290  *     using improved PFAs  *     using improved PFAs
3291  *     -------------------------------------------------  *     -------------------------------------------------
3292    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3293           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3294       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3295                            
# Line 3338  c            dedxtrk(nplanes-ip+1) = (de Line 3324  c            dedxtrk(nplanes-ip+1) = (de
3324              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)
3325              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)
3326                            
3327    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3328  *     -------------------------------------------------  *     -------------------------------------------------
3329  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3330  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3331  *     -------------------------------------------------  *     -------------------------------------------------
3332    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3333           else                             else                  
3334                                
3335              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3383  c            if(DEBUG)print*,'>>>> try t Line 3371  c            if(DEBUG)print*,'>>>> try t
3371                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3372                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3373                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3374       $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3375       $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3376         $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)
3377         $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
3378       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3379  *            *          
3380                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
# Line 3394  c     $              'ETA2','ETA2', Line 3384  c     $              'ETA2','ETA2',
3384       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest))
3385                                
3386                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3387                   distance = distance / RCHI2_STORE(ibest)!<<< MS
3388                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3389                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
3390       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3456  c            if(DEBUG)print*,'>>>> try t Line 3447  c            if(DEBUG)print*,'>>>> try t
3447                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
3448  *                                                !jump to the next couple  *                                                !jump to the next couple
3449  *----- try cluster x -----------------------------------------------  *----- try cluster x -----------------------------------------------
3450                 if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
3451                   if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
3452  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3453                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
3454  c     $              'ETA2','ETA2',  c     $              'ETA2','ETA2',
3455       $              PFA,PFA,       $              PFA,PFA,
3456       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3457                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3458  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
3459                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3460       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
3461                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3483  c                  dedxmm = dedx(icx) !( Line 3475  c                  dedxmm = dedx(icx) !(
3475                 endif                                   endif                  
3476  11881          continue  11881          continue
3477  *----- try cluster y -----------------------------------------------  *----- try cluster y -----------------------------------------------
3478                 if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
3479                   if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3480  *                                              !jump to the next couple  *                                              !jump to the next couple
3481                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
3482  c     $              'ETA2','ETA2',  c     $              'ETA2','ETA2',
3483       $              PFA,PFA,       $              PFA,PFA,
3484       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest))
3485                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3486                   distance = distance / RCHI2_STORE(ibest)!<<< MS
3487                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3488       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
3489                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3512  c                 dedxmm = dedx(icy)  !( Line 3506  c                 dedxmm = dedx(icy)  !(
3506  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------    
3507              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3508                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3509                 if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
3510                   if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
3511       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3512       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3513                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
# Line 3528  c     $                 'ETA2','ETA2', Line 3523  c     $                 'ETA2','ETA2',
3523                 endif                 endif
3524    
3525                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3526                   distance = distance / RCHI2_STORE(ibest)!<<< MS
3527                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3528       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
3529                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3557  c                  dedxmm = dedx(icl)   Line 3553  c                  dedxmm = dedx(icl)  
3553                                
3554                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3555  *              ----------------------------  *              ----------------------------
3556    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3557                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3558                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
3559                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
3560                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
3561       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
3562         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
3563         $                 ,', norm.dist.= ',distmin
3564         $                 ,', cut ',clinc,' )'
3565                 else                 else
3566                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
3567                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
3568                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
3569       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
3570         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
3571         $                 ,', norm.dist.= ', distmin
3572         $                 ,', cut ',clinc,' )'
3573                 endif                 endif
3574    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3575  *              ----------------------------  *              ----------------------------
3576                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
3577                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3596  c              dedxtrk(nplanes-ip+1) = d Line 3600  c              dedxtrk(nplanes-ip+1) = d
3600  *                                                 *  *                                                 *
3601  *                                                 *  *                                                 *
3602  **************************************************  **************************************************
3603    cccccc 12/08/2006 modified by elena ---> (1)
3604    *
3605        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3606    
3607        include 'commontracker.f'        include 'commontracker.f'
3608          include 'level1.f'
3609        include 'common_momanhough.f'        include 'common_momanhough.f'
3610        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3611          include 'level2.f'        !(1)
3612  c      include 'calib.f'  c      include 'calib.f'
3613  c      include 'level1.f'  c      include 'level1.f'
3614    
       logical DEBUG  
       common/dbg/DEBUG  
3615    
3616    
3617        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3617  c      include 'level1.f' Line 3622  c      include 'level1.f'
3622              if(id.ne.0)then              if(id.ne.0)then
3623                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3624                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3625                 cl_used(iclx)=1  !tag used clusters  c               cl_used(iclx)=1  !tag used clusters
3626                 cl_used(icly)=1  !tag used clusters  c               cl_used(icly)=1  !tag used clusters
3627                   cl_used(iclx)=ntrk  !tag used clusters !(1)
3628                   cl_used(icly)=ntrk  !tag used clusters !(1)
3629              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3630                 cl_used(icl)=1   !tag used clusters  c               cl_used(icl)=1   !tag used clusters
3631                   cl_used(icl)=ntrk   !tag used clusters !1)
3632              endif              endif
3633                            
3634  c               if(DEBUG)then  c               if(DEBUG)then
# Line 3779  c$$$ Line 3787  c$$$
3787    
3788        subroutine init_level2        subroutine init_level2
3789    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3790        include 'commontracker.f'        include 'commontracker.f'
3791          include 'level1.f'
3792        include 'common_momanhough.f'        include 'common_momanhough.f'
3793        include 'level2.f'        include 'level2.f'
3794        include 'level1.f'  c      include 'level1.f'
   
   
3795    
3796        good2 = 0!.false.        do i=1,nviews
3797  c$$$      nev2 = nev1           good2(i)=good1(i)
3798          enddo
3799    
 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*****************************************************  
3800    
3801        NTRK = 0        NTRK = 0
3802        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3803           IMAGE(IT)=0           IMAGE(IT)=0
3804           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
          BdL(IT) = 0.  
3805           do ip=1,nplanes           do ip=1,nplanes
3806              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3807              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3826  c*************************************** Line 3810  c***************************************
3810              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3811              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3812              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3813              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3814              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3815  c******************************************************              CLTRX(IP,IT) = 0
3816                CLTRY(IP,IT) = 0
3817           enddo           enddo
3818           do ipa=1,5           do ipa=1,5
3819              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3839  c*************************************** Line 3822  c***************************************
3822              enddo                                enddo                  
3823           enddo                             enddo                  
3824        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3825        nclsx=0        nclsx=0
3826        nclsy=0              nclsy=0      
3827        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3828          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3829          xs(1,ip)=0          xs(1,ip)=0
3830          xs(2,ip)=0          xs(2,ip)=0
3831          sgnlxs(ip)=0          sgnlxs(ip)=0
3832          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3833          ys(1,ip)=0          ys(1,ip)=0
3834          ys(2,ip)=0          ys(2,ip)=0
3835          sgnlys(ip)=0          sgnlys(ip)=0
3836        enddo        enddo
 c*******************************************************  
3837        end        end
3838    
3839    
# Line 3872  c*************************************** Line 3848  c***************************************
3848  ************************************************************  ************************************************************
3849    
3850    
3851          subroutine init_hough
3852    
3853          include 'commontracker.f'
3854          include 'level1.f'
3855          include 'common_momanhough.f'
3856          include 'common_hough.f'
3857          include 'level2.f'
3858    
3859          ntrpt_nt=0
3860          ndblt_nt=0
3861          NCLOUDS_XZ_nt=0
3862          NCLOUDS_YZ_nt=0
3863          do idb=1,ndblt_max_nt
3864             db_cloud_nt(idb)=0
3865             alfayz1_nt(idb)=0      
3866             alfayz2_nt(idb)=0      
3867          enddo
3868          do itr=1,ntrpl_max_nt
3869             tr_cloud_nt(itr)=0
3870             alfaxz1_nt(itr)=0      
3871             alfaxz2_nt(itr)=0      
3872             alfaxz3_nt(itr)=0      
3873          enddo
3874          do idb=1,ncloyz_max      
3875            ptcloud_yz_nt(idb)=0    
3876            alfayz1_av_nt(idb)=0    
3877            alfayz2_av_nt(idb)=0    
3878          enddo                    
3879          do itr=1,ncloxz_max      
3880            ptcloud_xz_nt(itr)=0    
3881            alfaxz1_av_nt(itr)=0    
3882            alfaxz2_av_nt(itr)=0    
3883            alfaxz3_av_nt(itr)=0    
3884          enddo                    
3885    
3886          ntrpt=0                  
3887          ndblt=0                  
3888          NCLOUDS_XZ=0              
3889          NCLOUDS_YZ=0              
3890          do idb=1,ndblt_max        
3891            db_cloud(idb)=0        
3892            cpyz1(idb)=0            
3893            cpyz2(idb)=0            
3894            alfayz1(idb)=0          
3895            alfayz2(idb)=0          
3896          enddo                    
3897          do itr=1,ntrpl_max        
3898            tr_cloud(itr)=0        
3899            cpxz1(itr)=0            
3900            cpxz2(itr)=0            
3901            cpxz3(itr)=0            
3902            alfaxz1(itr)=0          
3903            alfaxz2(itr)=0          
3904            alfaxz3(itr)=0          
3905          enddo                    
3906          do idb=1,ncloyz_max      
3907            ptcloud_yz(idb)=0      
3908            alfayz1_av(idb)=0      
3909            alfayz2_av(idb)=0      
3910            do idbb=1,ncouplemaxtot
3911              cpcloud_yz(idb,idbb)=0
3912            enddo                  
3913          enddo                    
3914          do itr=1,ncloxz_max      
3915            ptcloud_xz(itr)=0      
3916            alfaxz1_av(itr)=0      
3917            alfaxz2_av(itr)=0      
3918            alfaxz3_av(itr)=0      
3919            do itrr=1,ncouplemaxtot
3920              cpcloud_xz(itr,itrr)=0
3921            enddo                  
3922          enddo                    
3923          end
3924    ************************************************************
3925    *
3926    *
3927    *
3928    *
3929    *
3930    *
3931    *
3932    ************************************************************
3933    
3934    
3935        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3936    
3937  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3882  c*************************************** Line 3942  c***************************************
3942    
3943            
3944        include 'commontracker.f'        include 'commontracker.f'
3945    c      include 'level1.f'
3946          include 'level1.f'
3947          include 'common_momanhough.f'
3948        include 'level2.f'        include 'level2.f'
3949        include 'common_mini_2.f'        include 'common_mini_2.f'
3950        real sinth,phi,pig        !(4)        real sinth,phi,pig      
3951        pig=acos(-1.)        pig=acos(-1.)
3952    
       good2=1!.true.  
3953        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3954          nstep_nt(ntr)       = nstep
3955    
3956          phi   = al(4)          
3957          sinth = al(3)            
3958          if(sinth.lt.0)then      
3959             sinth = -sinth        
3960             phi = phi + pig      
3961          endif                    
3962          npig = aint(phi/(2*pig))
3963          phi = phi - npig*2*pig  
3964          if(phi.lt.0)            
3965         $     phi = phi + 2*pig  
3966          al(4) = phi              
3967          al(3) = sinth            
3968    
       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)  
 *****************************************************  
3969        do i=1,5        do i=1,5
3970           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3971           do j=1,5           do j=1,5
3972              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3973           enddo           enddo
 c     print*,al_nt(i,ntr)  
3974        enddo        enddo
3975                
3976        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3924  c     print*,al_nt(i,ntr) Line 3986  c     print*,al_nt(i,ntr)
3986           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3987           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3988           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))
 c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3989           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)
3990           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)                     dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)  
3991      
3992             id  = CP_STORE(ip,IDCAND)
3993             icl = CLS_STORE(ip,IDCAND)
3994             if(id.ne.0)then
3995                cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3996                cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3997    c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3998             elseif(icl.ne.0)then
3999                if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl
4000                if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl
4001    c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
4002             endif          
4003    
4004        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)  
4005    
4006    
4007        end        end
4008    
4009        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*****************************************************  
4010    
4011  *     -------------------------------------------------------  *     -------------------------------------------------------
4012  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3952  c*************************************** Line 4015  c***************************************
4015  *     -------------------------------------------------------  *     -------------------------------------------------------
4016    
4017        include 'commontracker.f'        include 'commontracker.f'
4018        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
4019        include 'calib.f'        include 'calib.f'
4020          include 'level1.f'
4021        include 'common_momanhough.f'        include 'common_momanhough.f'
4022          include 'level2.f'
4023        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
4024    
4025  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
4026        good2=1!.true.  c      good2=1!.true.
4027        nclsx = 0        nclsx = 0
4028        nclsy = 0        nclsy = 0
4029    
4030          do iv = 1,nviews
4031             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
4032          enddo
4033    
4034        do icl=1,nclstr1        do icl=1,nclstr1
4035           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
4036              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
# Line 3970  c*************************************** Line 4038  c***************************************
4038                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4039                 planex(nclsx) = ip                 planex(nclsx) = ip
4040                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
4041                   clsx(nclsx)   = icl
4042                 do is=1,2                 do is=1,2
4043  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4044                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
# Line 3984  c$$$               print*,'xs(2,nclsx)   Line 4053  c$$$               print*,'xs(2,nclsx)  
4053                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4054                 planey(nclsy) = ip                 planey(nclsy) = ip
4055                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
4056                   clsy(nclsy)   = icl
4057                 do is=1,2                 do is=1,2
4058  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4059                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
# Line 3997  c$$$               print*,'ys(2,nclsy)   Line 4067  c$$$               print*,'ys(2,nclsy)  
4067              endif              endif
4068           endif           endif
4069  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
4070    
4071    ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4072             whichtrack(icl) = cl_used(icl)
4073    
4074        enddo        enddo
4075        end        end
4076    
4077    ***************************************************
4078    *                                                 *
4079    *                                                 *
4080    *                                                 *
4081    *                                                 *
4082    *                                                 *
4083    *                                                 *
4084    **************************************************
4085    
4086          subroutine fill_hough
4087    
4088    *     -------------------------------------------------------
4089    *     This routine fills the  variables related to the hough
4090    *     transform, for the debig n-tuple
4091    *     -------------------------------------------------------
4092    
4093          include 'commontracker.f'
4094          include 'level1.f'
4095          include 'common_momanhough.f'
4096          include 'common_hough.f'
4097          include 'level2.f'
4098    
4099          if(.false.
4100         $     .or.ntrpt.gt.ntrpt_max_nt
4101         $     .or.ndblt.gt.ndblt_max_nt
4102         $     .or.NCLOUDS_XZ.gt.ncloxz_max
4103         $     .or.NCLOUDS_yZ.gt.ncloyz_max
4104         $     )then
4105             ntrpt_nt=0
4106             ndblt_nt=0
4107             NCLOUDS_XZ_nt=0
4108             NCLOUDS_YZ_nt=0
4109          else
4110             ndblt_nt=ndblt
4111             ntrpt_nt=ntrpt
4112             if(ndblt.ne.0)then
4113                do id=1,ndblt
4114                   alfayz1_nt(id)=alfayz1(id) !Y0
4115                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4116                enddo
4117             endif
4118             if(ndblt.ne.0)then
4119                do it=1,ntrpt
4120                   alfaxz1_nt(it)=alfaxz1(it) !X0
4121                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4122                   alfaxz3_nt(it)=alfaxz3(it) !1/r
4123                enddo
4124             endif
4125             nclouds_yz_nt=nclouds_yz
4126             nclouds_xz_nt=nclouds_xz
4127             if(nclouds_yz.ne.0)then
4128                nnn=0
4129                do iyz=1,nclouds_yz
4130                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4131                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4132                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4133                   nnn=nnn+ptcloud_yz(iyz)
4134                enddo
4135                do ipt=1,nnn
4136                   db_cloud_nt(ipt)=db_cloud(ipt)
4137                 enddo
4138             endif
4139             if(nclouds_xz.ne.0)then
4140                nnn=0
4141                do ixz=1,nclouds_xz
4142                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4143                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4144                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4145                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4146                   nnn=nnn+ptcloud_xz(ixz)              
4147                enddo
4148                do ipt=1,nnn
4149                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4150                 enddo
4151             endif
4152          endif
4153          end
4154          
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23