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

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

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

revision 1.1 by mocchiut, Fri May 19 13:15:54 2006 UTC revision 1.43 by pam-fi, Mon Dec 14 10:51:40 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                
       logical DEBUG  
       common/dbg/DEBUG  
   
56  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
57  *     STEP 1  *     STEP 1
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 44  Line 72 
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                
       logical DEBUG  
       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      iflag=0 Line 269  c      iflag=0
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
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           call mini_2(jstep,ifail)           iprint=0
376           if(ifail.ne.0) then  c         if(DEBUG.EQ.1)iprint=1
377              if(DEBUG)then           if(VERBOSE.EQ.1)iprint=1
378             if(DEBUG.EQ.1)iprint=2
379             call mini2(jstep,ifail,iprint)
380    cc         if(ifail.ne.0) then
381             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 *** (mini_2) '       $              '*** 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 263  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 279  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
473     117        continue
474                trackimage(i)=ncp   !number of matching couples
475                if(ncp>ncpmax)then
476                   ncpmax=ncp
477                   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              enddo
517              if(  iimage.ne.0.and.  
518  c     $           RCHI2_STORE(i).le.CHI2MAX.and.              if(DEBUG.EQ.1)then
519  c     $           RCHI2_STORE(i).gt.0.and.                 print*,'Track candidate ',iimage
      $           .true.)then  
                if(DEBUG)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    
 *     --- and store the results --------------------------------  
          ntrk = ntrk + 1                   !counter of found tracks  
          if(.not.FIMAGE  
      $        .and.iimage.eq.0) image(ntrk)= 0  
          if(.not.FIMAGE  
      $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next  
          if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous  
          call fill_level2_tracks(ntrk)     !==> good2=.true.  
 c         print*,'++++++++++ iimage,fimage,ntrk,image '  
 c     $        ,iimage,fimage,ntrk,image(ntrk)  
529    
530    *     --- and store the results --------------------------------
531           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
532              if(DEBUG)              if(verbose.eq.1)
533       $           print*,       $           print*,
534       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
535       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 318  c     $        ,iimage,fimage,ntrk,image Line 537  c     $        ,iimage,fimage,ntrk,image
537  cc            good2=.false.  cc            good2=.false.
538              goto 880            !fill ntp and go to next event              goto 880            !fill ntp and go to next event
539           endif           endif
540    
541             ntrk = ntrk + 1                   !counter of found tracks
542             if(.not.FIMAGE
543         $        .and.iimage.eq.0) image(ntrk)= 0
544             if(.not.FIMAGE
545         $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
546             if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
547             call fill_level2_tracks(ntrk)     !==> good2=.true.
548    
549    c$$$         if(ntrk.eq.NTRKMAX)then
550    c$$$            if(verbose.eq.1)
551    c$$$     $           print*,
552    c$$$     $           '** warning ** number of identified '//
553    c$$$     $           'tracks exceeds vector dimension '
554    c$$$     $           ,'( ',NTRKMAX,' )'
555    c$$$cc            good2=.false.
556    c$$$            goto 880            !fill ntp and go to next event
557    c$$$         endif
558           if(iimage.ne.0)then           if(iimage.ne.0)then
559              FIMAGE=.true.       !              FIMAGE=.true.       !
560              goto 1212           !>>> fit image-track              goto 1212           !>>> fit image-track
561           endif           endif
562    
 *     --- 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  
 *     **********************************************  
   
   
563    
564   880     return   880     return
565        end        end
566    
567    
   
   
 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$$$  
   
568                
569  ************************************************************  ************************************************************
570  ************************************************************  ************************************************************
# Line 556  c$$$ Line 589  c$$$
589  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
590  *     angx   - Projected angle in x  *     angx   - Projected angle in x
591  *     angy   - Projected angle in y  *     angy   - Projected angle in y
592    *     bfx    - x-component of magnetci field
593    *     bfy    - y-component of magnetic field
594  *  *
595  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
596  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 594  c$$$ Line 629  c$$$
629  *  *
630  *  *
631    
632        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
633    
 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*****************************************************  
634                
635        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
636        include 'level1.f'        include 'level1.f'
637          include 'calib.f'
638        include 'common_align.f'        include 'common_align.f'
639        include 'common_mech.f'        include 'common_mech.f'
640        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
       include 'common_resxy.f'  
   
       logical DEBUG  
       common/dbg/DEBUG  
641    
642        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
643        integer sensor        integer sensor
644        integer viewx,viewy        integer viewx,viewy
645        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
646        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
647          real angx,angy            !X-Y effective angle
648          real bfx,bfy              !X-Y b-field components
649    
650        real stripx,stripy        real stripx,stripy
651    
652          double precision xi,yi,zi
653          double precision xi_A,yi_A,zi_A
654          double precision xi_B,yi_B,zi_B
655        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
656        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
657        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  
658                
659    
660        parameter (ndivx=30)        parameter (ndivx=30)
661    
662    
663                
664        resxPAM = 0        resxPAM = 0
665        resyPAM = 0        resyPAM = 0
666    
667        xPAM = 0.        xPAM = 0.D0
668        yPAM = 0.        yPAM = 0.D0
669        zPAM = 0.        zPAM = 0.D0
670        xPAM_A = 0.        xPAM_A = 0.D0
671        yPAM_A = 0.        yPAM_A = 0.D0
672        zPAM_A = 0.        zPAM_A = 0.D0
673        xPAM_B = 0.        xPAM_B = 0.D0
674        yPAM_B = 0.        yPAM_B = 0.D0
675        zPAM_B = 0.        zPAM_B = 0.D0
676    
677          if(sensor.lt.1.or.sensor.gt.2)then
678             print*,'xyz_PAM   ***ERROR*** wrong input '
679             print*,'sensor ',sensor
680             icx=0
681             icy=0
682          endif
683    
684  *     -----------------  *     -----------------
685  *     CLUSTER X  *     CLUSTER X
686  *     -----------------  *     -----------------      
   
687        if(icx.ne.0)then        if(icx.ne.0)then
688           viewx = VIEW(icx)  
689           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
690           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
691           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
692            c         resxPAM = RESXAV
693           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
694           if(PFAx.eq.'COG1')then  !(1)  
695              stripx = stripx      !(1)           if(
696              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
697           elseif(PFAx.eq.'COG2')then       $        viewx.gt.12.or.
698              stripx = stripx + cog(2,icx)                   $        nldx.lt.1.or.
699              resxPAM = resxPAM*fbad_cog(2,icx)       $        nldx.gt.3.or.
700           elseif(PFAx.eq.'ETA2')then       $        stripx.lt.1.or.
701  c            cog2 = cog(2,icx)       $        stripx.gt.3072.or.
702  c            etacorr = pfa_eta2(cog2,viewx,nldx,angx)                   $        .false.)then
703  c            stripx = stripx + etacorr              print*,'xyz_PAM   ***ERROR*** wrong input '
704              stripx = stripx + pfa_eta2(icx,angx)            !(3)              print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
705              resxPAM = risx_eta2(angx)                       !   (4)              icx = 0
706              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 + pfa_eta3(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 + pfa_eta4(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 + pfa_eta(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  
707           endif           endif
708    
709    *        --------------------------
710    *        magnetic-field corrections
711    *        --------------------------
712             stripx  = stripx +  fieldcorr(viewx,bfy)      
713             angx    = effectiveangle(ax,viewx,bfy)
714            
715             call applypfa(PFAx,icx,angx,corr,res)
716             stripx  = stripx + corr
717             resxPAM = res
718    
719     10   continue
720        endif        endif
721              
722  *     -----------------  *     -----------------
723  *     CLUSTER Y  *     CLUSTER Y
724  *     -----------------  *     -----------------
725    
726        if(icy.ne.0)then        if(icy.ne.0)then
727    
728           viewy = VIEW(icy)           viewy = VIEW(icy)
729           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
730           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
731           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!  c         resyPAM = RESYAV
732             stripy = float(MAXS(icy))
733    
734             if(
735         $        viewy.lt.1.or.        
736         $        viewy.gt.12.or.
737         $        nldy.lt.1.or.
738         $        nldy.gt.3.or.
739         $        stripy.lt.1.or.
740         $        stripy.gt.3072.or.
741         $        .false.)then
742                print*,'xyz_PAM   ***ERROR*** wrong input '
743                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
744                icy = 0
745                goto 20
746             endif
747    
748           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
749              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG.EQ.1) then
750       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
751         $              ,icx,icy
752                endif
753              goto 100              goto 100
754           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 = pfa_eta2(cog2,viewy,nldy,angy)  
 c            stripy = stripy + etacorr  
             stripy = stripy + pfa_eta2(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 + pfa_eta3(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 + pfa_eta4(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 + pfa_eta(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  
755    
756    *        --------------------------
757    *        magnetic-field corrections
758    *        --------------------------
759             stripy  = stripy +  fieldcorr(viewy,bfx)      
760             angy    = effectiveangle(ay,viewy,bfx)
761            
762             call applypfa(PFAy,icy,angy,corr,res)
763             stripy  = stripy + corr
764             resyPAM = res
765    
766     20   continue
767        endif        endif
768    
769          
770  c===========================================================  c===========================================================
771  C     COUPLE  C     COUPLE
772  C===========================================================  C===========================================================
# Line 772  C======================================= Line 775  C=======================================
775  c------------------------------------------------------------------------  c------------------------------------------------------------------------
776  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
777  c------------------------------------------------------------------------  c------------------------------------------------------------------------
778           xi = acoordsi(stripx,viewx)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
779           yi = acoordsi(stripy,viewy)       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
780           zi = 0.              if(DEBUG.EQ.1) then
781                   print*,'xyz_PAM (couple):',
782         $              ' WARNING: false X strip: strip ',stripx
783                endif
784             endif
785             xi = dcoordsi(stripx,viewx)
786             yi = dcoordsi(stripy,viewy)
787             zi = 0.D0
788                    
   
789  c------------------------------------------------------------------------  c------------------------------------------------------------------------
790  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
791  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 810  c--------------------------------------- Line 819  c---------------------------------------
819           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
820           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
821    
822           xPAM_A = 0.           xPAM_A = 0.D0
823           yPAM_A = 0.           yPAM_A = 0.D0
824           zPAM_A = 0.           zPAM_A = 0.D0
825    
826           xPAM_B = 0.           xPAM_B = 0.D0
827           yPAM_B = 0.           yPAM_B = 0.D0
828           zPAM_B = 0.           zPAM_B = 0.D0
829    
830        elseif(        elseif(
831       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 836  C======================================= Line 845  C=======================================
845              nldx = nldy              nldx = nldy
846              viewx = viewy + 1              viewx = viewy + 1
847    
848              yi   = acoordsi(stripy,viewy)              xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
849                yi = dcoordsi(stripy,viewy)
850                zi = 0.D0
851    
852              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
853              yi_A = yi              yi_A = yi
# Line 846  C======================================= Line 857  C=======================================
857              yi_B = yi              yi_B = yi
858              zi_B = 0.              zi_B = 0.
859    
 c            print*,'Y-cl ',icy,stripy,' --> ',yi  
 c            print*,xi_A,' <--> ',xi_B  
860                            
861           elseif(icx.ne.0)then           elseif(icx.ne.0)then
862  c===========================================================  c===========================================================
# Line 858  C======================================= Line 867  C=======================================
867              nldy = nldx              nldy = nldx
868              viewy = viewx - 1              viewy = viewx - 1
869    
870              xi   = acoordsi(stripx,viewx)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
871         $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
872                   if(DEBUG.EQ.1) then
873                      print*,'xyz_PAM (X-singlet):',
874         $                 ' WARNING: false X strip: strip ',stripx
875                   endif
876                endif
877    
878                xi = dcoordsi(stripx,viewx)
879                yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
880                zi = 0.D0
881    
882              xi_A = xi              xi_A = xi
883              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 874  C======================================= Line 893  C=======================================
893                 yi_B = yi                 yi_B = yi
894              endif              endif
895    
 c            print*,'X-cl ',icx,stripx,' --> ',xi  
 c            print*,yi_A,' <--> ',yi_B  
896    
897           else           else
898                if(DEBUG.EQ.1) then
899              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
900              print *,'icx = ',icx                 print *,'icx = ',icx
901              print *,'icy = ',icy                 print *,'icy = ',icy
902                endif
903              goto 100              goto 100
904                            
905           endif           endif
# Line 920  c     N.B. I convert angles from microra Line 938  c     N.B. I convert angles from microra
938       $        + zi_B       $        + zi_B
939       $        + dz(nplx,nldx,sensor)       $        + dz(nplx,nldx,sensor)
940    
941    
942    
943             xrt = xi
944         $        - omega(nplx,nldx,sensor)*yi
945         $        + gamma(nplx,nldx,sensor)*zi
946         $        + dx(nplx,nldx,sensor)
947            
948             yrt = omega(nplx,nldx,sensor)*xi
949         $        + yi
950         $        - beta(nplx,nldx,sensor)*zi
951         $        + dy(nplx,nldx,sensor)
952            
953             zrt = -gamma(nplx,nldx,sensor)*xi
954         $        + beta(nplx,nldx,sensor)*yi
955         $        + zi
956         $        + dz(nplx,nldx,sensor)
957    
958    
959                    
960  c      xrt = xi  c      xrt = xi
961  c      yrt = yi  c      yrt = yi
# Line 930  c     (xPAM,yPAM,zPAM) = measured coordi Line 966  c     (xPAM,yPAM,zPAM) = measured coordi
966  c                        in PAMELA reference system  c                        in PAMELA reference system
967  c------------------------------------------------------------------------  c------------------------------------------------------------------------
968    
969           xPAM = 0.           xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
970           yPAM = 0.           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
971           zPAM = 0.           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
972    c$$$         xPAM = 0.D0
973    c$$$         yPAM = 0.D0
974    c$$$         zPAM = 0.D0
975    
976           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
977           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 943  c--------------------------------------- Line 982  c---------------------------------------
982           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
983                    
984    
 c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  
985    
986        else        else
987                       if(DEBUG.EQ.1) then
988           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
989           print *,'icx = ',icx              print *,'icx = ',icx
990           print *,'icy = ',icy              print *,'icy = ',icy
991                         endif
992        endif        endif
993                    
994    
995    
996   100  continue   100  continue
997        end        end
998    
999    ************************************************************************
1000    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1001    *     (done to be called from c/c++)
1002    ************************************************************************
1003    
1004          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1005    
1006          include 'commontracker.f'
1007          include 'level1.f'
1008          include 'common_mini_2.f'
1009          include 'common_xyzPAM.f'
1010          include 'common_mech.f'
1011          include 'calib.f'
1012          
1013    *     flag to chose PFA
1014    c$$$      character*10 PFA
1015    c$$$      common/FINALPFA/PFA
1016    
1017          integer icx,icy           !X-Y cluster ID
1018          integer sensor
1019          character*4 PFAx,PFAy     !PFA to be used
1020          real ax,ay                !X-Y geometric angle
1021          real bfx,bfy              !X-Y b-field components
1022    
1023          ipx=0
1024          ipy=0      
1025          
1026    c$$$      PFAx = 'COG4'!PFA
1027    c$$$      PFAy = 'COG4'!PFA
1028    
1029    
1030          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1031                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1032         $           ,' do not exists (n.clusters=',nclstr1,')'
1033                icx = -1*icx
1034                icy = -1*icy
1035                return
1036            
1037          endif
1038          
1039          call idtoc(pfaid,PFAx)
1040          call idtoc(pfaid,PFAy)
1041    
1042          
1043          if(icx.ne.0.and.icy.ne.0)then
1044    
1045             ipx=npl(VIEW(icx))
1046             ipy=npl(VIEW(icy))
1047            
1048             if( (nplanes-ipx+1).ne.ip )then
1049                print*,'xyzpam: ***WARNING*** cluster ',icx
1050         $           ,' does not belong to plane: ',ip
1051                icx = -1*icx
1052                return
1053             endif
1054             if( (nplanes-ipy+1).ne.ip )then
1055                print*,'xyzpam: ***WARNING*** cluster ',icy
1056         $           ,' does not belong to plane: ',ip
1057                icy = -1*icy
1058                return
1059             endif
1060    
1061             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1062    
1063             xgood(ip) = 1.
1064             ygood(ip) = 1.
1065             resx(ip)  = resxPAM
1066             resy(ip)  = resyPAM
1067    
1068             xm(ip) = xPAM
1069             ym(ip) = yPAM
1070             zm(ip) = zPAM
1071             xm_A(ip) = 0.D0
1072             ym_A(ip) = 0.D0
1073             xm_B(ip) = 0.D0
1074             ym_B(ip) = 0.D0
1075    
1076    c         zv(ip) = zPAM
1077    
1078          elseif(icx.eq.0.and.icy.ne.0)then
1079    
1080             ipy=npl(VIEW(icy))
1081             if( (nplanes-ipy+1).ne.ip )then
1082                print*,'xyzpam: ***WARNING*** cluster ',icy
1083         $           ,' does not belong to plane: ',ip
1084                icy = -1*icy
1085                return
1086             endif
1087    
1088             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1089            
1090             xgood(ip) = 0.
1091             ygood(ip) = 1.
1092             resx(ip)  = 1000.
1093             resy(ip)  = resyPAM
1094    
1095    c$$$         xm(ip) = -100.
1096    c$$$         ym(ip) = -100.
1097    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1098             xm(ip) = xPAM
1099             ym(ip) = yPAM
1100             zm(ip) = zPAM
1101             xm_A(ip) = xPAM_A
1102             ym_A(ip) = yPAM_A
1103             xm_B(ip) = xPAM_B
1104             ym_B(ip) = yPAM_B
1105    
1106    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1107            
1108          elseif(icx.ne.0.and.icy.eq.0)then
1109    
1110             ipx=npl(VIEW(icx))
1111    
1112             if( (nplanes-ipx+1).ne.ip )then
1113                print*,'xyzpam: ***WARNING*** cluster ',icx
1114         $           ,' does not belong to plane: ',ip
1115                icx = -1*icx
1116                return
1117             endif
1118    
1119             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1120          
1121             xgood(ip) = 1.
1122             ygood(ip) = 0.
1123             resx(ip)  = resxPAM
1124             resy(ip)  = 1000.
1125    
1126    c$$$         xm(ip) = -100.
1127    c$$$         ym(ip) = -100.
1128    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1129             xm(ip) = xPAM
1130             ym(ip) = yPAM
1131             zm(ip) = zPAM
1132             xm_A(ip) = xPAM_A
1133             ym_A(ip) = yPAM_A
1134             xm_B(ip) = xPAM_B
1135             ym_B(ip) = yPAM_B
1136            
1137    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1138    
1139          else
1140    
1141             il = 2
1142             if(lad.ne.0)il=lad
1143             is = 1
1144             if(sensor.ne.0)is=sensor
1145    
1146             xgood(ip) = 0.
1147             ygood(ip) = 0.
1148             resx(ip)  = 1000.
1149             resy(ip)  = 1000.
1150    
1151             xm(ip) = -100.
1152             ym(ip) = -100.          
1153             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1154             xm_A(ip) = 0.
1155             ym_A(ip) = 0.
1156             xm_B(ip) = 0.
1157             ym_B(ip) = 0.
1158    
1159    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1160    
1161          endif
1162    
1163          if(DEBUG.EQ.1)then
1164    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1165             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1166         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1167         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1168          endif
1169          end
1170    
1171  ********************************************************************************  ********************************************************************************
1172  ********************************************************************************  ********************************************************************************
# Line 977  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1188  c         print*,'A-(',xPAM_A,yPAM_A,')
1188  *      *    
1189  ********************************************************************************  ********************************************************************************
1190    
1191        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1192    
1193        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1194    
# Line 986  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1197  c         print*,'A-(',xPAM_A,yPAM_A,')
1197  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1198  *     -----------------------------------  *     -----------------------------------
1199    
1200          real rXPP,rYPP
1201          double precision XPP,YPP
1202        double precision distance,RE        double precision distance,RE
1203        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1204    
1205          XPP=DBLE(rXPP)
1206          YPP=DBLE(rYPP)
1207    
1208  *     ----------------------  *     ----------------------
1209        if (        if (
1210       +     xPAM.eq.0.and.  c$$$     +     xPAM.eq.0.and.
1211       +     yPAM.eq.0.and.  c$$$     +     yPAM.eq.0.and.
1212       +     zPAM.eq.0.and.  c$$$     +     zPAM.eq.0.and.
1213       +     xPAM_A.ne.0.and.       +     xPAM_A.ne.0.and.
1214       +     yPAM_A.ne.0.and.       +     yPAM_A.ne.0.and.
1215       +     zPAM_A.ne.0.and.       +     zPAM_A.ne.0.and.
# Line 1031  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1247  c         print*,'A-(',xPAM_A,yPAM_A,')
1247           endif                   endif        
1248    
1249           distance=           distance=
1250       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1251    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1252           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1253    
 c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',re  
1254    
1255                    
1256  *     ----------------------  *     ----------------------
1257        elseif(        elseif(
1258       +     xPAM.ne.0.and.  c$$$     +     xPAM.ne.0.and.
1259       +     yPAM.ne.0.and.  c$$$     +     yPAM.ne.0.and.
1260       +     zPAM.ne.0.and.  c$$$     +     zPAM.ne.0.and.
1261       +     xPAM_A.eq.0.and.       +     xPAM_A.eq.0.and.
1262       +     yPAM_A.eq.0.and.       +     yPAM_A.eq.0.and.
1263       +     zPAM_A.eq.0.and.       +     zPAM_A.eq.0.and.
# Line 1056  c$$$         print*,' resolution ',re Line 1270  c$$$         print*,' resolution ',re
1270  *     ----------------------  *     ----------------------
1271                    
1272           distance=           distance=
1273       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1274       $        +       $       +
1275       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1276    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1277    c$$$     $        +
1278    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1279           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1280    
 c$$$         print*,xPAM,yPAM,zPAM  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',resxPAM,resyPAM  
1281                    
1282        else        else
1283                    
          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  
1284        endif          endif  
1285    
1286        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1111  c$$$         print*,' resolution ',resxP Line 1320  c$$$         print*,' resolution ',resxP
1320        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1321        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1322        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1323        real*8 yvvv,xvvv        double precision yvvv,xvvv
1324        double precision xi,yi,zi        double precision xi,yi,zi
1325        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1326        real AA,BB        real AA,BB
1327        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1328    
1329  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1330        real ptoll        real ptoll
1331        data ptoll/150./          !um        data ptoll/150./          !um
1332    
1333        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1334    
1335        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1336        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1136  c$$$         print*,' resolution ',resxP Line 1345  c$$$         print*,' resolution ',resxP
1345  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1346  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1347  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1348                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1349                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1350                 zi = 0.                 zi = 0.D0
1351  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1352  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1353  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1163  c--------------------------------------- Line 1372  c---------------------------------------
1372    
1373                 yvv(iv)=sngl(yvvv)                 yvv(iv)=sngl(yvvv)
1374                 xvv(iv)=sngl(xvvv)                 xvv(iv)=sngl(xvvv)
 c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '  
 c     $              ,iv,xvv(iv),yvv(iv)  
1375              enddo               !end loop on sensor vertexes              enddo               !end loop on sensor vertexes
1376    
1377              dtot=0.              dtot=0.
# Line 1263  c     $              ,iv,xvv(iv),yvv(iv) Line 1470  c     $              ,iv,xvv(iv),yvv(iv)
1470  *     it returns the plane number  *     it returns the plane number
1471  *      *    
1472        include 'commontracker.f'        include 'commontracker.f'
1473          include 'level1.f'
1474  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1475        include 'common_momanhough.f'        include 'common_momanhough.f'
1476                
# Line 1288  c      include 'common_analysis.f' Line 1496  c      include 'common_analysis.f'
1496        is_cp=0        is_cp=0
1497        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1498        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
       if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  
1499    
1500        return        return
1501        end        end
# Line 1300  c      include 'common_analysis.f' Line 1507  c      include 'common_analysis.f'
1507  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1508  *      *    
1509        include 'commontracker.f'        include 'commontracker.f'
1510          include 'level1.f'
1511  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1512        include 'common_momanhough.f'        include 'common_momanhough.f'
1513                
# Line 1328  c      include 'common_analysis.f' Line 1536  c      include 'common_analysis.f'
1536  *     positive if sensor =2  *     positive if sensor =2
1537  *  *
1538        include 'commontracker.f'        include 'commontracker.f'
1539          include 'level1.f'
1540  c      include 'calib.f'  c      include 'calib.f'
1541  c      include 'level1.f'  c      include 'level1.f'
1542  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1357  c      include 'common_analysis.f' Line 1566  c      include 'common_analysis.f'
1566  *************************************************************************  *************************************************************************
1567  *************************************************************************  *************************************************************************
1568  *************************************************************************  *************************************************************************
 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  
1569                
1570    
1571  ***************************************************  ***************************************************
# Line 1639  c$$$      end Line 1580  c$$$      end
1580        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1581    
1582        include 'commontracker.f'        include 'commontracker.f'
1583          include 'level1.f'
1584        include 'common_momanhough.f'        include 'common_momanhough.f'
1585        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1586        include 'calib.f'        include 'calib.f'
1587        include 'level1.f'  c      include 'level1.f'
   
       logical DEBUG  
       common/dbg/DEBUG  
1588    
1589  *     output flag  *     output flag
1590  *     --------------  *     --------------
# Line 1654  c$$$      end Line 1593  c$$$      end
1593  *     --------------  *     --------------
1594        integer iflag        integer iflag
1595    
1596        integer badseed,badcl        integer badseed,badclx,badcly
1597    
1598          iflag = iflag
1599          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1600    
1601    cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
1602    
1603  *     init variables  *     init variables
       ncp_tot=0  
1604        do ip=1,nplanes        do ip=1,nplanes
1605           do ico=1,ncouplemax           do ico=1,ncouplemax
1606              clx(ip,ico)=0              clx(ip,ico)=0
# Line 1670  c$$$      end Line 1613  c$$$      end
1613           ncls(ip)=0           ncls(ip)=0
1614        enddo        enddo
1615        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1616           cl_single(icl)=1           cl_single(icl) = 1     !all are single
1617           cl_good(icl)=0           cl_good(icl)   = 0     !all are bad
1618          enddo
1619          do iv=1,nviews
1620             ncl_view(iv)  = 0
1621             mask_view(iv) = 0      !all included
1622        enddo        enddo
1623                
1624    *     count number of cluster per view
1625          do icl=1,nclstr1
1626             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1627          enddo
1628    *     mask views with too many clusters
1629          do iv=1,nviews
1630             if( ncl_view(iv).gt. nclusterlimit)then
1631    c            mask_view(iv) = 1
1632                mask_view(iv) = mask_view(iv) + 2**0
1633                if(DEBUG.EQ.1)
1634         $        print*,' * WARNING * cl_to_couple: n.clusters > '
1635         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1636             endif
1637          enddo
1638    
1639    
1640  *     start association  *     start association
1641        ncouples=0        ncouples=0
1642        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1643           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1644                    
1645             if(cl_used(icx).ne.0)goto 10
1646    
1647    *     ----------------------------------------------------
1648    *     jump masked views (X VIEW)
1649    *     ----------------------------------------------------
1650             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1651  *     ----------------------------------------------------  *     ----------------------------------------------------
1652  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1653           if(dedx(icx).lt.dedx_x_min)then  *     ----------------------------------------------------
1654             if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
1655              cl_single(icx)=0              cl_single(icx)=0
1656              goto 10              goto 10
1657           endif           endif
1658    *     ----------------------------------------------------
1659  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1660    *     ----------------------------------------------------
1661           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1662           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1663           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1693  c$$$      end Line 1665  c$$$      end
1665           else           else
1666              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1667           endif           endif
1668           badcl=badseed           badclx=badseed
1669           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1670              ibad=1              ibad=1
1671              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1703  c$$$      end Line 1675  c$$$      end
1675       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1676       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1677              endif              endif
1678              badcl=badcl*ibad              badclx=badclx*ibad
1679           enddo           enddo
1680    *     ----------------------------------------------------
1681    *     >>> eliminato il taglio sulle BAD <<<
1682    *     ----------------------------------------------------
1683  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1684  c     cl_single(icx)=0  c     cl_single(icx)=0
1685  c     goto 10  c     goto 10
# Line 1717  c     endif Line 1692  c     endif
1692                    
1693           do icy=1,nclstr1       !loop on cluster (Y)           do icy=1,nclstr1       !loop on cluster (Y)
1694              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1695    
1696                if(cl_used(icx).ne.0)goto 20
1697                            
1698  *     ----------------------------------------------------  *     ----------------------------------------------------
1699    *     jump masked views (Y VIEW)
1700    *     ----------------------------------------------------
1701                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1702    
1703    *     ----------------------------------------------------
1704  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1705              if(dedx(icy).lt.dedx_y_min)then  *     ----------------------------------------------------
1706                if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
1707                 cl_single(icy)=0                 cl_single(icy)=0
1708                 goto 20                 goto 20
1709              endif              endif
1710    *     ----------------------------------------------------
1711  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1712    *     ----------------------------------------------------
1713              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1714              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1715              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 1732  c     endif Line 1717  c     endif
1717              else              else
1718                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1719              endif              endif
1720              badcl=badseed              badcly=badseed
1721              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1722                 ibad=1                 ibad=1
1723                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 1741  c     endif Line 1726  c     endif
1726       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1727       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1728       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1729                 badcl=badcl*ibad                 badcly=badcly*ibad
1730              enddo              enddo
1731    *     ----------------------------------------------------
1732    *     >>> eliminato il taglio sulle BAD <<<
1733    *     ----------------------------------------------------
1734  c     if(badcl.eq.0)then  c     if(badcl.eq.0)then
1735  c     cl_single(icy)=0  c     cl_single(icy)=0
1736  c     goto 20  c     goto 20
1737  c     endif  c     endif
1738  *     ----------------------------------------------------  *     ----------------------------------------------------
1739                            
               
1740              cl_good(icy)=1                                cl_good(icy)=1                  
1741              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
1742              nldy=nld(MAXS(icy),VIEW(icy))              nldy=nld(MAXS(icy),VIEW(icy))
# Line 1760  c     endif Line 1747  c     endif
1747  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1748              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1749  *     charge correlation  *     charge correlation
1750                 ddd=(dedx(icy)  *     (modified to be applied only below saturation... obviously)
1751       $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
1752                 ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1753                 cut=chcut*sch(nplx,nldx)       $              .and.
1754                 if(abs(ddd).gt.cut)goto 20 !charge not consistent       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1755                       $              .and.
1756                       $              (badclx.eq.1.and.badcly.eq.1)
1757  *     ------------------> COUPLE <------------------       $              .and.
1758  *     check to do not overflow vector dimentions       $              .true.)then
1759                 if(ncp_plane(nplx).gt.ncouplemax)then  
1760                    if(DEBUG)print*,                    ddd=(sgnl(icy)
1761       $                    ' ** warning ** number of identified'//       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1762       $                    ' couples on plane ',nplx,                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1763       $                    ' exceeds vector dimention'//  
1764       $                    ' ( ',ncouplemax,' )'  c                  cut = chcut * sch(nplx,nldx)
1765  c     good2=.false.  
1766  c     goto 880   !fill ntp and go to next event                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1767                    iflag=1       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1768                    return                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1769                      cut = chcut * (16 + sss/50.)
1770    
1771                      if(abs(ddd).gt.cut)then
1772                         goto 20    !charge not consistent
1773                      endif
1774                 endif                 endif
1775                  
1776                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1777                    if(DEBUG)print*,                    if(verbose.eq.1)print*,
1778       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1779       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1780       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1781       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1782  c     good2=.false.  c                  mask_view(nviewx(nplx)) = 2
1783  c     goto 880   !fill ntp and go to next event                      c                  mask_view(nviewy(nply)) = 2
1784                    iflag=1                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1785                    return                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1786                      goto 10
1787                 endif                 endif
1788                                
1789    *     ------------------> COUPLE <------------------
1790                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1791                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1792                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1793                 cl_single(icx)=0                 cl_single(icx)=0
1794                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1795  *     ----------------------------------------------  *     ----------------------------------------------
1796    
1797                endif                              
1798    
1799   20         continue   20         continue
1800           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1801                    
1802   10      continue   10      continue
1803        enddo                     !end loop on clusters(X)        enddo                     !end loop on clusters(X)
1804                
         
1805        do icl=1,nclstr1        do icl=1,nclstr1
1806           if(cl_single(icl).eq.1)then           if(cl_single(icl).eq.1)then
1807              ip=npl(VIEW(icl))              ip=npl(VIEW(icl))
# Line 1815  c     goto 880   !fill ntp and go to nex Line 1809  c     goto 880   !fill ntp and go to nex
1809              cls(ip,ncls(ip))=icl              cls(ip,ncls(ip))=icl
1810           endif           endif
1811        enddo        enddo
1812    
1813    c 80   continue
1814          continue
1815                
1816                
1817        if(DEBUG)then        if(DEBUG.EQ.1)then
1818           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1819           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1820           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'used    ',(cl_used(i),i=1,nclstr1)
1821             print*,'singlets',(cl_single(i),i=1,nclstr1)
1822           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1823        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(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
1824    
1825        logical DEBUG    
1826        common/dbg/DEBUG        if(.not.RECOVER_SINGLETS)goto 81
1827    
1828  *     output flag  *     ////////////////////////////////////////////////
1829  *     --------------  *     PATCH to recover events with less than 3 couples
1830  *     0 = good event  *     ////////////////////////////////////////////////    
1831  *     1 = bad event  *     loop over singlet and create "fake" couples
1832  *     --------------  *     (with clx=0 or cly=0)
1833        integer iflag  *    
   
       integer badseed,badcl  
1834    
1835  *     init variables        if(DEBUG.EQ.1)
1836        ncp_tot=0       $     print*,'>>> Recover singlets '
1837        do ip=1,nplanes       $     ,'(creates fake couples) <<<'
1838           do ico=1,ncouplemax        do icl=1,nclstr1
1839              clx(ip,ico)=0           if(
1840              cly(ip,ico)=0       $        cl_single(icl).eq.1.and.
1841           enddo       $        cl_used(icl).eq.0.and.
1842           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  
           
1843  *     ----------------------------------------------------  *     ----------------------------------------------------
1844  *     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  
1845  *     ----------------------------------------------------  *     ----------------------------------------------------
1846                        if( mask_view(VIEW(icl)).ne.0 ) goto 21
1847           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  
               
1848  *     ----------------------------------------------------  *     ----------------------------------------------------
1849  *     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  
1850  *     ----------------------------------------------------  *     ----------------------------------------------------
1851                               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  
 *     ===========================================================  
1852                                
1853                   nplx=npl(VIEW(icl))
1854    *     ------------------> (FAKE) COUPLE <-----------------
1855                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1856                   clx(nplx,ncp_plane(nplx))=icl
1857                   cly(nplx,ncp_plane(nplx))=0
1858    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1859    *     ----------------------------------------------------
1860                                
1861  *     ------------------> COUPLE <------------------              else                !=== Y views
1862  *     check to do not overflow vector dimentions  *     ----------------------------------------------------
1863                 if(ncp_plane(nplx).gt.ncouplemax)then  *     cut on charge (Y VIEW)
1864                    if(DEBUG)print*,  *     ----------------------------------------------------
1865       $                    ' ** warning ** number of identified'//                 if(sgnl(icl).lt.dedx_y_min) goto 21
      $                    ' couples on plane ',nplx,  
      $                    ' exceeds vector dimention'//  
      $                    ' ( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
                   iflag=1  
                   return  
                endif  
1866                                
1867                 if(ncp_plane(nplx).eq.ncouplemax)then                 nply=npl(VIEW(icl))
1868                    if(DEBUG)print*,  *     ------------------> (FAKE) COUPLE <-----------------
1869       $                 '** warning ** number of identified '//                 ncp_plane(nply) = ncp_plane(nply) + 1
1870       $                 'couples on plane ',nplx,                 clx(nply,ncp_plane(nply))=0
1871       $                 'exceeds vector dimention '                 cly(nply,ncp_plane(nply))=icl
1872       $                 ,'( ',ncouplemax,' )'  c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1873  c     good2=.false.  *     ----------------------------------------------------
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
1874                                
1875                 ncp_plane(nplx) = ncp_plane(nplx) + 1              endif
1876                 clx(nplx,ncp_plane(nplx))=icx           endif                  !end "single" condition
1877                 cly(nply,ncp_plane(nplx))=icy   21      continue
1878                 cl_single(icx)=0        enddo                     !end loop over clusters
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
1879    
1880   20         continue        if(DEBUG.EQ.1)
1881           enddo                  !end loop on clusters(Y)       $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1882            
1883   10      continue  
1884        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  
1885                
1886        do ip=1,6        ncp_tot=0
1887           ncp_tot=ncp_tot+ncp_plane(ip)        do ip=1,NPLANES
1888             ncp_tot = ncp_tot + ncp_plane(ip)
1889        enddo        enddo
1890  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)        if(DEBUG.EQ.1)
1891               $     print*,'n.couple tot:      ',ncp_tot
       if(ncp_tot.gt.ncp_max)then  
          if(DEBUG)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1892                
1893        return        return
1894        end        end
   
 c$$$      subroutine cl_to_couples_2(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$      include 'level1.f'  
 c$$$  
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$*         print*,'icx ',icx,badcl  
 c$$$         if(badcl.eq.0)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$*            print*,'icy ',icy,badcl  
 c$$$            if(badcl.eq.0)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$  
 c$$$c$$$*     charge correlation  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$c               print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(DEBUG)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
1895                
1896  ***************************************************  ***************************************************
1897  *                                                 *  *                                                 *
# Line 2283  c$$$      end Line 1903  c$$$      end
1903  **************************************************  **************************************************
1904    
1905        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1906    
1907        include 'commontracker.f'        include 'commontracker.f'
1908          include 'level1.f'
1909        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1910        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1911        include 'common_mini_2.f'        include 'common_mini_2.f'
1912        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1913    
       logical DEBUG  
       common/dbg/DEBUG  
1914    
1915  *     output flag  *     output flag
1916  *     --------------  *     --------------
# Line 2328  c      double precision xm3,ym3,zm3 Line 1942  c      double precision xm3,ym3,zm3
1942        real xc,zc,radius        real xc,zc,radius
1943  *     -----------------------------  *     -----------------------------
1944    
1945          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1946    
1947    *     --------------------------------------------
1948    *     put a limit to the maximum number of couples
1949    *     per plane, in order to apply hough transform
1950    *     (couples recovered during track refinement)
1951    *     --------------------------------------------
1952          do ip=1,nplanes
1953             if(ncp_plane(ip).gt.ncouplelimit)then
1954                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1955                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1956             endif
1957          enddo
1958    
1959    
1960        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1961        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1962                
1963        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1964           do is1=1,2             !loop on sensors - COPPIA 1  c$$$         print*,'(1) ip ',ip1
1965                c$$$     $        ,mask_view(nviewx(ip1))
1966    c$$$     $        ,mask_view(nviewy(ip1))        
1967             if(  mask_view(nviewx(ip1)).ne.0 .or.
1968         $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1969             do is1=1,2             !loop on sensors - COPPIA 1            
1970              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1971                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1972                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1973                  
1974    c$$$               print*,'(1) ip ',ip1,' icp ',icp1
1975    
1976  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1977                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1978                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1979                 xm1=xPAM                 xm1=xPAM
1980                 ym1=yPAM                 ym1=yPAM
1981                 zm1=zPAM                                   zm1=zPAM                  
1982  c     print*,'***',is1,xm1,ym1,zm1  
1983                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1984                    do is2=1,2    !loop on sensors -ndblt COPPIA 2  c$$$                  print*,'(2) ip ',ip2
1985                        c$$$     $                 ,mask_view(nviewx(ip2))
1986    c$$$     $                 ,mask_view(nviewy(ip2))
1987                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1988         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1989                      do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
1990                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1991                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1992                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1993    
1994    c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
1995    
1996  c                        call xyz_PAM  c                        call xyz_PAM
1997  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1998    c                        call xyz_PAM
1999    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2000                          call xyz_PAM                          call xyz_PAM
2001       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2002                          xm2=xPAM                          xm2=xPAM
2003                          ym2=yPAM                          ym2=yPAM
2004                          zm2=zPAM                          zm2=zPAM
2005                                                    
2006    *                       ---------------------------------------------------
2007    *                       both couples must have a y-cluster
2008    *                       (condition necessary when in RECOVER_SINGLETS mode)
2009    *                       ---------------------------------------------------
2010                            if(icy1.eq.0.or.icy2.eq.0)goto 111
2011    
2012                            if(cl_used(icy1).ne.0)goto 111
2013                            if(cl_used(icy2).ne.0)goto 111
2014    
2015                            
2016  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2017  *     track parameters on Y VIEW  *     track parameters on Y VIEW
2018  *     (2 couples needed)  *     (2 couples needed)
2019  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2020                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2021                             if(DEBUG)print*,                             if(verbose.eq.1)print*,
2022       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2023       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2024       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2025  c                           good2=.false.  c     good2=.false.
2026  c                           goto 880 !fill ntp and go to next event  c     goto 880 !fill ntp and go to next event
2027                               do iv=1,12
2028    c     mask_view(iv) = 3
2029                                  mask_view(iv) = mask_view(iv)+ 2**2
2030                               enddo
2031                             iflag=1                             iflag=1
2032                             return                             return
2033                          endif                          endif
2034                            
2035                            
2036    ccc                        print*,'<doublet> ',icp1,icp2
2037    
2038                          ndblt = ndblt + 1                          ndblt = ndblt + 1
2039  *     store doublet info  *     store doublet info
2040                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
# Line 2381  c                           goto 880 !fi Line 2043  c                           goto 880 !fi
2043                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2044  *     y0 (cm)  *     y0 (cm)
2045                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
2046                                                      
2047  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2048  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
2049  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2050                            if(SECOND_SEARCH)goto 111
2051                          if(                          if(
2052       $                       abs(alfayz2(ndblt)).gt.alfyz2_max       $                       abs(alfayz2(ndblt)).gt.alfyz2_max
2053       $                       .or.       $                       .or.
2054       $                       abs(alfayz1(ndblt)).gt.alfyz1_max       $                       abs(alfayz1(ndblt)).gt.alfyz1_max
2055       $                       )ndblt = ndblt-1       $                       )ndblt = ndblt-1
2056                                                    
2057  c$$$      if(iev.eq.33)then                          
2058  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$$$  
2059  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2060  *     track parameters on Y VIEW - end  *     track parameters on Y VIEW - end
2061  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2062    
2063    
2064                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(icx1.ne.0)then
2065                               if(cl_used(icx1).ne.0)goto 31
2066                            endif
2067                            if(icx2.ne.0)then
2068                               if(cl_used(icx2).ne.0)goto 31
2069                            endif
2070    
2071                            if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2072    
2073                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2074    c$$$                           print*,'(3) ip ',ip3
2075    c$$$     $                          ,mask_view(nviewx(ip3))
2076    c$$$     $                          ,mask_view(nviewy(ip3))                          
2077                               if(  mask_view(nviewx(ip3)).ne.0 .or.
2078         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2079                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
2080                                                                
2081                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
2082                                   icx3=clx(ip3,icp3)                                   icx3=clx(ip3,icp3)
2083                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2084    
2085    c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
2086    
2087    *     ---------------------------------------------------
2088    *     all three couples must have a x-cluster
2089    *     (condition necessary when in RECOVER_SINGLETS mode)
2090    *     ---------------------------------------------------
2091                                     if(
2092         $                                icx1.eq.0.or.
2093         $                                icx2.eq.0.or.
2094         $                                icx3.eq.0.or.
2095         $                                .false.)goto 29
2096                                    
2097                                     if(cl_used(icx1).ne.0)goto 29
2098                                     if(cl_used(icx2).ne.0)goto 29
2099                                     if(cl_used(icx3).ne.0)goto 29
2100    
2101  c                                 call xyz_PAM  c                                 call xyz_PAM
2102  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2103    c                                 call xyz_PAM
2104    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2105                                   call xyz_PAM                                   call xyz_PAM
2106       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2107         $                                ,0.,0.,0.,0.)
2108                                   xm3=xPAM                                   xm3=xPAM
2109                                   ym3=yPAM                                   ym3=yPAM
2110                                   zm3=zPAM                                   zm3=zPAM
2111    
2112    
2113  *     find the circle passing through the three points  *     find the circle passing through the three points
2114                                     iflag_t = DEBUG
2115                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2116       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag_t)
 c     print*,xc,zc,radius  
2117  *     the circle must intersect the reference plane  *     the circle must intersect the reference plane
2118                                   if(  cc                                 if(iflag.ne.0)goto 29
2119  c     $                                 (xc.le.-1.*xclimit.or.                                   if(iflag_t.ne.0)then
2120  c     $                                 xc.ge.xclimit).and.  *     if tricircle fails, evaluate a straight line
2121       $                                radius**2.ge.(ZINI-zc)**2.and.                                      if(DEBUG.eq.1)
2122       $                                iflag.eq.0.and.       $                                   print*,'TRICIRCLE failure'
2123       $                                .true.)then       $                                   ,' >>> straight line'
2124                                                                        radius=0.
2125                                        xc=0.
2126                                        yc=0.
2127                                        
2128                                        SZZ=0.                  
2129                                        SZX=0.
2130                                        SSX=0.
2131                                        SZ=0.
2132                                        S1=0.
2133                                        X0=0.
2134                                        Ax=0.
2135                                        BX=0.
2136                                        DO I=1,3
2137                                           XX = XP(I)
2138                                           SZZ=SZZ+ZP(I)*ZP(I)
2139                                           SZX=SZX+ZP(I)*XX
2140                                           SSX=SSX+XX
2141                                           SZ=SZ+ZP(I)
2142                                           S1=S1+1.
2143                                        ENDDO
2144                                        DET=SZZ*S1-SZ*SZ
2145                                        AX=(SZX*S1-SZ*SSX)/DET
2146                                        BX=(SZZ*SSX-SZX*SZ)/DET
2147                                        X0  = AX*ZINI+BX
2148                                        
2149                                     endif
2150    
2151                                     if(  .not.SECOND_SEARCH.and.
2152         $                                radius**2.lt.(ZINI-zc)**2)goto 29
2153                                                                      
2154  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2155  *     track parameters on X VIEW  *     track parameters on X VIEW
2156  *     (3 couples needed)  *     (3 couples needed)
2157  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2158                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2159                                      if(DEBUG)print*,                                      if(verbose.eq.1)print*,
2160       $                     '** warning ** number of identified '//       $                                   '** warning **'//
2161       $                     'triplets exceeds vector dimention '       $                                   ' number of identified '//
2162       $                    ,'( ',ntrpt_max,' )'       $                                   'triplets exceeds'//
2163  c                                    good2=.false.       $                                   ' vector dimention '
2164  c                                    goto 880 !fill ntp and go to next event       $                                   ,'( ',ntrpt_max,' )'
2165    c     good2=.false.
2166    c     goto 880 !fill ntp and go to next event
2167                                        do iv=1,nviews
2168    c     mask_view(iv) = 4
2169                                           mask_view(iv) =
2170         $                                      mask_view(iv)+ 2**3
2171                                        enddo
2172                                      iflag=1                                      iflag=1
2173                                      return                                      return
2174                                   endif                                   endif
2175    
2176    ccc                                 print*,'<triplet> ',icp1,icp2,icp3
2177                                    
2178                                   ntrpt = ntrpt +1                                   ntrpt = ntrpt +1
2179  *     store triplet info  *     store triplet info
2180                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
2181                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
2182                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
2183                                                                    
2184                                   if(xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2185  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2186                alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
2187                alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2188                alfaxz3(ntrpt) = 1/radius               alfaxz3(ntrpt) = 1/radius
2189                                   else                                  else if(radius.ne.0.and.xc.ge.0)then
2190  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2191                alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)               alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
2192                alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2193                alfaxz3(ntrpt) = -1/radius               alfaxz3(ntrpt) = -1/radius
2194                                   endif                                  else if(radius.eq.0)then
2195                                    *************straight fit
2196                 alfaxz1(ntrpt) = X0
2197                 alfaxz2(ntrpt) = AX
2198                 alfaxz3(ntrpt) = 0.
2199                                    endif
2200    
2201    c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
2202    c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
2203    c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
2204                                        
2205  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2206  ****  reject non phisical triplets                   ****  ****  reject non phisical triplets                   ****
2207  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2208                                   if(                                  if(SECOND_SEARCH)goto 29
2209       $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max                                  if(
2210       $                                .or.       $                               abs(alfaxz2(ntrpt)).gt.
2211       $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max       $                               alfxz2_max
2212       $                                )ntrpt = ntrpt-1       $                               .or.
2213                                         $                               abs(alfaxz1(ntrpt)).gt.
2214                                         $                               alfxz1_max
2215  c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)       $                               )ntrpt = ntrpt-1
2216                                                                    
2217  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2218  *     track parameters on X VIEW - end  *     track parameters on X VIEW - end
2219  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2220                                endif                                
2221     29                           continue
2222                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2223                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2224     30                     continue
2225                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2226   30                  continue  
2227                         31                  continue                    
2228   1                enddo         !end loop on COPPIA 2  c 1                enddo         !end loop on COPPIA 2
2229                     enddo         !end loop on COPPIA 2
2230                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2231     20            continue
2232              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2233                
2234    c 11         continue
2235              continue
2236           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2237        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2238     10   continue
2239        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2240                
2241        if(DEBUG)then        if(DEBUG.EQ.1)then
2242           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2243           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2244        endif        endif
# Line 2514  c     goto 880               !ntp fill Line 2263  c     goto 880               !ntp fill
2263        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2264    
2265        include 'commontracker.f'        include 'commontracker.f'
2266          include 'level1.f'
2267        include 'common_momanhough.f'        include 'common_momanhough.f'
2268        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2269    
       logical DEBUG  
       common/dbg/DEBUG  
2270    
2271  *     output flag  *     output flag
2272  *     --------------  *     --------------
# Line 2537  c     goto 880               !ntp fill Line 2285  c     goto 880               !ntp fill
2285        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2286        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2287    
2288          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2289    
2290  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2291  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2550  c     goto 880               !ntp fill Line 2299  c     goto 880               !ntp fill
2299        distance=0        distance=0
2300        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2301        npt_tot=0        npt_tot=0
2302          nloop=0                  
2303     90   continue                  
2304        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2305           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
2306                            
 c     print*,'--------------'  
 c     print*,'** ',idb1,' **'  
               
2307           do icp=1,ncp_tot           do icp=1,ncp_tot
2308              cp_useds1(icp)=0    !init              cp_useds1(icp)=0    !init
2309              cp_useds2(icp)=0    !init              cp_useds2(icp)=0    !init
# Line 2591  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2339  ccccc if(db_used(idbref).eq.1)goto 1188
2339  *     doublet distance in parameter space  *     doublet distance in parameter space
2340                 distance=                 distance=
2341       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
2342       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2                     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
2343                 distance = sqrt(distance)                 distance = sqrt(distance)
2344                                
 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  
2345                 if(distance.lt.cutdistyz)then                 if(distance.lt.cutdistyz)then
2346    
 c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  
2347                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
2348                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
2349                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
# Line 2619  c     print*,idb1,idb2,distance,' cloud Line 2359  c     print*,idb1,idb2,distance,' cloud
2359    
2360                    temp1 = temp1 + alfayz1(idb2)                    temp1 = temp1 + alfayz1(idb2)
2361                    temp2 = temp2 + alfayz2(idb2)                    temp2 = temp2 + alfayz2(idb2)
 c     print*,'*   idbref,idb2 ',idbref,idb2  
2362                 endif                               endif              
2363                                
2364   1118          continue   1118          continue
2365              enddo               !end loop (2) on DOUBLETS              enddo               !end loop (2) on DOUBLETS
2366                            
2367   1188       continue  c 1188       continue
2368                continue
2369           enddo                  !end loop on... bo?           enddo                  !end loop on... bo?
2370                    
2371           nptloop=npv           nptloop=npv
# Line 2642  c     print*,'*   idbref,idb2 ',idbref,i Line 2382  c     print*,'*   idbref,idb2 ',idbref,i
2382           enddo           enddo
2383           ncpused=0           ncpused=0
2384           do icp=1,ncp_tot           do icp=1,ncp_tot
2385              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2386         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2387         $           .true.)then
2388                 ncpused=ncpused+1                 ncpused=ncpused+1
2389                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2390                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2652  c     print*,'*   idbref,idb2 ',idbref,i Line 2394  c     print*,'*   idbref,idb2 ',idbref,i
2394           do ip=1,nplanes           do ip=1,nplanes
2395              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2396           enddo           enddo
2397  c     print*,'>>>> ',ncpused,npt,nplused          
          if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  
          if(npt.lt.nptyz_min)goto 2228 !next doublet  
2398           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2399                    
2400  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2401  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2402    
2403           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2404              if(DEBUG)print*,              if(verbose.eq.1)print*,
2405       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2406       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2407       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2408  c               good2=.false.  c               good2=.false.
2409  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2410                do iv=1,nviews
2411    c               mask_view(iv) = 5
2412                   mask_view(iv) = mask_view(iv) + 2**4
2413                enddo
2414              iflag=1              iflag=1
2415              return              return
2416           endif           endif
# Line 2682  c     goto 880         !fill ntp and go Line 2426  c     goto 880         !fill ntp and go
2426  c     ptcloud_yz_nt(nclouds_yz)=npt  c     ptcloud_yz_nt(nclouds_yz)=npt
2427           do ipt=1,npt           do ipt=1,npt
2428              db_cloud(npt_tot+ipt) = db_all(ipt)              db_cloud(npt_tot+ipt) = db_all(ipt)
 c     print*,'>> ',ipt,db_all(ipt)  
2429           enddo             enddo  
2430           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2431           if(DEBUG)then           if(DEBUG.EQ.1)then
             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'  
2432              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2433              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2434              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2435              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2436              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2437              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2438  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2439  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)  
2440           endif           endif
2441  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2442  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2704  c$$$     $           ,(db_cloud(iii),iii Line 2444  c$$$     $           ,(db_cloud(iii),iii
2444        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2445                
2446                
2447        if(DEBUG)then        if(nloop.lt.nstepy)then      
2448           print*,'---------------------- '          cutdistyz = cutdistyz+cutystep
2449            nloop     = nloop+1          
2450            goto 90                
2451          endif                    
2452          
2453          if(DEBUG.EQ.1)then
2454           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
          print*,' '  
2455        endif        endif
2456                
2457                
# Line 2730  c$$$     $           ,(db_cloud(iii),iii Line 2474  c$$$     $           ,(db_cloud(iii),iii
2474        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2475    
2476        include 'commontracker.f'        include 'commontracker.f'
2477          include 'level1.f'
2478        include 'common_momanhough.f'        include 'common_momanhough.f'
2479        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2480    
       logical DEBUG  
       common/dbg/DEBUG  
2481    
2482  *     output flag  *     output flag
2483  *     --------------  *     --------------
# Line 2754  c$$$     $           ,(db_cloud(iii),iii Line 2497  c$$$     $           ,(db_cloud(iii),iii
2497        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2498        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2499    
2500          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2501    
2502  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2503  *     classification of TRIPLETS  *     classification of TRIPLETS
2504  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2765  c$$$     $           ,(db_cloud(iii),iii Line 2510  c$$$     $           ,(db_cloud(iii),iii
2510        distance=0        distance=0
2511        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2512        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2513          nloop=0                  
2514     91   continue                  
2515        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2516           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,' **'  
2517                    
2518           do icp=1,ncp_tot           do icp=1,ncp_tot
2519              cp_useds1(icp)=0              cp_useds1(icp)=0
# Line 2802  c         tr_temp(1)=itr1 Line 2547  c         tr_temp(1)=itr1
2547              do itr2=1,ntrpt     !loop (2) on TRIPLETS              do itr2=1,ntrpt     !loop (2) on TRIPLETS
2548                 if(itr2.eq.itr1)goto 11188       !next triplet                 if(itr2.eq.itr1)goto 11188       !next triplet
2549                 if(tr_used(itr2).eq.1)goto 11188 !next triplet                               if(tr_used(itr2).eq.1)goto 11188 !next triplet              
2550    
2551    
2552  *     triplet distance in parameter space  *     triplet distance in parameter space
2553  *     solo i due parametri spaziali per il momemnto  *     solo i due parametri spaziali per il momemnto
2554                 distance=                 distance=
2555       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2556       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2557                 distance = sqrt(distance)                 distance = sqrt(distance)
2558                  
2559                 if(distance.lt.cutdistxz)then  
2560  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  *     ------------------------------------------------------------------------
2561    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2562    *     ------------------------------------------------------------------------
2563    *     (added in august 2007)
2564                   istrimage=0
2565                   if(
2566         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2567         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2568         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2569         $              .true.)istrimage=1
2570    
2571                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2572                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2573                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2574                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
# Line 2829  c     print*,idb1,idb2,distance,' cloud Line 2587  c     print*,idb1,idb2,distance,' cloud
2587                    temp1 = temp1 + alfaxz1(itr2)                    temp1 = temp1 + alfaxz1(itr2)
2588                    temp2 = temp2 + alfaxz2(itr2)                    temp2 = temp2 + alfaxz2(itr2)
2589                    temp3 = temp3 + alfaxz3(itr2)                    temp3 = temp3 + alfaxz3(itr2)
 c     print*,'*   itrref,itr2 ',itrref,itr2,distance  
2590                 endif                               endif              
2591                                
2592  11188          continue  11188          continue
2593              enddo               !end loop (2) on TRIPLETS              enddo               !end loop (2) on TRIPLETS
2594                                                
2595  11888       continue  c11888       continue
2596                continue
2597           enddo                  !end loop on... bo?               enddo                  !end loop on... bo?    
2598                    
2599           nptloop=npv           nptloop=npv
# Line 2850  c     print*,'*   itrref,itr2 ',itrref,i Line 2608  c     print*,'*   itrref,itr2 ',itrref,i
2608  *     1bis)  *     1bis)
2609  *     2) it is not already stored  *     2) it is not already stored
2610  *     ------------------------------------------  *     ------------------------------------------
 c     print*,'check cp_used'  
2611           do ip=1,nplanes           do ip=1,nplanes
2612              hit_plane(ip)=0              hit_plane(ip)=0
2613           enddo           enddo
2614           ncpused=0           ncpused=0
2615           do icp=1,ncp_tot           do icp=1,ncp_tot
2616              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2617         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2618         $           .true.)then
2619                 ncpused=ncpused+1                 ncpused=ncpused+1
2620                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2621                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2866  c     print*,'check cp_used' Line 2625  c     print*,'check cp_used'
2625           do ip=1,nplanes           do ip=1,nplanes
2626              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2627           enddo           enddo
2628           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  
2629                    
2630  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2631  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2632           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2633              if(DEBUG)print*,              if(verbose.eq.1)print*,
2634       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2635       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2636       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2637  c     good2=.false.  c     good2=.false.
2638  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2639                do iv=1,nviews
2640    c               mask_view(iv) = 6
2641                   mask_view(iv) =  mask_view(iv) + 2**5
2642                enddo
2643              iflag=1              iflag=1
2644              return              return
2645           endif           endif
# Line 2896  c     goto 880         !fill ntp and go Line 2657  c     goto 880         !fill ntp and go
2657           enddo           enddo
2658           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2659                    
2660           if(DEBUG)then           if(DEBUG.EQ.1)then
2661              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
2662              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2663              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2664              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2665              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2666              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2667              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cpcloud_xz '
2668         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2669              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)  
2670           endif           endif
2671  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2672  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2673  22288    continue  22288    continue
2674        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2675          
2676        if(DEBUG)then         if(nloop.lt.nstepx)then      
2677           print*,'---------------------- '           cutdistxz=cutdistxz+cutxstep
2678             nloop=nloop+1          
2679             goto 91                
2680           endif                    
2681          
2682          if(DEBUG.EQ.1)then
2683           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
          print*,' '  
2684        endif        endif
2685                
2686                
# Line 2936  c$$$     $           ,(tr_cloud(iii),iii Line 2698  c$$$     $           ,(tr_cloud(iii),iii
2698  **************************************************  **************************************************
2699    
2700        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2701    
2702        include 'commontracker.f'        include 'commontracker.f'
2703          include 'level1.f'
2704        include 'common_momanhough.f'        include 'common_momanhough.f'
2705        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2706        include 'common_mini_2.f'        include 'common_mini_2.f'
2707        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
2708    
2709        logical DEBUG  
       common/dbg/DEBUG  
2710    
2711  *     output flag  *     output flag
2712  *     --------------  *     --------------
# Line 2964  c*************************************** Line 2722  c***************************************
2722  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2723  *     list of matching couples in the combination  *     list of matching couples in the combination
2724  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2725        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2726        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2727  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2728        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2729  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2730  *     variables for track fitting  *     variables for track fitting
2731        double precision AL_INI(5)        double precision AL_INI(5)
       double precision tath  
2732  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2733    
2734          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2735    
2736    
2737        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
2738                
2739        do iyz=1,nclouds_yz       !loop on YZ couds        do iyz=1,nclouds_yz       !loop on YZ clouds
2740           do ixz=1,nclouds_xz    !loop on XZ couds           do ixz=1,nclouds_xz    !loop on XZ clouds
2741                            
2742  *     --------------------------------------------------  *     --------------------------------------------------
2743  *     check of consistency of the clouds  *     check of consistency of the clouds
# Line 2989  c      real fitz(nplanes)        !z coor Line 2746  c      real fitz(nplanes)        !z coor
2746  *     of the two clouds  *     of the two clouds
2747  *     --------------------------------------------------  *     --------------------------------------------------
2748              do ip=1,nplanes              do ip=1,nplanes
2749                 hit_plane(ip)=0                 hit_plane(ip)=0  !n.matching couples (REAL couples, not SINGLETS)
2750                 ncp_match(ip)=0                 ncp_match(ip)=0  !n.matching couples per plane
2751                 do icpp=1,ncouplemax                 do icpp=1,ncouplemax
2752                    cp_match(ip,icpp)=0 !init couple list                    cp_match(ip,icpp)=0 !init couple list
2753                 enddo                 enddo
2754              enddo              enddo
2755              ncp_ok=0              ncp_ok=0            !count n.matching-couples
2756              do icp=1,ncp_tot    !loop on couples              ncpx_ok=0           !count n.matching-couples with x cluster
2757  *     get info on              ncpy_ok=0           !count n.matching-couples with y cluster
2758                 cpintersec(icp)=min(  
2759       $              cpcloud_yz(iyz,icp),  
2760       $              cpcloud_xz(ixz,icp))              do icp=1,ncp_tot    !loop over couples
2761                 if(  
2762       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.                 if(.not.RECOVER_SINGLETS)then
2763       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.  *     ------------------------------------------------------
2764       $              .false.)cpintersec(icp)=0  *     if NOT in RECOVER_SINGLETS mode, take the intersection
2765    *     between xz yz clouds
2766    *     ------------------------------------------------------
2767                      cpintersec(icp)=min(
2768         $                 cpcloud_yz(iyz,icp),
2769         $                 cpcloud_xz(ixz,icp))
2770    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2771    *     ------------------------------------------------------
2772    *     discard the couple if the sensor is in conflict
2773    *     ------------------------------------------------------
2774                      if(
2775         $       (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2776         $       (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2777         $                 .false.)cpintersec(icp)=0
2778                   else
2779    *     ------------------------------------------------------
2780    *     if RECOVER_SINGLETS take the union
2781    *     (otherwise the fake couples formed by singlets would be
2782    *     discarded)    
2783    *     ------------------------------------------------------
2784                      cpintersec(icp)=max(
2785         $                 cpcloud_yz(iyz,icp),
2786         $                 cpcloud_xz(ixz,icp))
2787    c$$$                  if(cpcloud_yz(iyz,icp).gt.0)
2788    c$$$     $                 cpintersec(icp)=cpcloud_yz(iyz,icp)
2789    *     cpintersec is >0 if either yz or xz cloud contains the couple icp
2790                   endif
2791    
2792    c$$$               print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
2793    
2794                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
                   ncp_ok=ncp_ok+1    
2795                                        
2796                    ip=ip_cp(icp)                    ip=ip_cp(icp)
2797                    hit_plane(ip)=1                    hit_plane(ip)=1
2798    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
2799    c$$$     $                 ncp_ok=ncp_ok+1  
2800    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
2801    c$$$     $                 ncpx_ok=ncpx_ok+1
2802    c$$$                  if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
2803    c$$$     $                 ncpy_ok=ncpy_ok+1
2804    
2805                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2806         $                 cpcloud_xz(ixz,icp).gt.0)
2807         $                 ncp_ok=ncp_ok+1  
2808                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2809         $                 cpcloud_xz(ixz,icp).eq.0)
2810         $                 ncpy_ok=ncpy_ok+1  
2811                      if(  cpcloud_yz(iyz,icp).eq.0.and.
2812         $                 cpcloud_xz(ixz,icp).gt.0)
2813         $                 ncpx_ok=ncpx_ok+1  
2814    
2815                    if(cpintersec(icp).eq.1)then                    if(cpintersec(icp).eq.1)then
2816  *     1) only the couple image in sensor 1 matches  *     1) only the couple image in sensor 1 matches
2817                       id=-icp                       id=-icp
# Line 3036  c      real fitz(nplanes)        !z coor Line 2838  c      real fitz(nplanes)        !z coor
2838              do ip=1,nplanes              do ip=1,nplanes
2839                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2840              enddo              enddo
2841                                        
2842              if(nplused.lt.nplxz_min)goto 888 !next doublet              if(nplused.lt.3)goto 888 !next combination
2843              if(ncp_ok.lt.ncpok)goto 888 !next cloud  ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
2844                ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
2845              if(DEBUG)then  *     -----------------------------------------------------------
2846                 print*,'Combination ',iyz,ixz  *     if in RECOVER_SINGLET mode, the two clouds must have
2847       $              ,' db ',ptcloud_yz(iyz)  *     at least ONE intersecting real couple
2848       $              ,' tr ',ptcloud_xz(ixz)  *     -----------------------------------------------------------
2849       $              ,'  -----> # matching couples ',ncp_ok              if(ncp_ok.lt.1)goto 888 !next combination
2850              endif  
2851  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'              if(DEBUG.EQ.1)then
2852  c$$$  print*,'Configurazione cluster XZ'                 print*,'////////////////////////////'
2853  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))                 print*,'Cloud combination (Y,X): ',iyz,ixz
2854  c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))                 print*,' db ',ptcloud_yz(iyz)
2855  c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))                 print*,' tr ',ptcloud_xz(ixz)
2856  c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))                 print*,'  -----> # matching couples ',ncp_ok
2857  c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))                 print*,'  -----> # fake couples (X)',ncpx_ok
2858  c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))                 print*,'  -----> # fake couples (Y)',ncpy_ok
2859  c$$$  print*,'Configurazione cluster YZ'                 do icp=1,ncp_tot
2860  c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))                    print*,'cp ',icp,' >'
2861  c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))       $                 ,' x ',cpcloud_xz(ixz,icp)
2862  c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))       $                 ,' y ',cpcloud_yz(iyz,icp)
2863  c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))       $                 ,' ==> ',cpintersec(icp)
2864  c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))                 enddo
2865  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))              endif
2866  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'                          
2867                            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  
2868                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2869                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2870                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 3114  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2897  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2897                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2898                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2899                                                                
2900                                                                if(DEBUG.eq.1)
2901         $                             print*,'combination: '
2902         $                             ,cp_match(6,icp1)
2903         $                             ,cp_match(5,icp2)
2904         $                             ,cp_match(4,icp3)
2905         $                             ,cp_match(3,icp4)
2906         $                             ,cp_match(2,icp5)
2907         $                             ,cp_match(1,icp6)
2908    
2909    
2910    *                             ---------------------------------------
2911    *                             check if this group of couples has been
2912    *                             already fitted    
2913    *                             ---------------------------------------
2914                                  do ica=1,ntracks
2915                                     isthesame=1
2916                                     do ip=1,NPLANES
2917                                        if(hit_plane(ip).ne.0)then
2918                                           if(  CP_STORE(nplanes-ip+1,ica)
2919         $                                      .ne.
2920         $                                      cp_match(ip,hit_plane(ip)) )
2921         $                                      isthesame=0
2922                                        else
2923                                           if(  CP_STORE(nplanes-ip+1,ica)
2924         $                                      .ne.
2925         $                                      0 )
2926         $                                      isthesame=0
2927                                        endif
2928                                     enddo
2929                                     if(isthesame.eq.1)then
2930                                        if(DEBUG.eq.1)
2931         $                                   print*,'(already fitted)'
2932                                        goto 666 !jump to next combination
2933                                     endif
2934                                  enddo
2935    
2936                                call track_init !init TRACK common                                call track_init !init TRACK common
2937    
2938                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2939                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2940                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2941                                      is=is_cp(id)                                      is=is_cp(id)
# Line 3131  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2949  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2949  *                                   *************************  *                                   *************************
2950  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2951  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2952    c                                    call xyz_PAM(icx,icy,is, !(1)
2953    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2954                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2955       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2956  *                                   *************************  *                                   *************************
2957  *                                   -----------------------------  *                                   -----------------------------
2958                                      xgood(nplanes-ip+1)=1.                                      if(icx.gt.0.and.icy.gt.0)then
2959                                      ygood(nplanes-ip+1)=1.                                         xgood(nplanes-ip+1)=1.
2960                                      xm(nplanes-ip+1)=xPAM                                         ygood(nplanes-ip+1)=1.
2961                                      ym(nplanes-ip+1)=yPAM                                         xm(nplanes-ip+1)=xPAM
2962                                      zm(nplanes-ip+1)=zPAM                                         ym(nplanes-ip+1)=yPAM
2963                                      resx(nplanes-ip+1)=resxPAM                                         zm(nplanes-ip+1)=zPAM
2964                                      resy(nplanes-ip+1)=resyPAM                                         resx(nplanes-ip+1)=resxPAM
2965                                           resy(nplanes-ip+1)=resyPAM
2966                                           if(DEBUG.EQ.1)print*,'(X,Y)'
2967         $                                      ,nplanes-ip+1,xPAM,yPAM
2968                                        else
2969                                           xm_A(nplanes-ip+1) = xPAM_A
2970                                           ym_A(nplanes-ip+1) = yPAM_A
2971                                           xm_B(nplanes-ip+1) = xPAM_B
2972                                           ym_B(nplanes-ip+1) = yPAM_B
2973                                           zm(nplanes-ip+1)
2974         $                                      = (zPAM_A+zPAM_B)/2.
2975                                           resx(nplanes-ip+1) = resxPAM
2976                                           resy(nplanes-ip+1) = resyPAM
2977                                           if(icx.eq.0.and.icy.gt.0)then
2978                                              xgood(nplanes-ip+1)=0.
2979                                              ygood(nplanes-ip+1)=1.
2980                                              resx(nplanes-ip+1) = 1000.
2981                                              if(DEBUG.EQ.1)print*,'(  Y)'
2982         $                                         ,nplanes-ip+1,xPAM,yPAM
2983                                           elseif(icx.gt.0.and.icy.eq.0)then
2984                                              xgood(nplanes-ip+1)=1.
2985                                              ygood(nplanes-ip+1)=0.
2986                                              if(DEBUG.EQ.1)print*,'(X  )'
2987         $                                         ,nplanes-ip+1,xPAM,yPAM
2988                                              resy(nplanes-ip+1) = 1000.
2989                                           else
2990                                              print*,'both icx=0 and icy=0'
2991         $                                         ,' ==> IMPOSSIBLE!!'
2992                                           endif
2993                                        endif
2994  *                                   -----------------------------  *                                   -----------------------------
2995                                   endif                                   endif
2996                                enddo !end loop on planes                                enddo !end loop on planes
2997  *     **********************************************************  *     **********************************************************
2998  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2999  *     **********************************************************  *     **********************************************************
3000    cccc  scommentare se si usa al_ini della nuvola
3001    c$$$                              do i=1,5
3002    c$$$                                 AL(i)=AL_INI(i)
3003    c$$$                              enddo
3004                                  call guess()
3005                                do i=1,5                                do i=1,5
3006                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
3007                                enddo                                enddo
3008                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3009                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3010                                call mini_2(jstep,ifail)                                iprint=0
3011    c                              if(DEBUG.EQ.1)iprint=1
3012                                  if(DEBUG.EQ.1)iprint=2
3013                                  call mini2(jstep,ifail,iprint)
3014                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3015                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
3016                                      print *,                                      print *,
3017       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3018       $                              //'(mini_2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
3019                                        print*,'initial guess: '
3020    
3021                                        print*,'AL_INI(1) = ',AL_INI(1)
3022                                        print*,'AL_INI(2) = ',AL_INI(2)
3023                                        print*,'AL_INI(3) = ',AL_INI(3)
3024                                        print*,'AL_INI(4) = ',AL_INI(4)
3025                                        print*,'AL_INI(5) = ',AL_INI(5)
3026                                   endif                                   endif
3027                                   chi2=-chi2  c                                 chi2=-chi2
3028                                endif                                endif
3029  *     **********************************************************  *     **********************************************************
3030  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
3031  *     **********************************************************  *     **********************************************************
3032    
3033                                if(chi2.le.0.)goto 666                                              if(chi2.le.0.)goto 666              
3034                                  if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3035                                  if(chi2.ne.chi2)goto 666  !OPTIMIZATION
3036    
3037  *     --------------------------  *     --------------------------
3038  *     STORE candidate TRACK INFO  *     STORE candidate TRACK INFO
3039  *     --------------------------  *     --------------------------
3040                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3041                                                                    
3042                                   if(DEBUG)print*,                                   if(verbose.eq.1)print*,
3043       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3044       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3045       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3046  c                                 good2=.false.  c                                 good2=.false.
3047  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3048                                     do iv=1,nviews
3049    c                                    mask_view(iv) = 7
3050                                        mask_view(iv) = mask_view(iv) + 2**6
3051                                     enddo
3052                                   iflag=1                                   iflag=1
3053                                   return                                   return
3054                                endif                                endif
3055                                                                
3056                                ntracks = ntracks + 1                                ntracks = ntracks + 1
3057                                                                
3058  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
3059                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
3060                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3061                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
3062                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))
3063                                   XM_STORE(ip,ntracks)=sngl(xm(ip))                                   XM_STORE(ip,ntracks)=sngl(xm(ip))
3064                                   YM_STORE(ip,ntracks)=sngl(ym(ip))                                   YM_STORE(ip,ntracks)=sngl(ym(ip))
3065                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))
# Line 3205  c$$$     $                               Line 3072  c$$$     $                              
3072                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
3073                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
3074                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
3075    *                                NB! hit_plane is defined from bottom to top
3076                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3077                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
3078       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
3079                                        SENSOR_STORE(nplanes-ip+1,ntracks)
3080         $                              = is_cp(cp_match(ip,hit_plane(ip)))
3081                                        
3082                                        icl=
3083         $                                   clx(ip,icp_cp(
3084         $                                   cp_match(ip,hit_plane(ip)
3085         $                                   )));
3086                                        if(icl.eq.0)
3087         $                                   icl=
3088         $                                   cly(ip,icp_cp(
3089         $                                   cp_match(ip,hit_plane(ip)
3090         $                                   )));
3091    
3092                                        LADDER_STORE(nplanes-ip+1,ntracks)
3093         $                                   = LADDER(icl);
3094                                   else                                   else
3095                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
3096                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
3097                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
3098                                   endif                                   endif
3099                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
3100                                     BY_STORE(ip,ntracks)=0!I dont need it now
3101                                     CLS_STORE(ip,ntracks)=0
3102                                   do i=1,5                                   do i=1,5
3103                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
3104                                   enddo                                   enddo
3105                                enddo                                enddo
3106                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
3107                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
3108                                                                
3109  *     --------------------------------  *     --------------------------------
# Line 3241  c$$$  rchi2=chi2/dble(ndof) Line 3124  c$$$  rchi2=chi2/dble(ndof)
3124                
3125        if(ntracks.eq.0)then        if(ntracks.eq.0)then
3126           iflag=1           iflag=1
3127           return  cc         return
3128        endif        endif
3129                
3130        if(DEBUG)then        if(DEBUG.EQ.1)then
3131           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
3132           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
3133           do i=1,ntracks          do i=1,ntracks
3134              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
3135       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
3136           enddo              ndof=ndof           !(1)
3137           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
3138         $           +int(ygood_store(ii,i)) !(1)
3139              enddo                 !(1)
3140              print*,i,' --- ',rchi2_store(i),' --- '
3141         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3142            enddo
3143            print*,'*****************************************'
3144        endif        endif
3145                
3146                
# Line 3270  c$$$  rchi2=chi2/dble(ndof) Line 3159  c$$$  rchi2=chi2/dble(ndof)
3159    
3160        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3161    
 c******************************************************  
 cccccc 06/10/2005 modified by elena vannuccini ---> (1)  
 cccccc 31/01/2006 modified by elena vannuccini ---> (2)  
 c******************************************************  
