/[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.4 by pam-fi, Thu Sep 28 14:04:40 2006 UTC revision 1.18 by pam-fi, Fri Feb 16 14:56:02 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    
23        include 'momanhough_init.f'  c      include 'momanhough_init.f'
24                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
25  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
26  *     STEP 1  *     STEP 1
27  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 47  c      common/dbg/DEBUG Line 44  c      common/dbg/DEBUG
44  c      iflag=0  c      iflag=0
45        call cl_to_couples(iflag)        call cl_to_couples(iflag)
46        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
47           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
48        endif        endif
49                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     selezione di tracce pulite per diagnostica  
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 c$$$         if(DEBUG)then  
 c$$$            do ip=1,nplanes  
 c$$$               if(ncp_plane(ip).ne.1)good2=.false.  
 c$$$            enddo  
 c$$$c            if(good2.eq.0)goto 100!next event  
 c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event  
 c$$$         endif  
           
   
   
50  *-----------------------------------------------------  *-----------------------------------------------------
51  *-----------------------------------------------------  *-----------------------------------------------------
52  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 94  c$$$         endif Line 78  c$$$         endif
78  c      iflag=0  c      iflag=0
79        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
80        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
81           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
82        endif        endif
83                
84                
# Line 123  c      iflag=0 Line 107  c      iflag=0
107  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    *     count number of hit planes
111          planehit=0                
112          do np=1,nplanes          
113            if(ncp_plane(np).ne.0)then  
114              planehit=planehit+1  
115            endif                  
116          enddo                    
117          if(planehit.lt.3) goto 880 ! exit              
118          
119          nptxz_min=x_min_start              
120          nplxz_min=x_min_start              
121                
122          nptyz_min=y_min_start              
123          nplyz_min=y_min_start              
124    
125  c      iflag=0        cutdistyz=cutystart      
126          cutdistxz=cutxstart      
127    
128     878  continue
129        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
130        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
131           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
132        endif        endif
133  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
134            if(cutdistyz.lt.maxcuty/2)then
135              cutdistyz=cutdistyz+cutystep
136            else
137              cutdistyz=cutdistyz+(3*cutystep)
138            endif
139            goto 878                
140          endif                    
141    
142          if(planehit.eq.3) goto 881          
143          
144     879  continue  
145        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
146        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
147           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
148        endif        endif
149                                    
150          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
151            cutdistxz=cutdistxz+cutxstep
152            goto 879                
153          endif                    
154    
155        
156     881  continue  
157    *     if there is at least three planes on the Y view decreases cuts on X view
158          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
159         $     nplxz_min.ne.y_min_start)then
160            nptxz_min=x_min_step    
161            nplxz_min=x_min_start-x_min_step              
162            goto 879                
163          endif                    
164            
165   880  return   880  return
166        end        end
167    
# Line 144  c      iflag=0 Line 171  c      iflag=0
171        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
172    
173        include 'commontracker.f'        include 'commontracker.f'
174          include 'level1.f'
175        include 'common_momanhough.f'        include 'common_momanhough.f'
176        include 'common_mech.f'        include 'common_mech.f'
177        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
178        include 'common_mini_2.f'        include 'common_mini_2.f'
179        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
180        include 'level2.f'        include 'level2.f'
181    
182        include 'momanhough_init.f'  c      include 'momanhough_init.f'
183                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
184        logical FIMAGE            !        logical FIMAGE            !
185          real*8 AL_GUESS(5)
186    
187  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
188  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 198  c         iflag=0 Line 223  c         iflag=0
223           ibest=0                !best track among candidates           ibest=0                !best track among candidates
224           iimage=0               !track image           iimage=0               !track image
225  *     ------------- select the best track -------------  *     ------------- select the best track -------------
226           rchi2best=1000000000.  c$$$         rchi2best=1000000000.
227    c$$$         do i=1,ntracks
228    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
229    c$$$     $         RCHI2_STORE(i).gt.0)then
230    c$$$               ibest=i
231    c$$$               rchi2best=RCHI2_STORE(i)
232    c$$$            endif
233    c$$$         enddo
234    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
235    
236    *     -------------------------------------------------------
237    *     order track-candidates according to:
238    *     1st) decreasing n.points
239    *     2nd) increasing chi**2
240    *     -------------------------------------------------------
241             rchi2best=1000000000.
242             ndofbest=0             !(1)
243           do i=1,ntracks           do i=1,ntracks
244              if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0               !(1)
245       $         RCHI2_STORE(i).gt.0)then             do ii=1,nplanes      !(1)
246                 ndof=ndof          !(1)
247         $            +int(xgood_store(ii,i)) !(1)
248         $            +int(ygood_store(ii,i)) !(1)
249               enddo                !(1)
250               if(ndof.gt.ndofbest)then !(1)
251                 ibest=i
252                 rchi2best=RCHI2_STORE(i)
253                 ndofbest=ndof      !(1)
254               elseif(ndof.eq.ndofbest)then !(1)
255                 if(RCHI2_STORE(i).lt.rchi2best.and.
256         $            RCHI2_STORE(i).gt.0)then
257                 ibest=i                 ibest=i
258                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
259              endif                 ndofbest=ndof    !(1)
260           enddo               endif              !(1)
261               endif
262             enddo
263    
264    c$$$         rchi2best=1000000000.
265    c$$$         ndofbest=0             !(1)
266    c$$$         do i=1,ntracks
267    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
268    c$$$     $          RCHI2_STORE(i).gt.0)then
269    c$$$             ndof=0             !(1)
270    c$$$             do ii=1,nplanes    !(1)
271    c$$$               ndof=ndof        !(1)
272    c$$$     $              +int(xgood_store(ii,i)) !(1)
273    c$$$     $              +int(ygood_store(ii,i)) !(1)
274    c$$$             enddo              !(1)
275    c$$$             if(ndof.ge.ndofbest)then !(1)
276    c$$$               ibest=i
277    c$$$               rchi2best=RCHI2_STORE(i)
278    c$$$               ndofbest=ndof    !(1)
279    c$$$             endif              !(1)
280    c$$$           endif
281    c$$$         enddo
282    
283           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
284  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
285  *     The best track candidate (ibest) is selected and a new fitting is performed.  *     The best track candidate (ibest) is selected and a new fitting is performed.
# Line 236  c         iflag=0 Line 310  c         iflag=0
310  *     **********************************************************  *     **********************************************************
311  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
312  *     **********************************************************  *     **********************************************************
313             call guess()
314             do i=1,5
315                AL_GUESS(i)=AL(i)
316             enddo
317    c         print*,'## guess: ',al
318    
319           do i=1,5           do i=1,5
320              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
321           enddo           enddo
322            
323           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
324           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
325           jstep=0                !# minimization steps           jstep=0                !# minimization steps
326    
327           call mini_2(jstep,ifail)           iprint=0
328    c         if(DEBUG)iprint=1
329             if(VERBOSE)iprint=1
330             if(DEBUG)iprint=2
331             call mini2(jstep,ifail,iprint)
332           if(ifail.ne.0) then           if(ifail.ne.0) then
333              if(DEBUG)then              if(VERBOSE)then
334                 print *,                 print *,
335       $              '*** MINIMIZATION FAILURE *** (mini_2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
336       $              ,iev       $              ,iev
337    
338    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
339    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
340    c$$$               print*,'result:   ',(al(i),i=1,5)
341    c$$$               print*,'xgood ',xgood
342    c$$$               print*,'ygood ',ygood
343    c$$$               print*,'----------------------------------------------'
344              endif              endif
345              chi2=-chi2  c            chi2=-chi2
346           endif           endif
347                    
348           if(DEBUG)then           if(DEBUG)then
# Line 311  c         print*,'++++++++++ iimage,fima Line 403  c         print*,'++++++++++ iimage,fima
403  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
404    
405           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
406              if(DEBUG)              if(verbose)
407       $           print*,       $           print*,
408       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
409       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 361  c     $        rchi2best.lt.15..and. Line 453  c     $        rchi2best.lt.15..and.
453        end        end
454    
455    
   
   
 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$$$  
   
