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

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

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

revision 1.6 by pam-fi, Mon Oct 2 13:12:48 2006 UTC revision 1.20 by pam-fi, Fri Apr 27 10:39:58 2007 UTC
# Line 12  Line 12 
12        subroutine track_finding(iflag)        subroutine track_finding(iflag)
13    
14        include 'commontracker.f'        include 'commontracker.f'
15          include 'level1.f'
16        include 'common_momanhough.f'        include 'common_momanhough.f'
17        include 'common_mech.f'        include 'common_mech.f'
18        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
19        include 'common_mini_2.f'        include 'common_mini_2.f'
20        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
       include 'momanhough_init.f'  
         
 c      logical DEBUG  
 c      common/dbg/DEBUG  
23    
24    
25    c      print*,'======================================================'
26    c$$$      do ic=1,NCLSTR1
27    c$$$         if(.false.
28    c$$$     $        .or.nsatstrips(ic).gt.0
29    c$$$c     $        .or.nbadstrips(0,ic).gt.0
30    c$$$c     $        .or.nbadstrips(4,ic).gt.0
31    c$$$c     $        .or.nbadstrips(3,ic).gt.0
32    c$$$     $        .or..false.)then
33    c$$$            print*,'--- cl-',ic,' ------------------------'
34    c$$$            istart = INDSTART(IC)
35    c$$$            istop  = TOTCLLENGTH
36    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
37    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
38    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
39    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
40    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
41    c$$$            print*,'view ',VIEW(ic)
42    c$$$            print*,'maxs ',MAXS(ic)
43    c$$$            print*,'COG4 ',cog(4,ic)
44    c$$$            ff = fbad_cog(4,ic)
45    c$$$            print*,'fbad ',ff
46    c$$$            print*,(CLBAD(i),i=istart,istop)
47    c$$$            bb=nbadstrips(0,ic)
48    c$$$            print*,'#BAD (tot)',bb
49    c$$$            bb=nbadstrips(4,ic)
50    c$$$            print*,'#BAD (4)',bb
51    c$$$            bb=nbadstrips(3,ic)
52    c$$$            print*,'#BAD (3)',bb
53    c$$$            ss=nsatstrips(ic)
54    c$$$            print*,'#saturated ',ss
55    c$$$         endif
56    c$$$      enddo
57          
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
59  *     STEP 1  *     STEP 1
60  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 44  c      common/dbg/DEBUG Line 74  c      common/dbg/DEBUG
74  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
75  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
76    
77  c      iflag=0  
78        call cl_to_couples(iflag)        call cl_to_couples(iflag)
79        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
80           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
81        endif        endif
82                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     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  
           
   
   
