/[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.8 by pam-fi, Wed Oct 25 16:18:41 2006 UTC revision 1.42 by pam-fi, Tue Nov 10 15:50:16 2009 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      print*,'======================================================'
24    c$$$      do ic=1,NCLSTR1
25    c$$$         if(.false.
26    c$$$     $        .or.nsatstrips(ic).gt.0
27    c$$$c     $        .or.nbadstrips(0,ic).gt.0
28    c$$$c     $        .or.nbadstrips(4,ic).gt.0
29    c$$$c     $        .or.nbadstrips(3,ic).gt.0
30    c$$$     $        .or..false.)then
31    c$$$            print*,'--- cl-',ic,' ------------------------'
32    c$$$            istart = INDSTART(IC)
33    c$$$            istop  = TOTCLLENGTH
34    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
35    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
36    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
37    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
38    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
39    c$$$            print*,'view ',VIEW(ic)
40    c$$$            print*,'maxs ',MAXS(ic)
41    c$$$            print*,'COG4 ',cog(4,ic)
42    c$$$            ff = fbad_cog(4,ic)
43    c$$$            print*,'fbad ',ff
44    c$$$            print*,(CLBAD(i),i=istart,istop)
45    c$$$            bb=nbadstrips(0,ic)
46    c$$$            print*,'#BAD (tot)',bb
47    c$$$            bb=nbadstrips(4,ic)
48    c$$$            print*,'#BAD (4)',bb
49    c$$$            bb=nbadstrips(3,ic)
50    c$$$            print*,'#BAD (3)',bb
51    c$$$            ss=nsatstrips(ic)
52    c$$$            print*,'#saturated ',ss
53    c$$$         endif
54    c$$$      enddo
55                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
56  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
57  *     STEP 1  *     STEP 1
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 44  c      common/dbg/DEBUG Line 72  c      common/dbg/DEBUG
72  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
73  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
74    
 c      iflag=0  
75        call cl_to_couples(iflag)        call cl_to_couples(iflag)
76        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
77           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
78        endif        endif
79                if(ncp_tot.eq.0)goto 880  !go to next event    
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     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  
           
   
   
80  *-----------------------------------------------------  *-----------------------------------------------------
81  *-----------------------------------------------------  *-----------------------------------------------------
82  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 91  c$$$         endif Line 105  c$$$         endif
105  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
106  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
107    
108  c      iflag=0  
109        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
110        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
111           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
112        endif        endif
113          if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event    
114                
115                
116  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 123  c      iflag=0 Line 138  c      iflag=0
138  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
139  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
140  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
141    *     count number of hit planes
142          planehit=0                
143          do np=1,nplanes          
144            if(ncp_plane(np).ne.0)then  
145              planehit=planehit+1  
146            endif                  
147          enddo                    
148          if(planehit.lt.3) goto 880 ! exit              
149          
150          nptxz_min=x_min_start              
151          nplxz_min=x_min_start              
152                
153          nptyz_min=y_min_start              
154          nplyz_min=y_min_start              
155    
156          cutdistyz=cutystart      
157          cutdistxz=cutxstart      
158    
159  c      iflag=0   878  continue
160        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
161        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
162           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
163          endif
164    *     ------------------------------------------------
165    *     first try to release the tolerance
166    *     ------------------------------------------------
167          if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
168             if(cutdistyz.lt.maxcuty/2)then
169                cutdistyz=cutdistyz+cutystep
170             else
171                cutdistyz=cutdistyz+(3*cutystep)
172             endif
173             if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
174             goto 878                
175          endif    
176    *     ------------------------------------------------
177    *     hence reduce the minimum number of plane
178    *     ------------------------------------------------
179          if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
180             nplyz_min=nplyz_min-1
181             if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
182             goto 878
183        endif        endif
184  c      iflag=0  
185          if(nclouds_yz.eq.0)goto 880 !go to next event    
186    
187    
188    ccc   if(planehit.eq.3) goto 881    
189          
190     879  continue  
191        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
192        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
193           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
194        endif        endif
195          *     ------------------------------------------------
196    *     first try to release the tolerance
197    *     ------------------------------------------------                          
198          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
199            cutdistxz=cutdistxz+cutxstep
200             if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz
201            goto 879                
202          endif                    
203    *     ------------------------------------------------
204    *     hence reduce the minimum number of plane
205    *     ------------------------------------------------
206          if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then
207             nplxz_min=nplxz_min-1
208             if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min
209             goto 879
210          endif
211    
212          if(nclouds_xz.eq.0)goto 880 !go to next event    
213    
214        
215    c$$$ 881  continue  
216    c$$$*     if there is at least three planes on the Y view decreases cuts on X view
217    c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
218    c$$$     $     nplxz_min.ne.y_min_start)then
219    c$$$        nptxz_min=x_min_step    
220    c$$$        nplxz_min=x_min_start-x_min_step              
221    c$$$        goto 879                
222    c$$$      endif                    
223            
224   880  return   880  return
225        end        end
226    
# Line 144  c      iflag=0 Line 230  c      iflag=0
230        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
231    
232        include 'commontracker.f'        include 'commontracker.f'
233          include 'level1.f'
234        include 'common_momanhough.f'        include 'common_momanhough.f'
235        include 'common_mech.f'        include 'common_mech.f'
236        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
237        include 'common_mini_2.f'        include 'common_mini_2.f'
238        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
239        include 'level2.f'        include 'level2.f'
240    
241        include 'momanhough_init.f'  c      include 'momanhough_init.f'
242                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
243        logical FIMAGE            !        logical FIMAGE            !
244          real trackimage(NTRACKSMAX)
245          real*8 AL_GUESS(5)
246    
247  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
248  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 184  c      common/dbg/DEBUG Line 269  c      common/dbg/DEBUG
269  *  *
270  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
271  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
272           ntrk=0                 !counter of identified physical tracks  ccc         ntrk=0                 !counter of identified physical tracks
273    
274  11111    continue               !<<<<<<< come here when performing a new search  c11111    continue               !<<<<<<< come here when performing a new search
275             continue                  !<<<<<<< come here when performing a new search
276    
277             if(nclouds_xz.eq.0)goto 880 !go to next event    
278             if(nclouds_yz.eq.0)goto 880 !go to next event    
279    
280  c         iflag=0  c         iflag=0
281           call clouds_to_ctrack(iflag)           call clouds_to_ctrack(iflag)
282           if(iflag.eq.1)then     !no candidate tracks found           if(iflag.eq.1)then     !no candidate tracks found
283              goto 880            !fill ntp and go to next event                goto 880            !fill ntp and go to next event  
284           endif           endif
285             if(ntracks.eq.0)goto 880 !go to next event    
286    
287           FIMAGE=.false.         !processing best track (not track image)           FIMAGE=.false.         !processing best track (not track image)
288           ibest=0                !best track among candidates           ibest=0                !best track among candidates
289           iimage=0               !track image           iimage=0               !track image
290  *     ------------- select the best track -------------  *     ------------- select the best track -------------
291    c$$$         rchi2best=1000000000.
292    c$$$         do i=1,ntracks
293    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
294    c$$$     $         RCHI2_STORE(i).gt.0)then
295    c$$$               ibest=i
296    c$$$               rchi2best=RCHI2_STORE(i)
297    c$$$            endif
298    c$$$         enddo
299    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
300    
301    *     -------------------------------------------------------
302    *     order track-candidates according to:
303    *     1st) decreasing n.points
304    *     2nd) increasing chi**2
305    *     -------------------------------------------------------
306           rchi2best=1000000000.           rchi2best=1000000000.
307             ndofbest=0            
308           do i=1,ntracks           do i=1,ntracks
309              if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0              
310       $         RCHI2_STORE(i).gt.0)then             do ii=1,nplanes    
311                 ndof=ndof        
312         $            +int(xgood_store(ii,i))
313         $            +int(ygood_store(ii,i))
314               enddo              
315               if(ndof.gt.ndofbest)then
316                 ibest=i
317                 rchi2best=RCHI2_STORE(i)
318                 ndofbest=ndof    
319               elseif(ndof.eq.ndofbest)then
320                 if(RCHI2_STORE(i).lt.rchi2best.and.
321         $            RCHI2_STORE(i).gt.0)then
322                 ibest=i                 ibest=i
323                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
324              endif                 ndofbest=ndof  
325                 endif            
326               endif
327           enddo           enddo
328    
329    
330           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
331  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
332  *     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 224  c         iflag=0 Line 345  c         iflag=0
345              iimage=0              iimage=0
346           endif           endif
347           if(icand.eq.0)then           if(icand.eq.0)then
348              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE.EQ.1)then
349       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
350         $              ,ibest,iimage
351                endif
352              return              return
353           endif           endif
354    
# Line 236  c         iflag=0 Line 359  c         iflag=0
359  *     **********************************************************  *     **********************************************************
360  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
361  *     **********************************************************  *     **********************************************************
362             call guess()
363             do i=1,5
364                AL_GUESS(i)=AL(i)
365             enddo
366    
367           do i=1,5           do i=1,5
368              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
369           enddo           enddo
370            
371           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
372           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
373           jstep=0                !# minimization steps           jstep=0                !# minimization steps
374    
375           iprint=0           iprint=0
376           if(DEBUG)iprint=1  c         if(DEBUG.EQ.1)iprint=1
377             if(VERBOSE.EQ.1)iprint=1
378             if(DEBUG.EQ.1)iprint=2
379           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
380           if(ifail.ne.0) then  cc         if(ifail.ne.0) then
381              if(DEBUG)then           if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
382                if(CHI2.ne.CHI2)CHI2=-9999. !new
383                if(VERBOSE.EQ.1)then
384                 print *,                 print *,
385       $              '*** MINIMIZATION FAILURE *** (mini2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
386       $              ,iev       $              ,iev
387              endif              endif
             chi2=-chi2  
388           endif           endif
389                    
390           if(DEBUG)then           if(DEBUG.EQ.1)then
391              print*,'----------------------------- improved track coord'              print*,'----------------------------- improved track coord'
392  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
393              do ip=1,6              do ip=1,6
# Line 265  c         iflag=0 Line 398  c         iflag=0
398           endif           endif
399    
400  c         rchi2=chi2/dble(ndof)  c         rchi2=chi2/dble(ndof)
401           if(DEBUG)then           if(DEBUG.EQ.1)then
402              print*,' '              print*,' '
403              print*,'****** SELECTED TRACK *************'              print*,'****** SELECTED TRACK *************'
404              print*,'#         R. chi2        RIG'              print*,'#         R. chi2        RIG'
# Line 281  c         rchi2=chi2/dble(ndof) Line 414  c         rchi2=chi2/dble(ndof)
414  *     ------------- search if the track has an IMAGE -------------  *     ------------- search if the track has an IMAGE -------------
415  *     ------------- (also this is stored )           -------------  *     ------------- (also this is stored )           -------------
416           if(FIMAGE)goto 122     !>>> jump! (this is already an image)           if(FIMAGE)goto 122     !>>> jump! (this is already an image)
417  *     now search for track-image, by comparing couples IDs  
418    *     -----------------------------------------------------
419    *     first check if the track is ambiguous
420    *     -----------------------------------------------------
421    *     (modified on august 2007 by ElenaV)
422             is1=0
423             do ip=1,NPLANES
424                if(SENSOR_STORE(ip,icand).ne.0)then
425                   is1=SENSOR_STORE(ip,icand)
426                   if(ip.eq.6)is1=3-is1 !last plane inverted
427                endif
428             enddo
429             if(is1.eq.0)then
430                if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
431                goto 122            !jump
432             endif
433             do ip=1,NPLANES
434                is2 = SENSOR_STORE(ip,icand) !sensor
435                if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
436                if(
437         $           (is1.ne.is2.and.is2.ne.0)
438         $           )goto 122      !jump (this track cannot have an image)
439             enddo
440             if(DEBUG.eq.1)print*,' >>> ambiguous track! '
441    *     ---------------------------------------------------------------
442    *     take the candidate with the greatest number of matching couples
443    *     if more than one satisfies the requirement,
444    *     choose the one with more points and lower chi2
445    *     ---------------------------------------------------------------
446    *     count the number of matching couples
447             ncpmax = 0
448             iimage   = 0           !id of candidate with better matching
449           do i=1,ntracks           do i=1,ntracks
450              iimage=i              ncp=0
451              do ip=1,nplanes              do ip=1,nplanes
452                 if(     CP_STORE(nplanes-ip+1,icand).ne.                 if(CP_STORE(nplanes-ip+1,icand).ne.0)then
453       $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0                    if(
454         $                 CP_STORE(nplanes-ip+1,i).ne.0
455         $                 .and.
456         $                 CP_STORE(nplanes-ip+1,icand).eq.
457         $                 -1*CP_STORE(nplanes-ip+1,i)
458         $                 )then
459                         ncp=ncp+1  !ok
460                      elseif(
461         $                    CP_STORE(nplanes-ip+1,i).ne.0
462         $                    .and.
463         $                    CP_STORE(nplanes-ip+1,icand).ne.
464         $                    -1*CP_STORE(nplanes-ip+1,i)
465         $                    )then
466                         ncp = 0
467                         goto 117   !it is not an image candidate
468                      else
469                        
470                      endif
471                   endif
472              enddo              enddo
473              if(  iimage.ne.0.and.   117        continue
474  c     $           RCHI2_STORE(i).le.CHI2MAX.and.              trackimage(i)=ncp   !number of matching couples
475  c     $           RCHI2_STORE(i).gt.0.and.              if(ncp>ncpmax)then
476       $           .true.)then                 ncpmax=ncp
477                 if(DEBUG)print*,'Track candidate ',iimage                 iimage=i
478                endif
479             enddo
480    *     check if there are more than one image candidates
481             ngood=0
482             do i=1,ntracks
483                if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
484             enddo
485             if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
486    *     if there are, choose the best one
487             if(ngood.gt.0)then
488    *     -------------------------------------------------------
489    *     order track-candidates according to:
490    *     1st) decreasing n.points
491    *     2nd) increasing chi**2
492    *     -------------------------------------------------------
493                rchi2best=1000000000.
494                ndofbest=0            
495                do i=1,ntracks
496                   if( trackimage(i).eq.ncpmax )then
497                      ndof=0              
498                      do ii=1,nplanes    
499                         ndof=ndof        
500         $                    +int(xgood_store(ii,i))
501         $                    +int(ygood_store(ii,i))
502                      enddo              
503                      if(ndof.gt.ndofbest)then
504                         iimage=i
505                         rchi2best=RCHI2_STORE(i)
506                         ndofbest=ndof    
507                      elseif(ndof.eq.ndofbest)then
508                         if(RCHI2_STORE(i).lt.rchi2best.and.
509         $                    RCHI2_STORE(i).gt.0)then
510                            iimage=i
511                            rchi2best=RCHI2_STORE(i)
512                            ndofbest=ndof  
513                         endif            
514                      endif
515                   endif
516                enddo
517    
518                if(DEBUG.EQ.1)then
519                   print*,'Track candidate ',iimage
520       $              ,' >>> TRACK IMAGE >>> of'       $              ,' >>> TRACK IMAGE >>> of'
521       $              ,ibest       $              ,ibest
                goto 122         !image track found  