3162    
3163        include 'commontracker.f'        include 'commontracker.f'
3164          include 'level1.f'
3165        include 'common_momanhough.f'        include 'common_momanhough.f'
3166        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3167        include 'common_mini_2.f'        include 'common_mini_2.f'
3168        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
       include 'level1.f'  
3169        include 'calib.f'        include 'calib.f'
3170    
       logical DEBUG  
       common/dbg/DEBUG  
   
3171  *     flag to chose PFA  *     flag to chose PFA
3172        character*10 PFA        character*10 PFA
3173        common/FINALPFA/PFA        common/FINALPFA/PFA
3174    
3175          real k(6)
3176          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3177    
3178          real xp,yp,zp
3179          real xyzp(3),bxyz(3)
3180          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3181    
3182          if(DEBUG.EQ.1)print*,'refine_track:'
3183  *     =================================================  *     =================================================
3184  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3185  *                          and  *                          and
# Line 3299  c*************************************** Line 3188  c***************************************
3188        call track_init        call track_init
3189        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3190    
3191             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3192    
3193             xP=XV_STORE(nplanes-ip+1,ibest)
3194             yP=YV_STORE(nplanes-ip+1,ibest)
3195             zP=ZV_STORE(nplanes-ip+1,ibest)
3196             call gufld(xyzp,bxyz)
3197             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3198             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3199    c$$$  bxyz(1)=0
3200    c$$$         bxyz(2)=0
3201    c$$$         bxyz(3)=0
3202    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3203  *     -------------------------------------------------  *     -------------------------------------------------
3204  *     If the plane has been already included, it just  *     If the plane has been already included, it just
3205  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
3206  *     using improved PFAs  *     using improved PFAs
3207  *     -------------------------------------------------  *     -------------------------------------------------
3208           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3209    c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3210    c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3211    c$$$            
3212    c$$$            id=CP_STORE(nplanes-ip+1,ibest)
3213    c$$$            
3214    c$$$            is=is_cp(id)
3215    c$$$            icp=icp_cp(id)
3216    c$$$            if(ip_cp(id).ne.ip)
3217    c$$$     $           print*,'OKKIO!!'
3218    c$$$     $           ,'id ',id,is,icp
3219    c$$$     $           ,ip_cp(id),ip
3220    c$$$            icx=clx(ip,icp)
3221    c$$$            icy=cly(ip,icp)
3222    c$$$c            call xyz_PAM(icx,icy,is,
3223    c$$$c     $           PFA,PFA,
3224    c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
3225    c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
3226    c$$$            call xyz_PAM(icx,icy,is,
3227    c$$$     $           PFA,PFA,
3228    c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
3229    c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
3230    c$$$     $           bxyz(1),
3231    c$$$     $           bxyz(2)
3232    c$$$     $           )
3233    c$$$
3234    c$$$            xm(nplanes-ip+1) = xPAM
3235    c$$$            ym(nplanes-ip+1) = yPAM
3236    c$$$            zm(nplanes-ip+1) = zPAM
3237    c$$$            xgood(nplanes-ip+1) = 1
3238    c$$$            ygood(nplanes-ip+1) = 1
3239    c$$$            resx(nplanes-ip+1) = resxPAM
3240    c$$$            resy(nplanes-ip+1) = resyPAM
3241    c$$$
3242    c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3243    c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3244             if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
3245       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3246                            
3247              id=CP_STORE(nplanes-ip+1,ibest)              id=CP_STORE(nplanes-ip+1,ibest)
# Line 3317  c*************************************** Line 3254  c***************************************
3254       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3255              icx=clx(ip,icp)              icx=clx(ip,icp)
3256              icy=cly(ip,icp)              icy=cly(ip,icp)
3257    c            call xyz_PAM(icx,icy,is,
3258    c     $           PFA,PFA,
3259    c     $           AXV_STORE(nplanes-ip+1,ibest),
3260    c     $           AYV_STORE(nplanes-ip+1,ibest))
3261              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3262       $           PFA,PFA,       $           PFA,PFA,
3263       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3264       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3265  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3266  c$$$  $              'COG2','COG2',       $           bxyz(2)
3267  c$$$  $              0.,       $           )
3268  c$$$  $              0.)  
3269              xm(nplanes-ip+1) = xPAM              if(icx.gt.0.and.icy.gt.0)then
3270              ym(nplanes-ip+1) = yPAM                 xm(nplanes-ip+1) = xPAM
3271              zm(nplanes-ip+1) = zPAM                 ym(nplanes-ip+1) = yPAM
3272              xgood(nplanes-ip+1) = 1                 zm(nplanes-ip+1) = zPAM
3273              ygood(nplanes-ip+1) = 1                 xm_A(nplanes-ip+1) = 0.
3274              resx(nplanes-ip+1) = resxPAM                 ym_A(nplanes-ip+1) = 0.
3275              resy(nplanes-ip+1) = resyPAM                 xm_B(nplanes-ip+1) = 0.
3276                   ym_B(nplanes-ip+1) = 0.
3277  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)                 xgood(nplanes-ip+1) = 1
3278              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                 ygood(nplanes-ip+1) = 1
3279              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                 resx(nplanes-ip+1) = resxPAM
3280                   resy(nplanes-ip+1) = resyPAM
3281                   dedxtrk_x(nplanes-ip+1)=
3282         $              sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3283                   dedxtrk_y(nplanes-ip+1)=
3284         $              sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3285                else
3286                   xm(nplanes-ip+1) = 0.
3287                   ym(nplanes-ip+1) = 0.
3288                   zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
3289                   xm_A(nplanes-ip+1) = xPAM_A
3290                   ym_A(nplanes-ip+1) = yPAM_A
3291                   xm_B(nplanes-ip+1) = xPAM_B
3292                   ym_B(nplanes-ip+1) = yPAM_B
3293                   xgood(nplanes-ip+1) = 0
3294                   ygood(nplanes-ip+1) = 0
3295                   resx(nplanes-ip+1) = 1000.!resxPAM
3296                   resy(nplanes-ip+1) = 1000.!resyPAM
3297                   dedxtrk_x(nplanes-ip+1)= 0
3298                   dedxtrk_y(nplanes-ip+1)= 0
3299                   if(icx.gt.0)then
3300                      xgood(nplanes-ip+1) = 1
3301                      resx(nplanes-ip+1) = resxPAM
3302                      dedxtrk_x(nplanes-ip+1)=
3303         $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3304                   elseif(icy.gt.0)then
3305                      ygood(nplanes-ip+1) = 1
3306                      resy(nplanes-ip+1) = resyPAM
3307                      dedxtrk_y(nplanes-ip+1)=
3308         $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3309                   endif
3310                endif
3311                            
3312    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3313  *     -------------------------------------------------  *     -------------------------------------------------
3314  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3315  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3316  *     -------------------------------------------------  *     -------------------------------------------------
3317    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3318           else                             else                  
3319                                
3320              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
3321              ygood(nplanes-ip+1)=0              ygood(nplanes-ip+1)=0
3322    
3323                CP_STORE(nplanes-ip+1,ibest)=0 !re-init
3324                CLS_STORE(nplanes-ip+1,ibest)=0
3325    
3326                                
3327  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3328  *     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)  
3329              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3330  *     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
3331              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3332    
3333                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3334                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3335  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3336    
3337              if(DEBUG)then              if(DEBUG.EQ.1)then
3338                 print*,                 print*,
3339       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3340       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3368  c            dedxtrk(nplanes-ip+1) = (de Line 3345  c            dedxtrk(nplanes-ip+1) = (de
3345  *     ===========================================  *     ===========================================
3346  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3347  *     ===========================================  *     ===========================================
3348  c            if(DEBUG)print*,'>>>> try to include a new couple'              distmin=100000000.
             distmin=1000000.  
3349              xmm = 0.              xmm = 0.
3350              ymm = 0.              ymm = 0.
3351              zmm = 0.              zmm = 0.
# Line 3382  c            if(DEBUG)print*,'>>>> try t Line 3358  c            if(DEBUG)print*,'>>>> try t
3358              do icp=1,ncp_plane(ip) !loop on couples on plane icp              do icp=1,ncp_plane(ip) !loop on couples on plane icp
3359                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3360                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3361                   if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
3362                 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
3363       $              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
3364       $              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
3365         $              cl_used(icx).ne.0.or. !or the X cluster is already used
3366         $              cl_used(icy).ne.0.or. !or the Y cluster is already used
3367       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3368  *            *          
3369                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3370       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3371       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3372       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3373         $              bxyz(1),
3374         $              bxyz(2)
3375         $              )
3376                                
3377                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3378    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3379                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3380                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)
3381       $              ,' ) normalized distance ',distance       $              print*,'( couple ',id
3382         $              ,' ) distance ',distance
3383                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3384                    xmm = xPAM                    xmm = xPAM
3385                    ymm = yPAM                    ymm = yPAM
# Line 3405  c     $              'ETA2','ETA2', Line 3388  c     $              'ETA2','ETA2',
3388                    rymm = resyPAM                    rymm = resyPAM
3389                    distmin = distance                    distmin = distance
3390                    idm = id                                      idm = id                  
3391  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3392                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3393                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    clincnewc=10*sqrt(rymm**2+rxmm**2
3394         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
3395                 endif                 endif
3396   1188          continue   1188          continue
3397              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3398              if(distmin.le.clinc)then                                if(distmin.le.clincnewc)then    
3399  *              -----------------------------------  *              -----------------------------------
3400                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3401                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3402                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3403                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3404                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3405                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3406                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3407  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3408                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3409  *              -----------------------------------  *              -----------------------------------
3410                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3411                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3412       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
3413                 goto 133         !next plane                 goto 133         !next plane
3414              endif              endif
3415  *     ================================================  *     ================================================
3416  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3417  *                     either from a couple or single  *                     either from a couple or single
3418  *     ================================================  *     ================================================
 c            if(DEBUG)print*,'>>>> try to include a new cluster'  
3419              distmin=1000000.              distmin=1000000.
3420              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3421              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3452  c            if(DEBUG)print*,'>>>> try t Line 3434  c            if(DEBUG)print*,'>>>> try t
3434              do icp=1,ncp_plane(ip) !loop on cluster inside couples              do icp=1,ncp_plane(ip) !loop on cluster inside couples
3435                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3436                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3437                   if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
3438                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3439                 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
3440  *                                                !jump to the next couple  *                                                !jump to the next couple
3441  *----- try cluster x -----------------------------------------------  *----- try cluster x -----------------------------------------------
3442                 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
3443                   if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
3444  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3445    c               call xyz_PAM(icx,0,ist,
3446    c     $              PFA,PFA,
3447    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3448                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3449       $              PFA,PFA,       $              PFA,PFA,
3450       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3451         $              bxyz(1),
3452         $              bxyz(2)
3453         $              )              
3454                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3455  c     if(DEBUG)print*,'normalized distance ',distance  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3456                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)
3457       $              ,' in cp ',id,' ) normalized distance ',distance       $              print*,'( cl-X ',icx
3458         $              ,' in cp ',id,' ) distance ',distance
3459                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3460                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3461                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3477  c     if(DEBUG)print*,'normalized distan Line 3467  c     if(DEBUG)print*,'normalized distan
3467                    rymm = resyPAM                    rymm = resyPAM
3468                    distmin = distance                    distmin = distance
3469                    iclm = icx                    iclm = icx
3470  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3471                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3472                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3473                 endif                                   endif                  
3474  11881          continue  11881          continue
3475  *----- try cluster y -----------------------------------------------  *----- try cluster y -----------------------------------------------
3476                 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
3477                   if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3478  *                                              !jump to the next couple  *                                              !jump to the next couple
3479    c               call xyz_PAM(0,icy,ist,
3480    c     $              PFA,PFA,
3481    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3482                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3483       $              PFA,PFA,       $              PFA,PFA,
3484       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3485         $              bxyz(1),
3486         $              bxyz(2)
3487         $              )
3488                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3489                 if(DEBUG)print*,'( cl-Y ',icy  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3490       $              ,' in cp ',id,' ) normalized distance ',distance                 if(DEBUG.EQ.1)
3491         $              print*,'( cl-Y ',icy
3492         $              ,' in cp ',id,' ) distance ',distance
3493                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3494                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3495                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3503  c     $              'ETA2','ETA2', Line 3501  c     $              'ETA2','ETA2',
3501                    rymm = resyPAM                    rymm = resyPAM
3502                    distmin = distance                    distmin = distance
3503                    iclm = icy                    iclm = icy
3504  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3505                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3506                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3507                 endif                                   endif                  
3508  11882          continue  11882          continue
3509              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3510  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3511              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3512                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3513                 if(cl_used(icl).eq.1.or.     !if the cluster is already used                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
3514       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3515       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3516                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3517                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3518       $                 PFA,PFA,       $                 PFA,PFA,
3519       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3520         $                 bxyz(1),
3521         $                 bxyz(2)
3522         $                 )
3523                 else                         !<---- Y view                 else                         !<---- Y view
3524                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3525       $                 PFA,PFA,       $                 PFA,PFA,
3526       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3527         $                 bxyz(1),
3528         $                 bxyz(2)
3529         $                 )
3530                 endif                 endif
3531    
3532                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3533                 if(DEBUG)print*,'( cl-s ',icl  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3534       $              ,' ) normalized distance ',distance                 if(DEBUG.EQ.1)
3535         $              print*,'( cl-s ',icl
3536         $              ,' ) distance ',distance
3537                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3538                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3539                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3541  c     $                 'ETA2','ETA2', Line 3545  c     $                 'ETA2','ETA2',
3545                    rymm = resyPAM                    rymm = resyPAM
3546                    distmin = distance                      distmin = distance  
3547                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3548                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3549                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3550                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3551                    else          !<---- Y view                    else          !<---- Y view
3552                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3553                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3554                    endif                    endif
3555                 endif                                   endif                  
3556  18882          continue  18882          continue
3557              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3558    
3559              if(distmin.le.clinc)then                                if(iclm.ne.0)then
                 
                CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<      
 *              ----------------------------  