456                
457  ************************************************************  ************************************************************
458  ************************************************************  ************************************************************
# Line 557  c$$$ Line 477  c$$$
477  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
478  *     angx   - Projected angle in x  *     angx   - Projected angle in x
479  *     angy   - Projected angle in y  *     angy   - Projected angle in y
480    *     bfx    - x-component of magnetci field
481    *     bfy    - y-component of magnetic field
482  *  *
483  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
484  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 595  c$$$ Line 517  c$$$
517  *  *
518  *  *
519    
520        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
521    
 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*****************************************************  
522                
523        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
524        include 'level1.f'        include 'level1.f'
525          include 'calib.f'
526    c      include 'level1.f'
527        include 'common_align.f'        include 'common_align.f'
528        include 'common_mech.f'        include 'common_mech.f'
529        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
530        include 'common_resxy.f'  c      include 'common_resxy.f'
531    
532  c      logical DEBUG  c      logical DEBUG
533  c      common/dbg/DEBUG  c      common/dbg/DEBUG
# Line 621  c      common/dbg/DEBUG Line 536  c      common/dbg/DEBUG
536        integer sensor        integer sensor
537        integer viewx,viewy        integer viewx,viewy
538        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
539        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
540          real angx,angy            !X-Y effective angle
541          real bfx,bfy              !X-Y b-field components
542    
543        real stripx,stripy        real stripx,stripy
544    
# Line 647  c      double precision xi_B,yi_B,zi_B Line 564  c      double precision xi_B,yi_B,zi_B
564        xPAM_B = 0.        xPAM_B = 0.
565        yPAM_B = 0.        yPAM_B = 0.
566        zPAM_B = 0.        zPAM_B = 0.
567    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
568  *     -----------------  *     -----------------
569  *     CLUSTER X  *     CLUSTER X
570  *     -----------------  *     -----------------
571    
572        if(icx.ne.0)then        if(icx.ne.0)then
573    
574           viewx = VIEW(icx)           viewx = VIEW(icx)
575           nldx = nld(MAXS(icx),VIEW(icx))           nldx = nld(MAXS(icx),VIEW(icx))
576           nplx = npl(VIEW(icx))           nplx = npl(VIEW(icx))
577           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resxPAM = RESXAV
           
578           stripx = float(MAXS(icx))           stripx = float(MAXS(icx))
579           if(PFAx.eq.'COG1')then  !(1)  *        --------------------------
580              stripx = stripx      !(1)  *        magnetic-field corrections
581              resxPAM = resxPAM    !(1)  *        --------------------------
582    c$$$         print*,nplx,ax,bfy/10.
583             angtemp  = ax
584             bfytemp  = bfy
585             if(nplx.eq.6) angtemp = -1. * ax
586             if(nplx.eq.6) bfytemp = -1. * bfy
587             tgtemp = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
588             angx = 180.*atan(tgtemp)/acos(-1.)
589             stripx = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
590    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
591    c$$$         print*,'========================'
592    *        --------------------------
593    
594             if(PFAx.eq.'COG1')then
595                stripx = stripx    
596                resxPAM = resxPAM  
597           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
598              stripx = stripx + cog(2,icx)                          stripx = stripx + cog(2,icx)            
599              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
600             elseif(PFAx.eq.'COG3')then
601                stripx = stripx + cog(3,icx)            
602                resxPAM = resxPAM*fbad_cog(3,icx)
603             elseif(PFAx.eq.'COG4')then
604    c            print*,'COG4'
605                stripx = stripx + cog(4,icx)            
606                resxPAM = resxPAM*fbad_cog(4,icx)
607           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
608  c            cog2 = cog(2,icx)              stripx = stripx + pfaeta2(icx,angx)          
609  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)                          resxPAM = risx_eta2(abs(angx))                    
 c            stripx = stripx + etacorr  
             stripx = stripx + pfa_eta2(icx,angx)            !(3)  
             resxPAM = risx_eta2(angx)                       !   (4)  
610              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
611       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
612              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
613           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
614              stripx = stripx + pfa_eta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)          
615              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(abs(angx))                      
616              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)            
617       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
618              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
619           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                        
620              stripx = stripx + pfa_eta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            
621              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(abs(angx))                      
622              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              
623       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
624              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
625           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then  
626              stripx = stripx + pfa_eta(icx,angx)             !(3)  c            print*,'ETA'
627              resxPAM = ris_eta(icx,angx)                     !   (4)              stripx = stripx + pfaeta(icx,angx)            
628              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              resxPAM = ris_eta(icx,angx)                    
629       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              
630  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
631              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              resxPAM = resxPAM*fbad_eta(icx,angx)            
632           elseif(PFAx.eq.'COG')then           !(2)           elseif(PFAx.eq.'COG')then          
633              stripx = stripx + cog(0,icx)     !(2)                      stripx = stripx + cog(0,icx)            
634              resxPAM = risx_cog(angx)                        !   (4)              resxPAM = risx_cog(abs(angx))                    
635              resxPAM = resxPAM*fbad_cog(0,icx)!(2)              resxPAM = resxPAM*fbad_cog(0,icx)
636           else           else
637              print*,'*** Non valid p.f.a. (x) --> ',PFAx              print*,'*** Non valid p.f.a. (x) --> ',PFAx
638           endif           endif
639    
640    c         print*,'%%%%%%%%%%%%'
641    
642        endif        endif
 c      if(icy.eq.0.and.icx.ne.0)  
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
643                
644  *     -----------------  *     -----------------
645  *     CLUSTER Y  *     CLUSTER Y
646  *     -----------------  *     -----------------
647    
648        if(icy.ne.0)then        if(icy.ne.0)then
649    
650           viewy = VIEW(icy)           viewy = VIEW(icy)
651           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
652           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
653           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
654             stripy = float(MAXS(icy))
655    
656           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
657              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
658       $           ,icx,icy       $           ,icx,icy
659              goto 100              goto 100
660           endif           endif
661    *        --------------------------
662    *        magnetic-field corrections
663    *        --------------------------
664             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
665             angy    = 180.*atan(tgtemp)/acos(-1.)
666             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
667    *        --------------------------
668                    
          stripy = float(MAXS(icy))  