522              endif              endif
523           enddo              
524             endif
525    
526    
527   122     continue   122     continue
528    
529    
530  *     --- and store the results --------------------------------  *     --- and store the results --------------------------------
531           ntrk = ntrk + 1                   !counter of found tracks           ntrk = ntrk + 1                   !counter of found tracks
532           if(.not.FIMAGE           if(.not.FIMAGE
# Line 308  c     $           RCHI2_STORE(i).gt.0.an Line 535  c     $           RCHI2_STORE(i).gt.0.an
535       $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next       $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
536           if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous           if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
537           call fill_level2_tracks(ntrk)     !==> good2=.true.           call fill_level2_tracks(ntrk)     !==> good2=.true.
 c         print*,'++++++++++ iimage,fimage,ntrk,image '  
 c     $        ,iimage,fimage,ntrk,image(ntrk)  
538    
539           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
540              if(verbose)              if(verbose.eq.1)
541       $           print*,       $           print*,
542       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
543       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 325  cc            good2=.false. Line 550  cc            good2=.false.
550              goto 1212           !>>> fit image-track              goto 1212           !>>> fit image-track
551           endif           endif
552    
 *     --- then remove selected clusters (ibest+iimage) from clouds ----  
          call clean_XYclouds(ibest,iflag)  
          if(iflag.eq.1)then     !bad event  
             goto 880            !fill ntp and go to next event              
          endif  
   
 *     **********************************************************  
 *     condition to start a new search  
 *     **********************************************************  
          ixznew=0  
          do ixz=1,nclouds_xz  
             if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1  
          enddo  
          iyznew=0  
          do iyz=1,nclouds_yz  
             if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1  
          enddo  
           
          if(ixznew.ne.0.and.  
      $      iyznew.ne.0.and.  
      $        rchi2best.le.CHI2MAX.and.  
 c     $        rchi2best.lt.15..and.  
      $        .true.)then  
             if(DEBUG)then  
                print*,'***** NEW SEARCH ****'  
             endif  
             goto 11111          !try new search  
               
          endif  
 *     **********************************************  
   
   
553    
554   880     return   880     return
555        end        end
556    
557    
   
   
 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$$$  
   
558                
559  ************************************************************  ************************************************************
560  ************************************************************  ************************************************************
# Line 558  c$$$ Line 579  c$$$
579  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
580  *     angx   - Projected angle in x  *     angx   - Projected angle in x
581  *     angy   - Projected angle in y  *     angy   - Projected angle in y
582    *     bfx    - x-component of magnetci field
583    *     bfy    - y-component of magnetic field
584  *  *
585  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
586  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 596  c$$$ Line 619  c$$$
619  *  *
620  *  *
621    
622        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
623    
 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*****************************************************  
624                
625        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
626        include 'level1.f'        include 'level1.f'
627          include 'calib.f'
628        include 'common_align.f'        include 'common_align.f'
629        include 'common_mech.f'        include 'common_mech.f'
630        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
       include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
631    
632        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
633        integer sensor        integer sensor
634        integer viewx,viewy        integer viewx,viewy
635        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
636        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
637          real angx,angy            !X-Y effective angle
638          real bfx,bfy              !X-Y b-field components
639    
640        real stripx,stripy        real stripx,stripy
641    
642          double precision xi,yi,zi
643          double precision xi_A,yi_A,zi_A
644          double precision xi_B,yi_B,zi_B
645        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
646        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
647        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
 c      double precision xi,yi,zi  
 c      double precision xi_A,yi_A,zi_A  
 c      double precision xi_B,yi_B,zi_B  
648                
649    
650        parameter (ndivx=30)        parameter (ndivx=30)
651    
652    
653                
654        resxPAM = 0        resxPAM = 0
655        resyPAM = 0        resyPAM = 0
656    
657        xPAM = 0.        xPAM = 0.D0
658        yPAM = 0.        yPAM = 0.D0
659        zPAM = 0.        zPAM = 0.D0
660        xPAM_A = 0.        xPAM_A = 0.D0
661        yPAM_A = 0.        yPAM_A = 0.D0
662        zPAM_A = 0.        zPAM_A = 0.D0
663        xPAM_B = 0.        xPAM_B = 0.D0
664        yPAM_B = 0.        yPAM_B = 0.D0
665        zPAM_B = 0.        zPAM_B = 0.D0
666    
667          if(sensor.lt.1.or.sensor.gt.2)then
668             print*,'xyz_PAM   ***ERROR*** wrong input '
669             print*,'sensor ',sensor
670             icx=0
671             icy=0
672          endif
673    
674  *     -----------------  *     -----------------
675  *     CLUSTER X  *     CLUSTER X
676  *     -----------------  *     -----------------      
   
677        if(icx.ne.0)then        if(icx.ne.0)then
678           viewx = VIEW(icx)  
679           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
680           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
681           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
682            c         resxPAM = RESXAV
683           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
684           if(PFAx.eq.'COG1')then  !(1)  
685              stripx = stripx      !(1)           if(
686              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
687           elseif(PFAx.eq.'COG2')then       $        viewx.gt.12.or.
688              stripx = stripx + cog(2,icx)                   $        nldx.lt.1.or.
689              resxPAM = resxPAM*fbad_cog(2,icx)       $        nldx.gt.3.or.
690           elseif(PFAx.eq.'ETA2')then       $        stripx.lt.1.or.
691  c            cog2 = cog(2,icx)       $        stripx.gt.3072.or.
692  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                   $        .false.)then
693  c            stripx = stripx + etacorr              print*,'xyz_PAM   ***ERROR*** wrong input '
694              stripx = stripx + pfaeta2(icx,angx)            !(3)              print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
695              resxPAM = risx_eta2(angx)                       !   (4)              icx = 0
696              if(DEBUG.and.fbad_cog(2,icx).ne.1)              goto 10
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
             resxPAM = resxPAM*fbad_cog(2,icx)  
          elseif(PFAx.eq.'ETA3')then                         !(3)  
             stripx = stripx + pfaeta3(icx,angx)            !(3)  
             resxPAM = risx_eta3(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)  
             resxPAM = resxPAM*fbad_cog(3,icx)               !(3)  
          elseif(PFAx.eq.'ETA4')then                         !(3)  
             stripx = stripx + pfaeta4(icx,angx)            !(3)  
             resxPAM = risx_eta4(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)  
             resxPAM = resxPAM*fbad_cog(4,icx)               !(3)  
          elseif(PFAx.eq.'ETA')then                          !(3)  
             stripx = stripx + pfaeta(icx,angx)             !(3)  
             resxPAM = ris_eta(icx,angx)                     !   (4)  
             if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)  
 c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
             resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)  
          elseif(PFAx.eq.'COG')then           !(2)  
             stripx = stripx + cog(0,icx)     !(2)          
             resxPAM = risx_cog(angx)                        !   (4)  
             resxPAM = resxPAM*fbad_cog(0,icx)!(2)  
          else  
             print*,'*** Non valid p.f.a. (x) --> ',PFAx  
697           endif           endif
698    
699    *        --------------------------
700    *        magnetic-field corrections
701    *        --------------------------
702             stripx  = stripx +  fieldcorr(viewx,bfy)      
703             angx    = effectiveangle(ax,viewx,bfy)
704            
705             call applypfa(PFAx,icx,angx,corr,res)
706             stripx  = stripx + corr
707             resxPAM = res
708    
709     10   continue
710        endif        endif
711  c      if(icy.eq.0.and.icx.ne.0)      
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
         
712  *     -----------------  *     -----------------
713  *     CLUSTER Y  *     CLUSTER Y
714  *     -----------------  *     -----------------
715    
716        if(icy.ne.0)then        if(icy.ne.0)then
717    
718           viewy = VIEW(icy)           viewy = VIEW(icy)
719           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
720           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
721           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!  c         resyPAM = RESYAV
722             stripy = float(MAXS(icy))
723    
724             if(
725         $        viewy.lt.1.or.        
726         $        viewy.gt.12.or.
727         $        nldy.lt.1.or.
728         $        nldy.gt.3.or.
729         $        stripy.lt.1.or.
730         $        stripy.gt.3072.or.
731         $        .false.)then
732                print*,'xyz_PAM   ***ERROR*** wrong input '
733                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
734                icy = 0
735                goto 20
736             endif
737    
738           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
739              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG.EQ.1) then
740       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
741         $              ,icx,icy
742                endif
743              goto 100              goto 100
744           endif           endif
           
          stripy = float(MAXS(icy))  
          if(PFAy.eq.'COG1')then !(1)  
             stripy = stripy     !(1)  
             resyPAM = resyPAM   !(1)  
          elseif(PFAy.eq.'COG2')then  
             stripy = stripy + cog(2,icy)  
             resyPAM = resyPAM*fbad_cog(2,icy)  
          elseif(PFAy.eq.'ETA2')then  
 c            cog2 = cog(2,icy)  
 c            etacorr = pfaeta2(cog2,viewy,nldy,angy)  
 c            stripy = stripy + etacorr  
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
             resyPAM = resyPAM*fbad_cog(2,icy)  
             if(DEBUG.and.fbad_cog(2,icy).ne.1)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
          elseif(PFAy.eq.'ETA3')then                         !(3)  
             stripy = stripy + pfaeta3(icy,angy)            !(3)  
             resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
             if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)  
          elseif(PFAy.eq.'ETA4')then                         !(3)  
             stripy = stripy + pfaeta4(icy,angy)            !(3)  
             resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
             if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
          elseif(PFAy.eq.'ETA')then                          !(3)  
             stripy = stripy + pfaeta(icy,angy)             !(3)  
             resyPAM = ris_eta(icy,angy)                     !   (4)  
 c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  
             resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
             if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
          elseif(PFAy.eq.'COG')then  
             stripy = stripy + cog(0,icy)              
             resyPAM = risy_cog(angy)                        !   (4)  
 c            resyPAM = ris_eta(icy,angy)                    !   (4)  
             resyPAM = resyPAM*fbad_cog(0,icy)  
          else  
             print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
745    
746    *        --------------------------
747    *        magnetic-field corrections
748    *        --------------------------
749             stripy  = stripy +  fieldcorr(viewy,bfx)      
750             angy    = effectiveangle(ay,viewy,bfx)
751            
752             call applypfa(PFAy,icy,angy,corr,res)
753             stripy  = stripy + corr
754             resyPAM = res
755    
756     20   continue
757        endif        endif
758    
759          
760  c===========================================================  c===========================================================
761  C     COUPLE  C     COUPLE
762  C===========================================================  C===========================================================
# Line 778  c     (xi,yi,zi) = mechanical coordinate Line 767  c     (xi,yi,zi) = mechanical coordinate
767  c------------------------------------------------------------------------  c------------------------------------------------------------------------
768           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
769       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
770              print*,'xyz_PAM (couple):',              if(DEBUG.EQ.1) then
771       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
772         $              ' WARNING: false X strip: strip ',stripx
773                endif
774           endif           endif
775           xi = acoordsi(stripx,viewx)           xi = dcoordsi(stripx,viewx)
776           yi = acoordsi(stripy,viewy)           yi = dcoordsi(stripy,viewy)
777           zi = 0.           zi = 0.D0
778                    
   