3560                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3561                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3562                    resx(nplanes-ip+1)=rxmm       $                 20*
3563                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3564       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'                 else if(mod(VIEW(iclm),2).ne.0)then
3565                 else                    clincnew=
3566                    YGOOD(nplanes-ip+1)=1.       $                 10*
3567                    resy(nplanes-ip+1)=rymm       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3568                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 endif
3569       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  
3570                   if(distmin.le.clincnew)then  
3571                      
3572                      CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3573    *     ----------------------------
3574                      if(mod(VIEW(iclm),2).eq.0)then
3575                         XGOOD(nplanes-ip+1)=1.
3576                         resx(nplanes-ip+1)=rxmm
3577                         if(DEBUG.EQ.1)
3578         $                    print*,'%%%% included X-cl ',iclm
3579         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3580         $                    ,', dist.= ',distmin
3581         $                    ,', cut ',clincnew,' )'
3582                      else
3583                         YGOOD(nplanes-ip+1)=1.
3584                         resy(nplanes-ip+1)=rymm
3585                         if(DEBUG.EQ.1)
3586         $                    print*,'%%%% included Y-cl ',iclm
3587         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3588         $                    ,', dist.= ', distmin
3589         $                    ,', cut ',clincnew,' )'
3590                      endif
3591    *     ----------------------------
3592                      xm_A(nplanes-ip+1) = xmm_A
3593                      ym_A(nplanes-ip+1) = ymm_A
3594                      xm_B(nplanes-ip+1) = xmm_B
3595                      ym_B(nplanes-ip+1) = ymm_B
3596                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3597                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3598                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3599    *     ----------------------------
3600                 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)  
 *              ----------------------------  