83  *-----------------------------------------------------  *-----------------------------------------------------
84  *-----------------------------------------------------  *-----------------------------------------------------
85  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 91  c$$$         endif Line 108  c$$$         endif
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    
111  c      iflag=0  
112        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
113        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
114           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
115        endif        endif
116                
117                
# Line 123  c      iflag=0 Line 140  c      iflag=0
140  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
141  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
142  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
143    *     count number of hit planes
144          planehit=0                
145          do np=1,nplanes          
146            if(ncp_plane(np).ne.0)then  
147              planehit=planehit+1  
148            endif                  
149          enddo                    
150          if(planehit.lt.3) goto 880 ! exit              
151          
152          nptxz_min=x_min_start              
153          nplxz_min=x_min_start              
154                
155          nptyz_min=y_min_start              
156          nplyz_min=y_min_start              
157    
158          cutdistyz=cutystart      
159          cutdistxz=cutxstart      
160    
161  c      iflag=0   878  continue
162        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
163        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
164           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
165        endif        endif
166  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
167            if(cutdistyz.lt.maxcuty/2)then
168              cutdistyz=cutdistyz+cutystep
169            else
170              cutdistyz=cutdistyz+(3*cutystep)
171            endif
172            goto 878                
173          endif                    
174    
175          if(planehit.eq.3) goto 881          
176          
177     879  continue  
178        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
179        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
180           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
181        endif        endif
182                                    
183          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
184            cutdistxz=cutdistxz+cutxstep
185            goto 879                
186          endif                    
187    
188        
189     881  continue  
190    *     if there is at least three planes on the Y view decreases cuts on X view
191          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
192         $     nplxz_min.ne.y_min_start)then
193            nptxz_min=x_min_step    
194            nplxz_min=x_min_start-x_min_step              
195            goto 879                
196          endif                    
197            
198   880  return   880  return
199        end        end
200    
# Line 144  c      iflag=0 Line 204  c      iflag=0
204        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
205    
206        include 'commontracker.f'        include 'commontracker.f'
207          include 'level1.f'
208        include 'common_momanhough.f'        include 'common_momanhough.f'
209        include 'common_mech.f'        include 'common_mech.f'
210        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
211        include 'common_mini_2.f'        include 'common_mini_2.f'
212        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
213        include 'level2.f'        include 'level2.f'
214    
215        include 'momanhough_init.f'  c      include 'momanhough_init.f'
216                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
217        logical FIMAGE            !        logical FIMAGE            !
218          real*8 AL_GUESS(5)
219    
220  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
221  *     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 256  c         iflag=0
256           ibest=0                !best track among candidates           ibest=0                !best track among candidates
257           iimage=0               !track image           iimage=0               !track image
258  *     ------------- select the best track -------------  *     ------------- select the best track -------------
259           rchi2best=1000000000.  c$$$         rchi2best=1000000000.
260    c$$$         do i=1,ntracks
261    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
262    c$$$     $         RCHI2_STORE(i).gt.0)then
263    c$$$               ibest=i
264    c$$$               rchi2best=RCHI2_STORE(i)
265    c$$$            endif
266    c$$$         enddo
267    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
268    
269    *     -------------------------------------------------------
270    *     order track-candidates according to:
271    *     1st) decreasing n.points
272    *     2nd) increasing chi**2
273    *     -------------------------------------------------------
274             rchi2best=1000000000.
275             ndofbest=0            
276           do i=1,ntracks           do i=1,ntracks
277              if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0              
278       $         RCHI2_STORE(i).gt.0)then             do ii=1,nplanes    
279                 ndof=ndof        
280         $            +int(xgood_store(ii,i))
281         $            +int(ygood_store(ii,i))
282               enddo              
283               if(ndof.gt.ndofbest)then
284                 ibest=i
285                 rchi2best=RCHI2_STORE(i)
286                 ndofbest=ndof    
287               elseif(ndof.eq.ndofbest)then
288                 if(RCHI2_STORE(i).lt.rchi2best.and.
289         $            RCHI2_STORE(i).gt.0)then
290                 ibest=i                 ibest=i
291                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
292              endif                 ndofbest=ndof  
293           enddo               endif            
294               endif
295             enddo
296    
297    c$$$         rchi2best=1000000000.
298    c$$$         ndofbest=0             !(1)
299    c$$$         do i=1,ntracks
300    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
301    c$$$     $          RCHI2_STORE(i).gt.0)then
302    c$$$             ndof=0             !(1)
303    c$$$             do ii=1,nplanes    !(1)
304    c$$$               ndof=ndof        !(1)
305    c$$$     $              +int(xgood_store(ii,i)) !(1)
306    c$$$     $              +int(ygood_store(ii,i)) !(1)
307    c$$$             enddo              !(1)
308    c$$$             if(ndof.ge.ndofbest)then !(1)
309    c$$$               ibest=i
310    c$$$               rchi2best=RCHI2_STORE(i)
311    c$$$               ndofbest=ndof    !(1)
312    c$$$             endif              !(1)
313    c$$$           endif
314    c$$$         enddo
315    
316           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
317  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
318  *     The best track candidate (ibest) is selected and a new fitting is performed.  *     The best track candidate (ibest) is selected and a new fitting is performed.
# Line 236  c         iflag=0 Line 343  c         iflag=0
343  *     **********************************************************  *     **********************************************************
344  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
345  *     **********************************************************  *     **********************************************************
346             call guess()
347             do i=1,5
348                AL_GUESS(i)=AL(i)
349             enddo
350    c         print*,'## guess: ',al
351    
352           do i=1,5           do i=1,5
353              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
354           enddo           enddo
355            
356           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
357           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
358           jstep=0                !# minimization steps           jstep=0                !# minimization steps
359    
360           call mini_2(jstep,ifail)           iprint=0
361    c         if(DEBUG)iprint=1
362             if(VERBOSE)iprint=1
363             if(DEBUG)iprint=2
364             call mini2(jstep,ifail,iprint)
365           if(ifail.ne.0) then           if(ifail.ne.0) then
366              if(DEBUG)then              if(VERBOSE)then
367                 print *,                 print *,
368       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
369       $              ,iev       $              ,iev
370    
371    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
372    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
373    c$$$               print*,'result:   ',(al(i),i=1,5)
374    c$$$               print*,'xgood ',xgood
375    c$$$               print*,'ygood ',ygood
376    c$$$               print*,'----------------------------------------------'
377              endif              endif
378              chi2=-chi2  c            chi2=-chi2
379           endif           endif
380                    
381           if(DEBUG)then           if(DEBUG)then
# Line 361  c     $        rchi2best.lt.15..and. Line 486  c     $        rchi2best.lt.15..and.
486        end        end
487    
488    
   
   
 c$$$************************************************************  
 c$$$  
 c$$$      subroutine readmipparam  
 c$$$              
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('trk-LADDER',i1,'-mip.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READMIPPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do iv=1,nviews  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            mip(int(pip),ilad)  
 c$$$c            print*,ilad,iv,pip,mip(int(pip),ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readchargeparam  
 c$$$        
 c$$$        
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('charge-l',i1,'.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READCHARGEPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do ip=1,nplanes  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)          
 c$$$c            print*,ilad,ip,pip,kch(ip,ilad),  
 c$$$c     $           cch(ip,ilad),sch(ip,ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readetaparam  
 c$$$*     -----------------------------------------  
 c$$$*     read eta2,3,4 calibration parameters  
 c$$$*     and fill variables:  
 c$$$*  
 c$$$*     eta2(netabin,nladders_view,nviews)  
 c$$$*     eta3(2*netabin,nladders_view,nviews)  
 c$$$*     eta4(2*netabin,nladders_view,nviews)  
 c$$$*  
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*40 fname_binning  
 c$$$      character*40 fname_param  
 c$$$c      character*120 cmd1  
 c$$$c      character*120 cmd2  
 c$$$  
 c$$$  
 c$$$******retrieve ANGULAR BINNING info  
 c$$$      fname_binning='binning.dat'  
 c$$$      print *,'Opening file: ',fname_binning  
 c$$$      open(10,  
 c$$$     $     FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning))  
 c$$$     $     ,STATUS='UNKNOWN'  
 c$$$     $     ,IOSTAT=iostat  
 c$$$     $     )  
 c$$$      if(iostat.ne.0)then  
 c$$$         print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$         return  
 c$$$      endif  
 c$$$      print*,'---- ANGULAR BINNING ----'  
 c$$$      print*,'Bin   -   angL   -   angR'  
 c$$$ 101  format(i2,'       ',f6.2,'     ',f6.2)  
 c$$$      do ibin=1,nangmax  
 c$$$         read(10,*  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )xnn,angL(ibin),angR(ibin)  
 c$$$         if(iostat.ne.0)goto 1000  
 c$$$         write(*,101)int(xnn),angL(ibin),angR(ibin)  
 c$$$      enddo          
 c$$$ 1000 nangbin=int(xnn)  
 c$$$      close(10)  
 c$$$      print*,'-------------------------'  
 c$$$        
 c$$$  
 c$$$  
 c$$$      do ieta=2,4               !loop on eta 2,3,4          
 c$$$******retrieve correction parameters  
 c$$$ 200     format(' Opening eta',i1,' files...')  
 c$$$         write(*,200)ieta  
 c$$$  
 c$$$ 201     format('eta',i1,'-bin',i1,'-l',i1,'.dat')  
 c$$$ 202     format('eta',i1,'-bin',i2,'-l',i1,'.dat')  
 c$$$         do iang=1,nangbin  
 c$$$            do ilad=1,nladders_view  
 c$$$               if(iang.lt.10)write(fname_param,201)ieta,iang,ilad  
 c$$$               if(iang.ge.10)write(fname_param,202)ieta,iang,ilad  
 c$$$c               print *,'Opening file: ',fname_param  
 c$$$               open(10,  
 c$$$     $             FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $              ,STATUS='UNKNOWN'  
 c$$$     $              ,IOSTAT=iostat  
 c$$$     $              )  
 c$$$               if(iostat.ne.0)then  
 c$$$                  print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$                  return  
 c$$$               endif  
 c$$$               do ival=1,netavalmax  
 c$$$                  if(ieta.eq.2)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta2(ival,iang),  
 c$$$     $                 (feta2(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.3)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta3(ival,iang),  
 c$$$     $                 (feta3(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.4)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta4(ival,iang),  
 c$$$     $                 (feta4(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(iostat.ne.0)then  
 c$$$                     netaval=ival-1  
 c$$$c$$$                     if(eta2(1,iang).ne.-eta2(netaval,iang))  
 c$$$c$$$     $                    print*,'**** ERROR on parameters !!! ****'  
 c$$$                     goto 2000  
 c$$$                  endif  
 c$$$               enddo  
 c$$$ 2000          close(10)  
 c$$$*               print*,'... done'  
 c$$$            enddo  
 c$$$         enddo  
 c$$$  
 c$$$      enddo                     !end loop on eta 2,3,4  
 c$$$  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
489                
490  ************************************************************  ************************************************************
491  ************************************************************  ************************************************************
# Line 557  c$$$ Line 510  c$$$
510  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
511  *     angx   - Projected angle in x  *     angx   - Projected angle in x
512  *     angy   - Projected angle in y  *     angy   - Projected angle in y
513    *     bfx    - x-component of magnetci field
514    *     bfy    - y-component of magnetic field
515  *  *
516  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
517  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 595  c$$$ Line 550  c$$$
550  *  *
551  *  *
552    
553        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
554    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c     01/02/2006 modified by elena vannuccini --> (2)  
 c     02/02/2006 modified by Elena Vannuccini --> (3)  
 c                (implemented new p.f.a.)  
 c     03/02/2006 modified by Elena Vannuccini --> (4)  
 c                (implemented variable resolution)  
 c*****************************************************  
555                
556        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
557        include 'level1.f'        include 'level1.f'
558          include 'calib.f'
559        include 'common_align.f'        include 'common_align.f'
560        include 'common_mech.f'        include 'common_mech.f'
561        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
       include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
562    
563        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
564        integer sensor        integer sensor
565        integer viewx,viewy        integer viewx,viewy
566        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
567        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
568          real angx,angy            !X-Y effective angle
569          real bfx,bfy              !X-Y b-field components
570    
571        real stripx,stripy        real stripx,stripy
572    
573        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
574        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
575        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
 c      double precision xi,yi,zi  
 c      double precision xi_A,yi_A,zi_A  
 c      double precision xi_B,yi_B,zi_B  
576                
577    
578        parameter (ndivx=30)        parameter (ndivx=30)
# Line 647  c      double precision xi_B,yi_B,zi_B Line 589  c      double precision xi_B,yi_B,zi_B
589        xPAM_B = 0.        xPAM_B = 0.
590        yPAM_B = 0.        yPAM_B = 0.
591        zPAM_B = 0.        zPAM_B = 0.
592    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
593    
594  *     -----------------  *     -----------------
595  *     CLUSTER X  *     CLUSTER X
596  *     -----------------  *     -----------------
597    
598        if(icx.ne.0)then        if(icx.ne.0)then
599           viewx = VIEW(icx)  
600           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
601           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
602           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
603                     resxPAM = RESXAV
604           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
605           if(PFAx.eq.'COG1')then  !(1)  *        --------------------------
606              stripx = stripx      !(1)  *        magnetic-field corrections
607              resxPAM = resxPAM    !(1)  *        --------------------------
608             angtemp  = ax
609             bfytemp  = bfy
610             if(nplx.eq.6) angtemp = -1. * ax
611             if(nplx.eq.6) bfytemp = -1. * bfy
612             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
613             angx     = 180.*atan(tgtemp)/acos(-1.)
614             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
615    c$$$         print*,nplx,ax,bfy/10.
616    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
617    c$$$         print*,'========================'
618    *        --------------------------
619    
620    c$$$         print*,'--- x-cl ---'
621    c$$$         istart = INDSTART(ICX)
622    c$$$         istop  = TOTCLLENGTH
623    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
624    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
625    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
626    c$$$         print*,INDMAX(icx)
627    c$$$         print*,cog(4,icx)
628    c$$$         print*,fbad_cog(4,icx)
629            
630    
631             if(PFAx.eq.'COG1')then
632    
633                stripx  = stripx
634                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
635    
636           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
637              stripx = stripx + cog(2,icx)              
638                stripx  = stripx + cog(2,icx)            
639                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
640              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
641    
642             elseif(PFAx.eq.'COG3')then
643    
644                stripx  = stripx + cog(3,icx)            
645                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
646                resxPAM = resxPAM*fbad_cog(3,icx)
647    
648             elseif(PFAx.eq.'COG4')then
649    
650                stripx  = stripx + cog(4,icx)            
651                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
652                resxPAM = resxPAM*fbad_cog(4,icx)
653    
654           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
655  c            cog2 = cog(2,icx)  
656  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)                          stripx  = stripx + pfaeta2(icx,angx)          
657  c            stripx = stripx + etacorr              resxPAM = risx_eta2(abs(angx))
658              stripx = stripx + pfa_eta2(icx,angx)            !(3)              resxPAM = resxPAM*fbad_cog(2,icx)
             resxPAM = risx_eta2(angx)                       !   (4)  
659              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
660       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
661              resxPAM = resxPAM*fbad_cog(2,icx)  
662           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
663              stripx = stripx + pfa_eta3(icx,angx)            !(3)  
664              resxPAM = risx_eta3(angx)                       !   (4)              stripx  = stripx + pfaeta3(icx,angx)          
665              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              resxPAM = risx_eta3(abs(angx))                      
666       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
667              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)            
668           elseif(PFAx.eq.'ETA4')then                         !(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
669              stripx = stripx + pfa_eta4(icx,angx)            !(3)  
670              resxPAM = risx_eta4(angx)                       !   (4)           elseif(PFAx.eq.'ETA4')then                        
671              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
672       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)              stripx  = stripx + pfaeta4(icx,angx)            
673              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = risx_eta4(abs(angx))                      
674           elseif(PFAx.eq.'ETA')then                          !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
675              stripx = stripx + pfa_eta(icx,angx)             !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              
676              resxPAM = ris_eta(icx,angx)                     !   (4)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
677              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
678       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)           elseif(PFAx.eq.'ETA')then  
679  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
680              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              stripx  = stripx + pfaeta(icx,angx)            
681           elseif(PFAx.eq.'COG')then           !(2)              resxPAM = ris_eta(icx,angx)                    
682              stripx = stripx + cog(0,icx)     !(2)                      resxPAM = resxPAM*fbad_eta(icx,angx)            
683              resxPAM = risx_cog(angx)                        !   (4)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              
684              resxPAM = resxPAM*fbad_cog(0,icx)!(2)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
685    
686             elseif(PFAx.eq.'COG')then          
687    
688                stripx  = stripx + cog(0,icx)            
689                resxPAM = risx_cog(abs(angx))                    
690                resxPAM = resxPAM*fbad_cog(0,icx)
691    
692           else           else
693              print*,'*** Non valid p.f.a. (x) --> ',PFAx              print*,'*** Non valid p.f.a. (x) --> ',PFAx
694           endif           endif
695    
696    
697    *     ======================================
698    *     temporary patch for saturated clusters
699    *     ======================================
700             if( nsatstrips(icx).gt.0 )then
701                stripx  = stripx + cog(4,icx)            
702                resxPAM = pitchX*1e-4/sqrt(12.)
703                cc=cog(4,icx)
704    c$$$            print*,icx,' *** ',cc
705    c$$$            print*,icx,' *** ',resxPAM
706             endif
707    
708        endif        endif
 c      if(icy.eq.0.and.icx.ne.0)  
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
709                
710  *     -----------------  *     -----------------
711  *     CLUSTER Y  *     CLUSTER Y
712  *     -----------------  *     -----------------
713    
714        if(icy.ne.0)then        if(icy.ne.0)then
715    
716           viewy = VIEW(icy)           viewy = VIEW(icy)
717           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
718           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
719           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
720             stripy = float(MAXS(icy))
721    
722           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
723              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
724       $           ,icx,icy       $           ,icx,icy
725              goto 100              goto 100
726           endif           endif
727            *        --------------------------
728           stripy = float(MAXS(icy))  *        magnetic-field corrections
729           if(PFAy.eq.'COG1')then !(1)  *        --------------------------
730              stripy = stripy     !(1)           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
731              resyPAM = resyPAM   !(1)           angy    = 180.*atan(tgtemp)/acos(-1.)
732             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
733    *        --------------------------
734            
735    c$$$         print*,'--- y-cl ---'
736    c$$$         istart = INDSTART(ICY)
737    c$$$         istop  = TOTCLLENGTH
738    c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
739    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
740    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
741    c$$$         print*,INDMAX(icy)
742    c$$$         print*,cog(4,icy)
743    c$$$         print*,fbad_cog(4,icy)
744    
745             if(PFAy.eq.'COG1')then
746    
747                stripy  = stripy    
748                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
749    
750           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
751              stripy = stripy + cog(2,icy)  
752                stripy  = stripy + cog(2,icy)
753                resyPAM = risy_cog(abs(angy))!TEMPORANEO
754              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
755    
756             elseif(PFAy.eq.'COG3')then
757    
758                stripy  = stripy + cog(3,icy)
759                resyPAM = risy_cog(abs(angy))!TEMPORANEO
760                resyPAM = resyPAM*fbad_cog(3,icy)
761    
762             elseif(PFAy.eq.'COG4')then
763    
764                stripy  = stripy + cog(4,icy)
765                resyPAM = risy_cog(abs(angy))!TEMPORANEO
766                resyPAM = resyPAM*fbad_cog(4,icy)
767    
768           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
769  c            cog2 = cog(2,icy)  
770  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
771  c            stripy = stripy + etacorr              resyPAM = risy_eta2(abs(angy))              
             stripy = stripy + pfa_eta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
772              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
773              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
774       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
775           elseif(PFAy.eq.'ETA3')then                         !(3)  
776              stripy = stripy + pfa_eta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
777              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
778              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
779       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
780           elseif(PFAy.eq.'ETA4')then                         !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)
781              stripy = stripy + pfa_eta4(icy,angy)            !(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
782              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
783              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
784       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
785           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
786              stripy = stripy + pfa_eta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
787              resyPAM = ris_eta(icy,angy)                     !   (4)              if(DEBUG.and.fbad_cog(4,icy).ne.1)
788  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
789              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
790              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
791       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
792                stripy  = stripy + pfaeta(icy,angy)
793                resyPAM = ris_eta(icy,angy)  
794                resyPAM = resyPAM*fbad_eta(icy,angy)
795                if(DEBUG.and.fbad_cog(2,icy).ne.1)
796         $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
797    
798           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
799              stripy = stripy + cog(0,icy)              
800              resyPAM = risy_cog(angy)                        !   (4)              stripy  = stripy + cog(0,icy)            
801  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
802              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
803    
804           else           else
805              print*,'*** Non valid p.f.a. (x) --> ',PFAx              print*,'*** Non valid p.f.a. (x) --> ',PFAx
806           endif           endif
807    
808    
809    *     ======================================
810    *     temporary patch for saturated clusters
811    *     ======================================
812             if( nsatstrips(icy).gt.0 )then
813                stripy  = stripy + cog(4,icy)            
814                resyPAM = pitchY*1e-4/sqrt(12.)
815                cc=cog(4,icy)
816    c$$$            print*,icy,' *** ',cc
817    c$$$            print*,icy,' *** ',resyPAM
818             endif
819    
820    
821        endif        endif
822    
823          c      print*,'## stripx,stripy ',stripx,stripy
824    
825  c===========================================================  c===========================================================
826  C     COUPLE  C     COUPLE
827  C===========================================================  C===========================================================
# Line 968  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1023  c         print*,'A-(',xPAM_A,yPAM_A,')
1023                            
1024        endif        endif
1025                    
1026    
1027    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1028    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1029    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1030    
1031   100  continue   100  continue
1032        end        end
1033    
# Line 1046  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1106  c         print*,'A-(',xPAM_A,yPAM_A,')
1106           endif                   endif        
1107    
1108           distance=           distance=
1109       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1110    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1111           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1112    
1113  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
# Line 1071  c$$$         print*,' resolution ',re Line 1132  c$$$         print*,' resolution ',re
1132  *     ----------------------  *     ----------------------
1133                    
1134           distance=           distance=
1135       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1136       $        +       $       +
1137       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1138    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1139    c$$$     $        +
1140    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1141           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1142    
1143  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1284  c     $              ,iv,xvv(iv),yvv(iv) Line 1348  c     $              ,iv,xvv(iv),yvv(iv)
1348  *     it returns the plane number  *     it returns the plane number
1349  *      *    
1350        include 'commontracker.f'        include 'commontracker.f'
1351          include 'level1.f'
1352  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1353        include 'common_momanhough.f'        include 'common_momanhough.f'
1354                
# Line 1321  c      include 'common_analysis.f' Line 1386  c      include 'common_analysis.f'
1386  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1387  *      *    
1388        include 'commontracker.f'        include 'commontracker.f'
1389          include 'level1.f'
1390  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1391        include 'common_momanhough.f'        include 'common_momanhough.f'
1392                
# Line 1349  c      include 'common_analysis.f' Line 1415  c      include 'common_analysis.f'
1415  *     positive if sensor =2  *     positive if sensor =2
1416  *  *
1417        include 'commontracker.f'        include 'commontracker.f'
1418          include 'level1.f'
1419  c      include 'calib.f'  c      include 'calib.f'
1420  c      include 'level1.f'  c      include 'level1.f'
1421  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1378  c      include 'common_analysis.f' Line 1445  c      include 'common_analysis.f'
1445  *************************************************************************  *************************************************************************
1446  *************************************************************************  *************************************************************************
1447  *************************************************************************  *************************************************************************
 c$$$      subroutine book_debug  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$  
 c$$$      character*35 block1,block2,block3!,block4  
 c$$$     $     ,block5!,block6        
 c$$$        
 c$$$* * * * * * * * * * * * * * * * * * * * * * * *  
 c$$$*     HOUGH TRANSFORM PARAMETERS  
 c$$$        
 c$$$      call HBOOK2(1003  
 c$$$     $     ,'y vs tg thyz'  
 c$$$     $     ,300,-1.,1.         !x          
 c$$$     $     ,3000,-70.,70.,0.)   !y  
 c$$$  
 c$$$      call HBOOK1(1004  
 c$$$     $     ,'Dy'  
 c$$$     $     ,100,0.,2.,0.)    
 c$$$  
 c$$$      call HBOOK1(1005  
 c$$$     $     ,'D thyz'  
 c$$$     $     ,100,0.,.05,0.)    
 c$$$  
 c$$$  
 c$$$  
 c$$$*     DEBUG ntuple:  
 c$$$      call HBNT(ntp_level2+1,'LEVEL2',' ')  
 c$$$        
 c$$$      call HBNAME(ntp_level2+1,'EVENT',good2_nt,  
 c$$$     $     'GOOD2:L,NEV2:I')  
 c$$$  
 c$$$ 411  format('NDBLT:I::[0,',I5,']')  
 c$$$      write(block1,411) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,  
 c$$$     $     block1//'  
 c$$$     $     ,ALFAYZ1(NDBLT):R  
 c$$$     $     ,ALFAYZ2(NDBLT):R  
 c$$$     $     ,DB_CLOUD(NDBLT):I  
 c$$$     $     ')    
 c$$$  
 c$$$ 412  format('NTRPT:I::[0,',I5,']')  
 c$$$      write(block2,412) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,  
 c$$$     $     block2//'  
 c$$$     $     ,ALFAXZ1(NTRPT):R  
 c$$$     $     ,ALFAXZ2(NTRPT):R    
 c$$$     $     ,ALFAXZ3(NTRPT):R  
 c$$$     $     ,TR_CLOUD(NTRPT):I  
 c$$$     $     ')  
 c$$$        
 c$$$    
 c$$$ 413  format('NCLOUDS_YZ:I::[0,',I4,']')  
 c$$$c$$$ 414  format('DB_CLOUD(',I4,'):I')  
 c$$$      write(block3,413) ncloyz_max  
 c$$$c$$$      write(block4,414) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,  
 c$$$     $     block3//'  
 c$$$     $     ,ALFAYZ1_AV(NCLOUDS_YZ):R  
 c$$$     $     ,ALFAYZ2_AV(NCLOUDS_YZ):R  
 c$$$     $     ,PTCLOUD_YZ(NCLOUDS_YZ):I'  
 c$$$c$$$     $     ,'//block4  
 c$$$     $     )  
 c$$$  
 c$$$ 415  format('NCLOUDS_XZ:I::[0,',I4,']')  
 c$$$c$$$ 416  format('TR_CLOUD(',I5,'):I')  
 c$$$      write(block5,415) ncloxz_max  
 c$$$c$$$      write(block6,416) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,  
 c$$$     $     block5//'  
 c$$$     $     ,ALFAXZ1_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ2_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ3_AV(NCLOUDS_XZ):R  
 c$$$     $     ,PTCLOUD_XZ(NCLOUDS_XZ):I'  
 c$$$c$$$     $     ,'//block6  
 c$$$     $     )  
 c$$$  
 c$$$        
 c$$$      return  
 c$$$      end  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 c$$$      subroutine book_level2  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2)  
 c$$$c*****************************************************  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      character*35 block1,block2  
 c$$$  
 c$$$c      print*,'__________ booking LEVEL2 n-tuple __________'  
 c$$$  
 c$$$*     LEVEL1 ntuple:  
 c$$$      call HBNT(ntp_level2,'LEVEL2',' ')  
 c$$$        
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I')  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$c     $     ,WHIC_CALIB:I,SWCODE:I')  
 c$$$      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$     $     ,WHICH_CALIB:I,SWCODE:I,CRC(12):L')  
 c$$$c*********************************************************  
 c$$$    
 c$$$  
 c$$$c$$$# ifndef TEST2003  
 c$$$c$$$  
 c$$$c$$$      call HBNAME(ntp_level2,'CPU',pkt_type  
 c$$$c$$$     $     ,'PKT_TYPE:I::[0,50]  
 c$$$c$$$     $     ,PKT_NUM:I  
 c$$$c$$$     $     ,OBT:I'//  
 c$$$c$$$c********************************************************  
 c$$$c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c$$$c     $     ,WHICH_CALIB:I::[0,50]')  
 c$$$c$$$     $     ',CPU_CRC:L')  
 c$$$c$$$c********************************************************  
 c$$$c$$$  
 c$$$c$$$# endif  
 c$$$  
 c$$$ 417  format('NTRK:I::[0,',I4,']')  
 c$$$ 418  format(',IMAGE(NTRK):I::[0,',I4,']')  
 c$$$      write(block1,417)NTRKMAX  
 c$$$      write(block2,418)NTRKMAX  
 c$$$      call HBNAME(ntp_level2,'TRACKS',NTRK,  
 c$$$     $     block1//  
 c$$$     $     block2//'  
 c$$$     $     ,XM(6,NTRK):R  
 c$$$     $     ,YM(6,NTRK):R  
 c$$$     $     ,ZM(6,NTRK):R  
 c$$$     $     ,RESX(6,NTRK):R  
 c$$$     $     ,RESY(6,NTRK):R  
 c$$$     $     ,AL(5,NTRK):R  
 c$$$     $     ,COVAL(5,5,NTRK):R  
 c$$$     $     ,CHI2(NTRK):R  
 c$$$     $     ,XGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,YGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,XV(6,NTRK):R  
 c$$$     $     ,YV(6,NTRK):R  
 c$$$     $     ,ZV(6,NTRK):R  
 c$$$     $     ,AXV(6,NTRK):R  
 c$$$     $     ,AYV(6,NTRK):R'//  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     ,DEDXP(6,NTRK):R'//  
 c$$$c     $     ')  
 c$$$     $     ',DEDX_X(6,NTRK):R  
 c$$$     $     ,DEDX_Y(6,NTRK):R'//  
 c$$$c****************************************************  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c     $     ,CRC(12):L  
 c$$$     $     ',BdL(NTRK):R'  
 c$$$     $     )  
 c$$$c****************************************************  
 c$$$  
 c$$$    
 c$$$      call HBNAME(ntp_level2,'SINGLETX',nclsx,  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     'NCLSX(6):I,NCLSY(6):I')  
 c$$$     $     'NCLSX:I::[0,500],PLANEX(NCLSX):I  
 c$$$     $     ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2)  
 c$$$c    $     ,XS(NCLSX):R,SGNLXS(NCLSX):R')   !(2)  
 c$$$      call HBNAME(ntp_level2,'SINGLETY',nclsy,  
 c$$$     $     'NCLSY:I::[0,500],PLANEY(NCLSY):I  
 c$$$     $     ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2)  
 c$$$c    $     ,YS(NCLSY):R,SGNLYS(NCLSY):R')   !(2)  
 c$$$      return  
 c$$$      end  
   
 c$$$      subroutine fill_level2_clouds  
 c$$$c*****************************************************  
 c$$$c     29/11/2005 created by elena vannuccini  
 c$$$c*****************************************************  
 c$$$  
 c$$$*     -------------------------------------------------------  
 c$$$*     This routine fills the  variables related to the hough  
 c$$$*     transform, for the debig n-tuple  
 c$$$*     -------------------------------------------------------  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      good2_nt=.true.!good2  
 c$$$c      nev2_nt=nev2  
 c$$$        
 c$$$      if(.false.  
 c$$$     $     .or.ntrpt.gt.ntrpt_max_nt  
 c$$$     $     .or.ndblt.gt.ndblt_max_nt  
 c$$$     $     .or.NCLOUDS_XZ.gt.ncloxz_max  
 c$$$     $     .or.NCLOUDS_yZ.gt.ncloyz_max  
 c$$$     $     )then  
 c$$$         good2_nt=.false.  
 c$$$         ntrpt_nt=0  
 c$$$         ndblt_nt=0  
 c$$$         NCLOUDS_XZ_nt=0  
 c$$$         NCLOUDS_YZ_nt=0  
 c$$$      else  
 c$$$         ndblt_nt=ndblt  
 c$$$         ntrpt_nt=ntrpt  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do id=1,ndblt  
 c$$$               alfayz1_nt(id)=alfayz1(id) !Y0  
 c$$$               alfayz2_nt(id)=alfayz2(id) !tg theta-yz  
 c$$$c               db_cloud_nt(id)=db_cloud(id)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do it=1,ntrpt  
 c$$$               alfaxz1_nt(it)=alfaxz1(it) !X0  
 c$$$               alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz  
 c$$$               alfaxz3_nt(it)=alfaxz3(it) !1/r  
 c$$$c               tr_cloud_nt(it)=tr_cloud(it)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         nclouds_yz_nt=nclouds_yz  
 c$$$         nclouds_xz_nt=nclouds_xz  
 c$$$         if(nclouds_yz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do iyz=1,nclouds_yz  
 c$$$               ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)  
 c$$$               alfayz1_av_nt(iyz)=alfayz1_av(iyz)  
 c$$$               alfayz2_av_nt(iyz)=alfayz2_av(iyz)  
 c$$$               nnn=nnn+ptcloud_yz(iyz)  
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               db_cloud_nt(ipt)=db_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_yz '  
 c$$$c     $           ,(ptcloud_yz(i),i=1,nclouds_yz)  
 c$$$c            print*,'#### ntupla #### db_cloud '  
 c$$$c     $           ,(db_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$         if(nclouds_xz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do ixz=1,nclouds_xz  
 c$$$               ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)  
 c$$$               alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)  
 c$$$               alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)  
 c$$$               alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)  
 c$$$               nnn=nnn+ptcloud_xz(ixz)                
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               tr_cloud_nt(ipt)=tr_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_xz '  
 c$$$c     $           ,(ptcloud_xz(i),i=1,nclouds_xz)  
 c$$$c            print*,'#### ntupla #### tr_cloud '  
 c$$$c     $           ,(tr_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$      endif  
 c$$$      end  
1448                
1449    
1450  ***************************************************  ***************************************************
# Line 1660  c$$$      end Line 1459  c$$$      end
1459        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1460    
1461        include 'commontracker.f'        include 'commontracker.f'
1462          include 'level1.f'
1463        include 'common_momanhough.f'        include 'common_momanhough.f'
1464        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1465        include 'calib.f'        include 'calib.f'
1466        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1467    
1468  *     output flag  *     output flag
1469  *     --------------  *     --------------
# Line 1691  c      common/dbg/DEBUG Line 1488  c      common/dbg/DEBUG
1488           ncls(ip)=0           ncls(ip)=0
1489        enddo        enddo
1490        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1491           cl_single(icl)=1           cl_single(icl) = 1
1492           cl_good(icl)=0           cl_good(icl)   = 0
1493          enddo
1494          do iv=1,nviews
1495             ncl_view(iv)  = 0
1496             mask_view(iv) = 0      !all included
1497        enddo        enddo
1498                
1499    *     count number of cluster per view
1500          do icl=1,nclstr1
1501             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1502          enddo
1503    *     mask views with too many clusters
1504          do iv=1,nviews
1505             if( ncl_view(iv).gt. nclusterlimit)then
1506                mask_view(iv) = 1
1507                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1508         $           ,nclusterlimit,' on view ', iv,' --> masked!'
1509             endif
1510          enddo
1511    
1512    
1513  *     start association  *     start association
1514        ncouples=0        ncouples=0
1515        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1516           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1517                    
1518  *     ----------------------------------------------------  *     ----------------------------------------------------
1519    *     jump masked views (X VIEW)
1520    *     ----------------------------------------------------
1521             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1522    *     ----------------------------------------------------
1523  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1524  *     ----------------------------------------------------  *     ----------------------------------------------------
1525           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1526              cl_single(icx)=0              cl_single(icx)=0
1527              goto 10              goto 10
1528           endif           endif
# Line 1746  c     endif Line 1565  c     endif
1565              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1566                            
1567  *     ----------------------------------------------------  *     ----------------------------------------------------
1568    *     jump masked views (Y VIEW)
1569    *     ----------------------------------------------------
1570                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1571    
1572    *     ----------------------------------------------------
1573  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1574  *     ----------------------------------------------------  *     ----------------------------------------------------
1575              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1576                 cl_single(icy)=0                 cl_single(icy)=0
1577                 goto 20                 goto 20
1578              endif              endif
# Line 1794  c     endif Line 1618  c     endif
1618  *     charge correlation  *     charge correlation
1619  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1620    
1621  *     -------------------------------------------------------------                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
 *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<  
 *     -------------------------------------------------------------  
                if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)  