779  c------------------------------------------------------------------------  c------------------------------------------------------------------------
780  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
781  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 819  c--------------------------------------- Line 809  c---------------------------------------
809           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
810           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
811    
812           xPAM_A = 0.           xPAM_A = 0.D0
813           yPAM_A = 0.           yPAM_A = 0.D0
814           zPAM_A = 0.           zPAM_A = 0.D0
815    
816           xPAM_B = 0.           xPAM_B = 0.D0
817           yPAM_B = 0.           yPAM_B = 0.D0
818           zPAM_B = 0.           zPAM_B = 0.D0
819    
820        elseif(        elseif(
821       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 845  C======================================= Line 835  C=======================================
835              nldx = nldy              nldx = nldy
836              viewx = viewy + 1              viewx = viewy + 1
837    
838              yi   = acoordsi(stripy,viewy)              xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
839                yi = dcoordsi(stripy,viewy)
840                zi = 0.D0
841    
842              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
843              yi_A = yi              yi_A = yi
# Line 855  C======================================= Line 847  C=======================================
847              yi_B = yi              yi_B = yi
848              zi_B = 0.              zi_B = 0.
849    
 c            print*,'Y-cl ',icy,stripy,' --> ',yi  
 c            print*,xi_A,' <--> ',xi_B  
850                            
851           elseif(icx.ne.0)then           elseif(icx.ne.0)then
852  c===========================================================  c===========================================================
# Line 867  C======================================= Line 857  C=======================================
857              nldy = nldx              nldy = nldx
858              viewy = viewx - 1              viewy = viewx - 1
859    
 c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx  
 c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  
860              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
861       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
862                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG.EQ.1) then
863       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
864         $                 ' WARNING: false X strip: strip ',stripx
865                   endif
866              endif              endif
867              xi   = acoordsi(stripx,viewx)  
868                xi = dcoordsi(stripx,viewx)
869                yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
870                zi = 0.D0
871    
872              xi_A = xi              xi_A = xi
873              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 890  c            if((stripx.le.3).or.(stripx Line 883  c            if((stripx.le.3).or.(stripx
883                 yi_B = yi                 yi_B = yi
884              endif              endif
885    
 c            print*,'X-cl ',icx,stripx,' --> ',xi  
 c            print*,yi_A,' <--> ',yi_B  
886    
887           else           else
888                if(DEBUG.EQ.1) then
889              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
890              print *,'icx = ',icx                 print *,'icx = ',icx
891              print *,'icy = ',icy                 print *,'icy = ',icy
892                endif
893              goto 100              goto 100
894                            
895           endif           endif
# Line 936  c     N.B. I convert angles from microra Line 928  c     N.B. I convert angles from microra
928       $        + zi_B       $        + zi_B
929       $        + dz(nplx,nldx,sensor)       $        + dz(nplx,nldx,sensor)
930    
931    
932    
933             xrt = xi
934         $        - omega(nplx,nldx,sensor)*yi
935         $        + gamma(nplx,nldx,sensor)*zi
936         $        + dx(nplx,nldx,sensor)
937            
938             yrt = omega(nplx,nldx,sensor)*xi
939         $        + yi
940         $        - beta(nplx,nldx,sensor)*zi
941         $        + dy(nplx,nldx,sensor)
942            
943             zrt = -gamma(nplx,nldx,sensor)*xi
944         $        + beta(nplx,nldx,sensor)*yi
945         $        + zi
946         $        + dz(nplx,nldx,sensor)
947    
948    
949                    
950  c      xrt = xi  c      xrt = xi
951  c      yrt = yi  c      yrt = yi
# Line 946  c     (xPAM,yPAM,zPAM) = measured coordi Line 956  c     (xPAM,yPAM,zPAM) = measured coordi
956  c                        in PAMELA reference system  c                        in PAMELA reference system
957  c------------------------------------------------------------------------  c------------------------------------------------------------------------
958    
959           xPAM = 0.           xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
960           yPAM = 0.           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
961           zPAM = 0.           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
962    c$$$         xPAM = 0.D0
963    c$$$         yPAM = 0.D0
964    c$$$         zPAM = 0.D0
965    
966           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
967           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 959  c--------------------------------------- Line 972  c---------------------------------------
972           zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4           zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
973                    
974    
 c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  
975    
976        else        else
977                       if(DEBUG.EQ.1) then
978           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
979           print *,'icx = ',icx              print *,'icx = ',icx
980           print *,'icy = ',icy              print *,'icy = ',icy
981                         endif
982        endif        endif
983                    
984    
985    
986   100  continue   100  continue
987        end        end
988    
989    ************************************************************************
990    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
991    *     (done to be called from c/c++)
992    ************************************************************************
993    
994          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
995    
996          include 'commontracker.f'
997          include 'level1.f'
998          include 'common_mini_2.f'
999          include 'common_xyzPAM.f'
1000          include 'common_mech.f'
1001          include 'calib.f'
1002          
1003    *     flag to chose PFA
1004    c$$$      character*10 PFA
1005    c$$$      common/FINALPFA/PFA
1006    
1007          integer icx,icy           !X-Y cluster ID
1008          integer sensor
1009          character*4 PFAx,PFAy     !PFA to be used
1010          real ax,ay                !X-Y geometric angle
1011          real bfx,bfy              !X-Y b-field components
1012    
1013          ipx=0
1014          ipy=0      
1015          
1016    c$$$      PFAx = 'COG4'!PFA
1017    c$$$      PFAy = 'COG4'!PFA
1018    
1019    
1020          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1021                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1022         $           ,' do not exists (n.clusters=',nclstr1,')'
1023                icx = -1*icx
1024                icy = -1*icy
1025                return
1026            
1027          endif
1028          
1029          call idtoc(pfaid,PFAx)
1030          call idtoc(pfaid,PFAy)
1031    
1032          
1033          if(icx.ne.0.and.icy.ne.0)then
1034    
1035             ipx=npl(VIEW(icx))
1036             ipy=npl(VIEW(icy))
1037            
1038             if( (nplanes-ipx+1).ne.ip )then
1039                print*,'xyzpam: ***WARNING*** cluster ',icx
1040         $           ,' does not belong to plane: ',ip
1041                icx = -1*icx
1042                return
1043             endif
1044             if( (nplanes-ipy+1).ne.ip )then
1045                print*,'xyzpam: ***WARNING*** cluster ',icy
1046         $           ,' does not belong to plane: ',ip
1047                icy = -1*icy
1048                return
1049             endif
1050    
1051             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1052    
1053             xgood(ip) = 1.
1054             ygood(ip) = 1.
1055             resx(ip)  = resxPAM
1056             resy(ip)  = resyPAM
1057    
1058             xm(ip) = xPAM
1059             ym(ip) = yPAM
1060             zm(ip) = zPAM
1061             xm_A(ip) = 0.D0
1062             ym_A(ip) = 0.D0
1063             xm_B(ip) = 0.D0
1064             ym_B(ip) = 0.D0
1065    
1066    c         zv(ip) = zPAM
1067    
1068          elseif(icx.eq.0.and.icy.ne.0)then
1069    
1070             ipy=npl(VIEW(icy))
1071             if( (nplanes-ipy+1).ne.ip )then
1072                print*,'xyzpam: ***WARNING*** cluster ',icy
1073         $           ,' does not belong to plane: ',ip
1074                icy = -1*icy
1075                return
1076             endif
1077    
1078             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1079            
1080             xgood(ip) = 0.
1081             ygood(ip) = 1.
1082             resx(ip)  = 1000.
1083             resy(ip)  = resyPAM
1084    
1085    c$$$         xm(ip) = -100.
1086    c$$$         ym(ip) = -100.
1087    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1088             xm(ip) = xPAM
1089             ym(ip) = yPAM
1090             zm(ip) = zPAM
1091             xm_A(ip) = xPAM_A
1092             ym_A(ip) = yPAM_A
1093             xm_B(ip) = xPAM_B
1094             ym_B(ip) = yPAM_B
1095    
1096    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1097            
1098          elseif(icx.ne.0.and.icy.eq.0)then
1099    
1100             ipx=npl(VIEW(icx))
1101    
1102             if( (nplanes-ipx+1).ne.ip )then
1103                print*,'xyzpam: ***WARNING*** cluster ',icx
1104         $           ,' does not belong to plane: ',ip
1105                icx = -1*icx
1106                return
1107             endif
1108    
1109             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1110          
1111             xgood(ip) = 1.
1112             ygood(ip) = 0.
1113             resx(ip)  = resxPAM
1114             resy(ip)  = 1000.
1115    
1116    c$$$         xm(ip) = -100.
1117    c$$$         ym(ip) = -100.
1118    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1119             xm(ip) = xPAM
1120             ym(ip) = yPAM
1121             zm(ip) = zPAM
1122             xm_A(ip) = xPAM_A
1123             ym_A(ip) = yPAM_A
1124             xm_B(ip) = xPAM_B
1125             ym_B(ip) = yPAM_B
1126            
1127    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1128    
1129          else
1130    
1131             il = 2
1132             if(lad.ne.0)il=lad
1133             is = 1
1134             if(sensor.ne.0)is=sensor
1135    
1136             xgood(ip) = 0.
1137             ygood(ip) = 0.
1138             resx(ip)  = 1000.
1139             resy(ip)  = 1000.
1140    
1141             xm(ip) = -100.
1142             ym(ip) = -100.          
1143             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1144             xm_A(ip) = 0.
1145             ym_A(ip) = 0.
1146             xm_B(ip) = 0.
1147             ym_B(ip) = 0.
1148    
1149    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1150    
1151          endif
1152    
1153          if(DEBUG.EQ.1)then
1154    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1155             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1156         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1157         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1158          endif
1159          end
1160    
1161  ********************************************************************************  ********************************************************************************
1162  ********************************************************************************  ********************************************************************************
# Line 993  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1178  c         print*,'A-(',xPAM_A,yPAM_A,')
1178  *      *    
1179  ********************************************************************************  ********************************************************************************
1180    
1181        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1182    
1183        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1184    
# Line 1002  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1187  c         print*,'A-(',xPAM_A,yPAM_A,')
1187  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1188  *     -----------------------------------  *     -----------------------------------
1189    
1190          real rXPP,rYPP
1191          double precision XPP,YPP
1192        double precision distance,RE        double precision distance,RE
1193        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1194    
1195          XPP=DBLE(rXPP)
1196          YPP=DBLE(rYPP)
1197    
1198  *     ----------------------  *     ----------------------
1199        if (        if (
1200       +     xPAM.eq.0.and.  c$$$     +     xPAM.eq.0.and.
1201       +     yPAM.eq.0.and.  c$$$     +     yPAM.eq.0.and.
1202       +     zPAM.eq.0.and.  c$$$     +     zPAM.eq.0.and.
1203       +     xPAM_A.ne.0.and.       +     xPAM_A.ne.0.and.
1204       +     yPAM_A.ne.0.and.       +     yPAM_A.ne.0.and.
1205       +     zPAM_A.ne.0.and.       +     zPAM_A.ne.0.and.
# Line 1047  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1237  c         print*,'A-(',xPAM_A,yPAM_A,')
1237           endif                   endif        
1238    
1239           distance=           distance=
1240       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1241    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1242           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1243    
 c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',re  
1244    
1245                    
1246  *     ----------------------  *     ----------------------
1247        elseif(        elseif(
1248       +     xPAM.ne.0.and.  c$$$     +     xPAM.ne.0.and.
1249       +     yPAM.ne.0.and.  c$$$     +     yPAM.ne.0.and.
1250       +     zPAM.ne.0.and.  c$$$     +     zPAM.ne.0.and.
1251       +     xPAM_A.eq.0.and.       +     xPAM_A.eq.0.and.
1252       +     yPAM_A.eq.0.and.       +     yPAM_A.eq.0.and.
1253       +     zPAM_A.eq.0.and.       +     zPAM_A.eq.0.and.
# Line 1072  c$$$         print*,' resolution ',re Line 1260  c$$$         print*,' resolution ',re
1260  *     ----------------------  *     ----------------------
1261                    
1262           distance=           distance=
1263       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1264       $        +       $       +
1265       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1266    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1267    c$$$     $        +
1268    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1269           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1270    
 c$$$         print*,xPAM,yPAM,zPAM  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',resxPAM,resyPAM  
1271                    
1272        else        else
1273                    
          print*  
      $        ,' function distance_to ---> wrong usage!!!'  
          print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  
          print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '  
      $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
1274        endif          endif  
1275    
1276        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1127  c$$$         print*,' resolution ',resxP Line 1310  c$$$         print*,' resolution ',resxP
1310        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1311        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1312        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1313        real*8 yvvv,xvvv        double precision yvvv,xvvv
1314        double precision xi,yi,zi        double precision xi,yi,zi
1315        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1316        real AA,BB        real AA,BB
1317        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1318    
1319  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1320        real ptoll        real ptoll
1321        data ptoll/150./          !um        data ptoll/150./          !um
1322    
1323        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1324    
1325        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1326        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1152  c$$$         print*,' resolution ',resxP Line 1335  c$$$         print*,' resolution ',resxP
1335  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1336  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1337  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1338                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 xi = dcoordsi(stripx,viewx)
1339       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...                 yi = dcoordsi(stripy,viewy)
1340  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...                 zi = 0.D0
                   print*,'whichsensor: ',  
      $                ' WARNING: false X strip: strip ',stripx  
                endif  
                xi = acoordsi(stripx,viewx)  
                yi = acoordsi(stripy,viewy)  
                zi = 0.  
1341  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1342  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1343  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1185  c--------------------------------------- Line 1362  c---------------------------------------
1362    
1363                 yvv(iv)=sngl(yvvv)                 yvv(iv)=sngl(yvvv)
1364                 xvv(iv)=sngl(xvvv)                 xvv(iv)=sngl(xvvv)
 c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '  
 c     $              ,iv,xvv(iv),yvv(iv)  
1365              enddo               !end loop on sensor vertexes              enddo               !end loop on sensor vertexes
1366    
1367              dtot=0.              dtot=0.
# Line 1285  c     $              ,iv,xvv(iv),yvv(iv) Line 1460  c     $              ,iv,xvv(iv),yvv(iv)
1460  *     it returns the plane number  *     it returns the plane number
1461  *      *    
1462        include 'commontracker.f'        include 'commontracker.f'
1463          include 'level1.f'
1464  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1465        include 'common_momanhough.f'        include 'common_momanhough.f'
1466                
# Line 1310  c      include 'common_analysis.f' Line 1486  c      include 'common_analysis.f'
1486        is_cp=0        is_cp=0
1487        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1488        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
       if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  
1489    
1490        return        return
1491        end        end
# Line 1322  c      include 'common_analysis.f' Line 1497  c      include 'common_analysis.f'
1497  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1498  *      *    
1499        include 'commontracker.f'        include 'commontracker.f'
1500          include 'level1.f'
1501  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1502        include 'common_momanhough.f'        include 'common_momanhough.f'
1503                
# Line 1350  c      include 'common_analysis.f' Line 1526  c      include 'common_analysis.f'
1526  *     positive if sensor =2  *     positive if sensor =2
1527  *  *
1528        include 'commontracker.f'        include 'commontracker.f'
1529          include 'level1.f'
1530  c      include 'calib.f'  c      include 'calib.f'
1531  c      include 'level1.f'  c      include 'level1.f'
1532  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1379  c      include 'common_analysis.f' Line 1556  c      include 'common_analysis.f'
1556  *************************************************************************  *************************************************************************
1557  *************************************************************************  *************************************************************************
1558  *************************************************************************  *************************************************************************
 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  
1559                
1560    
1561  ***************************************************  ***************************************************
# Line 1661  c$$$      end Line 1570  c$$$      end
1570        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1571    
1572        include 'commontracker.f'        include 'commontracker.f'
1573          include 'level1.f'
1574        include 'common_momanhough.f'        include 'common_momanhough.f'
1575        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1576        include 'calib.f'        include 'calib.f'
1577        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1578    
1579  *     output flag  *     output flag
1580  *     --------------  *     --------------
# Line 1678  c      common/dbg/DEBUG Line 1585  c      common/dbg/DEBUG
1585    
1586        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1587    
1588          iflag = iflag
1589          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1590    
1591    cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
1592    
1593  *     init variables  *     init variables
       ncp_tot=0  
1594        do ip=1,nplanes        do ip=1,nplanes
1595           do ico=1,ncouplemax           do ico=1,ncouplemax
1596              clx(ip,ico)=0              clx(ip,ico)=0
# Line 1692  c      common/dbg/DEBUG Line 1603  c      common/dbg/DEBUG
1603           ncls(ip)=0           ncls(ip)=0
1604        enddo        enddo
1605        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1606           cl_single(icl)=1           cl_single(icl) = 1     !all are single
1607           cl_good(icl)=0           cl_good(icl)   = 0     !all are bad
1608          enddo
1609          do iv=1,nviews
1610             ncl_view(iv)  = 0
1611             mask_view(iv) = 0      !all included
1612        enddo        enddo
1613                
1614    *     count number of cluster per view
1615          do icl=1,nclstr1
1616             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1617          enddo
1618    *     mask views with too many clusters
1619          do iv=1,nviews
1620             if( ncl_view(iv).gt. nclusterlimit)then
1621    c            mask_view(iv) = 1
1622                mask_view(iv) = mask_view(iv) + 2**0
1623                if(DEBUG.EQ.1)
1624         $        print*,' * WARNING * cl_to_couple: n.clusters > '
1625         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1626             endif
1627          enddo
1628    
1629    
1630  *     start association  *     start association
1631        ncouples=0        ncouples=0
1632        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1633           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1634                    
1635             if(cl_used(icx).ne.0)goto 10
1636    
1637  *     ----------------------------------------------------  *     ----------------------------------------------------
1638  *     cut on charge (X VIEW)  *     jump masked views (X VIEW)
1639  *     ----------------------------------------------------  *     ----------------------------------------------------
1640           if(dedx(icx).lt.dedx_x_min)then           if( mask_view(VIEW(icx)).ne.0 ) goto 10
             cl_single(icx)=0  
             goto 10  
          endif  
1641  *     ----------------------------------------------------  *     ----------------------------------------------------
1642  *     cut on multiplicity (X VIEW)  *     cut on charge (X VIEW)
1643  *     ----------------------------------------------------  *     ----------------------------------------------------
1644           if(mult(icx).ge.mult_x_max)then           if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
1645              cl_single(icx)=0              cl_single(icx)=0
1646              goto 10              goto 10
1647           endif           endif
# Line 1752  c     endif Line 1682  c     endif
1682                    
1683           do icy=1,nclstr1       !loop on cluster (Y)           do icy=1,nclstr1       !loop on cluster (Y)
1684              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1685    
1686                if(cl_used(icx).ne.0)goto 20
1687                            
1688  *     ----------------------------------------------------  *     ----------------------------------------------------
1689  *     cut on charge (Y VIEW)  *     jump masked views (Y VIEW)
1690  *     ----------------------------------------------------  *     ----------------------------------------------------
1691              if(dedx(icy).lt.dedx_y_min)then              if( mask_view(VIEW(icy)).ne.0 ) goto 20
1692                 cl_single(icy)=0  
                goto 20  
             endif  
1693  *     ----------------------------------------------------  *     ----------------------------------------------------
1694  *     cut on multiplicity (X VIEW)  *     cut on charge (Y VIEW)
1695  *     ----------------------------------------------------  *     ----------------------------------------------------
1696              if(mult(icy).ge.mult_y_max)then              if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
1697                 cl_single(icy)=0                 cl_single(icy)=0
1698                 goto 20                 goto 20
1699              endif              endif
# Line 1809  c     endif Line 1739  c     endif
1739  *     charge correlation  *     charge correlation
1740  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1741    
1742  *     -------------------------------------------------------------                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
 *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<  
 *     -------------------------------------------------------------  
                if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)  
1743       $              .and.       $              .and.
1744       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1745       $              .and.       $              .and.
1746       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1747       $              .and.       $              .and.
1748       $              .true.)then       $              .true.)then
1749    
1750                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1751       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1752                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1753    
1754  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1755    
1756                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1757       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1758                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1759                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1835  c                  cut = chcut * sch(npl Line 1762  c                  cut = chcut * sch(npl
1762                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1763                    endif                    endif
1764                 endif                 endif
1765                  
1766  *     ------------------> COUPLE <------------------                 if(ncp_plane(nplx).gt.ncouplemax)then
1767  *     check to do not overflow vector dimentions                    if(verbose.eq.1)print*,
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
1768       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1769       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1770       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1771       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1772  c     good2=.false.  c                  mask_view(nviewx(nplx)) = 2
1773  c     goto 880   !fill ntp and go to next event                      c                  mask_view(nviewy(nply)) = 2
1774                    iflag=1                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1775                    return                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1776                      goto 10
1777                 endif                 endif
1778                                
1779    *     ------------------> COUPLE <------------------
1780                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1781                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1782                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1783                 cl_single(icx)=0                 cl_single(icx)=0
1784                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1785  *     ----------------------------------------------  *     ----------------------------------------------
1786    
1787                endif                              
1788    
1789   20         continue   20         continue
1790           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1791                    
1792   10      continue   10      continue
1793        enddo                     !end loop on clusters(X)        enddo                     !end loop on clusters(X)
1794                
         
1795        do icl=1,nclstr1        do icl=1,nclstr1
1796           if(cl_single(icl).eq.1)then           if(cl_single(icl).eq.1)then
1797              ip=npl(VIEW(icl))              ip=npl(VIEW(icl))
# Line 1884  c     goto 880   !fill ntp and go to nex Line 1799  c     goto 880   !fill ntp and go to nex
1799              cls(ip,ncls(ip))=icl              cls(ip,ncls(ip))=icl
1800           endif           endif
1801        enddo        enddo
1802    
1803    c 80   continue
1804          continue
1805                
1806                
1807        if(DEBUG)then        if(DEBUG.EQ.1)then
1808           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1809           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1810           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'used    ',(cl_used(i),i=1,nclstr1)
1811             print*,'singlets',(cl_single(i),i=1,nclstr1)
1812           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1813        endif        endif
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
1814    
1815  c      logical DEBUG    
1816  c      common/dbg/DEBUG        if(.not.RECOVER_SINGLETS)goto 81
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
1817    
1818        integer badseed,badcl  *     ////////////////////////////////////////////////
1819    *     PATCH to recover events with less than 3 couples
1820    *     ////////////////////////////////////////////////    
1821    *     loop over singlet and create "fake" couples
1822    *     (with clx=0 or cly=0)
1823    *    
1824    
1825  *     init variables        if(DEBUG.EQ.1)
1826        ncp_tot=0       $     print*,'>>> Recover singlets '
1827        do ip=1,nplanes       $     ,'(creates fake couples) <<<'
1828           do ico=1,ncouplemax        do icl=1,nclstr1
1829              clx(ip,ico)=0           if(
1830              cly(ip,ico)=0       $        cl_single(icl).eq.1.and.
1831           enddo       $        cl_used(icl).eq.0.and.
1832           ncp_plane(ip)=0       $        .true.)then
          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  
           
1833  *     ----------------------------------------------------  *     ----------------------------------------------------
1834  *     cut on charge (X VIEW)  *     jump masked views
          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  
1835  *     ----------------------------------------------------  *     ----------------------------------------------------
1836                        if( mask_view(VIEW(icl)).ne.0 ) goto 21
1837           cl_good(icx)=1              if(mod(VIEW(icl),2).eq.0)then !=== X views
          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  
               
1838  *     ----------------------------------------------------  *     ----------------------------------------------------
1839  *     cut on charge (Y VIEW)  *     cut on charge (X 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  
1840  *     ----------------------------------------------------  *     ----------------------------------------------------
1841                               if(sgnl(icl).lt.dedx_x_min) goto 21
               
             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  
 *     ===========================================================  
1842                                
1843                   nplx=npl(VIEW(icl))
1844    *     ------------------> (FAKE) COUPLE <-----------------
1845                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1846                   clx(nplx,ncp_plane(nplx))=icl
1847                   cly(nplx,ncp_plane(nplx))=0
1848    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1849    *     ----------------------------------------------------
1850                                
1851  *     ------------------> COUPLE <------------------              else                !=== Y views
1852  *     check to do not overflow vector dimentions  *     ----------------------------------------------------
1853  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  *     cut on charge (Y VIEW)
1854  c$$$                  if(DEBUG)print*,  *     ----------------------------------------------------
1855  c$$$     $                    ' ** warning ** number of identified'//                 if(sgnl(icl).lt.dedx_y_min) goto 21
 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  
1856                                
1857                 if(ncp_plane(nplx).eq.ncouplemax)then                 nply=npl(VIEW(icl))
1858                    if(verbose)print*,  *     ------------------> (FAKE) COUPLE <-----------------
1859       $                 '** warning ** number of identified '//                 ncp_plane(nply) = ncp_plane(nply) + 1
1860       $                 'couples on plane ',nplx,                 clx(nply,ncp_plane(nply))=0
1861       $                 'exceeds vector dimention '                 cly(nply,ncp_plane(nply))=icl
1862       $                 ,'( ',ncouplemax,' )'  c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1863  c     good2=.false.  *     ----------------------------------------------------
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
1864                                
1865                 ncp_plane(nplx) = ncp_plane(nplx) + 1              endif
1866                 clx(nplx,ncp_plane(nplx))=icx           endif                  !end "single" condition
1867                 cly(nply,ncp_plane(nplx))=icy   21      continue
1868                 cl_single(icx)=0        enddo                     !end loop over clusters
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
1869    
1870   20         continue        if(DEBUG.EQ.1)
1871           enddo                  !end loop on clusters(Y)       $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1872            
1873   10      continue  
1874        enddo                     !end loop on clusters(X)   81   continue
         
         
       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  
1875                
1876        do ip=1,6        ncp_tot=0
1877           ncp_tot=ncp_tot+ncp_plane(ip)        do ip=1,NPLANES
1878             ncp_tot = ncp_tot + ncp_plane(ip)
1879        enddo        enddo
1880  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)        if(DEBUG.EQ.1)
1881               $     print*,'n.couple tot:      ',ncp_tot
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1882                
1883        return        return
1884        end        end
   
1885                
1886  ***************************************************  ***************************************************
1887  *                                                 *  *                                                 *
# Line 2143  c     goto 880       !fill ntp and go to Line 1893  c     goto 880       !fill ntp and go to
1893  **************************************************  **************************************************
1894    
1895        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1896    
1897        include 'commontracker.f'        include 'commontracker.f'
1898          include 'level1.f'
1899        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1900        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1901        include 'common_mini_2.f'        include 'common_mini_2.f'
1902        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1903    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1904    
1905  *     output flag  *     output flag
1906  *     --------------  *     --------------
# Line 2188  c      double precision xm3,ym3,zm3 Line 1932  c      double precision xm3,ym3,zm3
1932        real xc,zc,radius        real xc,zc,radius
1933  *     -----------------------------  *     -----------------------------
1934    
1935          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1936    
1937    *     --------------------------------------------
1938    *     put a limit to the maximum number of couples
1939    *     per plane, in order to apply hough transform
1940    *     (couples recovered during track refinement)
1941    *     --------------------------------------------
1942          do ip=1,nplanes
1943             if(ncp_plane(ip).gt.ncouplelimit)then
1944                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1945                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1946             endif
1947          enddo
1948    
1949    
1950        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1951        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1952                
1953        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1954           do is1=1,2             !loop on sensors - COPPIA 1  c$$$         print*,'(1) ip ',ip1
1955                c$$$     $        ,mask_view(nviewx(ip1))
1956    c$$$     $        ,mask_view(nviewy(ip1))        
1957             if(  mask_view(nviewx(ip1)).ne.0 .or.
1958         $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1959             do is1=1,2             !loop on sensors - COPPIA 1            
1960              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1961                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1962                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1963                  
1964    c$$$               print*,'(1) ip ',ip1,' icp ',icp1
1965    
1966  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1967                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1968                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1969                 xm1=xPAM                 xm1=xPAM
1970                 ym1=yPAM                 ym1=yPAM
1971                 zm1=zPAM                                   zm1=zPAM                  
1972  c     print*,'***',is1,xm1,ym1,zm1  
1973                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1974                    do is2=1,2    !loop on sensors -ndblt COPPIA 2  c$$$                  print*,'(2) ip ',ip2
1975                        c$$$     $                 ,mask_view(nviewx(ip2))
1976    c$$$     $                 ,mask_view(nviewy(ip2))
1977                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1978         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1979                      do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
1980                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1981                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1982                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1983    
1984    c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
1985    
1986  c                        call xyz_PAM  c                        call xyz_PAM
1987  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1988    c                        call xyz_PAM
1989    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1990                          call xyz_PAM                          call xyz_PAM
1991       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1992                          xm2=xPAM                          xm2=xPAM
1993                          ym2=yPAM                          ym2=yPAM
1994                          zm2=zPAM                          zm2=zPAM
1995                                                    
1996    *                       ---------------------------------------------------
1997    *                       both couples must have a y-cluster
1998    *                       (condition necessary when in RECOVER_SINGLETS mode)
1999    *                       ---------------------------------------------------
2000                            if(icy1.eq.0.or.icy2.eq.0)goto 111
2001    
2002                            if(cl_used(icy1).ne.0)goto 111
2003                            if(cl_used(icy2).ne.0)goto 111
2004    
2005                            
2006  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2007  *     track parameters on Y VIEW  *     track parameters on Y VIEW
2008  *     (2 couples needed)  *     (2 couples needed)
2009  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2010                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2011                             if(verbose)print*,                             if(verbose.eq.1)print*,
2012       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2013       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2014       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2015  c                           good2=.false.  c     good2=.false.
2016  c                           goto 880 !fill ntp and go to next event  c     goto 880 !fill ntp and go to next event
2017                               do iv=1,12
2018    c     mask_view(iv) = 3
2019                                  mask_view(iv) = mask_view(iv)+ 2**2
2020                               enddo
2021                             iflag=1                             iflag=1
2022                             return                             return
2023                          endif                          endif
2024                            
2025                            
2026    ccc                        print*,'<doublet> ',icp1,icp2
2027    
2028                          ndblt = ndblt + 1                          ndblt = ndblt + 1
2029  *     store doublet info  *     store doublet info
2030                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
# Line 2241  c                           goto 880 !fi Line 2033  c                           goto 880 !fi
2033                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2034  *     y0 (cm)  *     y0 (cm)
2035                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
2036                                                      
2037  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2038  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
2039  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2040                            if(SECOND_SEARCH)goto 111
2041                          if(                          if(
2042       $                       abs(alfayz2(ndblt)).gt.alfyz2_max       $                       abs(alfayz2(ndblt)).gt.alfyz2_max
2043       $                       .or.       $                       .or.
2044       $                       abs(alfayz1(ndblt)).gt.alfyz1_max       $                       abs(alfayz1(ndblt)).gt.alfyz1_max
2045       $                       )ndblt = ndblt-1       $                       )ndblt = ndblt-1
2046                                                    
2047  c$$$      if(iev.eq.33)then                          
2048  c$$$      print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2   111                    continue
 c$$$     $        ,' || ',icx1,icy1,icx2,icy2  
 c$$$     $        ,' || ',xm1,ym1,xm2,ym2  
 c$$$     $        ,' || ',alfayz2(ndblt),alfayz1(ndblt)  
 c$$$      endif  
 c$$$  
2049  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2050  *     track parameters on Y VIEW - end  *     track parameters on Y VIEW - end
2051  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2052    
2053    
2054                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(icx1.ne.0)then
2055                               if(cl_used(icx1).ne.0)goto 31
2056                            endif
2057                            if(icx2.ne.0)then
2058                               if(cl_used(icx2).ne.0)goto 31
2059                            endif
2060    
2061                            if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2062    
2063                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2064    c$$$                           print*,'(3) ip ',ip3
2065    c$$$     $                          ,mask_view(nviewx(ip3))
2066    c$$$     $                          ,mask_view(nviewy(ip3))                          
2067                               if(  mask_view(nviewx(ip3)).ne.0 .or.
2068         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2069                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
2070                                                                
2071                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
2072                                   icx3=clx(ip3,icp3)                                   icx3=clx(ip3,icp3)
2073                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2074    
2075    c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
2076    
2077    *     ---------------------------------------------------
2078    *     all three couples must have a x-cluster
2079    *     (condition necessary when in RECOVER_SINGLETS mode)
2080    *     ---------------------------------------------------
2081                                     if(
2082         $                                icx1.eq.0.or.
2083         $                                icx2.eq.0.or.
2084         $                                icx3.eq.0.or.
2085         $                                .false.)goto 29
2086                                    
2087                                     if(cl_used(icx1).ne.0)goto 29
2088                                     if(cl_used(icx2).ne.0)goto 29
2089                                     if(cl_used(icx3).ne.0)goto 29
2090    
2091  c                                 call xyz_PAM  c                                 call xyz_PAM
2092  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2093    c                                 call xyz_PAM
2094    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2095                                   call xyz_PAM                                   call xyz_PAM
2096       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2097         $                                ,0.,0.,0.,0.)
2098                                   xm3=xPAM                                   xm3=xPAM
2099                                   ym3=yPAM                                   ym3=yPAM
2100                                   zm3=zPAM                                   zm3=zPAM
2101    
2102    
2103  *     find the circle passing through the three points  *     find the circle passing through the three points
2104                                     iflag_t = DEBUG
2105                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2106       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag_t)
 c     print*,xc,zc,radius  
2107  *     the circle must intersect the reference plane  *     the circle must intersect the reference plane
2108                                   if(  cc                                 if(iflag.ne.0)goto 29
2109  c     $                                 (xc.le.-1.*xclimit.or.                                   if(iflag_t.ne.0)then
2110  c     $                                 xc.ge.xclimit).and.  *     if tricircle fails, evaluate a straight line
2111       $                                radius**2.ge.(ZINI-zc)**2.and.                                      if(DEBUG.eq.1)
2112       $                                iflag.eq.0.and.       $                                   print*,'TRICIRCLE failure'
2113       $                                .true.)then       $                                   ,' >>> straight line'
2114                                                                        radius=0.
2115                                        xc=0.
2116                                        yc=0.
2117                                        
2118                                        SZZ=0.                  
2119                                        SZX=0.
2120                                        SSX=0.
2121                                        SZ=0.
2122                                        S1=0.
2123                                        X0=0.
2124                                        Ax=0.
2125                                        BX=0.
2126                                        DO I=1,3
2127                                           XX = XP(I)
2128                                           SZZ=SZZ+ZP(I)*ZP(I)
2129                                           SZX=SZX+ZP(I)*XX
2130                                           SSX=SSX+XX
2131                                           SZ=SZ+ZP(I)
2132                                           S1=S1+1.
2133                                        ENDDO
2134                                        DET=SZZ*S1-SZ*SZ
2135                                        AX=(SZX*S1-SZ*SSX)/DET
2136                                        BX=(SZZ*SSX-SZX*SZ)/DET
2137                                        X0  = AX*ZINI+BX
2138                                        
2139                                     endif
2140    
2141                                     if(  .not.SECOND_SEARCH.and.
2142         $                                radius**2.lt.(ZINI-zc)**2)goto 29
2143                                                                      
2144  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2145  *     track parameters on X VIEW  *     track parameters on X VIEW
2146  *     (3 couples needed)  *     (3 couples needed)
2147  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2148                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2149                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2150       $                     '** warning ** number of identified '//       $                                   '** warning **'//
2151       $                     'triplets exceeds vector dimention '       $                                   ' number of identified '//
2152       $                    ,'( ',ntrpt_max,' )'       $                                   'triplets exceeds'//
2153  c                                    good2=.false.       $                                   ' vector dimention '
2154  c                                    goto 880 !fill ntp and go to next event       $                                   ,'( ',ntrpt_max,' )'
2155    c     good2=.false.
2156    c     goto 880 !fill ntp and go to next event
2157                                        do iv=1,nviews
2158    c     mask_view(iv) = 4
2159                                           mask_view(iv) =
2160         $                                      mask_view(iv)+ 2**3
2161                                        enddo
2162                                      iflag=1                                      iflag=1
2163                                      return                                      return
2164                                   endif                                   endif
2165    
2166    ccc                                 print*,'<triplet> ',icp1,icp2,icp3
2167                                    
2168                                   ntrpt = ntrpt +1                                   ntrpt = ntrpt +1
2169  *     store triplet info  *     store triplet info
2170                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
2171                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
2172                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
2173                                                                    
2174                                   if(xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2175  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2176                alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
2177                alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2178                alfaxz3(ntrpt) = 1/radius               alfaxz3(ntrpt) = 1/radius
2179                                   else                                  else if(radius.ne.0.and.xc.ge.0)then
2180  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2181                alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)               alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
2182                alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2183                alfaxz3(ntrpt) = -1/radius               alfaxz3(ntrpt) = -1/radius
2184                                   endif                                  else if(radius.eq.0)then
2185                                    *************straight fit
2186                 alfaxz1(ntrpt) = X0
2187                 alfaxz2(ntrpt) = AX
2188                 alfaxz3(ntrpt) = 0.
2189                                    endif
2190    
2191    c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
2192    c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
2193    c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
2194                                        
2195  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2196  ****  reject non phisical triplets                   ****  ****  reject non phisical triplets                   ****
2197  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2198                                   if(                                  if(SECOND_SEARCH)goto 29
2199       $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max                                  if(
2200       $                                .or.       $                               abs(alfaxz2(ntrpt)).gt.
2201       $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max       $                               alfxz2_max
2202       $                                )ntrpt = ntrpt-1       $                               .or.
2203                                         $                               abs(alfaxz1(ntrpt)).gt.
2204                                         $                               alfxz1_max
2205  c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)       $                               )ntrpt = ntrpt-1
2206                                                                    
2207  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2208  *     track parameters on X VIEW - end  *     track parameters on X VIEW - end
2209  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2210                                endif                                
2211     29                           continue
2212                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2213                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2214     30                     continue
2215                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2216   30                  continue  
2217                         31                  continue                    
2218   1                enddo         !end loop on COPPIA 2  c 1                enddo         !end loop on COPPIA 2
2219                     enddo         !end loop on COPPIA 2
2220                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2221     20            continue
2222              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2223                
2224    c 11         continue
2225              continue
2226           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2227        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2228     10   continue
2229        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2230                
2231        if(DEBUG)then        if(DEBUG.EQ.1)then
2232           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2233           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2234        endif        endif
# Line 2374  c     goto 880               !ntp fill Line 2253  c     goto 880               !ntp fill
2253        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2254    
2255        include 'commontracker.f'        include 'commontracker.f'
2256          include 'level1.f'
2257        include 'common_momanhough.f'        include 'common_momanhough.f'
2258        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2259    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2260    
2261  *     output flag  *     output flag
2262  *     --------------  *     --------------
# Line 2397  c      common/dbg/DEBUG Line 2275  c      common/dbg/DEBUG
2275        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2276        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2277    
2278          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2279    
2280  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2281  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2410  c      common/dbg/DEBUG Line 2289  c      common/dbg/DEBUG
2289        distance=0        distance=0
2290        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2291        npt_tot=0        npt_tot=0
2292          nloop=0                  
2293     90   continue                  
2294        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2295           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
2296                            
 c     print*,'--------------'  
 c     print*,'** ',idb1,' **'  
               
2297           do icp=1,ncp_tot           do icp=1,ncp_tot
2298              cp_useds1(icp)=0    !init              cp_useds1(icp)=0    !init
2299              cp_useds2(icp)=0    !init              cp_useds2(icp)=0    !init
# Line 2451  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2329  ccccc if(db_used(idbref).eq.1)goto 1188
2329  *     doublet distance in parameter space  *     doublet distance in parameter space
2330                 distance=                 distance=
2331       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
2332       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2                     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
2333                 distance = sqrt(distance)                 distance = sqrt(distance)
2334                                
 c$$$      if(iev.eq.33)then  
 c$$$      if(distance.lt.100)  
 c$$$     $ print*,'********* ',idb1,idbref,idb2,distance  
 c$$$      if(distance.lt.100)  
 c$$$     $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)  
 c$$$     $                    ,alfayz2(idbref),alfayz2(idb2)  
 c$$$      endif  
2335                 if(distance.lt.cutdistyz)then                 if(distance.lt.cutdistyz)then
2336    
 c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  
2337                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
2338                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
2339                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
# Line 2479  c     print*,idb1,idb2,distance,' cloud Line 2349  c     print*,idb1,idb2,distance,' cloud
2349    
2350                    temp1 = temp1 + alfayz1(idb2)                    temp1 = temp1 + alfayz1(idb2)
2351                    temp2 = temp2 + alfayz2(idb2)                    temp2 = temp2 + alfayz2(idb2)
 c     print*,'*   idbref,idb2 ',idbref,idb2  
2352                 endif                               endif              
2353                                
2354   1118          continue   1118          continue
2355              enddo               !end loop (2) on DOUBLETS              enddo               !end loop (2) on DOUBLETS
2356                            
2357   1188       continue  c 1188       continue
2358                continue
2359           enddo                  !end loop on... bo?           enddo                  !end loop on... bo?
2360                    
2361           nptloop=npv           nptloop=npv
# Line 2502  c     print*,'*   idbref,idb2 ',idbref,i Line 2372  c     print*,'*   idbref,idb2 ',idbref,i
2372           enddo           enddo
2373           ncpused=0           ncpused=0
2374           do icp=1,ncp_tot           do icp=1,ncp_tot
2375              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2376         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2377         $           .true.)then
2378                 ncpused=ncpused+1                 ncpused=ncpused+1
2379                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2380                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2512  c     print*,'*   idbref,idb2 ',idbref,i Line 2384  c     print*,'*   idbref,idb2 ',idbref,i
2384           do ip=1,nplanes           do ip=1,nplanes
2385              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2386           enddo           enddo
2387  c     print*,'>>>> ',ncpused,npt,nplused          
          if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  
          if(npt.lt.nptyz_min)goto 2228 !next doublet  
2388           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2389                    
2390  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2391  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2392    
2393           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2394              if(verbose)print*,              if(verbose.eq.1)print*,
2395       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2396       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2397       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2398  c               good2=.false.  c               good2=.false.
2399  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2400                do iv=1,nviews
2401    c               mask_view(iv) = 5
2402                   mask_view(iv) = mask_view(iv) + 2**4
2403                enddo
2404              iflag=1              iflag=1
2405              return              return
2406           endif           endif
# Line 2542  c     goto 880         !fill ntp and go Line 2416  c     goto 880         !fill ntp and go
2416  c     ptcloud_yz_nt(nclouds_yz)=npt  c     ptcloud_yz_nt(nclouds_yz)=npt
2417           do ipt=1,npt           do ipt=1,npt
2418              db_cloud(npt_tot+ipt) = db_all(ipt)              db_cloud(npt_tot+ipt) = db_all(ipt)
 c     print*,'>> ',ipt,db_all(ipt)  
2419           enddo             enddo  
2420           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2421           if(DEBUG)then           if(DEBUG.EQ.1)then
             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'  
2422              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2423              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2424              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2425              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2426              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2427              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2428  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2429  c$$$     $           ,ptcloud_yz(nclouds_yz)              print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
 c$$$            print*,'nt-uple: db_cloud(...) = '  
 c$$$     $           ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)  
2430           endif           endif
2431  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2432  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2564  c$$$     $           ,(db_cloud(iii),iii Line 2434  c$$$     $           ,(db_cloud(iii),iii
2434        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2435                
2436                
2437        if(DEBUG)then        if(nloop.lt.nstepy)then      
2438           print*,'---------------------- '          cutdistyz = cutdistyz+cutystep
2439            nloop     = nloop+1          
2440            goto 90                
2441          endif                    
2442          
2443          if(DEBUG.EQ.1)then
2444           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
          print*,' '  
2445        endif        endif
2446                
2447                
# Line 2590  c$$$     $           ,(db_cloud(iii),iii Line 2464  c$$$     $           ,(db_cloud(iii),iii
2464        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2465    
2466        include 'commontracker.f'        include 'commontracker.f'
2467          include 'level1.f'
2468        include 'common_momanhough.f'        include 'common_momanhough.f'
2469        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2470    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2471    
2472  *     output flag  *     output flag
2473  *     --------------  *     --------------
# Line 2614  c      common/dbg/DEBUG Line 2487  c      common/dbg/DEBUG
2487        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2488        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2489    
2490          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2491    
2492  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2493  *     classification of TRIPLETS  *     classification of TRIPLETS
2494  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2625  c      common/dbg/DEBUG Line 2500  c      common/dbg/DEBUG
2500        distance=0        distance=0
2501        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2502        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2503          nloop=0                  
2504     91   continue                  
2505        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2506           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
 c     print*,'--------------'  
 c     print*,'** ',itr1,' **'  
2507                    
2508           do icp=1,ncp_tot           do icp=1,ncp_tot
2509              cp_useds1(icp)=0              cp_useds1(icp)=0
# Line 2662  c         tr_temp(1)=itr1 Line 2537  c         tr_temp(1)=itr1
2537              do itr2=1,ntrpt     !loop (2) on TRIPLETS              do itr2=1,ntrpt     !loop (2) on TRIPLETS
2538                 if(itr2.eq.itr1)goto 11188       !next triplet                 if(itr2.eq.itr1)goto 11188       !next triplet
2539                 if(tr_used(itr2).eq.1)goto 11188 !next triplet                               if(tr_used(itr2).eq.1)goto 11188 !next triplet              
2540    
2541    
2542  *     triplet distance in parameter space  *     triplet distance in parameter space
2543  *     solo i due parametri spaziali per il momemnto  *     solo i due parametri spaziali per il momemnto
2544                 distance=                 distance=
2545       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2546       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2547                 distance = sqrt(distance)                 distance = sqrt(distance)
2548                  
2549                 if(distance.lt.cutdistxz)then  
2550  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  *     ------------------------------------------------------------------------
2551    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2552    *     ------------------------------------------------------------------------
2553    *     (added in august 2007)
2554                   istrimage=0
2555                   if(
2556         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2557         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2558         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2559         $              .true.)istrimage=1
2560    
2561                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2562                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2563                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2564                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
# Line 2689  c     print*,idb1,idb2,distance,' cloud Line 2577  c     print*,idb1,idb2,distance,' cloud
2577                    temp1 = temp1 + alfaxz1(itr2)                    temp1 = temp1 + alfaxz1(itr2)
2578                    temp2 = temp2 + alfaxz2(itr2)                    temp2 = temp2 + alfaxz2(itr2)
2579                    temp3 = temp3 + alfaxz3(itr2)                    temp3 = temp3 + alfaxz3(itr2)
 c     print*,'*   itrref,itr2 ',itrref,itr2,distance  
2580                 endif                               endif              
2581                                
2582  11188          continue  11188          continue
2583              enddo               !end loop (2) on TRIPLETS              enddo               !end loop (2) on TRIPLETS
2584                                                
2585  11888       continue  c11888       continue
2586                continue
2587           enddo                  !end loop on... bo?               enddo                  !end loop on... bo?    
2588                    
2589           nptloop=npv           nptloop=npv
# Line 2710  c     print*,'*   itrref,itr2 ',itrref,i Line 2598  c     print*,'*   itrref,itr2 ',itrref,i
2598  *     1bis)  *     1bis)
2599  *     2) it is not already stored  *     2) it is not already stored
2600  *     ------------------------------------------  *     ------------------------------------------
 c     print*,'check cp_used'  
2601           do ip=1,nplanes           do ip=1,nplanes
2602              hit_plane(ip)=0              hit_plane(ip)=0
2603           enddo           enddo
2604           ncpused=0           ncpused=0
2605           do icp=1,ncp_tot           do icp=1,ncp_tot
2606              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2607         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2608         $           .true.)then
2609                 ncpused=ncpused+1                 ncpused=ncpused+1
2610                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2611                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2726  c     print*,'check cp_used' Line 2615  c     print*,'check cp_used'
2615           do ip=1,nplanes           do ip=1,nplanes
2616              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2617           enddo           enddo
2618           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
          if(npt.lt.nptxz_min)goto 22288     !next triplet  
          if(nplused.lt.nplxz_min)goto 22288 !next doublet  
2619                    
2620  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2621  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2622           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2623              if(verbose)print*,              if(verbose.eq.1)print*,
2624       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2625       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2626       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2627  c     good2=.false.  c     good2=.false.
2628  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2629                do iv=1,nviews
2630    c               mask_view(iv) = 6
2631                   mask_view(iv) =  mask_view(iv) + 2**5
2632                enddo
2633              iflag=1              iflag=1
2634              return              return
2635           endif           endif
# Line 2756  c     goto 880         !fill ntp and go Line 2647  c     goto 880         !fill ntp and go
2647           enddo           enddo
2648           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2649                    
2650           if(DEBUG)then           if(DEBUG.EQ.1)then
2651              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
2652              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2653              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2654              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2655              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2656              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2657              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cpcloud_xz '
2658         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2659              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
 c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  
 c$$$     $           ,ptcloud_xz(nclouds_xz)  
 c$$$            print*,'nt-uple: tr_cloud(...) = '  
 c$$$     $           ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot)  