3601              endif              endif
3602           endif           endif
3603   133     continue   133     continue
# Line 3588  c              dedxtrk(nplanes-ip+1) = d Line 3608  c              dedxtrk(nplanes-ip+1) = d
3608        return        return
3609        end        end
3610    
3611    
3612  ***************************************************  ***************************************************
3613  *                                                 *  *                                                 *
3614  *                                                 *  *                                                 *
# Line 3596  c              dedxtrk(nplanes-ip+1) = d Line 3617  c              dedxtrk(nplanes-ip+1) = d
3617  *                                                 *  *                                                 *
3618  *                                                 *  *                                                 *
3619  **************************************************  **************************************************
3620    *
       subroutine clean_XYclouds(ibest,iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
       logical DEBUG  
       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))  
                cl_used(iclx)=1  !tag used clusters  
                cl_used(icly)=1  !tag used clusters  
             elseif(icl.ne.0)then  
                cl_used(icl)=1   !tag used clusters  
             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$$$  
3621    
3622    
3623    
# Line 3779  c$$$ Line 3625  c$$$
3625    
3626        subroutine init_level2        subroutine init_level2
3627    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3628        include 'commontracker.f'        include 'commontracker.f'
3629          include 'level1.f'
3630        include 'common_momanhough.f'        include 'common_momanhough.f'
3631        include 'level2.f'        include 'level2.f'
       include 'level1.f'  
   
   
   
       good2 = 0!.false.  
 c$$$      nev2 = nev1  
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
3632    
3633    *     ---------------------------------
3634    *     variables initialized from level1
3635    *     ---------------------------------
3636          do i=1,nviews
3637             good2(i)=good1(i)
3638             do j=1,nva1_view
3639                vkflag(i,j)=1
3640                if(cnnev(i,j).le.0)then
3641                   vkflag(i,j)=cnnev(i,j)
3642                endif
3643             enddo
3644          enddo
3645    *     ----------------
3646    *     level2 variables
3647    *     ----------------
3648        NTRK = 0        NTRK = 0
3649        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3650           IMAGE(IT)=0           IMAGE(IT)=0
3651           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
          BdL(IT) = 0.  