1622       $              .and.       $              .and.
1623       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1624       $              .and.       $              .and.
1625       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1626       $              .and.       $              .and.
1627       $              .true.)then       $              .true.)then
1628    
1629                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1630       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1631                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1632    
1633  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1634    
1635                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1636       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1637                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1638                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1820  c                  cut = chcut * sch(npl Line 1641  c                  cut = chcut * sch(npl
1641                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1642                    endif                    endif
1643                 endif                 endif
1644                  
1645  *     ------------------> COUPLE <------------------                 if(ncp_plane(nplx).gt.ncouplemax)then
 *     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
1646                    if(verbose)print*,                    if(verbose)print*,
1647       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1648       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1649       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1650       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1651  c     good2=.false.                    mask_view(nviewx(nplx)) = 2
1652  c     goto 880   !fill ntp and go to next event                                        mask_view(nviewy(nply)) = 2
1653                    iflag=1                    goto 10
                   return  
1654                 endif                 endif
1655                                
1656    *     ------------------> COUPLE <------------------
1657                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1658                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1659                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1660                 cl_single(icx)=0                 cl_single(icx)=0
1661                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1662  *     ----------------------------------------------  *     ----------------------------------------------
1663    
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     init variables  
       ncp_tot=0  
       do ip=1,nplanes  
          do ico=1,ncouplemax  
             clx(ip,ico)=0  
             cly(ip,ico)=0  
          enddo  
          ncp_plane(ip)=0  
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
 *     start association  
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
 *     ----------------------------------------------------  
 *     cut on charge (X VIEW)  
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     cut BAD (X VIEW)              
          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
          ifirst=INDSTART(icx)  
          if(icx.ne.nclstr1) then  
             ilast=INDSTART(icx+1)-1  
          else  
             ilast=TOTCLLENGTH  
          endif  
          badcl=badseed  
          do igood=-ngoodstr,ngoodstr  
             ibad=1  
             if((INDMAX(icx)+igood).gt.ifirst.and.  
      $           (INDMAX(icx)+igood).lt.ilast.and.  
      $           .true.)then  
                ibad=BAD(VIEW(icx),  
      $              nvk(MAXS(icx)+igood),  
      $              nst(MAXS(icx)+igood))  
             endif  
             badcl=badcl*ibad  
          enddo  
          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
             goto 10             !<<<<<<<<<<<<<< BAD cut  
          endif                  !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
           
          cl_good(icx)=1  
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
 *     ----------------------------------------------------  
 *     cut on charge (Y VIEW)  
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     cut BAD (Y VIEW)              
             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
             ifirst=INDSTART(icy)  
             if(icy.ne.nclstr1) then  
                ilast=INDSTART(icy+1)-1  
             else  
                ilast=TOTCLLENGTH  
             endif  
             badcl=badseed  
             do igood=-ngoodstr,ngoodstr  
                ibad=1  
                if((INDMAX(icy)+igood).gt.ifirst.and.  
      $              (INDMAX(icy)+igood).lt.ilast.and.  
      $              .true.)  
      $              ibad=BAD(VIEW(icy),  
      $              nvk(MAXS(icy)+igood),  
      $              nst(MAXS(icy)+igood))  
                badcl=badcl*ibad  
             enddo  
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
               
               
             cl_good(icy)=1                    
             nply=npl(VIEW(icy))  
             nldy=nld(MAXS(icy),VIEW(icy))  
               
 *     ----------------------------------------------  
 *     CONDITION TO FORM A COUPLE  
 *     ----------------------------------------------  
 *     geometrical consistency (same plane and ladder)  
             if(nply.eq.nplx.and.nldy.eq.nldx)then  
 *     charge correlation  
 *     ===========================================================  
 *     this version of the subroutine is used for the calibration  
 *     thus charge-correlation selection is obviously removed  
 *     ===========================================================  
 c$$$               ddd=(dedx(icy)  
 c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 *     ===========================================================  
                 
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
                 
                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  
1664              endif                                            endif                              
 *     ----------------------------------------------  
1665    
1666   20         continue   20         continue
1667           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
# Line 2099  c     goto 880   !fill ntp and go to nex Line 1687  c     goto 880   !fill ntp and go to nex
1687        endif        endif
1688                
1689        do ip=1,6        do ip=1,6
1690           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1691        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1692                
1693        return        return
1694        end        end
   
1695                
1696  ***************************************************  ***************************************************
1697  *                                                 *  *                                                 *
# Line 2128  c     goto 880       !fill ntp and go to Line 1703  c     goto 880       !fill ntp and go to
1703  **************************************************  **************************************************
1704    
1705        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1706    
1707        include 'commontracker.f'        include 'commontracker.f'
1708          include 'level1.f'
1709        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1710        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1711        include 'common_mini_2.f'        include 'common_mini_2.f'
1712        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1713    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1714    
1715  *     output flag  *     output flag
1716  *     --------------  *     --------------
# Line 2174  c      double precision xm3,ym3,zm3 Line 1743  c      double precision xm3,ym3,zm3
1743  *     -----------------------------  *     -----------------------------
1744    
1745    
1746    *     --------------------------------------------
1747    *     put a limit to the maximum number of couples
1748    *     per plane, in order to apply hough transform
1749    *     (couples recovered during track refinement)
1750    *     --------------------------------------------
1751          do ip=1,nplanes
1752             if(ncp_plane(ip).gt.ncouplelimit)then
1753                mask_view(nviewx(ip)) = 8
1754                mask_view(nviewy(ip)) = 8
1755             endif
1756          enddo
1757    
1758    
1759        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1760        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1761                
1762        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1763           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1764                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1765             do is1=1,2             !loop on sensors - COPPIA 1            
1766              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1767                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1768                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1769  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1770                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1771                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1772                 xm1=xPAM                 xm1=xPAM
1773                 ym1=yPAM                 ym1=yPAM
1774                 zm1=zPAM                                   zm1=zPAM                  
1775  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1776    
1777                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1778                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1779         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1780                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1781                                            
1782                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2198  c     print*,'***',is1,xm1,ym1,zm1 Line 1784  c     print*,'***',is1,xm1,ym1,zm1
1784                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1785  c                        call xyz_PAM  c                        call xyz_PAM
1786  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1787    c                        call xyz_PAM
1788    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1789                          call xyz_PAM                          call xyz_PAM
1790       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1791                          xm2=xPAM                          xm2=xPAM
1792                          ym2=yPAM                          ym2=yPAM
1793                          zm2=zPAM                          zm2=zPAM
# Line 2215  c     $                       (icx2,icy2 Line 1803  c     $                       (icx2,icy2
1803       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1804  c                           good2=.false.  c                           good2=.false.
1805  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1806                               do iv=1,12
1807                                  mask_view(iv) = 3
1808                               enddo
1809                             iflag=1                             iflag=1
1810                             return                             return
1811                          endif                          endif
# Line 2248  c$$$ Line 1839  c$$$
1839  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1840    
1841    
1842                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1843    
1844                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1845                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1846         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1847                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1848                                                                
1849                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2257  c$$$ Line 1851  c$$$
1851                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
1852  c                                 call xyz_PAM  c                                 call xyz_PAM
1853  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
1854    c                                 call xyz_PAM
1855    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
1856                                   call xyz_PAM                                   call xyz_PAM
1857       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
1858         $                                ,0.,0.,0.,0.)
1859                                   xm3=xPAM                                   xm3=xPAM
1860                                   ym3=yPAM                                   ym3=yPAM
1861                                   zm3=zPAM                                   zm3=zPAM
# Line 2285  c     $                                 Line 1882  c     $                                
1882       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
1883  c                                    good2=.false.  c                                    good2=.false.
1884  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
1885                                        do iv=1,nviews
1886                                           mask_view(iv) = 4
1887                                        enddo
1888                                      iflag=1                                      iflag=1
1889                                      return                                      return
1890                                   endif                                   endif
# Line 2323  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1923  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1923                                endif                                endif
1924                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1925                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1926     30                     continue
1927                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1928   30                  continue   31                  continue
1929                                            
1930   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1931                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1932     20            continue
1933              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1934                            
1935           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1936        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1937     10   continue
1938        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1939                
1940        if(DEBUG)then        if(DEBUG)then
# Line 2359  c     goto 880               !ntp fill Line 1962  c     goto 880               !ntp fill
1962        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
1963    
1964        include 'commontracker.f'        include 'commontracker.f'
1965          include 'level1.f'
1966        include 'common_momanhough.f'        include 'common_momanhough.f'
1967        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1968    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1969    
1970  *     output flag  *     output flag
1971  *     --------------  *     --------------
# Line 2395  c      common/dbg/DEBUG Line 1997  c      common/dbg/DEBUG
1997        distance=0        distance=0
1998        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
1999        npt_tot=0        npt_tot=0
2000          nloop=0                  
2001     90   continue                  
2002        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2003           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
2004                            
# Line 2498  c     print*,'*   idbref,idb2 ',idbref,i Line 2102  c     print*,'*   idbref,idb2 ',idbref,i
2102              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2103           enddo           enddo
2104  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
2105           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
2106           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
2107           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2108                    
# Line 2512  c     print*,'>>>> ',ncpused,npt,nplused Line 2116  c     print*,'>>>> ',ncpused,npt,nplused
2116       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2117  c               good2=.false.  c               good2=.false.
2118  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2119                do iv=1,nviews
2120                   mask_view(iv) = 5
2121                enddo
2122              iflag=1              iflag=1
2123              return              return
2124           endif           endif
# Line 2549  c$$$     $           ,(db_cloud(iii),iii Line 2156  c$$$     $           ,(db_cloud(iii),iii
2156        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2157                
2158                
2159          if(nloop.lt.nstepy)then      
2160            cutdistyz = cutdistyz+cutystep
2161            nloop     = nloop+1          
2162            goto 90                
2163          endif                    
2164          
2165        if(DEBUG)then        if(DEBUG)then
2166           print*,'---------------------- '           print*,'---------------------- '
2167           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2575  c$$$     $           ,(db_cloud(iii),iii Line 2188  c$$$     $           ,(db_cloud(iii),iii
2188        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2189    
2190        include 'commontracker.f'        include 'commontracker.f'
2191          include 'level1.f'
2192        include 'common_momanhough.f'        include 'common_momanhough.f'
2193        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2194    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2195    
2196  *     output flag  *     output flag
2197  *     --------------  *     --------------
# Line 2610  c      common/dbg/DEBUG Line 2222  c      common/dbg/DEBUG
2222        distance=0        distance=0
2223        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2224        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2225          nloop=0                  
2226     91   continue                  
2227        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2228           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
2229  c     print*,'--------------'  c     print*,'--------------'
# Line 2711  c     print*,'check cp_used' Line 2325  c     print*,'check cp_used'
2325           do ip=1,nplanes           do ip=1,nplanes
2326              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2327           enddo           enddo
2328           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2329           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2330           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2331                    
2332  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2333  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
# Line 2724  c     print*,'check cp_used' Line 2338  c     print*,'check cp_used'
2338       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2339  c     good2=.false.  c     good2=.false.
2340  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2341                do iv=1,nviews
2342                   mask_view(iv) = 6
2343                enddo
2344              iflag=1              iflag=1
2345              return              return
2346           endif           endif
# Line 2759  c$$$     $           ,(tr_cloud(iii),iii Line 2376  c$$$     $           ,(tr_cloud(iii),iii
2376  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2377  22288    continue  22288    continue
2378        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2379          
2380           if(nloop.lt.nstepx)then      
2381             cutdistxz=cutdistxz+cutxstep
2382             nloop=nloop+1          
2383             goto 91                
2384           endif                    
2385          
2386        if(DEBUG)then        if(DEBUG)then
2387           print*,'---------------------- '           print*,'---------------------- '
2388           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2781  c$$$     $           ,(tr_cloud(iii),iii Line 2404  c$$$     $           ,(tr_cloud(iii),iii
2404  **************************************************  **************************************************
2405    
2406        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2407    
2408        include 'commontracker.f'        include 'commontracker.f'
2409          include 'level1.f'
2410        include 'common_momanhough.f'        include 'common_momanhough.f'
2411        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2412        include 'common_mini_2.f'        include 'common_mini_2.f'
2413        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
2414    
2415  c      logical DEBUG  
 c      common/dbg/DEBUG  
2416    
2417  *     output flag  *     output flag
2418  *     --------------  *     --------------
# Line 2809  c      common/dbg/DEBUG Line 2428  c      common/dbg/DEBUG
2428  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2429  *     list of matching couples in the combination  *     list of matching couples in the combination
2430  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2431        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2432        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2433  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2434        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2435  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2436  *     variables for track fitting  *     variables for track fitting
2437        double precision AL_INI(5)        double precision AL_INI(5)
2438        double precision tath  c      double precision tath
2439  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2440  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2441    
# Line 2882  c      real fitz(nplanes)        !z coor Line 2501  c      real fitz(nplanes)        !z coor
2501                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2502              enddo              enddo
2503                            
2504              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2505                if(nplused.lt.nplyz_min)goto 888 !next doublet
2506              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2507                            
2508              if(DEBUG)then              if(DEBUG)then
# Line 2909  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2529  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2529  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2530                            
2531  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2532              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2533              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2534              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2535       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2536              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2537              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2538              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2539                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2540  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2541              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2542                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2543  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2544  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2545                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2546  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2547  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2548                c$$$            ELSE
2549    c$$$               AL_INI(4) = acos(-1.)/2
2550    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2551    c$$$            ENDIF
2552    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2553    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2554    c$$$            
2555    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2556    c$$$            
2557    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2558                            
2559              if(DEBUG)then              if(DEBUG)then
2560                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2561                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 2976  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2606  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2606  *                                   *************************  *                                   *************************
2607  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2608  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2609    c                                    call xyz_PAM(icx,icy,is, !(1)
2610    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2611                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2612       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2613  *                                   *************************  *                                   *************************
2614  *                                   -----------------------------  *                                   -----------------------------
2615                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 2993  c     $                                 Line 2625  c     $                                
2625  *     **********************************************************  *     **********************************************************
2626  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2627  *     **********************************************************  *     **********************************************************
2628    cccc  scommentare se si usa al_ini della nuvola
2629    c$$$                              do i=1,5
2630    c$$$                                 AL(i)=AL_INI(i)
2631    c$$$                              enddo
2632                                  call guess()
2633                                do i=1,5                                do i=1,5
2634                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2635                                enddo                                enddo
2636                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2637                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2638                                call mini_2(jstep,ifail)                                iprint=0
2639    c                              if(DEBUG)iprint=1
2640                                  if(DEBUG)iprint=2
2641                                  call mini2(jstep,ifail,iprint)
2642                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2643                                   if(DEBUG)then                                   if(DEBUG)then
2644                                      print *,                                      print *,
2645       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2646       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2647                                        print*,'initial guess: '
2648    
2649                                        print*,'AL_INI(1) = ',AL_INI(1)
2650                                        print*,'AL_INI(2) = ',AL_INI(2)
2651                                        print*,'AL_INI(3) = ',AL_INI(3)
2652                                        print*,'AL_INI(4) = ',AL_INI(4)
2653                                        print*,'AL_INI(5) = ',AL_INI(5)
2654                                   endif                                   endif
2655                                   chi2=-chi2  c                                 chi2=-chi2
2656                                endif                                endif
2657  *     **********************************************************  *     **********************************************************
2658  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3024  c     $                                 Line 2671  c     $                                
2671       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2672  c                                 good2=.false.  c                                 good2=.false.
2673  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2674                                     do iv=1,nviews
2675                                        mask_view(iv) = 7
2676                                     enddo
2677                                   iflag=1                                   iflag=1
2678                                   return                                   return
2679                                endif                                endif
2680                                                                
2681                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2682                                                                
 c$$$                              ndof=0                                  
2683                                do ip=1,nplanes                                do ip=1,nplanes
2684  c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2685                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2686                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2687                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 3053  c$$$     $                               Line 2700  c$$$     $                              
2700                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2701                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2702       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2703                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2704         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2705                                        LADDER_STORE(nplanes-ip+1,ntracks)
2706         $                                   = LADDER(
2707         $                                   clx(ip,icp_cp(
2708         $                                   cp_match(ip,hit_plane(ip)
2709         $                                   ))));
2710                                   else                                   else
2711                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2712                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2713                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2714                                   endif                                   endif
2715                                     BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2716                                     BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2717                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(nplanes-ip+1,ntracks)=0
2718                                   do i=1,5                                   do i=1,5
2719                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2720                                   enddo                                   enddo
2721                                enddo                                enddo
2722                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2723                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2724                                                                
2725  *     --------------------------------  *     --------------------------------
# Line 3089  c$$$  rchi2=chi2/dble(ndof) Line 2743  c$$$  rchi2=chi2/dble(ndof)
2743           return           return
2744        endif        endif
2745                
2746    c$$$      if(DEBUG)then
2747    c$$$         print*,'****** TRACK CANDIDATES ***********'
2748    c$$$         print*,'#         R. chi2        RIG'
2749    c$$$         do i=1,ntracks
2750    c$$$            print*,i,' --- ',rchi2_store(i),' --- '
2751    c$$$     $           ,1./abs(AL_STORE(5,i))
2752    c$$$         enddo
2753    c$$$         print*,'***********************************'
2754    c$$$      endif
2755        if(DEBUG)then        if(DEBUG)then
2756           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
2757           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
2758           do i=1,ntracks          do i=1,ntracks
2759              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
2760       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
2761           enddo              ndof=ndof           !(1)
2762           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
2763         $           +int(ygood_store(ii,i)) !(1)
2764              enddo                 !(1)
2765              print*,i,' --- ',rchi2_store(i),' --- '
2766         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
2767            enddo
2768            print*,'*****************************************'
2769        endif        endif
2770                
2771                
# Line 3115  c$$$  rchi2=chi2/dble(ndof) Line 2784  c$$$  rchi2=chi2/dble(ndof)
2784    
2785        subroutine refine_track(ibest)        subroutine refine_track(ibest)
2786    
 c******************************************************  
 cccccc 06/10/2005 modified by elena vannuccini ---> (1)  
 cccccc 31/01/2006 modified by elena vannuccini ---> (2)  
 cccccc 12/08/2006 modified by elena vannucicni ---> (3)  
 c******************************************************  
2787    
2788        include 'commontracker.f'        include 'commontracker.f'
2789          include 'level1.f'
2790        include 'common_momanhough.f'        include 'common_momanhough.f'
2791        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2792        include 'common_mini_2.f'        include 'common_mini_2.f'
2793        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
       include 'level1.f'  
2794        include 'calib.f'        include 'calib.f'
2795    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
2796  *     flag to chose PFA  *     flag to chose PFA
2797        character*10 PFA        character*10 PFA
2798        common/FINALPFA/PFA        common/FINALPFA/PFA
2799    
2800          real k(6)
2801          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
2802    
2803          real xp,yp,zp
2804          real xyzp(3),bxyz(3)
2805          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
2806    
2807  *     =================================================  *     =================================================
2808  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
2809  *                          and  *                          and
# Line 3145  c      common/dbg/DEBUG Line 2812  c      common/dbg/DEBUG
2812        call track_init        call track_init
2813        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2814    
2815             xP=XV_STORE(nplanes-ip+1,ibest)
2816             yP=YV_STORE(nplanes-ip+1,ibest)
2817             zP=ZV_STORE(nplanes-ip+1,ibest)
2818             call gufld(xyzp,bxyz)
2819             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
2820             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
2821    c$$$  bxyz(1)=0
2822    c$$$         bxyz(2)=0
2823    c$$$         bxyz(3)=0
2824    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2825  *     -------------------------------------------------  *     -------------------------------------------------
2826  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2827  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2828  *     using improved PFAs  *     using improved PFAs
2829  *     -------------------------------------------------  *     -------------------------------------------------
2830    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2831           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2832       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2833                            
# Line 3163  c      common/dbg/DEBUG Line 2841  c      common/dbg/DEBUG
2841       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
2842              icx=clx(ip,icp)              icx=clx(ip,icp)
2843              icy=cly(ip,icp)              icy=cly(ip,icp)
2844    c            call xyz_PAM(icx,icy,is,
2845    c     $           PFA,PFA,
2846    c     $           AXV_STORE(nplanes-ip+1,ibest),
2847    c     $           AYV_STORE(nplanes-ip+1,ibest))
2848              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
2849       $           PFA,PFA,       $           PFA,PFA,
2850       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
2851       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
2852  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
2853  c$$$  $              'COG2','COG2',       $           bxyz(2)
2854  c$$$  $              0.,       $           )
2855  c$$$  $              0.)  
2856              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
2857              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
2858              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 3180  c$$$  $              0.) Line 2861  c$$$  $              0.)
2861              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
2862              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
2863    
2864  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
2865              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
             dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  
2866                            
2867    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2868  *     -------------------------------------------------  *     -------------------------------------------------
2869  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
2870  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
2871  *     -------------------------------------------------  *     -------------------------------------------------
2872    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2873           else                             else                  
2874                                
2875              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3195  c            dedxtrk(nplanes-ip+1) = (de Line 2877  c            dedxtrk(nplanes-ip+1) = (de
2877                                
2878  *     --------------------------------------------------------------  *     --------------------------------------------------------------
2879  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
             xP=XV_STORE(nplanes-ip+1,ibest)  
             yP=YV_STORE(nplanes-ip+1,ibest)  
             zP=ZV_STORE(nplanes-ip+1,ibest)  
2880              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
2881  *     if the track hit the plane in a dead area, go to the next plane  *     if the track hit the plane in a dead area, go to the next plane
2882              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
2883    
2884                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
2885                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
2886  *     --------------------------------------------------------------  *     --------------------------------------------------------------
2887    
2888              if(DEBUG)then              if(DEBUG)then
# Line 3237  c     $              cl_used(icy).eq.1.o Line 2919  c     $              cl_used(icy).eq.1.o
2919  *            *          
2920                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
2921       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
2922       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
2923       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
2924         $              bxyz(1),
2925         $              bxyz(2)
2926         $              )
2927                                
2928                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2929    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
2930                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2931                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2932       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
2933                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
2934                    xmm = xPAM                    xmm = xPAM
2935                    ymm = yPAM                    ymm = yPAM
# Line 3253  c     $              'ETA2','ETA2', Line 2938  c     $              'ETA2','ETA2',
2938                    rymm = resyPAM                    rymm = resyPAM
2939                    distmin = distance                    distmin = distance
2940                    idm = id                                      idm = id                  
2941  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2942                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2943                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
2944                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
2945         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
2946                 endif                 endif
2947   1188          continue   1188          continue
2948              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
2949              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
2950                if(distmin.le.clincnewc)then     !QUIQUI              
2951  *              -----------------------------------  *              -----------------------------------
2952                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
2953                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
2954                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
2955                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
2956                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
2957                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
2958                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
2959  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
2960                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
2961  *              -----------------------------------  *              -----------------------------------
2962                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
2963                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG)print*,'%%%% included couple ',idm
2964       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
2965                 goto 133         !next plane                 goto 133         !next plane
2966              endif              endif
2967  *     ================================================  *     ================================================
# Line 3307  c            if(DEBUG)print*,'>>>> try t Line 2994  c            if(DEBUG)print*,'>>>> try t
2994  c               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
2995                 if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)                 if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
2996  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
2997    c               call xyz_PAM(icx,0,ist,
2998    c     $              PFA,PFA,
2999    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3000                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3001       $              PFA,PFA,       $              PFA,PFA,
3002       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3003         $              bxyz(1),
3004         $              bxyz(2)
3005         $              )              
3006                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3007  c     if(DEBUG)print*,'normalized distance ',distance  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3008                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3009       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3010                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3011                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3012                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3326  c     if(DEBUG)print*,'normalized distan Line 3018  c     if(DEBUG)print*,'normalized distan
3018                    rymm = resyPAM                    rymm = resyPAM
3019                    distmin = distance                    distmin = distance
3020                    iclm = icx                    iclm = icx
3021  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3022                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3023                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3024                 endif                                   endif                  
3025  11881          continue  11881          continue
# Line 3335  c                  dedxmm = dedx(icx) !( Line 3027  c                  dedxmm = dedx(icx) !(
3027  c               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
3028                 if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)                 if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3029  *                                              !jump to the next couple  *                                              !jump to the next couple
3030    c               call xyz_PAM(0,icy,ist,
3031    c     $              PFA,PFA,
3032    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3033                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3034       $              PFA,PFA,       $              PFA,PFA,
3035       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3036         $              bxyz(1),
3037         $              bxyz(2)
3038         $              )
3039                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3040    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3041                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3042       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3043                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3044                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3045                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3353  c     $              'ETA2','ETA2', Line 3051  c     $              'ETA2','ETA2',
3051                    rymm = resyPAM                    rymm = resyPAM
3052                    distmin = distance                    distmin = distance
3053                    iclm = icy                    iclm = icy
3054  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3055                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3056                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3057                 endif                                   endif                  
3058  11882          continue  11882          continue
3059              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3060  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3061    c            print*,'## ncls(',ip,') ',ncls(ip)
3062              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3063                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3064  c               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
# Line 3368  c               if(cl_used(icl).eq.1.or. Line 3067  c               if(cl_used(icl).eq.1.or.
3067       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3068                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3069                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3070       $                 PFA,PFA,       $                 PFA,PFA,
3071       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3072         $                 bxyz(1),
3073         $                 bxyz(2)
3074         $                 )
3075                 else                         !<---- Y view                 else                         !<---- Y view
3076                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3077       $                 PFA,PFA,       $                 PFA,PFA,
3078       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3079         $                 bxyz(1),
3080         $                 bxyz(2)
3081         $                 )
3082                 endif                 endif
3083    
3084                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3085    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3086                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3087       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance,'<',distmin,' ?'
3088                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3089                      if(DEBUG)print*,'YES'
3090                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3091                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3092                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3392  c     $                 'ETA2','ETA2', Line 3097  c     $                 'ETA2','ETA2',
3097                    rymm = resyPAM                    rymm = resyPAM
3098                    distmin = distance                      distmin = distance  
3099                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3100                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3101                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3102                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3103                    else          !<---- Y view                    else          !<---- Y view
3104                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3105                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3106                    endif                    endif
3107                 endif                                   endif                  
3108  18882          continue  18882          continue
3109              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3110    c            print*,'## distmin ', distmin,' clinc ',clinc
3111    
3112              if(distmin.le.clinc)then                    c     QUIQUI------------
3113                  c     anche qui: non ci vuole clinc???
3114                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
3115                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3116                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3117                    resx(nplanes-ip+1)=rxmm       $                 20*
3118                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3119       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'                 else if(mod(VIEW(iclm),2).ne.0)then
3120                 else                    clincnew=
3121                    YGOOD(nplanes-ip+1)=1.       $                 10*
3122                    resy(nplanes-ip+1)=rymm       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3123                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 endif
3124       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c     QUIQUI------------
3125                  
3126                   if(distmin.le.clincnew)then   !QUIQUI
3127    c     if(distmin.le.clinc)then          !QUIQUI          
3128                      
3129                      CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3130    *     ----------------------------
3131    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3132                      if(mod(VIEW(iclm),2).eq.0)then
3133                         XGOOD(nplanes-ip+1)=1.
3134                         resx(nplanes-ip+1)=rxmm
3135                         if(DEBUG)print*,'%%%% included X-cl ',iclm
3136         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3137         $                    ,', dist.= ',distmin
3138         $                    ,', cut ',clinc,' )'
3139                      else
3140                         YGOOD(nplanes-ip+1)=1.
3141                         resy(nplanes-ip+1)=rymm
3142                         if(DEBUG)print*,'%%%% included Y-cl ',iclm
3143         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3144         $                    ,', dist.= ', distmin
3145         $                    ,', cut ',clinc,' )'
3146                      endif
3147    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3148    *     ----------------------------
3149                      xm_A(nplanes-ip+1) = xmm_A
3150                      ym_A(nplanes-ip+1) = ymm_A
3151                      xm_B(nplanes-ip+1) = xmm_B
3152                      ym_B(nplanes-ip+1) = ymm_B
3153                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3154                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3155                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3156    *     ----------------------------
3157                 endif                 endif
 *              ----------------------------  
                xm_A(nplanes-ip+1) = xmm_A  
                ym_A(nplanes-ip+1) = ymm_A  
                xm_B(nplanes-ip+1) = xmm_B  
                ym_B(nplanes-ip+1) = ymm_B  
                zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.  
 c              dedxtrk(nplanes-ip+1) = dedxmm !<<<    !(1)  
                dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)  
                dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)  
 *              ----------------------------  