2660           endif           endif
2661  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2662  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2663  22288    continue  22288    continue
2664        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2665          
2666        if(DEBUG)then         if(nloop.lt.nstepx)then      
2667           print*,'---------------------- '           cutdistxz=cutdistxz+cutxstep
2668             nloop=nloop+1          
2669             goto 91                
2670           endif                    
2671          
2672          if(DEBUG.EQ.1)then
2673           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
          print*,' '  
2674        endif        endif
2675                
2676                
# Line 2796  c$$$     $           ,(tr_cloud(iii),iii Line 2688  c$$$     $           ,(tr_cloud(iii),iii
2688  **************************************************  **************************************************
2689    
2690        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2691    
2692        include 'commontracker.f'        include 'commontracker.f'
2693          include 'level1.f'
2694        include 'common_momanhough.f'        include 'common_momanhough.f'
2695        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2696        include 'common_mini_2.f'        include 'common_mini_2.f'
2697        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
2698    
2699  c      logical DEBUG  
 c      common/dbg/DEBUG  
2700    
2701  *     output flag  *     output flag
2702  *     --------------  *     --------------
# Line 2824  c      common/dbg/DEBUG Line 2712  c      common/dbg/DEBUG
2712  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2713  *     list of matching couples in the combination  *     list of matching couples in the combination
2714  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2715        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2716        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2717  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2718        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2719  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2720  *     variables for track fitting  *     variables for track fitting
2721        double precision AL_INI(5)        double precision AL_INI(5)
       double precision tath  