3652           do ip=1,nplanes           do ip=1,nplanes
3653              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3654              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
3655              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3656              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3657              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3658                TAILX_nt(IP,IT) = 0
3659                TAILY_nt(IP,IT) = 0
3660                XBAD(IP,IT) = 0
3661                YBAD(IP,IT) = 0
3662              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3663              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3664  c*****************************************************              LS(IP,IT) = 0
 cccccc 11/9/2005 modified by david fedele  
3665              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3666              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3667  c******************************************************              CLTRX(IP,IT) = 0
3668                CLTRY(IP,IT) = 0
3669                multmaxx(ip,it) = 0
3670                seedx(ip,it)    = 0
3671                xpu(ip,it)      = 0
3672                multmaxy(ip,it) = 0
3673                seedy(ip,it)    = 0
3674                ypu(ip,it)      = 0
3675           enddo           enddo
3676           do ipa=1,5           do ipa=1,5
3677              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3839  c*************************************** Line 3680  c***************************************
3680              enddo                                enddo                  
3681           enddo                             enddo                  
3682        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3683        nclsx=0        nclsx=0
3684        nclsy=0              nclsy=0      
3685        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3686          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3687          xs(1,ip)=0          xs(1,ip)=0
3688          xs(2,ip)=0          xs(2,ip)=0
3689          sgnlxs(ip)=0          sgnlxs(ip)=0
3690          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3691          ys(1,ip)=0          ys(1,ip)=0
3692          ys(2,ip)=0          ys(2,ip)=0
3693          sgnlys(ip)=0          sgnlys(ip)=0
3694            sxbad(ip)=0
3695            sybad(ip)=0
3696            multmaxsx(ip)=0
3697            multmaxsy(ip)=0
3698        enddo        enddo
 c*******************************************************  