669           if(PFAy.eq.'COG1')then !(1)           if(PFAy.eq.'COG1')then !(1)
670              stripy = stripy     !(1)              stripy = stripy     !(1)
671              resyPAM = resyPAM   !(1)              resyPAM = resyPAM   !(1)
672           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
673              stripy = stripy + cog(2,icy)              stripy = stripy + cog(2,icy)
674              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
675             elseif(PFAy.eq.'COG3')then
676                stripy = stripy + cog(3,icy)
677                resyPAM = resyPAM*fbad_cog(3,icy)
678             elseif(PFAy.eq.'COG4')then
679                stripy = stripy + cog(4,icy)
680                resyPAM = resyPAM*fbad_cog(4,icy)
681           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
682  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
683  c            etacorr = pfa_eta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
684  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
685              stripy = stripy + pfa_eta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
686              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(abs(angy))                       !   (4)
687              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
688              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
689       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
690           elseif(PFAy.eq.'ETA3')then                         !(3)           elseif(PFAy.eq.'ETA3')then                         !(3)
691              stripy = stripy + pfa_eta3(icy,angy)            !(3)              stripy = stripy + pfaeta3(icy,angy)            !(3)
692              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)
693              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)
694       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
695           elseif(PFAy.eq.'ETA4')then                         !(3)           elseif(PFAy.eq.'ETA4')then                         !(3)
696              stripy = stripy + pfa_eta4(icy,angy)            !(3)              stripy = stripy + pfaeta4(icy,angy)            !(3)
697              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)
698              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)
699       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
700           elseif(PFAy.eq.'ETA')then                          !(3)           elseif(PFAy.eq.'ETA')then                          !(3)
701              stripy = stripy + pfa_eta(icy,angy)             !(3)              stripy = stripy + pfaeta(icy,angy)             !(3)
702              resyPAM = ris_eta(icy,angy)                     !   (4)              resyPAM = ris_eta(icy,angy)                     !   (4)
703  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO
704              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)
# Line 757  c            resyPAM = resyPAM*fbad_cog( Line 706  c            resyPAM = resyPAM*fbad_cog(
706       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)
707           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
708              stripy = stripy + cog(0,icy)                          stripy = stripy + cog(0,icy)            
709              resyPAM = risy_cog(angy)                        !   (4)              resyPAM = risy_cog(abs(angy))                        !   (4)
710  c            resyPAM = ris_eta(icy,angy)                    !   (4)  c            resyPAM = ris_eta(icy,angy)                    !   (4)
711              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
712           else           else
# Line 766  c            resyPAM = ris_eta(icy,angy) Line 715  c            resyPAM = ris_eta(icy,angy)
715    
716        endif        endif
717    
718          c      print*,'## stripx,stripy ',stripx,stripy
719    
720  c===========================================================  c===========================================================
721  C     COUPLE  C     COUPLE
722  C===========================================================  C===========================================================
# Line 968  c         print*,'A-(',xPAM_A,yPAM_A,') Line 918  c         print*,'A-(',xPAM_A,yPAM_A,')
918                            
919        endif        endif
920                    
921    
922    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
923    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
924    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
925    
926   100  continue   100  continue
927        end        end
928    
# Line 1284  c     $              ,iv,xvv(iv),yvv(iv) Line 1239  c     $              ,iv,xvv(iv),yvv(iv)
1239  *     it returns the plane number  *     it returns the plane number
1240  *      *    
1241        include 'commontracker.f'        include 'commontracker.f'
1242          include 'level1.f'
1243  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1244        include 'common_momanhough.f'        include 'common_momanhough.f'
1245                
# Line 1321  c      include 'common_analysis.f' Line 1277  c      include 'common_analysis.f'
1277  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1278  *      *    
1279        include 'commontracker.f'        include 'commontracker.f'
1280          include 'level1.f'
1281  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1282        include 'common_momanhough.f'        include 'common_momanhough.f'
1283                
# Line 1349  c      include 'common_analysis.f' Line 1306  c      include 'common_analysis.f'
1306  *     positive if sensor =2  *     positive if sensor =2
1307  *  *
1308        include 'commontracker.f'        include 'commontracker.f'
1309          include 'level1.f'
1310  c      include 'calib.f'  c      include 'calib.f'
1311  c      include 'level1.f'  c      include 'level1.f'
1312  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1378  c      include 'common_analysis.f' Line 1336  c      include 'common_analysis.f'
1336  *************************************************************************  *************************************************************************
1337  *************************************************************************  *************************************************************************
1338  *************************************************************************  *************************************************************************
 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  