2722  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2723    
2724          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2725    
2726    
2727        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
2728                
2729        do iyz=1,nclouds_yz       !loop on YZ couds        do iyz=1,nclouds_yz       !loop on YZ clouds
2730           do ixz=1,nclouds_xz    !loop on XZ couds           do ixz=1,nclouds_xz    !loop on XZ clouds
2731                            
2732  *     --------------------------------------------------  *     --------------------------------------------------
2733  *     check of consistency of the clouds  *     check of consistency of the clouds
# Line 2849  c      real fitz(nplanes)        !z coor Line 2736  c      real fitz(nplanes)        !z coor
2736  *     of the two clouds  *     of the two clouds
2737  *     --------------------------------------------------  *     --------------------------------------------------
2738              do ip=1,nplanes              do ip=1,nplanes
2739                 hit_plane(ip)=0                 hit_plane(ip)=0  !n.matching couples (REAL couples, not SINGLETS)
2740                 ncp_match(ip)=0                 ncp_match(ip)=0  !n.matching couples per plane
2741                 do icpp=1,ncouplemax                 do icpp=1,ncouplemax
2742                    cp_match(ip,icpp)=0 !init couple list                    cp_match(ip,icpp)=0 !init couple list
2743                 enddo                 enddo
2744              enddo              enddo
2745              ncp_ok=0              ncp_ok=0            !count n.matching-couples
2746              do icp=1,ncp_tot    !loop on couples              ncpx_ok=0           !count n.matching-couples with x cluster
2747  *     get info on              ncpy_ok=0           !count n.matching-couples with y cluster
2748                 cpintersec(icp)=min(  
2749       $              cpcloud_yz(iyz,icp),  
2750       $              cpcloud_xz(ixz,icp))              do icp=1,ncp_tot    !loop over couples
2751                 if(  
2752       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.                 if(.not.RECOVER_SINGLETS)then
2753       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.  *     ------------------------------------------------------
2754       $              .false.)cpintersec(icp)=0  *     if NOT in RECOVER_SINGLETS mode, take the intersection
2755    *     between xz yz clouds
2756    *     ------------------------------------------------------
2757                      cpintersec(icp)=min(
2758         $                 cpcloud_yz(iyz,icp),
2759         $                 cpcloud_xz(ixz,icp))
2760    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2761    *     ------------------------------------------------------
2762    *     discard the couple if the sensor is in conflict
2763    *     ------------------------------------------------------
2764                      if(
2765         $       (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2766         $       (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2767         $                 .false.)cpintersec(icp)=0
2768                   else
2769    *     ------------------------------------------------------
2770    *     if RECOVER_SINGLETS take the union
2771    *     (otherwise the fake couples formed by singlets would be
2772    *     discarded)    
2773    *     ------------------------------------------------------
2774                      cpintersec(icp)=max(
2775         $                 cpcloud_yz(iyz,icp),
2776         $                 cpcloud_xz(ixz,icp))
2777    c$$$                  if(cpcloud_yz(iyz,icp).gt.0)
2778    c$$$     $                 cpintersec(icp)=cpcloud_yz(iyz,icp)
2779    *     cpintersec is >0 if either yz or xz cloud contains the couple icp
2780                   endif
2781    
2782    c$$$               print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
2783    
2784                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
                   ncp_ok=ncp_ok+1    
2785                                        
2786                    ip=ip_cp(icp)                    ip=ip_cp(icp)
2787                    hit_plane(ip)=1                    hit_plane(ip)=1
2788    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
2789    c$$$     $                 ncp_ok=ncp_ok+1  
2790    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
2791    c$$$     $                 ncpx_ok=ncpx_ok+1
2792    c$$$                  if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
2793    c$$$     $                 ncpy_ok=ncpy_ok+1
2794    
2795                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2796         $                 cpcloud_xz(ixz,icp).gt.0)
2797         $                 ncp_ok=ncp_ok+1  
2798                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2799         $                 cpcloud_xz(ixz,icp).eq.0)
2800         $                 ncpy_ok=ncpy_ok+1  
2801                      if(  cpcloud_yz(iyz,icp).eq.0.and.
2802         $                 cpcloud_xz(ixz,icp).gt.0)
2803         $                 ncpx_ok=ncpx_ok+1  
2804    
2805                    if(cpintersec(icp).eq.1)then                    if(cpintersec(icp).eq.1)then
2806  *     1) only the couple image in sensor 1 matches  *     1) only the couple image in sensor 1 matches
2807                       id=-icp                       id=-icp
# Line 2896  c      real fitz(nplanes)        !z coor Line 2828  c      real fitz(nplanes)        !z coor
2828              do ip=1,nplanes              do ip=1,nplanes
2829                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2830              enddo              enddo
2831                                        
2832              if(nplused.lt.nplxz_min)goto 888 !next doublet              if(nplused.lt.3)goto 888 !next combination
2833              if(ncp_ok.lt.ncpok)goto 888 !next cloud  ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
2834                ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
2835              if(DEBUG)then  *     -----------------------------------------------------------
2836                 print*,'Combination ',iyz,ixz  *     if in RECOVER_SINGLET mode, the two clouds must have
2837       $              ,' db ',ptcloud_yz(iyz)  *     at least ONE intersecting real couple
2838       $              ,' tr ',ptcloud_xz(ixz)  *     -----------------------------------------------------------
2839       $              ,'  -----> # matching couples ',ncp_ok              if(ncp_ok.lt.1)goto 888 !next combination
2840              endif  
2841  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'              if(DEBUG.EQ.1)then
2842  c$$$  print*,'Configurazione cluster XZ'                 print*,'////////////////////////////'
2843  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))                 print*,'Cloud combination (Y,X): ',iyz,ixz
2844  c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))                 print*,' db ',ptcloud_yz(iyz)
2845  c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))                 print*,' tr ',ptcloud_xz(ixz)
2846  c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))                 print*,'  -----> # matching couples ',ncp_ok
2847  c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))                 print*,'  -----> # fake couples (X)',ncpx_ok
2848  c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))                 print*,'  -----> # fake couples (Y)',ncpy_ok
2849  c$$$  print*,'Configurazione cluster YZ'                 do icp=1,ncp_tot
2850  c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))                    print*,'cp ',icp,' >'
2851  c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))       $                 ,' x ',cpcloud_xz(ixz,icp)
2852  c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))       $                 ,' y ',cpcloud_yz(iyz,icp)
2853  c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))       $                 ,' ==> ',cpintersec(icp)
2854  c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))                 enddo
2855  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))              endif
2856  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'                          
2857                            if(DEBUG.EQ.1)then
 *     -------> INITIAL GUESS <-------  
             AL_INI(1)=dreal(alfaxz1_av(ixz))  
             AL_INI(2)=dreal(alfayz1_av(iyz))  
             AL_INI(4)=datan(dreal(alfayz2_av(iyz))  
      $           /dreal(alfaxz2_av(ixz)))  
             tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  
             AL_INI(3)=tath/sqrt(1+tath**2)  
             AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  
               
 c     print*,'*******',AL_INI(5)  
             if(AL_INI(5).gt.defmax)goto 888 !next cloud  
               
 c     print*,'alfaxz2, alfayz2 '  
 c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  
               
 *     -------> INITIAL GUESS <-------  
 c     print*,'AL_INI ',(al_ini(i),i=1,5)  
               
             if(DEBUG)then  
