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

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

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

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

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.37

  ViewVC Help
Powered by ViewVC 1.1.23