3158              endif              endif
3159           endif           endif
3160   133     continue   133     continue
# Line 3447  c              dedxtrk(nplanes-ip+1) = d Line 3173  c              dedxtrk(nplanes-ip+1) = d
3173  *                                                 *  *                                                 *
3174  *                                                 *  *                                                 *
3175  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3176  *  *
3177        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3178    
3179        include 'commontracker.f'        include 'commontracker.f'
3180          include 'level1.f'
3181        include 'common_momanhough.f'        include 'common_momanhough.f'
3182        include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
3183    
3184        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3185    
# Line 3470  c      common/dbg/DEBUG Line 3189  c      common/dbg/DEBUG
3189              if(id.ne.0)then              if(id.ne.0)then
3190                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3191                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3192  c               cl_used(iclx)=1  !tag used clusters                 cl_used(iclx)=ntrk  !tag used clusters
3193  c               cl_used(icly)=1  !tag used clusters                 cl_used(icly)=ntrk  !tag used clusters
                cl_used(iclx)=ntrk  !tag used clusters !(1)  
                cl_used(icly)=ntrk  !tag used clusters !(1)  
3194              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3195  c               cl_used(icl)=1   !tag used clusters                 cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3196              endif              endif
3197                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3198  *     -----------------------------  *     -----------------------------
3199  *     remove the couple from clouds  *     remove the couple from clouds
3200  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3532  c               endif Line 3245  c               endif
3245    
3246    
3247    
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      real function fbad_cog(ncog,ic)  
 c$$$  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$*     --> signal of the central strip  
 c$$$      sc = CLSIGNAL(INDMAX(ic)) !center  
 c$$$  
 c$$$*     signal of adjacent strips  
 c$$$*     --> left  
 c$$$      sl1 = 0                  !left 1  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
 c$$$  
 c$$$      sl2 = 0                  !left 2  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
 c$$$  
 c$$$*     --> right  
 c$$$      sr1 = 0                  !right 1  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
 c$$$  
 c$$$      sr2 = 0                  !right 2  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
 c$$$  
 c$$$  
 c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  
 c$$$         f  = 4.  
 c$$$         si = 8.4  
 c$$$      else                              !X-view  
 c$$$         f  = 6.  
 c$$$         si = 3.9  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = 1.  
 c$$$      f0 = 1  
 c$$$      f1 = 1  
 c$$$      f2 = 1  
 c$$$      f3 = 1    
 c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then  
 c$$$          
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$          
 c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then  
 c$$$  
 c$$$  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sr1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = sqrt(fbad_cog)  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