3699        end        end
3700    
3701    
# Line 3872  c*************************************** Line 3710  c***************************************
3710  ************************************************************  ************************************************************
3711    
3712    
3713          subroutine init_hough
3714    
3715          include 'commontracker.f'
3716          include 'level1.f'
3717          include 'common_momanhough.f'
3718          include 'common_hough.f'
3719          include 'level2.f'
3720    
3721          ntrpt_nt=0
3722          ndblt_nt=0
3723          NCLOUDS_XZ_nt=0
3724          NCLOUDS_YZ_nt=0
3725          do idb=1,ndblt_max_nt
3726             db_cloud_nt(idb)=0
3727             alfayz1_nt(idb)=0      
3728             alfayz2_nt(idb)=0      
3729          enddo
3730          do itr=1,ntrpt_max_nt
3731             tr_cloud_nt(itr)=0
3732             alfaxz1_nt(itr)=0      
3733             alfaxz2_nt(itr)=0      
3734             alfaxz3_nt(itr)=0      
3735          enddo
3736          do idb=1,ncloyz_max      
3737            ptcloud_yz_nt(idb)=0    
3738            alfayz1_av_nt(idb)=0    
3739            alfayz2_av_nt(idb)=0    
3740          enddo                    
3741          do itr=1,ncloxz_max      
3742            ptcloud_xz_nt(itr)=0    
3743            alfaxz1_av_nt(itr)=0    
3744            alfaxz2_av_nt(itr)=0    
3745            alfaxz3_av_nt(itr)=0    
3746          enddo                    
3747    
3748          ntrpt=0                  
3749          ndblt=0                  
3750          NCLOUDS_XZ=0              
3751          NCLOUDS_YZ=0              
3752          do idb=1,ndblt_max        
3753            db_cloud(idb)=0        
3754            cpyz1(idb)=0            
3755            cpyz2(idb)=0            
3756            alfayz1(idb)=0          
3757            alfayz2(idb)=0          
3758          enddo                    
3759          do itr=1,ntrpt_max        
3760            tr_cloud(itr)=0        
3761            cpxz1(itr)=0            
3762            cpxz2(itr)=0            
3763            cpxz3(itr)=0            
3764            alfaxz1(itr)=0          
3765            alfaxz2(itr)=0          
3766            alfaxz3(itr)=0          
3767          enddo                    
3768          do idb=1,ncloyz_max      
3769            ptcloud_yz(idb)=0      
3770            alfayz1_av(idb)=0      
3771            alfayz2_av(idb)=0      
3772            do idbb=1,ncouplemaxtot
3773              cpcloud_yz(idb,idbb)=0
3774            enddo                  
3775          enddo                    
3776          do itr=1,ncloxz_max      
3777            ptcloud_xz(itr)=0      
3778            alfaxz1_av(itr)=0      
3779            alfaxz2_av(itr)=0      
3780            alfaxz3_av(itr)=0      
3781            do itrr=1,ncouplemaxtot
3782              cpcloud_xz(itr,itrr)=0
3783            enddo                  
3784          enddo                    
3785          end
3786    ************************************************************
3787    *
3788    *
3789    *
3790    *
3791    *
3792    *
3793    *
3794    ************************************************************
3795    
3796    
3797        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3798    
3799  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3882  c*************************************** Line 3804  c***************************************
3804    
3805            
3806        include 'commontracker.f'        include 'commontracker.f'
3807          include 'level1.f'
3808          include 'common_momanhough.f'
3809        include 'level2.f'        include 'level2.f'
3810        include 'common_mini_2.f'        include 'common_mini_2.f'
3811        real sinth,phi,pig        !(4)        include 'calib.f'
3812    
3813          character*10 PFA
3814          common/FINALPFA/PFA
3815    
3816          real sinth,phi,pig
3817          integer ssensor,sladder
3818        pig=acos(-1.)        pig=acos(-1.)
3819    
       good2=1!.true.  
       chi2_nt(ntr)        = sngl(chi2)  