2858                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2859                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2860                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2974  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2887  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2887                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2888                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2889                                                                
2890                                                                if(DEBUG.eq.1)
2891         $                             print*,'combination: '
2892         $                             ,cp_match(6,icp1)
2893         $                             ,cp_match(5,icp2)
2894         $                             ,cp_match(4,icp3)
2895         $                             ,cp_match(3,icp4)
2896         $                             ,cp_match(2,icp5)
2897         $                             ,cp_match(1,icp6)
2898    
2899    
2900    *                             ---------------------------------------
2901    *                             check if this group of couples has been
2902    *                             already fitted    
2903    *                             ---------------------------------------
2904                                  do ica=1,ntracks
2905                                     isthesame=1
2906                                     do ip=1,NPLANES
2907                                        if(hit_plane(ip).ne.0)then
2908                                           if(  CP_STORE(nplanes-ip+1,ica)
2909         $                                      .ne.
2910         $                                      cp_match(ip,hit_plane(ip)) )
2911         $                                      isthesame=0
2912                                        else
2913                                           if(  CP_STORE(nplanes-ip+1,ica)
2914         $                                      .ne.
2915         $                                      0 )
2916         $                                      isthesame=0
2917                                        endif
2918                                     enddo
2919                                     if(isthesame.eq.1)then
2920                                        if(DEBUG.eq.1)
2921         $                                   print*,'(already fitted)'
2922                                        goto 666 !jump to next combination
2923                                     endif
2924                                  enddo
2925    
2926                                call track_init !init TRACK common                                call track_init !init TRACK common
2927    
2928                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2929                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2930                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2931                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2991  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2939  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2939  *                                   *************************  *                                   *************************
2940  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2941  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2942    c                                    call xyz_PAM(icx,icy,is, !(1)
2943    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2944                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2945       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2946  *                                   *************************  *                                   *************************
2947  *                                   -----------------------------  *                                   -----------------------------
2948                                      xgood(nplanes-ip+1)=1.                                      if(icx.gt.0.and.icy.gt.0)then
2949                                      ygood(nplanes-ip+1)=1.                                         xgood(nplanes-ip+1)=1.
2950                                      xm(nplanes-ip+1)=xPAM                                         ygood(nplanes-ip+1)=1.
2951                                      ym(nplanes-ip+1)=yPAM                                         xm(nplanes-ip+1)=xPAM
2952                                      zm(nplanes-ip+1)=zPAM                                         ym(nplanes-ip+1)=yPAM
2953                                      resx(nplanes-ip+1)=resxPAM                                         zm(nplanes-ip+1)=zPAM
2954                                      resy(nplanes-ip+1)=resyPAM                                         resx(nplanes-ip+1)=resxPAM
2955                                           resy(nplanes-ip+1)=resyPAM
2956                                           if(DEBUG.EQ.1)print*,'(X,Y)'
2957         $                                      ,nplanes-ip+1,xPAM,yPAM
2958                                        else
2959                                           xm_A(nplanes-ip+1) = xPAM_A
2960                                           ym_A(nplanes-ip+1) = yPAM_A
2961                                           xm_B(nplanes-ip+1) = xPAM_B
2962                                           ym_B(nplanes-ip+1) = yPAM_B
2963                                           zm(nplanes-ip+1)
2964         $                                      = (zPAM_A+zPAM_B)/2.
2965                                           resx(nplanes-ip+1) = resxPAM
2966                                           resy(nplanes-ip+1) = resyPAM
2967                                           if(icx.eq.0.and.icy.gt.0)then
2968                                              xgood(nplanes-ip+1)=0.
2969                                              ygood(nplanes-ip+1)=1.
2970                                              resx(nplanes-ip+1) = 1000.
2971                                              if(DEBUG.EQ.1)print*,'(  Y)'
2972         $                                         ,nplanes-ip+1,xPAM,yPAM
2973                                           elseif(icx.gt.0.and.icy.eq.0)then
2974                                              xgood(nplanes-ip+1)=1.
2975                                              ygood(nplanes-ip+1)=0.
2976                                              if(DEBUG.EQ.1)print*,'(X  )'
2977         $                                         ,nplanes-ip+1,xPAM,yPAM
2978                                              resy(nplanes-ip+1) = 1000.
2979                                           else
2980                                              print*,'both icx=0 and icy=0'
2981         $                                         ,' ==> IMPOSSIBLE!!'
2982                                           endif
2983                                        endif
2984  *                                   -----------------------------  *                                   -----------------------------
2985                                   endif                                   endif
2986                                enddo !end loop on planes                                enddo !end loop on planes
2987  *     **********************************************************  *     **********************************************************
2988  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2989  *     **********************************************************  *     **********************************************************
2990    cccc  scommentare se si usa al_ini della nuvola
2991    c$$$                              do i=1,5
2992    c$$$                                 AL(i)=AL_INI(i)
2993    c$$$                              enddo
2994                                  call guess()
2995                                do i=1,5                                do i=1,5
2996                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2997                                enddo                                enddo
2998                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2999                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3000                                iprint=0                                iprint=0
3001                                if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
3002                                  if(DEBUG.EQ.1)iprint=2
3003                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
3004                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3005                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
3006                                      print *,                                      print *,
3007       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3008       $                              //'(mini2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
3009                                        print*,'initial guess: '
3010    
3011                                        print*,'AL_INI(1) = ',AL_INI(1)
3012                                        print*,'AL_INI(2) = ',AL_INI(2)
3013                                        print*,'AL_INI(3) = ',AL_INI(3)
3014                                        print*,'AL_INI(4) = ',AL_INI(4)
3015                                        print*,'AL_INI(5) = ',AL_INI(5)
3016                                   endif                                   endif
3017                                   chi2=-chi2  c                                 chi2=-chi2
3018                                endif                                endif
3019  *     **********************************************************  *     **********************************************************
3020  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
3021  *     **********************************************************  *     **********************************************************
3022    
3023                                if(chi2.le.0.)goto 666                                              if(chi2.le.0.)goto 666              
3024                                  if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3025                                  if(chi2.ne.chi2)goto 666  !OPTIMIZATION
3026    
3027  *     --------------------------  *     --------------------------
3028  *     STORE candidate TRACK INFO  *     STORE candidate TRACK INFO
3029  *     --------------------------  *     --------------------------
3030                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3031                                                                    
3032                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
3033       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3034       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3035       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3036  c                                 good2=.false.  c                                 good2=.false.
3037  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3038                                     do iv=1,nviews
3039    c                                    mask_view(iv) = 7
3040                                        mask_view(iv) = mask_view(iv) + 2**6
3041                                     enddo
3042                                   iflag=1                                   iflag=1
3043                                   return                                   return
3044                                endif                                endif
3045                                                                
3046                                ntracks = ntracks + 1                                ntracks = ntracks + 1
3047                                                                
3048  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
3049                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
3050                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3051                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
3052                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))
3053                                   XM_STORE(ip,ntracks)=sngl(xm(ip))                                   XM_STORE(ip,ntracks)=sngl(xm(ip))
3054                                   YM_STORE(ip,ntracks)=sngl(ym(ip))                                   YM_STORE(ip,ntracks)=sngl(ym(ip))
3055                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))
# Line 3067  c$$$     $                               Line 3062  c$$$     $                              
3062                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
3063                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
3064                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
3065    *                                NB! hit_plane is defined from bottom to top
3066                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3067                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
3068       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
3069                                        SENSOR_STORE(nplanes-ip+1,ntracks)
3070         $                              = is_cp(cp_match(ip,hit_plane(ip)))
3071                                        
3072                                        icl=
3073         $                                   clx(ip,icp_cp(
3074         $                                   cp_match(ip,hit_plane(ip)
3075         $                                   )));
3076                                        if(icl.eq.0)
3077         $                                   icl=
3078         $                                   cly(ip,icp_cp(
3079         $                                   cp_match(ip,hit_plane(ip)
3080         $                                   )));
3081    
3082                                        LADDER_STORE(nplanes-ip+1,ntracks)
3083         $                                   = LADDER(icl);
3084                                   else                                   else
3085                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
3086                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
3087                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
3088                                   endif                                   endif
3089                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
3090                                     BY_STORE(ip,ntracks)=0!I dont need it now
3091                                     CLS_STORE(ip,ntracks)=0
3092                                   do i=1,5                                   do i=1,5
3093                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
3094                                   enddo                                   enddo
3095                                enddo                                enddo
3096                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
3097                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
3098                                                                
3099  *     --------------------------------  *     --------------------------------
# Line 3103  c$$$  rchi2=chi2/dble(ndof) Line 3114  c$$$  rchi2=chi2/dble(ndof)
3114                
3115        if(ntracks.eq.0)then        if(ntracks.eq.0)then
3116           iflag=1           iflag=1
3117           return  cc         return
3118        endif        endif
3119                
3120        if(DEBUG)then        if(DEBUG.EQ.1)then
3121           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
3122           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
3123           do i=1,ntracks          do i=1,ntracks
3124              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
3125       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
3126           enddo              ndof=ndof           !(1)
3127           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
3128         $           +int(ygood_store(ii,i)) !(1)
3129              enddo                 !(1)
3130              print*,i,' --- ',rchi2_store(i),' --- '
3131         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3132            enddo
3133            print*,'*****************************************'
3134        endif        endif
3135                
3136                
# Line 3132  c$$$  rchi2=chi2/dble(ndof) Line 3149  c$$$  rchi2=chi2/dble(ndof)
3149    
3150        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3151    
 c******************************************************  
 cccccc 06/10/2005 modified by elena vannuccini ---> (1)  
 cccccc 31/01/2006 modified by elena vannuccini ---> (2)  
 cccccc 12/08/2006 modified by elena vannucicni ---> (3)  
 c******************************************************  
3152    
3153        include 'commontracker.f'        include 'commontracker.f'
3154          include 'level1.f'
3155        include 'common_momanhough.f'        include 'common_momanhough.f'
3156        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3157        include 'common_mini_2.f'        include 'common_mini_2.f'
3158        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
       include 'level1.f'  
3159        include 'calib.f'        include 'calib.f'
3160    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
3161  *     flag to chose PFA  *     flag to chose PFA
3162        character*10 PFA        character*10 PFA
3163        common/FINALPFA/PFA        common/FINALPFA/PFA
3164    
3165          real k(6)
3166          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3167    
3168          real xp,yp,zp
3169          real xyzp(3),bxyz(3)
3170          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3171    
3172          if(DEBUG.EQ.1)print*,'refine_track:'
3173  *     =================================================  *     =================================================
3174  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3175  *                          and  *                          and
# Line 3162  c      common/dbg/DEBUG Line 3178  c      common/dbg/DEBUG
3178        call track_init        call track_init
3179        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3180    
3181             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3182    
3183             xP=XV_STORE(nplanes-ip+1,ibest)
3184             yP=YV_STORE(nplanes-ip+1,ibest)
3185             zP=ZV_STORE(nplanes-ip+1,ibest)
3186             call gufld(xyzp,bxyz)
3187             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3188             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3189    c$$$  bxyz(1)=0
3190    c$$$         bxyz(2)=0
3191    c$$$         bxyz(3)=0
3192    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3193  *     -------------------------------------------------  *     -------------------------------------------------
3194  *     If the plane has been already included, it just  *     If the plane has been already included, it just
3195  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
3196  *     using improved PFAs  *     using improved PFAs
3197  *     -------------------------------------------------  *     -------------------------------------------------
3198           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3199    c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3200    c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3201    c$$$            
3202    c$$$            id=CP_STORE(nplanes-ip+1,ibest)
3203    c$$$            
3204    c$$$            is=is_cp(id)
3205    c$$$            icp=icp_cp(id)
3206    c$$$            if(ip_cp(id).ne.ip)
3207    c$$$     $           print*,'OKKIO!!'
3208    c$$$     $           ,'id ',id,is,icp
3209    c$$$     $           ,ip_cp(id),ip
3210    c$$$            icx=clx(ip,icp)
3211    c$$$            icy=cly(ip,icp)
3212    c$$$c            call xyz_PAM(icx,icy,is,
3213    c$$$c     $           PFA,PFA,
3214    c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
3215    c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
3216    c$$$            call xyz_PAM(icx,icy,is,
3217    c$$$     $           PFA,PFA,
3218    c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
3219    c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
3220    c$$$     $           bxyz(1),
3221    c$$$     $           bxyz(2)
3222    c$$$     $           )
3223    c$$$
3224    c$$$            xm(nplanes-ip+1) = xPAM
3225    c$$$            ym(nplanes-ip+1) = yPAM
3226    c$$$            zm(nplanes-ip+1) = zPAM
3227    c$$$            xgood(nplanes-ip+1) = 1
3228    c$$$            ygood(nplanes-ip+1) = 1
3229    c$$$            resx(nplanes-ip+1) = resxPAM
3230    c$$$            resy(nplanes-ip+1) = resyPAM
3231    c$$$
3232    c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3233    c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3234             if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
3235       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3236                            
3237              id=CP_STORE(nplanes-ip+1,ibest)              id=CP_STORE(nplanes-ip+1,ibest)
# Line 3180  c      common/dbg/DEBUG Line 3244  c      common/dbg/DEBUG
3244       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3245              icx=clx(ip,icp)              icx=clx(ip,icp)
3246              icy=cly(ip,icp)              icy=cly(ip,icp)
3247    c            call xyz_PAM(icx,icy,is,
3248    c     $           PFA,PFA,
3249    c     $           AXV_STORE(nplanes-ip+1,ibest),
3250    c     $           AYV_STORE(nplanes-ip+1,ibest))
3251              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3252       $           PFA,PFA,       $           PFA,PFA,
3253       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3254       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3255  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3256  c$$$  $              'COG2','COG2',       $           bxyz(2)
3257  c$$$  $              0.,       $           )
3258  c$$$  $              0.)  
3259              xm(nplanes-ip+1) = xPAM              if(icx.gt.0.and.icy.gt.0)then
3260              ym(nplanes-ip+1) = yPAM                 xm(nplanes-ip+1) = xPAM
3261              zm(nplanes-ip+1) = zPAM                 ym(nplanes-ip+1) = yPAM
3262              xgood(nplanes-ip+1) = 1                 zm(nplanes-ip+1) = zPAM
3263              ygood(nplanes-ip+1) = 1                 xm_A(nplanes-ip+1) = 0.
3264              resx(nplanes-ip+1) = resxPAM                 ym_A(nplanes-ip+1) = 0.
3265              resy(nplanes-ip+1) = resyPAM                 xm_B(nplanes-ip+1) = 0.
3266                   ym_B(nplanes-ip+1) = 0.
3267  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)                 xgood(nplanes-ip+1) = 1
3268              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                 ygood(nplanes-ip+1) = 1
3269              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                 resx(nplanes-ip+1) = resxPAM
3270                   resy(nplanes-ip+1) = resyPAM
3271                   dedxtrk_x(nplanes-ip+1)=
3272         $              sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3273                   dedxtrk_y(nplanes-ip+1)=
3274         $              sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3275                else
3276                   xm(nplanes-ip+1) = 0.
3277                   ym(nplanes-ip+1) = 0.
3278                   zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
3279                   xm_A(nplanes-ip+1) = xPAM_A
3280                   ym_A(nplanes-ip+1) = yPAM_A
3281                   xm_B(nplanes-ip+1) = xPAM_B
3282                   ym_B(nplanes-ip+1) = yPAM_B
3283                   xgood(nplanes-ip+1) = 0
3284                   ygood(nplanes-ip+1) = 0
3285                   resx(nplanes-ip+1) = 1000.!resxPAM
3286                   resy(nplanes-ip+1) = 1000.!resyPAM
3287                   dedxtrk_x(nplanes-ip+1)= 0
3288                   dedxtrk_y(nplanes-ip+1)= 0
3289                   if(icx.gt.0)then
3290                      xgood(nplanes-ip+1) = 1
3291                      resx(nplanes-ip+1) = resxPAM
3292                      dedxtrk_x(nplanes-ip+1)=
3293         $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3294                   elseif(icy.gt.0)then
3295                      ygood(nplanes-ip+1) = 1
3296                      resy(nplanes-ip+1) = resyPAM
3297                      dedxtrk_y(nplanes-ip+1)=
3298         $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3299                   endif
3300                endif
3301                            
3302    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3303  *     -------------------------------------------------  *     -------------------------------------------------
3304  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3305  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3306  *     -------------------------------------------------  *     -------------------------------------------------
3307    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3308           else                             else                  
3309                                
3310              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
3311              ygood(nplanes-ip+1)=0              ygood(nplanes-ip+1)=0
3312    
3313                CP_STORE(nplanes-ip+1,ibest)=0 !re-init
3314                CLS_STORE(nplanes-ip+1,ibest)=0
3315    
3316                                
3317  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3318  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
             xP=XV_STORE(nplanes-ip+1,ibest)  
             yP=YV_STORE(nplanes-ip+1,ibest)  
             zP=ZV_STORE(nplanes-ip+1,ibest)  