3248    
3249    
3250  *     ****************************************************  *     ****************************************************
3251    
3252        subroutine init_level2        subroutine init_level2
3253    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3254        include 'commontracker.f'        include 'commontracker.f'
3255          include 'level1.f'
3256        include 'common_momanhough.f'        include 'common_momanhough.f'
3257        include 'level2.f'        include 'level2.f'
       include 'level1.f'  
3258    
3259    *     ---------------------------------
3260    *     variables initialized from level1
3261    *     ---------------------------------
3262        do i=1,nviews        do i=1,nviews
3263           good2(i)=good1(i)           good2(i)=good1(i)
3264             do j=1,nva1_view
3265                vkflag(i,j)=1
3266                if(cnnev(i,j).le.0)then
3267                   vkflag(i,j)=cnnev(i,j)
3268                endif
3269             enddo
3270        enddo        enddo
3271    *     ----------------
3272  c      good2 = 0!.false.  *     level2 variables
3273  c$$$      nev2 = nev1  *     ----------------
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
   
3274        NTRK = 0        NTRK = 0
3275        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3276           IMAGE(IT)=0           IMAGE(IT)=0
3277           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3278           do ip=1,nplanes           do ip=1,nplanes
3279              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3280              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
3281              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3282              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3283              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3284                TAILX_nt(IP,IT) = 0
3285                TAILY_nt(IP,IT) = 0
3286                XBAD(IP,IT) = 0
3287                YBAD(IP,IT) = 0
3288              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3289              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3290  c*****************************************************              LS(IP,IT) = 0
 cccccc 11/9/2005 modified by david fedele  