3820    
3821        phi   = al(4)             !(4)  
3822        sinth = al(3)             !(4)  *     -------------------------------------
3823        if(sinth.lt.0)then        !(4)        chi2_nt(ntr)        = sngl(chi2)
3824           sinth = -sinth         !(4)        nstep_nt(ntr)       = nstep
3825           phi = phi + pig        !(4)  *     -------------------------------------
3826        endif                     !(4)        phi   = al(4)          
3827        npig = aint(phi/(2*pig))  !(4)        sinth = al(3)            
3828        phi = phi - npig*2*pig    !(4)        if(sinth.lt.0)then      
3829        if(phi.lt.0)              !(4)           sinth = -sinth        
3830       $     phi = phi + 2*pig    !(4)           phi = phi + pig      
3831        al(4) = phi               !(4)        endif                    
3832        al(3) = sinth             !(4)        npig = aint(phi/(2*pig))
3833  *****************************************************        phi = phi - npig*2*pig  
3834          if(phi.lt.0)            
3835         $     phi = phi + 2*pig  
3836          al(4) = phi              
3837          al(3) = sinth            
3838        do i=1,5        do i=1,5
3839           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3840           do j=1,5           do j=1,5
3841              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3842           enddo           enddo
 c     print*,al_nt(i,ntr)  