3319              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3320  *     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
3321              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3322    
3323                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3324                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3325  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3326    
3327              if(DEBUG)then              if(DEBUG.EQ.1)then
3328                 print*,                 print*,
3329       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3330       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3231  c            dedxtrk(nplanes-ip+1) = (de Line 3335  c            dedxtrk(nplanes-ip+1) = (de
3335  *     ===========================================  *     ===========================================
3336  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3337  *     ===========================================  *     ===========================================
3338  c            if(DEBUG)print*,'>>>> try to include a new couple'              distmin=100000000.
             distmin=1000000.  
3339              xmm = 0.              xmm = 0.
3340              ymm = 0.              ymm = 0.
3341              zmm = 0.              zmm = 0.
# Line 3245  c            if(DEBUG)print*,'>>>> try t Line 3348  c            if(DEBUG)print*,'>>>> try t
3348              do icp=1,ncp_plane(ip) !loop on couples on plane icp              do icp=1,ncp_plane(ip) !loop on couples on plane icp
3349                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3350                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3351                   if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
3352                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3353  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3354  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3355       $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)       $              cl_used(icx).ne.0.or. !or the X cluster is already used
3356       $              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
3357       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3358  *            *          
3359                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3360       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3361       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3362       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3363         $              bxyz(1),
3364         $              bxyz(2)
3365         $              )
3366                                
3367                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3368    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3369                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3370                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)
3371       $              ,' ) normalized distance ',distance       $              print*,'( couple ',id
3372         $              ,' ) distance ',distance
3373                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3374                    xmm = xPAM                    xmm = xPAM
3375                    ymm = yPAM                    ymm = yPAM
# Line 3270  c     $              'ETA2','ETA2', Line 3378  c     $              'ETA2','ETA2',
3378                    rymm = resyPAM                    rymm = resyPAM
3379                    distmin = distance                    distmin = distance
3380                    idm = id                                      idm = id                  
3381  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3382                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3383                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    clincnewc=10*sqrt(rymm**2+rxmm**2
3384         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
3385                 endif                 endif
3386   1188          continue   1188          continue
3387              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3388              if(distmin.le.clinc)then                                if(distmin.le.clincnewc)then    
3389  *              -----------------------------------  *              -----------------------------------
3390                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3391                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3392                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3393                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3394                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3395                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3396                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3397  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3398                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3399  *              -----------------------------------  *              -----------------------------------
3400                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3401                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3402       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
3403                 goto 133         !next plane                 goto 133         !next plane
3404              endif              endif
3405  *     ================================================  *     ================================================
3406  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3407  *                     either from a couple or single  *                     either from a couple or single
3408  *     ================================================  *     ================================================
 c            if(DEBUG)print*,'>>>> try to include a new cluster'  
3409              distmin=1000000.              distmin=1000000.
3410              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3411              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3317  c            if(DEBUG)print*,'>>>> try t Line 3424  c            if(DEBUG)print*,'>>>> try t
3424              do icp=1,ncp_plane(ip) !loop on cluster inside couples              do icp=1,ncp_plane(ip) !loop on cluster inside couples
3425                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3426                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3427                   if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
3428                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3429                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
3430  *                                                !jump to the next couple  *                                                !jump to the next couple
# Line 3324  c            if(DEBUG)print*,'>>>> try t Line 3432  c            if(DEBUG)print*,'>>>> try t
3432  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
3433                 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)
3434  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3435    c               call xyz_PAM(icx,0,ist,
3436    c     $              PFA,PFA,
3437    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3438                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3439       $              PFA,PFA,       $              PFA,PFA,
3440       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3441         $              bxyz(1),
3442         $              bxyz(2)
3443         $              )              
3444                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3445  c     if(DEBUG)print*,'normalized distance ',distance  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3446                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)
3447       $              ,' in cp ',id,' ) normalized distance ',distance       $              print*,'( cl-X ',icx
3448         $              ,' in cp ',id,' ) distance ',distance
3449                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3450                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3451                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3343  c     if(DEBUG)print*,'normalized distan Line 3457  c     if(DEBUG)print*,'normalized distan
3457                    rymm = resyPAM                    rymm = resyPAM
3458                    distmin = distance                    distmin = distance
3459                    iclm = icx                    iclm = icx
3460  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3461                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3462                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3463                 endif                                   endif                  
3464  11881          continue  11881          continue
# Line 3352  c                  dedxmm = dedx(icx) !( Line 3466  c                  dedxmm = dedx(icx) !(
3466  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
3467                 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)
3468  *                                              !jump to the next couple  *                                              !jump to the next couple
3469    c               call xyz_PAM(0,icy,ist,
3470    c     $              PFA,PFA,
3471    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3472                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3473       $              PFA,PFA,       $              PFA,PFA,
3474       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3475         $              bxyz(1),
3476         $              bxyz(2)
3477         $              )
3478                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3479                 if(DEBUG)print*,'( cl-Y ',icy  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3480       $              ,' in cp ',id,' ) normalized distance ',distance                 if(DEBUG.EQ.1)
3481         $              print*,'( cl-Y ',icy
3482         $              ,' in cp ',id,' ) distance ',distance
3483                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3484                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3485                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3370  c     $              'ETA2','ETA2', Line 3491  c     $              'ETA2','ETA2',
3491                    rymm = resyPAM                    rymm = resyPAM
3492                    distmin = distance                    distmin = distance
3493                    iclm = icy                    iclm = icy
3494  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3495                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3496                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3497                 endif                                   endif                  
3498  11882          continue  11882          continue
3499              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3500  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3501              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3502                 icl=cls(ip,ic)                 icl=cls(ip,ic)
 c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  
3503                 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)
3504       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3505       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3506                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3507                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3508       $                 PFA,PFA,       $                 PFA,PFA,
3509       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3510         $                 bxyz(1),
3511         $                 bxyz(2)
3512         $                 )
3513                 else                         !<---- Y view                 else                         !<---- Y view
3514                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3515       $                 PFA,PFA,       $                 PFA,PFA,
3516       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3517         $                 bxyz(1),
3518         $                 bxyz(2)
3519         $                 )
3520                 endif                 endif
3521    
3522                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3523                 if(DEBUG)print*,'( cl-s ',icl  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3524       $              ,' ) normalized distance ',distance                 if(DEBUG.EQ.1)
3525         $              print*,'( cl-s ',icl
3526         $              ,' ) distance ',distance
3527                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3528                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3529                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3409  c     $                 'ETA2','ETA2', Line 3535  c     $                 'ETA2','ETA2',
3535                    rymm = resyPAM                    rymm = resyPAM
3536                    distmin = distance                      distmin = distance  
3537                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3538                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3539                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3540                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3541                    else          !<---- Y view                    else          !<---- Y view
3542                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3543                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3544                    endif                    endif
3545                 endif                                   endif                  
3546  18882          continue  18882          continue
3547              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3548    
3549              if(distmin.le.clinc)then                                if(iclm.ne.0)then
                 
                CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<      
 *              ----------------------------  
3550                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3551                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3552                    resx(nplanes-ip+1)=rxmm       $                 20*
3553                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3554       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'                 else if(mod(VIEW(iclm),2).ne.0)then
3555                 else                    clincnew=
3556                    YGOOD(nplanes-ip+1)=1.       $                 10*
3557                    resy(nplanes-ip+1)=rymm       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3558                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 endif
3559       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  
3560                   if(distmin.le.clincnew)then  
3561                      
3562                      CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3563    *     ----------------------------
3564                      if(mod(VIEW(iclm),2).eq.0)then
3565                         XGOOD(nplanes-ip+1)=1.
3566                         resx(nplanes-ip+1)=rxmm
3567                         if(DEBUG.EQ.1)
3568         $                    print*,'%%%% included X-cl ',iclm
3569         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3570         $                    ,', dist.= ',distmin
3571         $                    ,', cut ',clincnew,' )'
3572                      else
3573                         YGOOD(nplanes-ip+1)=1.
3574                         resy(nplanes-ip+1)=rymm
3575                         if(DEBUG.EQ.1)
3576         $                    print*,'%%%% included Y-cl ',iclm
3577         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3578         $                    ,', dist.= ', distmin
3579         $                    ,', cut ',clincnew,' )'
3580                      endif
3581    *     ----------------------------
3582                      xm_A(nplanes-ip+1) = xmm_A
3583                      ym_A(nplanes-ip+1) = ymm_A
3584                      xm_B(nplanes-ip+1) = xmm_B
3585                      ym_B(nplanes-ip+1) = ymm_B
3586                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3587                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3588                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3589    *     ----------------------------
3590                 endif                 endif
 *              ----------------------------  
                xm_A(nplanes-ip+1) = xmm_A  
                ym_A(nplanes-ip+1) = ymm_A  
                xm_B(nplanes-ip+1) = xmm_B  
                ym_B(nplanes-ip+1) = ymm_B  
                zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.  
 c              dedxtrk(nplanes-ip+1) = dedxmm !<<<    !(1)  
                dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)  
                dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)  
 *              ----------------------------  
3591              endif              endif
3592           endif           endif
3593   133     continue   133     continue
# Line 3456  c              dedxtrk(nplanes-ip+1) = d Line 3598  c              dedxtrk(nplanes-ip+1) = d
3598        return        return
3599        end        end
3600    
3601    
3602  ***************************************************  ***************************************************
3603  *                                                 *  *                                                 *
3604  *                                                 *  *                                                 *
# Line 3464  c              dedxtrk(nplanes-ip+1) = d Line 3607  c              dedxtrk(nplanes-ip+1) = d
3607  *                                                 *  *                                                 *
3608  *                                                 *  *                                                 *
3609  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3610  *  *
       subroutine clean_XYclouds(ibest,iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
   
       do ip=1,nplanes           !loop on planes  
   
          id=CP_STORE(nplanes-ip+1,ibest)  
          icl=CLS_STORE(nplanes-ip+1,ibest)  
          if(id.ne.0.or.icl.ne.0)then                
             if(id.ne.0)then  
                iclx=clx(ip,icp_cp(id))  
                icly=cly(ip,icp_cp(id))  
 c               cl_used(iclx)=1  !tag used clusters  
 c               cl_used(icly)=1  !tag used clusters  
                cl_used(iclx)=ntrk  !tag used clusters !(1)  
                cl_used(icly)=ntrk  !tag used clusters !(1)  
             elseif(icl.ne.0)then  
 c               cl_used(icl)=1   !tag used clusters  
                cl_used(icl)=ntrk   !tag used clusters !1)  
             endif  
               
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
 *     -----------------------------  
 *     remove the couple from clouds  
 *     remove also vitual couples containing the  
 *     selected clusters  
 *     -----------------------------  
             do icp=1,ncp_plane(ip)  
                if(  
      $              clx(ip,icp).eq.iclx  
      $              .or.  
      $              clx(ip,icp).eq.icl  
      $              .or.  
      $              cly(ip,icp).eq.icly  
      $              .or.  
      $              cly(ip,icp).eq.icl  
      $              )then  
                   id=id_cp(ip,icp,1)  
                   if(DEBUG)then  
                      print*,ip,' <<< cp ',id  
      $                    ,' ( cl-x '  
      $                    ,clx(ip,icp)  
      $                    ,' cl-y '  
      $                    ,cly(ip,icp),' ) --> removed'  
                   endif  
 *     -----------------------------  
 *     remove the couple from clouds  
                   do iyz=1,nclouds_yz  
                      if(cpcloud_yz(iyz,abs(id)).ne.0)then  
                         ptcloud_yz(iyz)=ptcloud_yz(iyz)-1  
                         cpcloud_yz(iyz,abs(id))=0  
                      endif  
                   enddo  
                   do ixz=1,nclouds_xz  
                      if(cpcloud_xz(ixz,abs(id)).ne.0)then  
                         ptcloud_xz(ixz)=ptcloud_xz(ixz)-1  
                         cpcloud_xz(ixz,abs(id))=0  
                      endif  
                   enddo                      
 *     -----------------------------  
                endif  
             enddo  
               
          endif                
       enddo                     !end loop on planes  
         
       return  
       end  
   
   
   
   
 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$$$  
3611    
3612    
3613    
# Line 3652  c$$$ Line 3615  c$$$
3615    
3616        subroutine init_level2        subroutine init_level2
3617    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3618        include 'commontracker.f'        include 'commontracker.f'
3619          include 'level1.f'
3620        include 'common_momanhough.f'        include 'common_momanhough.f'
3621        include 'level2.f'        include 'level2.f'
       include 'level1.f'  
3622    
3623    *     ---------------------------------
3624    *     variables initialized from level1
3625    *     ---------------------------------
3626        do i=1,nviews        do i=1,nviews
3627           good2(i)=good1(i)           good2(i)=good1(i)
3628             do j=1,nva1_view
3629                vkflag(i,j)=1
3630                if(cnnev(i,j).le.0)then
3631                   vkflag(i,j)=cnnev(i,j)
3632                endif
3633             enddo
3634        enddo        enddo
3635    *     ----------------
3636  c      good2 = 0!.false.  *     level2 variables
3637  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*****************************************************  
   
3638        NTRK = 0        NTRK = 0
3639        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3640           IMAGE(IT)=0           IMAGE(IT)=0
3641           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3642           do ip=1,nplanes           do ip=1,nplanes
3643              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3644              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
3645              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3646              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3647              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3648                TAILX_nt(IP,IT) = 0
3649                TAILY_nt(IP,IT) = 0
3650                XBAD(IP,IT) = 0
3651                YBAD(IP,IT) = 0
3652              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3653              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3654  c*****************************************************              LS(IP,IT) = 0
 cccccc 11/9/2005 modified by david fedele  