3291              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3292              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3293              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3294              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3295           enddo           enddo
# Line 3700  cccccc 17/8/2006 modified by elena Line 3300  cccccc 17/8/2006 modified by elena
3300              enddo                                enddo                  
3301           enddo                             enddo                  
3302        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3303        nclsx=0        nclsx=0
3304        nclsy=0              nclsy=0      
3305        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3306          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3307          xs(1,ip)=0          xs(1,ip)=0
3308          xs(2,ip)=0          xs(2,ip)=0
3309          sgnlxs(ip)=0          sgnlxs(ip)=0
3310          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3311          ys(1,ip)=0          ys(1,ip)=0
3312          ys(2,ip)=0          ys(2,ip)=0
3313          sgnlys(ip)=0          sgnlys(ip)=0
3314        enddo        enddo
 c*******************************************************  
3315        end        end
3316    
3317    
# Line 3733  c*************************************** Line 3326  c***************************************
3326  ************************************************************  ************************************************************
3327    
3328    
3329          subroutine init_hough
3330    
3331          include 'commontracker.f'
3332          include 'level1.f'
3333          include 'common_momanhough.f'
3334          include 'common_hough.f'
3335          include 'level2.f'
3336    
3337          ntrpt_nt=0
3338          ndblt_nt=0
3339          NCLOUDS_XZ_nt=0
3340          NCLOUDS_YZ_nt=0
3341          do idb=1,ndblt_max_nt
3342             db_cloud_nt(idb)=0
3343             alfayz1_nt(idb)=0      
3344             alfayz2_nt(idb)=0      
3345          enddo
3346          do itr=1,ntrpt_max_nt
3347             tr_cloud_nt(itr)=0
3348             alfaxz1_nt(itr)=0      
3349             alfaxz2_nt(itr)=0      
3350             alfaxz3_nt(itr)=0      
3351          enddo
3352          do idb=1,ncloyz_max      
3353            ptcloud_yz_nt(idb)=0    
3354            alfayz1_av_nt(idb)=0    
3355            alfayz2_av_nt(idb)=0    
3356          enddo                    
3357          do itr=1,ncloxz_max      
3358            ptcloud_xz_nt(itr)=0    
3359            alfaxz1_av_nt(itr)=0    
3360            alfaxz2_av_nt(itr)=0    
3361            alfaxz3_av_nt(itr)=0    
3362          enddo                    
3363    
3364          ntrpt=0                  
3365          ndblt=0                  
3366          NCLOUDS_XZ=0              
3367          NCLOUDS_YZ=0              
3368          do idb=1,ndblt_max        
3369            db_cloud(idb)=0        
3370            cpyz1(idb)=0            
3371            cpyz2(idb)=0            
3372            alfayz1(idb)=0          
3373            alfayz2(idb)=0          
3374          enddo                    
3375          do itr=1,ntrpt_max        
3376            tr_cloud(itr)=0        
3377            cpxz1(itr)=0            
3378            cpxz2(itr)=0            
3379            cpxz3(itr)=0            
3380            alfaxz1(itr)=0          
3381            alfaxz2(itr)=0          
3382            alfaxz3(itr)=0          
3383          enddo                    
3384          do idb=1,ncloyz_max      
3385            ptcloud_yz(idb)=0      
3386            alfayz1_av(idb)=0      
3387            alfayz2_av(idb)=0      
3388            do idbb=1,ncouplemaxtot
3389              cpcloud_yz(idb,idbb)=0
3390            enddo                  
3391          enddo                    
3392          do itr=1,ncloxz_max      
3393            ptcloud_xz(itr)=0      
3394            alfaxz1_av(itr)=0      
3395            alfaxz2_av(itr)=0      
3396            alfaxz3_av(itr)=0      
3397            do itrr=1,ncouplemaxtot
3398              cpcloud_xz(itr,itrr)=0
3399            enddo                  
3400          enddo                    
3401          end
3402    ************************************************************
3403    *
3404    *
3405    *
3406    *
3407    *
3408    *
3409    *
3410    ************************************************************
3411    
3412    
3413        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3414    
3415  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3744  c*************************************** Line 3421  c***************************************
3421            
3422        include 'commontracker.f'        include 'commontracker.f'
3423        include 'level1.f'        include 'level1.f'
3424          include 'common_momanhough.f'
3425        include 'level2.f'        include 'level2.f'
3426        include 'common_mini_2.f'        include 'common_mini_2.f'
3427        include 'common_momanhough.f'        include 'calib.f'
3428        real sinth,phi,pig        !(4)  
3429          character*10 PFA
3430          common/FINALPFA/PFA
3431    
3432          real sinth,phi,pig
3433          integer ssensor,sladder
3434        pig=acos(-1.)        pig=acos(-1.)
3435    
3436  c      good2=1!.true.  *     -------------------------------------
3437        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3438        nstep_nt(ntr)       = 0!nstep        nstep_nt(ntr)       = nstep
3439    *     -------------------------------------
3440        phi   = al(4)             !(4)        phi   = al(4)          
3441        sinth = al(3)             !(4)        sinth = al(3)            
3442        if(sinth.lt.0)then        !(4)        if(sinth.lt.0)then      
3443           sinth = -sinth         !(4)           sinth = -sinth        
3444           phi = phi + pig        !(4)           phi = phi + pig      
3445        endif                     !(4)        endif                    
3446        npig = aint(phi/(2*pig))  !(4)        npig = aint(phi/(2*pig))
3447        phi = phi - npig*2*pig    !(4)        phi = phi - npig*2*pig  
3448        if(phi.lt.0)              !(4)        if(phi.lt.0)            
3449       $     phi = phi + 2*pig    !(4)       $     phi = phi + 2*pig  
3450        al(4) = phi               !(4)        al(4) = phi              
3451        al(3) = sinth             !(4)        al(3) = sinth            
 *****************************************************  