1339                
1340    
1341  ***************************************************  ***************************************************
# Line 1660  c$$$      end Line 1350  c$$$      end
1350        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1351    
1352        include 'commontracker.f'        include 'commontracker.f'
1353          include 'level1.f'
1354        include 'common_momanhough.f'        include 'common_momanhough.f'
1355        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1356        include 'calib.f'        include 'calib.f'
1357        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1358    
1359  *     output flag  *     output flag
1360  *     --------------  *     --------------
# Line 1675  c      common/dbg/DEBUG Line 1363  c      common/dbg/DEBUG
1363  *     --------------  *     --------------
1364        integer iflag        integer iflag
1365    
1366        integer badseed,badcl        integer badseed,badclx,badcly
1367    
1368  *     init variables  *     init variables
1369        ncp_tot=0        ncp_tot=0
# Line 1691  c      common/dbg/DEBUG Line 1379  c      common/dbg/DEBUG
1379           ncls(ip)=0           ncls(ip)=0
1380        enddo        enddo
1381        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1382           cl_single(icl)=1           cl_single(icl) = 1
1383           cl_good(icl)=0           cl_good(icl)   = 0
1384          enddo
1385          do iv=1,nviews
1386             ncl_view(iv)  = 0
1387             mask_view(iv) = 0      !all included
1388        enddo        enddo
1389                
1390    *     count number of cluster per view
1391          do icl=1,nclstr1
1392             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1393          enddo
1394    *     mask views with too many clusters
1395          do iv=1,nviews
1396             if( ncl_view(iv).gt. nclusterlimit)then
1397                mask_view(iv) = 1
1398                if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1399         $           ,nclusterlimit,' on view ', iv,' --> masked!'
1400             endif
1401          enddo
1402    
1403    
1404  *     start association  *     start association
1405        ncouples=0        ncouples=0
1406        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1407           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1408                    
1409  *     ----------------------------------------------------  *     ----------------------------------------------------
1410    *     jump masked views (X VIEW)
1411    *     ----------------------------------------------------
1412             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1413    *     ----------------------------------------------------
1414  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1415  *     ----------------------------------------------------  *     ----------------------------------------------------
1416           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1417              cl_single(icx)=0              cl_single(icx)=0
1418              goto 10              goto 10
1419           endif           endif
# Line 1717  c      common/dbg/DEBUG Line 1427  c      common/dbg/DEBUG
1427           else           else
1428              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1429           endif           endif
1430           badcl=badseed           badclx=badseed
1431           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1432              ibad=1              ibad=1
1433              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1727  c      common/dbg/DEBUG Line 1437  c      common/dbg/DEBUG
1437       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1438       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1439              endif              endif
1440              badcl=badcl*ibad              badclx=badclx*ibad
1441           enddo           enddo
1442  *     ----------------------------------------------------  *     ----------------------------------------------------
1443  *     >>> eliminato il taglio sulle BAD <<<  *     >>> eliminato il taglio sulle BAD <<<
# Line 1746  c     endif Line 1456  c     endif
1456              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1457                            
1458  *     ----------------------------------------------------  *     ----------------------------------------------------
1459    *     jump masked views (Y VIEW)
1460    *     ----------------------------------------------------
1461                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1462    
1463    *     ----------------------------------------------------
1464  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1465  *     ----------------------------------------------------  *     ----------------------------------------------------
1466              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1467                 cl_single(icy)=0                 cl_single(icy)=0
1468                 goto 20                 goto 20
1469              endif              endif
# Line 1762  c     endif Line 1477  c     endif
1477              else              else
1478                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1479              endif              endif
1480              badcl=badseed              badcly=badseed
1481              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1482                 ibad=1                 ibad=1
1483                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1771  c     endif Line 1486  c     endif
1486       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1487       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1488       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1489                 badcl=badcl*ibad                 badcly=badcly*ibad
1490              enddo              enddo
1491  *     ----------------------------------------------------  *     ----------------------------------------------------
1492  *     >>> eliminato il taglio sulle BAD <<<  *     >>> eliminato il taglio sulle BAD <<<
# Line 1794  c     endif Line 1509  c     endif
1509  *     charge correlation  *     charge correlation
1510  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1511    
1512  *     -------------------------------------------------------------                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1513  *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<       $              .and.
1514  *     -------------------------------------------------------------       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1515  c$$$               if(dedx(icy).lt.chsaty.or.dedx(icx).lt.chsatx)then       $              .and.
1516  c$$$                  ddd=(dedx(icy)       $              (badclx.eq.1.and.badcly.eq.1)
1517  c$$$     $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $              .and.
1518  c$$$                  ddd=ddd/sqrt(kch(nplx,nldx)**2+1)       $              .true.)then
1519  c$$$                  cut=chcut*sch(nplx,nldx)  
1520  c$$$                  if(abs(ddd).gt.cut)goto 20 !charge not consistent                    ddd=(sgnl(icy)
1521  c$$$               endif       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1522                                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1523  *     ------------------> COUPLE <------------------  
1524  *     check to do not overflow vector dimentions  c                  cut = chcut * sch(nplx,nldx)
1525                 if(ncp_plane(nplx).gt.ncouplemax)then  
1526                    if(DEBUG)print*,                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1527       $                    ' ** warning ** number of identified'//       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1528       $                    ' couples on plane ',nplx,                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1529       $                    ' exceeds vector dimention'//                    cut = chcut * (16 + sss/50.)
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
                 
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
   
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       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  
1530    
1531  *     output flag                    if(abs(ddd).gt.cut)then
1532  *     --------------                       goto 20    !charge not consistent
1533  *     0 = good event                    endif
1534  *     1 = bad event                 endif
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
1535    
 *     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  
1536                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1537                    if(DEBUG)print*,                    if(verbose)print*,
1538       $                    ' ** warning ** number of identified'//       $                 '** warning ** number of identified '//
1539       $                    ' couples on plane ',nplx,       $                 'couples on plane ',nplx,
1540       $                    ' exceeds vector dimention'//       $                 'exceeds vector dimention '
1541       $                    ' ( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1542  c     good2=.false.                    mask_view(nviewx(nplx)) = 2
1543  c     goto 880   !fill ntp and go to next event                    mask_view(nviewy(nply)) = 2
1544                    iflag=1                    goto 10
                   return  
1545                 endif                 endif
1546                                
1547  c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  *     ------------------> COUPLE <------------------
 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  
                 
1548                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1549                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1550                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1551                 cl_single(icx)=0                 cl_single(icx)=0
1552                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1553  *     ----------------------------------------------  *     ----------------------------------------------
1554    
1555                endif                              
1556    
1557   20         continue   20         continue
1558           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1559                    
# Line 2083  c$$$               endif Line 1578  c$$$               endif
1578        endif        endif
1579                
1580        do ip=1,6        do ip=1,6
1581           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1582        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1583                
1584        return        return
1585        end        end
   
 c$$$      subroutine cl_to_couples_2(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$      include 'level1.f'  
 c$$$  
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$*         print*,'icx ',icx,badcl  
 c$$$         if(badcl.eq.0)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$*            print*,'icy ',icy,badcl  
 c$$$            if(badcl.eq.0)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$  
 c$$$c$$$*     charge correlation  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(DEBUG)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
1586                
1587  ***************************************************  ***************************************************
1588  *                                                 *  *                                                 *
# Line 2321  c$$$      end Line 1594  c$$$      end
1594  **************************************************  **************************************************
1595    
1596        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1597    
1598        include 'commontracker.f'        include 'commontracker.f'
1599          include 'level1.f'
1600        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1601        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1602        include 'common_mini_2.f'        include 'common_mini_2.f'
1603        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1604    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1605    
1606  *     output flag  *     output flag
1607  *     --------------  *     --------------
# Line 2367  c      double precision xm3,ym3,zm3 Line 1634  c      double precision xm3,ym3,zm3
1634  *     -----------------------------  *     -----------------------------
1635    
1636    
1637    *     --------------------------------------------
1638    *     put a limit to the maximum number of couples
1639    *     per plane, in order to apply hough transform
1640    *     (couples recovered during track refinement)
1641    *     --------------------------------------------
1642          do ip=1,nplanes
1643             if(ncp_plane(ip).gt.ncouplelimit)then
1644                mask_view(nviewx(ip)) = 8
1645                mask_view(nviewy(ip)) = 8
1646             endif
1647          enddo
1648    
1649    
1650        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1651        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1652                
1653        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1654           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1655                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1656             do is1=1,2             !loop on sensors - COPPIA 1            
1657              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1658                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1659                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1660  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1661                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1662                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1663                 xm1=xPAM                 xm1=xPAM
1664                 ym1=yPAM                 ym1=yPAM
1665                 zm1=zPAM                                   zm1=zPAM                  
1666  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1667    
1668                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1669                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1670         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1671                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1672                                            
1673                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2391  c     print*,'***',is1,xm1,ym1,zm1 Line 1675  c     print*,'***',is1,xm1,ym1,zm1
1675                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1676  c                        call xyz_PAM  c                        call xyz_PAM
1677  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1678    c                        call xyz_PAM
1679    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1680                          call xyz_PAM                          call xyz_PAM
1681       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1682                          xm2=xPAM                          xm2=xPAM
1683                          ym2=yPAM                          ym2=yPAM
1684                          zm2=zPAM                          zm2=zPAM
# Line 2402  c     $                       (icx2,icy2 Line 1688  c     $                       (icx2,icy2
1688  *     (2 couples needed)  *     (2 couples needed)
1689  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1690                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1691                             if(DEBUG)print*,                             if(verbose)print*,
1692       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1693       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1694       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1695  c                           good2=.false.  c                           good2=.false.
1696  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1697                               do iv=1,12
1698                                  mask_view(iv) = 3
1699                               enddo
1700                             iflag=1                             iflag=1
1701                             return                             return
1702                          endif                          endif
# Line 2441  c$$$ Line 1730  c$$$
1730  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1731    
1732    
1733                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1734    
1735                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1736                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1737         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1738                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1739                                                                
1740                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2450  c$$$ Line 1742  c$$$
1742                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
1743  c                                 call xyz_PAM  c                                 call xyz_PAM
1744  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
1745    c                                 call xyz_PAM
1746    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
1747                                   call xyz_PAM                                   call xyz_PAM
1748       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
1749         $                                ,0.,0.,0.,0.)
1750                                   xm3=xPAM                                   xm3=xPAM
1751                                   ym3=yPAM                                   ym3=yPAM
1752                                   zm3=zPAM                                   zm3=zPAM
# Line 2472  c     $                                 Line 1767  c     $                                
1767  *     (3 couples needed)  *     (3 couples needed)
1768  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1769                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
1770                                      if(DEBUG)print*,                                      if(verbose)print*,
1771       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
1772       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
1773       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
1774  c                                    good2=.false.  c                                    good2=.false.
1775  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
1776                                        do iv=1,nviews
1777                                           mask_view(iv) = 4
1778                                        enddo
1779                                      iflag=1                                      iflag=1
1780                                      return                                      return
1781                                   endif                                   endif
# Line 2516  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1814  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1814                                endif                                endif
1815                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1816                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1817     30                     continue
1818                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1819   30                  continue   31                  continue
1820                                            
1821   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1822                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1823     20            continue
1824              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1825                            
1826           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1827        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1828     10   continue
1829        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1830                
1831        if(DEBUG)then        if(DEBUG)then
# Line 2552  c     goto 880               !ntp fill Line 1853  c     goto 880               !ntp fill
1853        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
1854    
1855        include 'commontracker.f'        include 'commontracker.f'
1856          include 'level1.f'
1857        include 'common_momanhough.f'        include 'common_momanhough.f'
1858        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1859    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1860    
1861  *     output flag  *     output flag
1862  *     --------------  *     --------------
# Line 2588  c      common/dbg/DEBUG Line 1888  c      common/dbg/DEBUG
1888        distance=0        distance=0
1889        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
1890        npt_tot=0        npt_tot=0
1891          nloop=0                  
1892     90   continue                  
1893        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
1894           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
1895                            
# Line 2691  c     print*,'*   idbref,idb2 ',idbref,i Line 1993  c     print*,'*   idbref,idb2 ',idbref,i
1993              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
1994           enddo           enddo
1995  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
1996           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
1997           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
1998           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
1999                    
# Line 2699  c     print*,'>>>> ',ncpused,npt,nplused Line 2001  c     print*,'>>>> ',ncpused,npt,nplused
2001  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2002    
2003           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2004              if(DEBUG)print*,              if(verbose)print*,
2005       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2006       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2007       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2008  c               good2=.false.  c               good2=.false.
2009  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2010                do iv=1,nviews
2011                   mask_view(iv) = 5
2012                enddo
2013              iflag=1              iflag=1
2014              return              return
2015           endif           endif
# Line 2742  c$$$     $           ,(db_cloud(iii),iii Line 2047  c$$$     $           ,(db_cloud(iii),iii
2047        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2048                
2049                
2050          if(nloop.lt.nstepy)then      
2051            cutdistyz = cutdistyz+cutystep
2052            nloop     = nloop+1          
2053            goto 90                
2054          endif                    
2055          
2056        if(DEBUG)then        if(DEBUG)then
2057           print*,'---------------------- '           print*,'---------------------- '
2058           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
# Line 2768  c$$$     $           ,(db_cloud(iii),iii Line 2079  c$$$     $           ,(db_cloud(iii),iii
2079        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2080    
2081        include 'commontracker.f'        include 'commontracker.f'
2082          include 'level1.f'
2083        include 'common_momanhough.f'        include 'common_momanhough.f'
2084        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2085    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2086    
2087  *     output flag  *     output flag
2088  *     --------------  *     --------------
# Line 2803  c      common/dbg/DEBUG Line 2113  c      common/dbg/DEBUG
2113        distance=0        distance=0
2114        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2115        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2116          nloop=0                  
2117     91   continue                  
2118        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2119           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
2120  c     print*,'--------------'  c     print*,'--------------'
# Line 2904  c     print*,'check cp_used' Line 2216  c     print*,'check cp_used'
2216           do ip=1,nplanes           do ip=1,nplanes
2217              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2218           enddo           enddo
2219           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2220           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2221           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next doublet
2222                    
2223  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2224  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2225           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2226              if(DEBUG)print*,              if(verbose)print*,
2227       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2228       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2229       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2230  c     good2=.false.  c     good2=.false.
2231  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2232                do iv=1,nviews
2233                   mask_view(iv) = 6
2234                enddo
2235              iflag=1              iflag=1
2236              return              return
2237           endif           endif
# Line 2952  c$$$     $           ,(tr_cloud(iii),iii Line 2267  c$$$     $           ,(tr_cloud(iii),iii
2267  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2268  22288    continue  22288    continue
2269        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2270          
2271           if(nloop.lt.nstepx)then      
2272             cutdistxz=cutdistxz+cutxstep
2273             nloop=nloop+1          
2274             goto 91                
2275           endif                    
2276          
2277        if(DEBUG)then        if(DEBUG)then
2278           print*,'---------------------- '           print*,'---------------------- '
2279           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
# Line 2979  c     02/02/2006 modified by Elena Vannu Line 2300  c     02/02/2006 modified by Elena Vannu
2300  c*****************************************************  c*****************************************************
2301    
2302        include 'commontracker.f'        include 'commontracker.f'
2303          include 'level1.f'
2304        include 'common_momanhough.f'        include 'common_momanhough.f'
2305        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2306        include 'common_mini_2.f'        include 'common_mini_2.f'
2307        include 'common_mech.f'        include 'common_mech.f'
2308        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2309    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2310    
2311  *     output flag  *     output flag
2312  *     --------------  *     --------------
# Line 3002  c      common/dbg/DEBUG Line 2322  c      common/dbg/DEBUG
2322  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2323  *     list of matching couples in the combination  *     list of matching couples in the combination
2324  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2325        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2326        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2327  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2328        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2329  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2330  *     variables for track fitting  *     variables for track fitting
2331        double precision AL_INI(5)        double precision AL_INI(5)
2332        double precision tath  c      double precision tath
2333  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2334  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2335    
# Line 3075  c      real fitz(nplanes)        !z coor Line 2395  c      real fitz(nplanes)        !z coor
2395                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2396              enddo              enddo
2397                            
2398              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2399                if(nplused.lt.nplyz_min)goto 888 !next doublet
2400              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2401                            
2402              if(DEBUG)then              if(DEBUG)then
# Line 3102  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2423  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2423  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2424                            
2425  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2426              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2427              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2428              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2429       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2430              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2431              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2432              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2433                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2434  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2435              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2436                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2437  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2438  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2439                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2440  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2441  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2442                c$$$            ELSE
2443    c$$$               AL_INI(4) = acos(-1.)/2
2444    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2445    c$$$            ENDIF
2446    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2447    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2448    c$$$            
2449    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2450    c$$$            
2451    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2452                            
2453              if(DEBUG)then              if(DEBUG)then
2454                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2455                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3169  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2500  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2500  *                                   *************************  *                                   *************************
2501  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2502  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2503    c                                    call xyz_PAM(icx,icy,is, !(1)
2504    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2505                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2506       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2507  *                                   *************************  *                                   *************************
2508  *                                   -----------------------------  *                                   -----------------------------
2509                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 3186  c     $                                 Line 2519  c     $                                
2519  *     **********************************************************  *     **********************************************************
2520  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2521  *     **********************************************************  *     **********************************************************
2522    cccc  scommentare se si usa al_ini della nuvola
2523    c$$$                              do i=1,5
2524    c$$$                                 AL(i)=AL_INI(i)
2525    c$$$                              enddo
2526                                  call guess()
2527                                do i=1,5                                do i=1,5
2528                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2529                                enddo                                enddo
2530                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2531                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2532                                call mini_2(jstep,ifail)                                iprint=0
2533    c                              if(DEBUG)iprint=1
2534                                  if(DEBUG)iprint=2
2535                                  call mini2(jstep,ifail,iprint)
2536                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2537                                   if(DEBUG)then                                   if(DEBUG)then
2538                                      print *,                                      print *,
2539       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2540       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2541                                        print*,'initial guess: '
2542    
2543                                        print*,'AL_INI(1) = ',AL_INI(1)
2544                                        print*,'AL_INI(2) = ',AL_INI(2)
2545                                        print*,'AL_INI(3) = ',AL_INI(3)
2546                                        print*,'AL_INI(4) = ',AL_INI(4)
2547                                        print*,'AL_INI(5) = ',AL_INI(5)
2548                                   endif                                   endif
2549                                   chi2=-chi2  c                                 chi2=-chi2
2550                                endif                                endif
2551  *     **********************************************************  *     **********************************************************
2552  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3211  c     $                                 Line 2559  c     $                                
2559  *     --------------------------  *     --------------------------
2560                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2561                                                                    
2562                                   if(DEBUG)print*,                                   if(verbose)print*,
2563       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2564       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2565       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2566  c                                 good2=.false.  c                                 good2=.false.
2567  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2568                                     do iv=1,nviews
2569                                        mask_view(iv) = 7
2570                                     enddo
2571                                   iflag=1                                   iflag=1
2572                                   return                                   return
2573                                endif                                endif
# Line 3315  cccccc 12/08/2006 modified by elena vann Line 2666  cccccc 12/08/2006 modified by elena vann
2666  c******************************************************  c******************************************************
2667    
2668        include 'commontracker.f'        include 'commontracker.f'
2669          include 'level1.f'
2670        include 'common_momanhough.f'        include 'common_momanhough.f'
2671        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2672        include 'common_mini_2.f'        include 'common_mini_2.f'
2673        include 'common_mech.f'        include 'common_mech.f'
2674        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2675        include 'level1.f'  c      include 'level1.f'
2676        include 'calib.f'        include 'calib.f'
2677    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
2678  *     flag to chose PFA  *     flag to chose PFA
2679        character*10 PFA        character*10 PFA
2680        common/FINALPFA/PFA        common/FINALPFA/PFA
2681    
2682          real xp,yp,zp
2683          real xyzp(3),bxyz(3)
2684          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
2685    
2686  *     =================================================  *     =================================================
2687  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
2688  *                          and  *                          and
# Line 3338  c      common/dbg/DEBUG Line 2691  c      common/dbg/DEBUG
2691        call track_init        call track_init
2692        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2693    
2694             xP=XV_STORE(nplanes-ip+1,ibest)
2695             yP=YV_STORE(nplanes-ip+1,ibest)
2696             zP=ZV_STORE(nplanes-ip+1,ibest)
2697             call gufld(xyzp,bxyz)
2698    c$$$         bxyz(1)=0
2699    c$$$         bxyz(2)=0
2700    c$$$         bxyz(3)=0
2701    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2702  *     -------------------------------------------------  *     -------------------------------------------------
2703  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2704  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2705  *     using improved PFAs  *     using improved PFAs
2706  *     -------------------------------------------------  *     -------------------------------------------------
2707    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2708           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2709       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2710                            
# Line 3356  c      common/dbg/DEBUG Line 2718  c      common/dbg/DEBUG
2718       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
2719              icx=clx(ip,icp)              icx=clx(ip,icp)
2720              icy=cly(ip,icp)              icy=cly(ip,icp)
2721    c            call xyz_PAM(icx,icy,is,
2722    c     $           PFA,PFA,
2723    c     $           AXV_STORE(nplanes-ip+1,ibest),
2724    c     $           AYV_STORE(nplanes-ip+1,ibest))
2725              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
2726       $           PFA,PFA,       $           PFA,PFA,
2727       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
2728       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
2729         $           bxyz(1),
2730         $           bxyz(2)
2731         $           )
2732  c$$$  call xyz_PAM(icx,icy,is,  c$$$  call xyz_PAM(icx,icy,is,
2733  c$$$  $              'COG2','COG2',  c$$$  $              'COG2','COG2',
2734  c$$$  $              0.,  c$$$  $              0.,
# Line 3373  c$$$  $              0.) Line 2741  c$$$  $              0.)
2741              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
2742              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
2743    
2744  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)  c            dedxtrk(nplanes-ip+1) = (sgnl(icx)+sgnl(icy))/2. !(1)
2745              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2746              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2747                            
2748    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2749  *     -------------------------------------------------  *     -------------------------------------------------
2750  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
2751  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
2752  *     -------------------------------------------------  *     -------------------------------------------------
2753    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2754           else                             else                  
2755                                
2756              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3388  c            dedxtrk(nplanes-ip+1) = (de Line 2758  c            dedxtrk(nplanes-ip+1) = (de
2758                                
2759  *     --------------------------------------------------------------  *     --------------------------------------------------------------
2760  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
2761              xP=XV_STORE(nplanes-ip+1,ibest)  c$$$            xP=XV_STORE(nplanes-ip+1,ibest)
2762              yP=YV_STORE(nplanes-ip+1,ibest)  c$$$            yP=YV_STORE(nplanes-ip+1,ibest)
2763              zP=ZV_STORE(nplanes-ip+1,ibest)  c$$$            zP=ZV_STORE(nplanes-ip+1,ibest)
2764              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
2765  *     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
2766              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
# Line 3428  c     $              cl_used(icy).eq.1.o Line 2798  c     $              cl_used(icy).eq.1.o
2798       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
2799       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
2800  *            *          
2801    c               call xyz_PAM(icx,icy,ist,
2802    c     $              PFA,PFA,
2803    c     $              AXV_STORE(nplanes-ip+1,ibest),
2804    c     $              AYV_STORE(nplanes-ip+1,ibest))
2805                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
2806       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
2807       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
2808       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
2809         $              bxyz(1),
2810         $              bxyz(2)
2811         $              )
2812                                
2813                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2814                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2815                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2816                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2817       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3446  c     $              'ETA2','ETA2', Line 2823  c     $              'ETA2','ETA2',
2823                    rymm = resyPAM                    rymm = resyPAM
2824                    distmin = distance                    distmin = distance
2825                    idm = id                                      idm = id                  
2826  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)  c                 dedxmm = (sgnl(icx)+sgnl(icy))/2. !(1)
2827                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2828                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2829                 endif                 endif
2830   1188          continue   1188          continue
2831              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
# Line 3500  c            if(DEBUG)print*,'>>>> try t Line 2877  c            if(DEBUG)print*,'>>>> try t
2877  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
2878                 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)
2879  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
2880    c               call xyz_PAM(icx,0,ist,
2881    c     $              PFA,PFA,
2882    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
2883                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
2884       $              PFA,PFA,       $              PFA,PFA,
2885       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
2886         $              bxyz(1),
2887         $              bxyz(2)
2888         $              )              
2889                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2890  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2891                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
2892       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2893                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3519  c     if(DEBUG)print*,'normalized distan Line 2901  c     if(DEBUG)print*,'normalized distan
2901                    rymm = resyPAM                    rymm = resyPAM
2902                    distmin = distance                    distmin = distance
2903                    iclm = icx                    iclm = icx
2904  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
2905                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2906                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
2907                 endif                                   endif                  
2908  11881          continue  11881          continue
# Line 3528  c                  dedxmm = dedx(icx) !( Line 2910  c                  dedxmm = dedx(icx) !(
2910  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
2911                 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)
2912  *                                              !jump to the next couple  *                                              !jump to the next couple
2913    c               call xyz_PAM(0,icy,ist,
2914    c     $              PFA,PFA,
2915    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
2916                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
2917       $              PFA,PFA,       $              PFA,PFA,
2918       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
2919         $              bxyz(1),
2920         $              bxyz(2)
2921         $              )
2922                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2923                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2924                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
2925       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2926                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3546  c     $              'ETA2','ETA2', Line 2934  c     $              'ETA2','ETA2',
2934                    rymm = resyPAM                    rymm = resyPAM
2935                    distmin = distance                    distmin = distance
2936                    iclm = icy                    iclm = icy
2937  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
2938                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
2939                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2940                 endif                                   endif                  
2941  11882          continue  11882          continue
2942              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
2943  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
2944    c            print*,'## ncls(',ip,') ',ncls(ip)
2945              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
2946                 icl=cls(ip,ic)                 icl=cls(ip,ic)
2947    c              print*,'## ic ',ic,' ist ',ist
2948  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
2949                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
2950       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
2951       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
2952                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
2953    c                  call xyz_PAM(icl,0,ist,
2954    c     $                 PFA,PFA,
2955    c     $                 AXV_STORE(nplanes-ip+1,ibest),0.)
2956                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
2957       $                 PFA,PFA,       $                 PFA,PFA,
2958       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
2959         $                 bxyz(1),
2960         $                 bxyz(2)
2961         $                 )
2962                 else                         !<---- Y view                 else                         !<---- Y view
2963    c                  call xyz_PAM(0,icl,ist,
2964    c     $                 PFA,PFA,
2965    c     $                 0.,AYV_STORE(nplanes-ip+1,ibest))
2966                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
2967       $                 PFA,PFA,       $                 PFA,PFA,
2968       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
2969         $                 bxyz(1),
2970         $                 bxyz(2)
2971         $                 )
2972                 endif                 endif
2973    
2974                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2975                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2976                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
2977       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance,'<',distmin,' ?'
2978                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
2979                      if(DEBUG)print*,'YES'
2980                    xmm_A = xPAM_A                    xmm_A = xPAM_A
2981                    ymm_A = yPAM_A                    ymm_A = yPAM_A
2982                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3585  c     $                 'ETA2','ETA2', Line 2987  c     $                 'ETA2','ETA2',
2987                    rymm = resyPAM                    rymm = resyPAM
2988                    distmin = distance                      distmin = distance  
2989                    iclm = icl                    iclm = icl
2990  c                  dedxmm = dedx(icl)                   !(1)  c                  dedxmm = sgnl(icl)                   !(1)
2991                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
2992                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
2993                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                       !(1)
2994                    else          !<---- Y view                    else          !<---- Y view
2995                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                       !(1)
2996                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
2997                    endif                    endif
2998                 endif                                   endif                  
2999  18882          continue  18882          continue
3000              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3001    c            print*,'## distmin ', distmin,' clinc ',clinc
3002              if(distmin.le.clinc)then                                if(distmin.le.clinc)then                  
3003                                
3004                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3005  *              ----------------------------  *              ----------------------------
3006    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3007                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3008                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
3009                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
3010                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
3011       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
3012         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
3013         $                 ,', norm.dist.= ',distmin
3014         $                 ,', cut ',clinc,' )'
3015                 else                 else
3016                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
3017                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
3018                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
3019       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
3020         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
3021         $                 ,', norm.dist.= ', distmin
3022         $                 ,', cut ',clinc,' )'
3023                 endif                 endif
3024    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3025  *              ----------------------------  *              ----------------------------
3026                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
3027                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3645  cccccc 12/08/2006 modified by elena ---> Line 3055  cccccc 12/08/2006 modified by elena --->
3055        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3056    
3057        include 'commontracker.f'        include 'commontracker.f'
3058          include 'level1.f'
3059        include 'common_momanhough.f'        include 'common_momanhough.f'
3060        include 'momanhough_init.f'  c      include 'momanhough_init.f'
3061        include 'level2.f'        !(1)        include 'level2.f'        !(1)
3062  c      include 'calib.f'  c      include 'calib.f'
3063  c      include 'level1.f'  c      include 'level1.f'
3064    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
3065    
3066    
3067        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
# Line 3725  c               endif Line 3134  c               endif
3134    
3135    
3136    
 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$$$  
   