3843        enddo        enddo
3844          *     -------------------------------------      
3845        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3846           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3847           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3919  c     print*,al_nt(i,ntr) Line 3850  c     print*,al_nt(i,ntr)
3850           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3851           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3852           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3853             TAILX_nt(IP,ntr) = 0.
3854             TAILY_nt(IP,ntr) = 0.
3855           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3856           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3857           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3858           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3859           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3860  c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3861           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           factor = sqrt(
3862           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)                 $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3863        enddo       $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3864  c      call CalcBdL(100,xxxx,IFAIL)       $        1. )
3865  c      if(ifps(xxxx).eq.1)BdL(ntr) = xxxx  
3866  c$$$      print*,'xgood(ip,ntr) ',(xgood_nt(ip,ntr),ip=1,6)           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3867  c$$$      print*,'ygood(ip,ntr) ',(ygood_nt(ip,ntr),ip=1,6)           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3868  c$$$      print*,'dedx_x(ip,ntr) ',(dedx_x(ip,ntr),ip=1,6)    
3869  c$$$      print*,'dedx_y(ip,ntr) ',(dedx_y(ip,ntr),ip=1,6)  
3870    ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3871    
3872             ax   = axv_nt(ip,ntr)
3873             ay   = ayv_nt(ip,ntr)
3874             bfx  = BX_STORE(ip,IDCAND)
3875             bfy  = BY_STORE(ip,IDCAND)
3876    c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3877    c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3878    c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3879    c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3880    c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3881    c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3882    
3883             angx = effectiveangle(ax,2*ip,bfy)
3884             angy = effectiveangle(ay,2*ip-1,bfx)
3885            
3886            
3887    
3888             id  = CP_STORE(ip,IDCAND) ! couple id
3889             icl = CLS_STORE(ip,IDCAND)
3890             ssensor = -1
3891             sladder = -1
3892             ssensor = SENSOR_STORE(ip,IDCAND)
3893             sladder = LADDER_STORE(ip,IDCAND)
3894             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3895             LS(IP,ntr)      = ssensor+10*sladder
3896    
3897    c         if(id.ne.0)then
3898    CCCCCC(10 novembre 2009) PATCH X NUCLEI
3899    C     non ho capito perche', ma durante il ritracciamento dei nuclei
3900    C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3901    C     che non viene reinizializzata correttamente e i cluster esclusi
3902    C     dal fit risultano ancora inclusi...
3903    C
3904             cltrx(ip,ntr)   = 0
3905             cltry(ip,ntr)   = 0
3906             if(
3907         $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3908         $        .and.
3909         $        id.ne.0)then
3910    
3911    c           >>> is a couple
3912                cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3913                cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3914    
3915                if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
3916    
3917                   cl_used(cltrx(ip,ntr)) = 1 !tag used clusters          
3918    
3919                   xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3920    
3921                   if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3922         $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3923                  
3924                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3925         $              +10000*mult(cltrx(ip,ntr))
3926                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3927         $              /clsigma(indmax(cltrx(ip,ntr)))
3928                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3929                   xpu(ip,ntr)      = corr
3930    
3931                endif
3932                if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then
3933    
3934                   cl_used(cltry(ip,ntr)) = 1 !tag used clusters          
3935    
3936                   ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3937    
3938                   if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3939         $              dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3940                  
3941                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3942         $              +10000*mult(cltry(ip,ntr))
3943                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3944         $              /clsigma(indmax(cltry(ip,ntr)))
3945                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3946                   ypu(ip,ntr)      = corr
3947                endif
3948    
3949             elseif(icl.ne.0)then
3950    
3951                cl_used(icl) = 1    !tag used clusters          
3952    
3953                if(mod(VIEW(icl),2).eq.0)then
3954                   cltrx(ip,ntr)=icl
3955                   xbad(ip,ntr) = nbadstrips(4,icl)
3956    
3957                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3958    
3959                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3960         $                         +10000*mult(cltrx(ip,ntr))
3961                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3962         $           /clsigma(indmax(cltrx(ip,ntr)))
3963                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3964                   xpu(ip,ntr)      = corr
3965    
3966                elseif(mod(VIEW(icl),2).eq.1)then
3967                   cltry(ip,ntr)=icl
3968                   ybad(ip,ntr) = nbadstrips(4,icl)
3969    
3970                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3971    
3972                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3973         $                         +10000*mult(cltry(ip,ntr))
3974                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3975         $           /clsigma(indmax(cltry(ip,ntr)))
3976                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3977                   ypu(ip,ntr)      = corr
3978                  
3979                endif
3980    
3981             endif          
3982    
3983          enddo
3984    
3985          if(DEBUG.eq.1)then
3986             print*,'> STORING TRACK ',ntr
3987             print*,'clusters: '
3988             do ip=1,6
3989                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3990             enddo
3991             print*,'dedx: '
3992             do ip=1,6
3993                print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3994             enddo
3995          endif
3996    
3997        end        end
3998    
3999        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*****************************************************  
4000    
4001  *     -------------------------------------------------------  *     -------------------------------------------------------
4002  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3952  c*************************************** Line 4005  c***************************************
4005  *     -------------------------------------------------------  *     -------------------------------------------------------
4006    
4007        include 'commontracker.f'        include 'commontracker.f'
       include 'level1.f'  
       include 'level2.f'  
4008        include 'calib.f'        include 'calib.f'
4009          include 'level1.f'
4010        include 'common_momanhough.f'        include 'common_momanhough.f'
4011          include 'level2.f'
4012        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
4013    
4014  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
       good2=1!.true.  
4015        nclsx = 0        nclsx = 0
4016        nclsy = 0        nclsy = 0
4017    
4018          do iv = 1,nviews
4019    c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
4020             good2(iv) = good2(iv) + mask_view(iv)*2**8
4021          enddo
4022    
4023          if(DEBUG.eq.1)then
4024             print*,'> STORING SINGLETS '
4025          endif
4026    
4027        do icl=1,nclstr1        do icl=1,nclstr1
4028    
4029             ip=nplanes-npl(VIEW(icl))+1            
4030            
4031           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
4032              ip=nplanes-npl(VIEW(icl))+1              
4033              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
4034    
4035                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4036                 planex(nclsx) = ip                 planex(nclsx) = ip
4037                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4038                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
4039                   clsx(nclsx)   = icl
4040                   sxbad(nclsx)  = nbadstrips(1,icl)
4041                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
4042                  
4043    
4044                 do is=1,2                 do is=1,2
4045  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4046                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4047                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
4048                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
4049                 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)  
4050              else                          !=== Y views              else                          !=== Y views
4051                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4052                 planey(nclsy) = ip                 planey(nclsy) = ip
4053                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4054                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
4055                   clsy(nclsy)   = icl
4056                   sybad(nclsy)  = nbadstrips(1,icl)
4057                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
4058    
4059    
4060                 do is=1,2                 do is=1,2
4061  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4062                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4063                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
4064                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
4065                 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)  
4066              endif              endif
4067           endif           endif
4068  c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
4069    ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4070             whichtrack(icl) = cl_used(icl)
4071    *     --------------------------------------------------
4072    *     per non perdere la testa...
4073    *     whichtrack(icl) e` una variabile del common level1
4074    *     che serve solo per sapere quali cluster sono stati
4075    *     associati ad una traccia, e permettere di salvare
4076    *     solo questi nell'albero di uscita
4077    *     --------------------------------------------------
4078                    
4079        enddo        enddo
4080        end        end
4081    
4082    ***************************************************
4083    *                                                 *
4084    *                                                 *
4085    *                                                 *
4086    *                                                 *
4087    *                                                 *
4088    *                                                 *
4089    **************************************************
4090    
4091          subroutine fill_hough
4092    
4093    *     -------------------------------------------------------
4094    *     This routine fills the  variables related to the hough
4095    *     transform, for the debig n-tuple
4096    *     -------------------------------------------------------
4097    
4098          include 'commontracker.f'
4099          include 'level1.f'
4100          include 'common_momanhough.f'
4101          include 'common_hough.f'
4102          include 'level2.f'
4103    
4104          if(.false.
4105         $     .or.ntrpt.gt.ntrpt_max_nt
4106         $     .or.ndblt.gt.ndblt_max_nt
4107         $     .or.NCLOUDS_XZ.gt.ncloxz_max
4108         $     .or.NCLOUDS_yZ.gt.ncloyz_max
4109         $     )then
4110             ntrpt_nt=0
4111             ndblt_nt=0
4112             NCLOUDS_XZ_nt=0
4113             NCLOUDS_YZ_nt=0
4114          else
4115             ndblt_nt=ndblt
4116             ntrpt_nt=ntrpt
4117             if(ndblt.ne.0)then
4118                do id=1,ndblt
4119                   alfayz1_nt(id)=alfayz1(id) !Y0
4120                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4121                enddo
4122             endif
4123             if(ndblt.ne.0)then
4124                do it=1,ntrpt
4125                   alfaxz1_nt(it)=alfaxz1(it) !X0
4126                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4127                   alfaxz3_nt(it)=alfaxz3(it) !1/r
4128                enddo
4129             endif
4130             nclouds_yz_nt=nclouds_yz
4131             nclouds_xz_nt=nclouds_xz
4132             if(nclouds_yz.ne.0)then
4133                nnn=0
4134                do iyz=1,nclouds_yz
4135                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4136                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4137                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4138                   nnn=nnn+ptcloud_yz(iyz)
4139                enddo
4140                do ipt=1,min(ndblt_max_nt,nnn)
4141                   db_cloud_nt(ipt)=db_cloud(ipt)
4142                 enddo
4143             endif
4144             if(nclouds_xz.ne.0)then
4145                nnn=0
4146                do ixz=1,nclouds_xz
4147                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4148                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4149                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4150                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4151                   nnn=nnn+ptcloud_xz(ixz)              
4152                enddo
4153                do ipt=1,min(ntrpt_max_nt,nnn)
4154                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4155                 enddo
4156             endif
4157          endif
4158          end
4159          

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.43

  ViewVC Help
Powered by ViewVC 1.1.23