3452        do i=1,5        do i=1,5
3453           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3454           do j=1,5           do j=1,5
3455              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3456           enddo           enddo
 c     print*,al_nt(i,ntr)  
3457        enddo        enddo
3458          *     -------------------------------------      
3459        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3460           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3461           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3783  c     print*,al_nt(i,ntr) Line 3464  c     print*,al_nt(i,ntr)
3464           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3465           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3466           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3467             TAILX_nt(IP,ntr) = 0.
3468             TAILY_nt(IP,ntr) = 0.
3469           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3470           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3471           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3472           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3473           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3474  c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3475           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           factor = sqrt(
3476           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3477         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3478         $        1. )
3479    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3480             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3481             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3482        
3483           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3484             ay   = ayv_nt(ip,ntr)
3485             bfx  = BX_STORE(ip,IDCAND)
3486             bfy  = BY_STORE(ip,IDCAND)
3487             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3488             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3489             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3490             angx     = 180.*atan(tgtemp)/acos(-1.)
3491             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3492             angy    = 180.*atan(tgtemp)/acos(-1.)
3493            
3494    c         print*,'* ',ip,bfx,bfy,angx,angy
3495    
3496             id  = CP_STORE(ip,IDCAND) ! couple id
3497           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3498             ssensor = -1
3499             sladder = -1
3500             ssensor = SENSOR_STORE(ip,IDCAND)
3501             sladder = LADDER_STORE(ip,IDCAND)
3502             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3503             LS(IP,ntr)      = ssensor+10*sladder
3504    
3505           if(id.ne.0)then           if(id.ne.0)then
3506    c           >>> is a couple
3507              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3508              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3509  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3510    c$$$            if(is_cp(id).ne.ssensor)
3511    c$$$     $           print*,'ERROR is sensor assignment (couple)'
3512    c$$$     $           ,is_cp(id),ssensor
3513    c$$$            if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder)
3514    c$$$     $           print*,'ERROR is ladder assignment (couple)'
3515    c$$$     $           ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder
3516                
3517                nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3518                nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3519                xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3520                ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3521    
3522                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3523         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3524                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3525         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3526    
3527           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3528              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  c           >>> is a singlet
3529              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl  c$$$            if(LADDER(icl).ne.sladder)
3530  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  c$$$     $           print*,'ERROR is ladder assignment (single)'
3531    c$$$     $           ,LADDER(icl),sladder
3532                if(mod(VIEW(icl),2).eq.0)then
3533                   cltrx(ip,ntr)=icl
3534                   nnnnn = npfastrips(icl,PFA,angx)
3535                   xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3536                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3537                elseif(mod(VIEW(icl),2).eq.1)then
3538                   cltry(ip,ntr)=icl
3539                   nnnnn = npfastrips(icl,PFA,angy)
3540                   ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3541                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3542                endif
3543           endif                     endif          
3544    
3545        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)  
3546    
3547    
3548    c$$$      print*,(xgood(i),i=1,6)
3549    c$$$      print*,(ygood(i),i=1,6)
3550    c$$$      print*,(ls(i,ntr),i=1,6)
3551    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3552    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3553    c$$$      print*,'-----------------------'
3554    
3555        end        end
3556    
3557        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*****************************************************  
3558    
3559  *     -------------------------------------------------------  *     -------------------------------------------------------
3560  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3829  c*************************************** Line 3563  c***************************************
3563  *     -------------------------------------------------------  *     -------------------------------------------------------
3564    
3565        include 'commontracker.f'        include 'commontracker.f'
       include 'level1.f'  
       include 'level2.f'  