3137    
3138    
3139  *     ****************************************************  *     ****************************************************
3140    
3141        subroutine init_level2        subroutine init_level2
3142    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3143        include 'commontracker.f'        include 'commontracker.f'
3144          include 'level1.f'
3145        include 'common_momanhough.f'        include 'common_momanhough.f'
3146        include 'level2.f'        include 'level2.f'
3147        include 'level1.f'  c      include 'level1.f'
3148    
3149        do i=1,nviews        do i=1,nviews
3150           good2(i)=good1(i)           good2(i)=good1(i)
3151        enddo        enddo
3152    
 c      good2 = 0!.false.  
 c$$$      nev2 = nev1  
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
3153    
3154        NTRK = 0        NTRK = 0
3155        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3156           IMAGE(IT)=0           IMAGE(IT)=0
3157           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3158           do ip=1,nplanes           do ip=1,nplanes
3159              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3160              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
# Line 3877  c         BdL(IT) = 0. Line 3163  c         BdL(IT) = 0.
3163              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3164              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3165              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3166              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3167              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3168              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3169              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3170           enddo           enddo
# Line 3893  cccccc 17/8/2006 modified by elena Line 3175  cccccc 17/8/2006 modified by elena
3175              enddo                                enddo                  
3176           enddo                             enddo                  
3177        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3178        nclsx=0        nclsx=0
3179        nclsy=0              nclsy=0      
3180        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3181          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3182          xs(1,ip)=0          xs(1,ip)=0
3183          xs(2,ip)=0          xs(2,ip)=0
3184          sgnlxs(ip)=0          sgnlxs(ip)=0
3185          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3186          ys(1,ip)=0          ys(1,ip)=0
3187          ys(2,ip)=0          ys(2,ip)=0
3188          sgnlys(ip)=0          sgnlys(ip)=0
3189        enddo        enddo
 c*******************************************************  