3655              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3656              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3657              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3658              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3659                multmaxx(ip,it) = 0
3660                seedx(ip,it)    = 0
3661                xpu(ip,it)      = 0
3662                multmaxy(ip,it) = 0
3663                seedy(ip,it)    = 0
3664                ypu(ip,it)      = 0
3665           enddo           enddo
3666           do ipa=1,5           do ipa=1,5
3667              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3717  cccccc 17/8/2006 modified by elena Line 3670  cccccc 17/8/2006 modified by elena
3670              enddo                                enddo                  
3671           enddo                             enddo                  
3672        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3673        nclsx=0        nclsx=0
3674        nclsy=0              nclsy=0      
3675        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3676          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3677          xs(1,ip)=0          xs(1,ip)=0
3678          xs(2,ip)=0          xs(2,ip)=0
3679          sgnlxs(ip)=0          sgnlxs(ip)=0
3680          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3681          ys(1,ip)=0          ys(1,ip)=0
3682          ys(2,ip)=0          ys(2,ip)=0
3683          sgnlys(ip)=0          sgnlys(ip)=0
3684            sxbad(ip)=0
3685            sybad(ip)=0
3686            multmaxsx(ip)=0
3687            multmaxsy(ip)=0
3688        enddo        enddo
 c*******************************************************  
3689        end        end
3690    
3691    
# Line 3750  c*************************************** Line 3700  c***************************************
3700  ************************************************************  ************************************************************
3701    
3702    
3703          subroutine init_hough
3704    
3705          include 'commontracker.f'
3706          include 'level1.f'
3707          include 'common_momanhough.f'
3708          include 'common_hough.f'
3709          include 'level2.f'
3710    
3711          ntrpt_nt=0
3712          ndblt_nt=0
3713          NCLOUDS_XZ_nt=0
3714          NCLOUDS_YZ_nt=0
3715          do idb=1,ndblt_max_nt
3716             db_cloud_nt(idb)=0
3717             alfayz1_nt(idb)=0      
3718             alfayz2_nt(idb)=0      
3719          enddo
3720          do itr=1,ntrpt_max_nt
3721             tr_cloud_nt(itr)=0
3722             alfaxz1_nt(itr)=0      
3723             alfaxz2_nt(itr)=0      
3724             alfaxz3_nt(itr)=0      
3725          enddo
3726          do idb=1,ncloyz_max      
3727            ptcloud_yz_nt(idb)=0    
3728            alfayz1_av_nt(idb)=0    
3729            alfayz2_av_nt(idb)=0    
3730          enddo                    
3731          do itr=1,ncloxz_max      
3732            ptcloud_xz_nt(itr)=0    
3733            alfaxz1_av_nt(itr)=0    
3734            alfaxz2_av_nt(itr)=0    
3735            alfaxz3_av_nt(itr)=0    
3736          enddo                    
3737    
3738          ntrpt=0                  
3739          ndblt=0                  
3740          NCLOUDS_XZ=0              
3741          NCLOUDS_YZ=0              
3742          do idb=1,ndblt_max        
3743            db_cloud(idb)=0        
3744            cpyz1(idb)=0            
3745            cpyz2(idb)=0            
3746            alfayz1(idb)=0          
3747            alfayz2(idb)=0          
3748          enddo                    
3749          do itr=1,ntrpt_max        
3750            tr_cloud(itr)=0        
3751            cpxz1(itr)=0            
3752            cpxz2(itr)=0            
3753            cpxz3(itr)=0            
3754            alfaxz1(itr)=0          
3755            alfaxz2(itr)=0          
3756            alfaxz3(itr)=0          
3757          enddo                    
3758          do idb=1,ncloyz_max      
3759            ptcloud_yz(idb)=0      
3760            alfayz1_av(idb)=0      
3761            alfayz2_av(idb)=0      
3762            do idbb=1,ncouplemaxtot
3763              cpcloud_yz(idb,idbb)=0
3764            enddo                  
3765          enddo                    
3766          do itr=1,ncloxz_max      
3767            ptcloud_xz(itr)=0      
3768            alfaxz1_av(itr)=0      
3769            alfaxz2_av(itr)=0      
3770            alfaxz3_av(itr)=0      
3771            do itrr=1,ncouplemaxtot
3772              cpcloud_xz(itr,itrr)=0
3773            enddo                  
3774          enddo                    
3775          end
3776    ************************************************************
3777    *
3778    *
3779    *
3780    *
3781    *
3782    *
3783    *
3784    ************************************************************
3785    
3786    
3787        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3788    
3789  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3761  c*************************************** Line 3795  c***************************************
3795            
3796        include 'commontracker.f'        include 'commontracker.f'
3797        include 'level1.f'        include 'level1.f'
3798          include 'common_momanhough.f'
3799        include 'level2.f'        include 'level2.f'
3800        include 'common_mini_2.f'        include 'common_mini_2.f'
3801        include 'common_momanhough.f'        include 'calib.f'
3802        real sinth,phi,pig        !(4)  
3803          character*10 PFA
3804          common/FINALPFA/PFA
3805    
3806          real sinth,phi,pig
3807          integer ssensor,sladder
3808        pig=acos(-1.)        pig=acos(-1.)
3809    
3810  c      good2=1!.true.  
3811    
3812    *     -------------------------------------
3813        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3814        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3815    *     -------------------------------------
3816        phi   = al(4)             !(4)        phi   = al(4)          
3817        sinth = al(3)             !(4)        sinth = al(3)            
3818        if(sinth.lt.0)then        !(4)        if(sinth.lt.0)then      
3819           sinth = -sinth         !(4)           sinth = -sinth        
3820           phi = phi + pig        !(4)           phi = phi + pig      
3821        endif                     !(4)        endif                    
3822        npig = aint(phi/(2*pig))  !(4)        npig = aint(phi/(2*pig))
3823        phi = phi - npig*2*pig    !(4)        phi = phi - npig*2*pig  
3824        if(phi.lt.0)              !(4)        if(phi.lt.0)            
3825       $     phi = phi + 2*pig    !(4)       $     phi = phi + 2*pig  
3826        al(4) = phi               !(4)        al(4) = phi              
3827        al(3) = sinth             !(4)        al(3) = sinth            
 *****************************************************  
3828        do i=1,5        do i=1,5
3829           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3830           do j=1,5           do j=1,5
3831              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3832           enddo           enddo
 c     print*,al_nt(i,ntr)  
3833        enddo        enddo
3834          *     -------------------------------------      
3835        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3836           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3837           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3800  c     print*,al_nt(i,ntr) Line 3840  c     print*,al_nt(i,ntr)
3840           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3841           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3842           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3843             TAILX_nt(IP,ntr) = 0.
3844             TAILY_nt(IP,ntr) = 0.
3845           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3846           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3847           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3848           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3849           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3850  c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3851           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           factor = sqrt(
3852           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3853         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3854         $        1. )
3855    
3856             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3857             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3858        
3859           id  = CP_STORE(ip,IDCAND)  
3860    ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3861    
3862             ax   = axv_nt(ip,ntr)
3863             ay   = ayv_nt(ip,ntr)
3864             bfx  = BX_STORE(ip,IDCAND)
3865             bfy  = BY_STORE(ip,IDCAND)
3866    c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3867    c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3868    c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3869    c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3870    c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3871    c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3872    
3873             angx = effectiveangle(ax,2*ip,bfy)
3874             angy = effectiveangle(ay,2*ip-1,bfx)
3875            
3876            
3877    
3878             id  = CP_STORE(ip,IDCAND) ! couple id
3879           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3880           if(id.ne.0)then           ssensor = -1
3881             sladder = -1
3882             ssensor = SENSOR_STORE(ip,IDCAND)
3883             sladder = LADDER_STORE(ip,IDCAND)
3884             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3885             LS(IP,ntr)      = ssensor+10*sladder
3886    
3887    c         if(id.ne.0)then
3888    CCCCCC(10 novembre 2009) PATCH X NUCLEI
3889    C     non ho capito perche', ma durante il ritracciamento dei nuclei
3890    C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3891    C     che non viene reinizializzata correttamente e i cluster esclusi
3892    C     dal fit risultano ancora inclusi...
3893    C
3894             cltrx(ip,ntr)   = 0
3895             cltry(ip,ntr)   = 0
3896             if(
3897         $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3898         $        .and.
3899         $        id.ne.0)then
3900    
3901    c           >>> is a couple
3902              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3903              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3904  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3905                if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
3906    
3907                   cl_used(cltrx(ip,ntr)) = 1 !tag used clusters          
3908    
3909                   xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3910    
3911                   if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3912         $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3913                  
3914                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3915         $              +10000*mult(cltrx(ip,ntr))
3916                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3917         $              /clsigma(indmax(cltrx(ip,ntr)))
3918                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3919                   xpu(ip,ntr)      = corr
3920    
3921                endif
3922                if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then
3923    
3924                   cl_used(cltry(ip,ntr)) = 1 !tag used clusters          
3925    
3926                   ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3927    
3928                   if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3929         $              dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3930                  
3931                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3932         $              +10000*mult(cltry(ip,ntr))
3933                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3934         $              /clsigma(indmax(cltry(ip,ntr)))
3935                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3936                   ypu(ip,ntr)      = corr
3937                endif
3938    
3939           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3940              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3941              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl              cl_used(icl) = 1    !tag used clusters          
3942  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3943                if(mod(VIEW(icl),2).eq.0)then
3944                   cltrx(ip,ntr)=icl
3945                   xbad(ip,ntr) = nbadstrips(4,icl)
3946    
3947                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3948    
3949                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3950         $                         +10000*mult(cltrx(ip,ntr))
3951                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3952         $           /clsigma(indmax(cltrx(ip,ntr)))
3953                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3954                   xpu(ip,ntr)      = corr
3955    
3956                elseif(mod(VIEW(icl),2).eq.1)then
3957                   cltry(ip,ntr)=icl
3958                   ybad(ip,ntr) = nbadstrips(4,icl)
3959    
3960                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3961    
3962                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3963         $                         +10000*mult(cltry(ip,ntr))
3964                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3965         $           /clsigma(indmax(cltry(ip,ntr)))
3966                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3967                   ypu(ip,ntr)      = corr
3968                  
3969                endif
3970    
3971           endif                     endif          
3972    
3973        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)  
3974    
3975          if(DEBUG.eq.1)then
3976             print*,'> STORING TRACK ',ntr
3977             print*,'clusters: '
3978             do ip=1,6
3979                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3980             enddo
3981             print*,'dedx: '
3982             do ip=1,6
3983                print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3984             enddo
3985          endif
3986    
3987        end        end
3988    
3989        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*****************************************************  
3990    
3991  *     -------------------------------------------------------  *     -------------------------------------------------------
3992  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3846  c*************************************** Line 3995  c***************************************
3995  *     -------------------------------------------------------  *     -------------------------------------------------------
3996    
3997        include 'commontracker.f'        include 'commontracker.f'
       include 'level1.f'  
       include 'level2.f'  
3998        include 'calib.f'        include 'calib.f'
3999          include 'level1.f'
4000        include 'common_momanhough.f'        include 'common_momanhough.f'
4001          include 'level2.f'
4002        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
4003    
4004  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
4005        nclsx = 0        nclsx = 0
4006        nclsy = 0        nclsy = 0
4007    
4008          do iv = 1,nviews
4009    c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
4010             good2(iv) = good2(iv) + mask_view(iv)*2**8
4011          enddo
4012    
4013          if(DEBUG.eq.1)then
4014             print*,'> STORING SINGLETS '
4015          endif
4016    
4017        do icl=1,nclstr1        do icl=1,nclstr1
4018    
4019             ip=nplanes-npl(VIEW(icl))+1            
4020            
4021           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
4022              ip=nplanes-npl(VIEW(icl))+1              
4023              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
4024    
4025                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4026                 planex(nclsx) = ip                 planex(nclsx) = ip
4027                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4028                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
4029                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
4030                   sxbad(nclsx)  = nbadstrips(1,icl)
4031                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
4032                  
4033    
4034                 do is=1,2                 do is=1,2
4035  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4036                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4037                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
4038                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
4039                 enddo                 enddo
 c$$$               print*,'nclsx         ',nclsx  
 c$$$               print*,'planex(nclsx) ',planex(nclsx)  
 c$$$               print*,'sgnlxs(nclsx) ',sgnlxs(nclsx)  
 c$$$               print*,'xs(1,nclsx)   ',xs(1,nclsx)  
 c$$$               print*,'xs(2,nclsx)   ',xs(2,nclsx)  
4040              else                          !=== Y views              else                          !=== Y views
4041                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4042                 planey(nclsy) = ip                 planey(nclsy) = ip
4043                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4044                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
4045                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
4046                   sybad(nclsy)  = nbadstrips(1,icl)
4047                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
4048    
4049    
4050                 do is=1,2                 do is=1,2
4051  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4052                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4053                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
4054                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
4055                 enddo                 enddo
 c$$$               print*,'nclsy         ',nclsy  
 c$$$               print*,'planey(nclsy) ',planey(nclsy)  
 c$$$               print*,'sgnlys(nclsy) ',sgnlys(nclsy)  
 c$$$               print*,'ys(1,nclsy)   ',ys(1,nclsy)  
 c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  
4056              endif              endif
4057           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
4058    
4059  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4060           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
4061    *     --------------------------------------------------
4062    *     per non perdere la testa...
4063    *     whichtrack(icl) e` una variabile del common level1
4064    *     che serve solo per sapere quali cluster sono stati
4065    *     associati ad una traccia, e permettere di salvare
4066    *     solo questi nell'albero di uscita
4067    *     --------------------------------------------------
4068                    
4069        enddo        enddo
4070        end        end
4071    
4072    ***************************************************
4073    *                                                 *
4074    *                                                 *
4075    *                                                 *
4076    *                                                 *
4077    *                                                 *
4078    *                                                 *
4079    **************************************************
4080    
4081          subroutine fill_hough
4082    
4083    *     -------------------------------------------------------
4084    *     This routine fills the  variables related to the hough
4085    *     transform, for the debig n-tuple
4086    *     -------------------------------------------------------
4087    
4088          include 'commontracker.f'
4089          include 'level1.f'
4090          include 'common_momanhough.f'
4091          include 'common_hough.f'
4092          include 'level2.f'
4093    
4094          if(.false.
4095         $     .or.ntrpt.gt.ntrpt_max_nt
4096         $     .or.ndblt.gt.ndblt_max_nt
4097         $     .or.NCLOUDS_XZ.gt.ncloxz_max
4098         $     .or.NCLOUDS_yZ.gt.ncloyz_max
4099         $     )then
4100             ntrpt_nt=0
4101             ndblt_nt=0
4102             NCLOUDS_XZ_nt=0
4103             NCLOUDS_YZ_nt=0
4104          else
4105             ndblt_nt=ndblt
4106             ntrpt_nt=ntrpt
4107             if(ndblt.ne.0)then
4108                do id=1,ndblt
4109                   alfayz1_nt(id)=alfayz1(id) !Y0
4110                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4111                enddo
4112             endif
4113             if(ndblt.ne.0)then
4114                do it=1,ntrpt
4115                   alfaxz1_nt(it)=alfaxz1(it) !X0
4116                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4117                   alfaxz3_nt(it)=alfaxz3(it) !1/r
4118                enddo
4119             endif
4120             nclouds_yz_nt=nclouds_yz
4121             nclouds_xz_nt=nclouds_xz
4122             if(nclouds_yz.ne.0)then
4123                nnn=0
4124                do iyz=1,nclouds_yz
4125                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4126                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4127                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4128                   nnn=nnn+ptcloud_yz(iyz)
4129                enddo
4130                do ipt=1,min(ndblt_max_nt,nnn)
4131                   db_cloud_nt(ipt)=db_cloud(ipt)
4132                 enddo
4133             endif
4134             if(nclouds_xz.ne.0)then
4135                nnn=0
4136                do ixz=1,nclouds_xz
4137                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4138                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4139                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4140                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4141                   nnn=nnn+ptcloud_xz(ixz)              
4142                enddo
4143                do ipt=1,min(ntrpt_max_nt,nnn)
4144                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4145                 enddo
4146             endif
4147          endif
4148          end
4149          

Legend:
Removed from v.1.8  
changed lines
  Added in v.1.42

  ViewVC Help
Powered by ViewVC 1.1.23