3566        include 'calib.f'        include 'calib.f'
3567          include 'level1.f'
3568        include 'common_momanhough.f'        include 'common_momanhough.f'
3569          include 'level2.f'
3570        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3571    
3572  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3573        nclsx = 0        nclsx = 0
3574        nclsy = 0        nclsy = 0
3575    
3576          do iv = 1,nviews
3577             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3578          enddo
3579    
3580        do icl=1,nclstr1        do icl=1,nclstr1
3581           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
3582              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
3583              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3584                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3585                 planex(nclsx) = ip                 planex(nclsx) = ip
3586                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3587                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3588                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3589                 do is=1,2                 do is=1,2
3590  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3591                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3592                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3593                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3594                 enddo                 enddo
3595  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3861  c$$$               print*,'xs(2,nclsx)   Line 3600  c$$$               print*,'xs(2,nclsx)  
3600              else                          !=== Y views              else                          !=== Y views
3601                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3602                 planey(nclsy) = ip                 planey(nclsy) = ip
3603                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3604                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3605                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3606                 do is=1,2                 do is=1,2
3607  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3608                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3609                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3610                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3611                 enddo                 enddo
3612  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3875  c$$$               print*,'ys(1,nclsy)   Line 3616  c$$$               print*,'ys(1,nclsy)  
3616  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3617              endif              endif
3618           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3619    
3620  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3621           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
# Line 3883  c      print*,icl,cl_used(icl),cl_good(i Line 3623  c      print*,icl,cl_used(icl),cl_good(i
3623        enddo        enddo
3624        end        end
3625    
3626    ***************************************************
3627    *                                                 *
3628    *                                                 *
3629    *                                                 *
3630    *                                                 *
3631    *                                                 *
3632    *                                                 *
3633    **************************************************
3634    
3635          subroutine fill_hough
3636    
3637    *     -------------------------------------------------------
3638    *     This routine fills the  variables related to the hough
3639    *     transform, for the debig n-tuple
3640    *     -------------------------------------------------------
3641    
3642          include 'commontracker.f'
3643          include 'level1.f'
3644          include 'common_momanhough.f'
3645          include 'common_hough.f'
3646          include 'level2.f'
3647    
3648          if(.false.
3649         $     .or.ntrpt.gt.ntrpt_max_nt
3650         $     .or.ndblt.gt.ndblt_max_nt
3651         $     .or.NCLOUDS_XZ.gt.ncloxz_max
3652         $     .or.NCLOUDS_yZ.gt.ncloyz_max
3653         $     )then
3654             ntrpt_nt=0
3655             ndblt_nt=0
3656             NCLOUDS_XZ_nt=0
3657             NCLOUDS_YZ_nt=0
3658          else
3659             ndblt_nt=ndblt
3660             ntrpt_nt=ntrpt
3661             if(ndblt.ne.0)then
3662                do id=1,ndblt
3663                   alfayz1_nt(id)=alfayz1(id) !Y0
3664                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3665                enddo
3666             endif
3667             if(ndblt.ne.0)then
3668                do it=1,ntrpt
3669                   alfaxz1_nt(it)=alfaxz1(it) !X0
3670                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3671                   alfaxz3_nt(it)=alfaxz3(it) !1/r
3672                enddo
3673             endif
3674             nclouds_yz_nt=nclouds_yz
3675             nclouds_xz_nt=nclouds_xz
3676             if(nclouds_yz.ne.0)then
3677                nnn=0
3678                do iyz=1,nclouds_yz
3679                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3680                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3681                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3682                   nnn=nnn+ptcloud_yz(iyz)
3683                enddo
3684                do ipt=1,nnn
3685                   db_cloud_nt(ipt)=db_cloud(ipt)
3686                 enddo
3687             endif
3688             if(nclouds_xz.ne.0)then
3689                nnn=0
3690                do ixz=1,nclouds_xz
3691                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
3692                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
3693                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
3694                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
3695                   nnn=nnn+ptcloud_xz(ixz)              
3696                enddo
3697                do ipt=1,nnn
3698                  tr_cloud_nt(ipt)=tr_cloud(ipt)
3699                 enddo
3700             endif
3701          endif
3702          end
3703          
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23