3190        end        end
3191    
3192    
# Line 3926  c*************************************** Line 3201  c***************************************
3201  ************************************************************  ************************************************************
3202    
3203    
3204          subroutine init_hough
3205    
3206          include 'commontracker.f'
3207          include 'level1.f'
3208          include 'common_momanhough.f'
3209          include 'common_hough.f'
3210          include 'level2.f'
3211    
3212          ntrpt_nt=0
3213          ndblt_nt=0
3214          NCLOUDS_XZ_nt=0
3215          NCLOUDS_YZ_nt=0
3216          do idb=1,ndblt_max_nt
3217             db_cloud_nt(idb)=0
3218             alfayz1_nt(idb)=0      
3219             alfayz2_nt(idb)=0      
3220          enddo
3221          do itr=1,ntrpt_max_nt
3222             tr_cloud_nt(itr)=0
3223             alfaxz1_nt(itr)=0      
3224             alfaxz2_nt(itr)=0      
3225             alfaxz3_nt(itr)=0      
3226          enddo
3227          do idb=1,ncloyz_max      
3228            ptcloud_yz_nt(idb)=0    
3229            alfayz1_av_nt(idb)=0    
3230            alfayz2_av_nt(idb)=0    
3231          enddo                    
3232          do itr=1,ncloxz_max      
3233            ptcloud_xz_nt(itr)=0    
3234            alfaxz1_av_nt(itr)=0    
3235            alfaxz2_av_nt(itr)=0    
3236            alfaxz3_av_nt(itr)=0    
3237          enddo                    
3238    
3239          ntrpt=0                  
3240          ndblt=0                  
3241          NCLOUDS_XZ=0              
3242          NCLOUDS_YZ=0              
3243          do idb=1,ndblt_max        
3244            db_cloud(idb)=0        
3245            cpyz1(idb)=0            
3246            cpyz2(idb)=0            
3247            alfayz1(idb)=0          
3248            alfayz2(idb)=0          
3249          enddo                    
3250          do itr=1,ntrpt_max        
3251            tr_cloud(itr)=0        
3252            cpxz1(itr)=0            
3253            cpxz2(itr)=0            
3254            cpxz3(itr)=0            
3255            alfaxz1(itr)=0          
3256            alfaxz2(itr)=0          
3257            alfaxz3(itr)=0          
3258          enddo                    
3259          do idb=1,ncloyz_max      
3260            ptcloud_yz(idb)=0      
3261            alfayz1_av(idb)=0      
3262            alfayz2_av(idb)=0      
3263            do idbb=1,ncouplemaxtot
3264              cpcloud_yz(idb,idbb)=0
3265            enddo                  
3266          enddo                    
3267          do itr=1,ncloxz_max      
3268            ptcloud_xz(itr)=0      
3269            alfaxz1_av(itr)=0      
3270            alfaxz2_av(itr)=0      
3271            alfaxz3_av(itr)=0      
3272            do itrr=1,ncouplemaxtot
3273              cpcloud_xz(itr,itrr)=0
3274            enddo                  
3275          enddo                    
3276          end
3277    ************************************************************
3278    *
3279    *
3280    *
3281    *
3282    *
3283    *
3284    *
3285    ************************************************************
3286    
3287    
3288        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3289    
3290  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3936  c*************************************** Line 3295  c***************************************
3295    
3296            
3297        include 'commontracker.f'        include 'commontracker.f'
3298    c      include 'level1.f'
3299        include 'level1.f'        include 'level1.f'
3300          include 'common_momanhough.f'
3301        include 'level2.f'        include 'level2.f'
3302        include 'common_mini_2.f'        include 'common_mini_2.f'
3303        include 'common_momanhough.f'        real sinth,phi,pig      
       real sinth,phi,pig        !(4)  
3304        pig=acos(-1.)        pig=acos(-1.)
3305    
 c      good2=1!.true.  
3306        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3307        nstep_nt(ntr)       = 0!nstep        nstep_nt(ntr)       = nstep
3308    
3309          phi   = al(4)          
3310          sinth = al(3)            
3311          if(sinth.lt.0)then      
3312             sinth = -sinth        
3313             phi = phi + pig      
3314          endif                    
3315          npig = aint(phi/(2*pig))
3316          phi = phi - npig*2*pig  
3317          if(phi.lt.0)            
3318         $     phi = phi + 2*pig  
3319          al(4) = phi              
3320          al(3) = sinth            
3321    
       phi   = al(4)             !(4)  
       sinth = al(3)             !(4)  
       if(sinth.lt.0)then        !(4)  
          sinth = -sinth         !(4)  
          phi = phi + pig        !(4)  
       endif                     !(4)  
       npig = aint(phi/(2*pig))  !(4)  
       phi = phi - npig*2*pig    !(4)  
       if(phi.lt.0)              !(4)  
      $     phi = phi + 2*pig    !(4)  
       al(4) = phi               !(4)  
       al(3) = sinth             !(4)  
 *****************************************************  
3322        do i=1,5        do i=1,5
3323           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3324           do j=1,5           do j=1,5
3325              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3326           enddo           enddo
 c     print*,al_nt(i,ntr)  
3327        enddo        enddo
3328                
3329        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
# Line 3980  c     print*,al_nt(i,ntr) Line 3338  c     print*,al_nt(i,ntr)
3338           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3339           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3340           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3341           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3342  c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3343           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           factor = sqrt(
3344           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)         $        sin( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3345         $        sin( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3346         $        1. )
3347    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3348             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3349             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3350        
3351           id  = CP_STORE(ip,IDCAND)           id  = CP_STORE(ip,IDCAND)
3352           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
# Line 3998  c            print*,ip,' ',cltrx(ip,ntr) Line 3361  c            print*,ip,' ',cltrx(ip,ntr)
3361           endif                     endif          
3362    
3363        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)  
3364    
3365    
3366        end        end
3367    
3368        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*****************************************************  
3369    
3370  *     -------------------------------------------------------  *     -------------------------------------------------------
3371  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 4022  c*************************************** Line 3374  c***************************************
3374  *     -------------------------------------------------------  *     -------------------------------------------------------
3375    
3376        include 'commontracker.f'        include 'commontracker.f'
3377        include 'level1.f'  c      include 'level1.f'
       include 'level2.f'  
3378        include 'calib.f'        include 'calib.f'
3379          include 'level1.f'
3380        include 'common_momanhough.f'        include 'common_momanhough.f'
3381          include 'level2.f'
3382        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3383    
3384  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
# Line 4033  c      good2=1!.true. Line 3386  c      good2=1!.true.
3386        nclsx = 0        nclsx = 0
3387        nclsy = 0        nclsy = 0
3388    
3389          do iv = 1,nviews
3390             if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3391          enddo
3392    
3393        do icl=1,nclstr1        do icl=1,nclstr1
3394           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
3395              ip=nplanes-npl(VIEW(icl))+1                          ip=nplanes-npl(VIEW(icl))+1            
3396              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3397                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3398                 planex(nclsx) = ip                 planex(nclsx) = ip
3399                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3400                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3401                 do is=1,2                 do is=1,2
3402  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3403                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3404                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3405                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3406                 enddo                 enddo
3407  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 4054  c$$$               print*,'xs(2,nclsx)   Line 3412  c$$$               print*,'xs(2,nclsx)  
3412              else                          !=== Y views              else                          !=== Y views
3413                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3414                 planey(nclsy) = ip                 planey(nclsy) = ip
3415                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3416                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3417                 do is=1,2                 do is=1,2
3418  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3419                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3420                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3421                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3422                 enddo                 enddo
3423  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 4076  c      print*,icl,cl_used(icl),cl_good(i Line 3435  c      print*,icl,cl_used(icl),cl_good(i
3435        enddo        enddo
3436        end        end
3437    
3438    ***************************************************
3439    *                                                 *
3440    *                                                 *
3441    *                                                 *
3442    *                                                 *
3443    *                                                 *
3444    *                                                 *
3445    **************************************************
3446    
3447          subroutine fill_hough
3448    
3449    *     -------------------------------------------------------
3450    *     This routine fills the  variables related to the hough
3451    *     transform, for the debig n-tuple
3452    *     -------------------------------------------------------
3453    
3454          include 'commontracker.f'
3455          include 'level1.f'
3456          include 'common_momanhough.f'
3457          include 'common_hough.f'
3458          include 'level2.f'
3459    
3460          if(.false.
3461         $     .or.ntrpt.gt.ntrpt_max_nt
3462         $     .or.ndblt.gt.ndblt_max_nt
3463         $     .or.NCLOUDS_XZ.gt.ncloxz_max
3464         $     .or.NCLOUDS_yZ.gt.ncloyz_max
3465         $     )then
3466             ntrpt_nt=0
3467             ndblt_nt=0
3468             NCLOUDS_XZ_nt=0
3469             NCLOUDS_YZ_nt=0
3470          else
3471             ndblt_nt=ndblt
3472             ntrpt_nt=ntrpt
3473             if(ndblt.ne.0)then
3474                do id=1,ndblt
3475                   alfayz1_nt(id)=alfayz1(id) !Y0
3476                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3477                enddo
3478             endif
3479             if(ndblt.ne.0)then
3480                do it=1,ntrpt
3481                   alfaxz1_nt(it)=alfaxz1(it) !X0
3482                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3483                   alfaxz3_nt(it)=alfaxz3(it) !1/r
3484                enddo
3485             endif
3486             nclouds_yz_nt=nclouds_yz
3487             nclouds_xz_nt=nclouds_xz
3488             if(nclouds_yz.ne.0)then
3489                nnn=0
3490                do iyz=1,nclouds_yz
3491                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3492                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3493                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3494                   nnn=nnn+ptcloud_yz(iyz)
3495                enddo
3496                do ipt=1,nnn
3497                   db_cloud_nt(ipt)=db_cloud(ipt)
3498                 enddo
3499             endif
3500             if(nclouds_xz.ne.0)then
3501                nnn=0
3502                do ixz=1,nclouds_xz
3503                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
3504                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
3505                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
3506                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
3507                   nnn=nnn+ptcloud_xz(ixz)              
3508                enddo
3509                do ipt=1,nnn
3510                  tr_cloud_nt(ipt)=tr_cloud(ipt)
3511                 enddo
3512             endif
3513          endif
3514          end
3515          
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.23